Starting with the Klein Gordon in position space,
\begin{align*}
\left(\frac{\partial^2}{\partial t^2} – \nabla^2+m^2\right)\phi(\mathbf{x},t) = 0
\end{align*}
And using the Fourier Transform: $\displaystyle\phi(\mathbf{x},t) = \int \frac{d^3p}{(2\pi)^3}e^{i \mathbf{p} \cdot\mathbf{x}}\phi(\mathbf{p},t)$:
\begin{align*}
\int \frac{d^3p}{(2\pi)^3}\left(\frac{\partial^2}{\partial t^2} – \nabla^2+m^2\right)e^{i \mathbf{p} \cdot\mathbf{x}}\phi(\mathbf{p},t)&=0 \\
\int \frac{d^3p}{(2\pi)^3}e^{i \mathbf{p} \cdot\mathbf{x}}\left(\frac{\partial^2}{\partial t^2} +|\mathbf{p}|^2+m^2\right)\phi(\mathbf{p},t)&=0
\end{align*}
Now I don't understand why we are able to get rid of the integral, to be left with \begin{align*}
\left(\frac{\partial^2}{\partial t^2} +|\mathbf{p}|^2+m^2\right)\phi(\mathbf{p},t)=0
\end{align*}
Fourier Transform – Applying Fourier Transformation to the Klein Gordon Equation
fourier transformklein-gordon-equation
Related Solutions
If at every time $t$, $\phi(\mathbf{x},t)$ is a nice enough function that it has a Fourier transform, then $$\phi(\mathbf{x},t)=\int_{-\infty}^{\infty}\frac{d^{3}k}{(2\pi)^{3}}\widetilde{\phi}(\mathbf{k},t)e^{i\mathbf{k}\cdot\mathbf{x}},$$ where $\widetilde{\phi}(\mathbf{k},t)$ is just the Fourier coefficient at that time $t.$
But now you ask that the whole function $\phi(\mathbf{x},t)$ (at every time) be a solution to the Klein-Gordon equation. Which means the function at different times need to have the right derivatives. If at every time there is a Fourier transform, then there are Fourier coefficients at every time. So if the wave evolves a certain way in time, then the Fourier coefficients need have values at different times need to be just right to let he wave have the right temporal derivatives.
OK, so the Klein-Gordon equation is second order so we can find the initial function $\phi(\mathbf{x},t=0)$ and its Fourier transform, call it $\theta(\mathbf{k})$ and we can take the initial temporal derivative $\partial_t \phi(\mathbf{x},t)\big\vert_{t=0}$ and it's Fourier transform, call it $\omega(\mathbf{k}).$ Then we know the initial Fourier coefficients and we know their derivatives and the second derivative is enforced to make it satisfy the Klein-Gordon equation, so
$\widetilde{\phi}(\mathbf{k},t)=\theta(\mathbf{k})\cos(E_kt) + \frac{1}{E_k}\omega(\mathbf{k})\sin(E_kt).$
Why? Because it has the right initial values and the right initial temporal derivative and it satisfies the Klein-Gordon equation when $E_k=\sqrt{m^2+k^2}.$
So at each time we have Fourier coefficients of our wave. They are designed so that the functions satisfies the Klein-Gordon equation when the Fourier coefficients evolve in time according to a second order equation.
I may have misunderstood your question. But the idea is that from the initial wave and the initial temporal derivative of the wave you get enough initial conditions to know the initial Fourier coefficients and the initial time rate of change of the Fourier coefficients which is all the freedom you have. The rest is determined that the Fourier coefficients have to evolve temporally a certain way so the wave evolves a certain way temporally.
edit to respond to comments
At each time you get a Fourier transform. And you then ask yourself how those Fourier coefficients depend (in time) on each other. In order for the wave to evolve in time by a second order equation, the Fourier transform needs to evolve in time by a second order equation. But the transform gets to do a pointwise evolution that is simpler and that's why we do it.
When you take the Fourier transform at each time, you get a Fourier transform at each time, so you get a $\widetilde{\phi}(\mathbf{k},t)$ that is itself a function of time. So it has a partial derivative $\partial_t \widetilde{\phi}(\mathbf{k},t)\neq 0.$
edit to further respond to comments
Let's say we want to solve $$\left(\partial_t^2-\nabla^2+m^2\right)\phi=.0$$ We start by noting that $\left(A_ke^{i\mathbf{k}\cdot\mathbf{x}}+B_ke^{-i\mathbf{k}\cdot\mathbf{x}}\right)\cos(t\sqrt{m^2+k^2})$ and $\frac{\left(C_ke^{i\mathbf{k}\cdot\mathbf{x}}+D_ke^{-i\mathbf{k}\cdot\mathbf{x}}\right)}{\sqrt{m^2+k^2}}\sin(t\sqrt{m^2+k^2})$ are both solutions for any $\mathbf{k},$ any $A_k, B_k$ and any $C_k, D_k$.
Great.
And the first solution has zero time rate of change at $t=0$ and the second has zero value at $t=0.$ So if our initial wave was $\phi(\mathbf{x},t=0)=3\cos(3x)$ and had an initial time rate of change equal to $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}=4\cos(2y)$ then we know exactly what the solution to $\left(\partial_t^2-\nabla^2+m^2\right)\phi=0$ is: $$\phi(\mathbf{x},t)=3\cos(3x)\cos(t\sqrt{m^2+3^2})+\frac{4}{\sqrt{m^2+2^2}}\cos(2y)\sin(t\sqrt{m^2+2^2}).$$
And any finite linear combination of $e^{\pm i\mathbf k \cdot \mathbf x}$ for the initial value $\phi(\mathbf{x},t=0)$ is as easily manageable by having a finite linear combination of solutions like $e^{\pm i\mathbf{k}\cdot\mathbf{x}}\cos(t\sqrt{m^2+k^2}).$ And similarly if the initial time rate of change $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}$ is a finite linear combination of $e^{\pm i\mathbf k \cdot \mathbf x}$ then we add terms that are a finite linear combination of $\frac{e^{\pm i\mathbf{k}\cdot\mathbf{x}}}{\sqrt{m^2+k^2}}\sin(t\sqrt{m^2+k^2}).$
All we are doing is taking solutions and adding them up in combinations that give us the right initial values and right initial time rate of changes. And it is super easy if the initial conditions happen to be a finite linear combination of sines and cosines.
But wait. What if instead of being a finite linear combination of terms like $e^{\pm i\mathbf k \cdot \mathbf x}$ the initial condition was just a function that has a Fourier Transform? Then you can try the same thing. Write your initial value $\phi(\mathbf{x},t=0)$ and your initial time rate of change $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}$ as inverse Fourier transforms of sines and cosine spatially. Then replace a spatial $\sin(\mathbf k \cdot \mathbf x)$ with the function $\sin(\mathbf k \cdot \mathbf x)\cos(t\sqrt{m^2+k^2})$ and replace the spatial function $\cos(\mathbf k \cdot \mathbf x)$ with the function $\cos(\mathbf k \cdot \mathbf x)\cos(t\sqrt{m^2+k^2}).$ Why? Because each of those satisfies the Klein-Gordon equation. And so for $t=0$ the inverse Fourier transform of those will be the initial value of the wave $\phi(\mathbf{x},t=0).$ So you are taking a $t=0$ spatial Fourier transform of $\phi(\mathbf{x},t=0)$ then replacing every spatial Fourier component $e^{\pm i \mathbf k \cdot \mathbf x}$ with $e^{\pm i \mathbf k \cdot \mathbf x}\cos(t\sqrt{m^2+k^2}).$ This solves the Klein-Gordon equation, has the right initial values and has $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}=0.$
Next, take your initial time rate of change $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}$ as an inverse Fourier transforms of sines and cosine spatially. Then replace a spatial $\sin(\mathbf k \cdot \mathbf x)$ with the function $\frac{\sin(\mathbf k \cdot \mathbf x)}{\sqrt{m^2+k^2}}\sin(t\sqrt{m^2+k^2})$ and replace the spatial function $\cos(\mathbf k \cdot \mathbf x)$ with the function $\frac{\cos(\mathbf k \cdot \mathbf x)}{\sqrt{m^2+k^2}}\sin(t\sqrt{m^2+k^2}).$ Why? Because each of those satisfies the Klein-Gordon equation. And yet for $t=0$ if you take the inverse spatial Fourier transform you get a function that has value zero at $t=0$ and has an initial time rate of change that equals $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}$ of the wave $\phi(\mathbf{x},t=0).$
Why have two solutions? Because this second one also solves the Klein-Gordon equation and has an initial value of zero but has an initial time rate of change that equals $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}.$ And the first one solves the Klein-Gordon equation, has an initial value that equals $\phi(\mathbf{x},t=0)$ and has a zero initial time rate of change.
So if you add those two solutions together you get a function that (1) solves the Klein-Gordon equation (2) has the right initial value and (3) has the right initial time rate of change. That's what you wanted all along.
If you understand that when the initial wave was $\phi(\mathbf{x},t=0)=3\cos(3x)$ and had an initial time rate of change equal to $\partial_t\phi(\mathbf{x},t)\big\vert_{t=0}=4\cos(2y)$ then he solution was $\left(\partial_t^2-\nabla^2+m^2\right)\phi=0$ is: $$\phi(\mathbf{x},t)=3\cos(3x)\cos(t\sqrt{m^2+3^2})+\frac{4}{\sqrt{m^2+2^2}}\cos(2y)\sin(t\sqrt{m^2+2^2}).$$ If you understand that, then everything else is the same idea applied to finite and "infinite" linear combinations of $e^{\pm i \mathbf k \cdot \mathbf x}$ for the initial conditions.
Notation: $x=(t,\boldsymbol x)$; $k=(k_0,\boldsymbol k)$; $kx=k_0t-\boldsymbol k\cdot\boldsymbol x$; $\mathrm dx=\mathrm dt\;\mathrm d^3\boldsymbol x$; etc.
You can in principle perform the Fourier decomposition on both space and time variables, but to do so you'll need several properties of the Dirac's delta funciton:
The first one is: let $\xi\in\mathbb R$; then $$ \delta(f(\xi))=\sum_{f(\xi_i)=0} \frac{\delta(\xi-\xi_i)}{|f'(\xi_i)|} \tag{1} $$ where the sum is over every $\xi_i$ such that $f(\xi_i)=0$, ie, over the roots of $f(\xi)$.
The second one is that, given $g(\xi)$ a known function, the distributional solution of $g(\xi)f(\xi)=0$ is $f(\xi)=h(\xi)\delta(g(\xi))$ for an arbitrary function $h(\xi)$. If you believe these, then the Fourier decomposition is as follows:
Let $\phi(x)$ be the solution of $$ (\partial^2+m^2)\phi(x)=0 $$
Take the Fourier transform of the equation to find $$ (k^2-m^2)\phi(k)=0 \tag{2} $$ where $$ \phi(k)=\int \mathrm dx\ \mathrm e^{ikx} \phi(x) $$
As $\phi(x)$ is a distribution, the solution of $(2)$ is $\phi(k)=h(k)\delta(k^2-m^2)$ for an arbitrary function $h(k)$. Inverting the Fourier Transform, we find $$ \phi(x)=\int\mathrm dk\ \mathrm e^{-ikx}h(k)\delta(k^2-m^2) $$
Next, use $(1)$ to expand the delta over the roots of $k^2-m^2$. These roots are easily found to be $k_0=\pm \omega(\boldsymbol k)$, where $\omega(\boldsymbol k)=+(\boldsymbol k^2+m^2)^{1/2}$. Therefore, it is immediate to get $$ \phi(x)=\int\mathrm dk\ \mathrm e^{-ikx}h(k)\frac{1}{2\omega}\left[\delta(k_0-\omega)+\delta(k_0+\omega)\right] $$ and, after integrating over $\mathrm dk_0$ using the deltas, we find $$ \phi(x)=\int\frac{\mathrm d \boldsymbol k}{2\omega}\ \left[\mathrm e^{-i\omega t} \mathrm e^{i\boldsymbol k\cdot\boldsymbol x}h(\omega,\boldsymbol k)+\mathrm e^{+i\omega t} \mathrm e^{i\boldsymbol k\cdot\boldsymbol x}h(-\omega,\boldsymbol k)\right] $$
Finally, make the change of variable $\boldsymbol k\to-\boldsymbol k$ in the second term, which yeilds the usual expansion $$ \phi(x)=\int\frac{\mathrm d \boldsymbol k}{2\omega}\ \left[\mathrm e^{-ikx}a(\boldsymbol a)+\mathrm e^{+ikx}b^\dagger(\boldsymbol k)\right] $$ where I have defined $a(\boldsymbol k)=h(\omega,\boldsymbol k)$ and $b^\dagger(\boldsymbol k)=h(-\omega,-\boldsymbol k)$.
As you can see, the solution is the same as yours (modulo some irrelevant prefactor that can be reabsorbed into the definition of $h(k)$), though the algebraic procedure to find it is a bit harder.
Best Answer
The functions $e^{i \bf p \cdot \bf x}$ as functions of $\bf x$ are linearly independent for different $\bf p$'s, hence every coefficient in the linear superposition (that is, in the integral) must be zero.