Fourier transform notation and Klein-Gordon equation The discussion of the Klein-Gordon equation (KGE) assumes a relativistic context. For convenience the KGE reads:
$$(\square - m^2)\phi(x^{\mu})= 0 \tag{1}.$$
From partial differential equation solving (PDES) perspective we can simply apply a Fourier transform (FT) on the scalar field of $(1)$:
$$\phi(x^{\mu}) = \mathcal{N}\int^{+\infty}_{-\infty}\phi(p^{\mu})e^{ip_{\nu}x^{\nu}}d^{4}p \equiv \mathcal{N}\int^{+\infty}_{-\infty}\phi(p^{\mu})e^{ip_{\nu}x^{\nu}}dp^{\mu} \tag{2},$$
where $\mathcal{N}$ is a normalization constant.
Now, $(2)$ is, in my opinion, the most natural way to think about solving $(1)$ using FT method. Because, you have a covariant equation, $(1)$, and a integral with a covariant differential $dp^{\mu} \equiv d^{4}p$, $(2)$. But it seems that, as soon as you want to merge physical intuition with PDES methods, the notations start to become confusing; I will explain.
Let the Klein-Gordon equation be written as:
$$\Big(-\frac{\partial^{2}}{\partial t^{2}} + \nabla^{2} - m^2\Big)\phi(t,\vec{x})= 0 \tag{3}.$$
Then, there are two annoying notations:
$$\phi(t,\vec{x}) = \mathcal{N}\int^{+\infty}_{-\infty}\phi(t,\vec{p})e^{i\vec{p}\cdot \vec{x}}d^{3}p \equiv \mathcal{N}\int^{+\infty}_{-\infty}\phi(t,\vec{p})e^{i\vec{p}\cdot \vec{x}}d\vec{p}, \tag{4}$$
$$\phi(t,\vec{x}) = \mathcal{N}\int^{+\infty}_{-\infty}a_{\vec{p}}(t)e^{i\vec{p}\cdot \vec{x}}d^{3}p. \tag{5}$$
So, I would like to ask:
Since you know you're dealing with a 4-dimension solution, you want to work with covariant expressions. Why in Earth do text books use a 3D FT( $(4), (5)$) and not the (simpler) 4D FT $(2)$ in the first place?

I would like to point out that I know the difference between a 4D formalism and $3+1$-formalism. I mean, the $3+1$-FT is totally different from $(4), (5)$:
$$\phi(x^{\mu}) = \mathcal{N}\int^{+\infty}_{-\infty}\phi(p^{\mu})e^{ip_{\nu}x^{\nu}}d^{4}p  = \phi(t,\vec{x}) =  \mathcal{NC}\int d\omega \int \phi(t,\vec{p})e^{\mp i(\omega t-\vec{p}\cdot\vec{x})}d^{3}p \neq \mathcal{N}\int^{+\infty}_{-\infty}\phi(t,\vec{p})e^{i\vec{p}\cdot \vec{x}}d^{3}p.$$
Where, $\mathcal{N} = \frac{1}{(2\pi)^{3}}$ and $\mathcal{C} = \frac{1}{2\pi}.$

 A: When you perform a Fourier transform in QFT, you actually do use a four-dimensional Fourier transform. However, as Ryder notes, due to the so-called 'on-shell' condition $p_0=\sqrt{\mathbf{p}^2+m^2}$, the fourth component is not independent of the other three. You can give your free particle whatever momentum you'd like, but if you want to be doing physics, it had better obey the energy-momentum dispersion relation above. Thus, due to the constraints, you have to append additional factors when you transform:
\begin{align*}
\phi(x)=\int\frac{\text{d}^4p}{(2\pi)^4}\theta(p_0)\delta(p_0^2-m^2-\mathbf{p}^2)\tilde{\phi}(p)e^{-ipx}.
\end{align*}
Here, $\theta(x)$ is Heaviside's step function. These terms force the particles to remain physical (technical parlance is 'remain on the mass shell' or, more simply, 'remaining on-shell'). The delta function constrains the energy to obey the energy-momentum dispersion relation, and the step function makes sure the energy is positive. This can be shown to simplify to the usual three-dimensional Lorentz-invariant measure (e.g. see here); I encourage you to try and work it out for yourself before looking at the solution there.
A: This is because the temporal Fourier component $p_0$ is not an independent component. You are trying to express your solution as a sum/integral of exponentials. But there can be no exponential $e^{i p^{\mu}x_{\mu}}$ in the sum such that $p^{\mu}p_{\mu}\neq m^2$. Such exponentials do not satisfy the Klein Gordon equation. You can plug an exponential into the KG equation to see that it requires $p^{\mu}p_{\mu}=m^2$ to be a solution.
For each $\vec{p}$, we get two fixed values of the fourth Fourier component $+\omega$ and $-\omega$ if we want to satisfy $p^{\mu}p_{\mu}=m^2$. This is why textbook solutions have an integral over $\vec {p}$ and and also a sum inside the integral specifying the two possible values of $p_0$
The KG equation is $(\Box +m^2)\phi=0$. Forget about Fourier transforms for a moment. Just plug in $e^{ip^{\mu}x_{\mu}}$ into this:
$$(\Box +m^2) e^{ip^{\mu}x_{\mu}} =0$$
$$(-p_{\mu}p^{\mu}+m^2)e^{ip^{\mu}x_{\mu}}=0$$
$$p^{\mu}p_{\mu}=m^2$$
The last equation is a constraint that the KG equation puts on the $p^{\mu}$ that are allowed in the solutions to the KG equation. You can have the three components of $\vec{p}$ vary freely from $-\infty$ to $\infty$ in the general solution. But the  fourth component $p_0$ is constrained to only two values $\sqrt{p_1^2+p_2^2+p_3^2+m^2}$ and $-\sqrt{p_1^2+p_2^2+p_3^2+m^2}$
