A question on using Fourier decomposition to solve the Klein Gordon equation Given the Klein Gordon equation $$\left(\Box +m^{2}\right)\phi(t,\mathbf{x})=0$$ it is possible to find a solution $\phi(t,\mathbf{x})$ by carrying out a Fourier decomposition of the scalar field $\phi$ at a given instant in time $t$, such that $$\phi(t,\mathbf{x})=\int\frac{d^{3}x}{(2\pi)^{3}}\tilde{\phi}\left(t,\mathbf{k}\right)e^{i\mathbf{k}\cdot\mathbf{x}}$$ where $\tilde{\phi}\left(t,\mathbf{k}\right)$ are the Fourier modes of the corresponding field $\phi(t,\mathbf{x})$.
From this we can calculate the required evolution of the Fourier modes $\tilde{\phi}\left(t,\mathbf{k}\right)$ such that at each instant in time $t$, $\phi(t,\mathbf{x})$ is a solution to the Klein Gordon equation. This can be done, following on from the above, as follows: $$\left(\Box +m^{2}\right)\phi(t,\mathbf{x})=\left(\Box +m^{2}\right)\int\frac{d^{3}x}{(2\pi)^{3}}\tilde{\phi}\left(t,\mathbf{k}\right)e^{i\mathbf{k}\cdot\mathbf{x}}\qquad\qquad\qquad\qquad\qquad\qquad\;\;\,\\ =\int\frac{d^{3}x}{(2\pi)^{3}}\left[\left(\partial^{2}_{t}+\mathbf{k}^{2}+m^{2}\right)\tilde{\phi}\left(t,\mathbf{k}\right)\right]e^{i\mathbf{k}\cdot\mathbf{x}} =0\\ \Rightarrow \left(\partial^{2}_{t}+\mathbf{k}^{2}+m^{2}\right)\tilde{\phi}\left(t,\mathbf{k}\right)=0 \qquad\qquad\qquad$$
Question: This is all well and good, but why is it that in this case we only perform a  Fourier decomposition of the spatial part only, whereas in other cases, such as for finding solutions for propagators (Green's functions), we perform a Fourier decomposition over all 4 spacetime coordinates? [e.g. $$G(x-y)=\int\frac{d^{4}x}{(2\pi)^{4}}\tilde{G}\left(t,\mathbf{k}\right)e^{ik\cdot x}$$ (where in this case $k\cdot x\equiv k_{\mu}x^{\mu}$).]
Is it simply because when we construct the appropriate QFT for a scalar field we do so in the Heisenberg picture, or is there something else to it?
Apologies if this is a really dumb question but it's really been bugging me for a while and I want to get the reasoning straight in my mind!
 A: 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 k)+\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.
