It is noted in Peskin and Schroeder's QFT text that the momenta used in the evaluation of the field operator $\phi(x)$ are "on mass-shell": $p^2=m^2$. Specifically, this is in relation to the expression:
$$\phi(x)=\int{d^3p\over (2\pi)^3}{1\over\sqrt{2E_p}}\bigg\{\bigg(\mathbf{a_p}+{i\over\sqrt{2E_p}}\tilde{j}(p)\bigg)e^{-ip\cdot x}+\;\text{h.c.}\bigg\}\quad \tag{2.64}.$$
Why is this the case? My educated guess is that the particle momenta must be "on shell" because without this condition, $\phi(x)$ will not be guaranteed to have Lorentz invariance.
Particle Creation by a Classical Source (on-shell mass momenta)
causalityfourier transformoperatorsquantum-field-theoryspecial-relativity
Related Solutions
I think the usual way we approach this issue is to first obtain the field equations, in the classical manner, and obtain classical solutions by free wave expansion or Green's functions. For example in equation (2.55) and (2.60), you can see that the propagators, which are derived from quantum mechanics, have "classical" solutions, that is they satisfy the classical field equation (Klein-Gordon). So, your equation of motion will be eqn. (2.61), \begin{align} (\partial^2 + m^2)\phi(x) & = j(x) \end{align} Take the Fourier Transform of this to momentum space, and you obtain, \begin{align} (-p^2 + m^2) \tilde{\phi}(p) = \tilde{j}(p) \end{align} To solve this, you solve for the free Green's function, \begin{align} (-p^2 + m^2) \tilde{\phi}(p) = -i \end{align} Then after you obtain a solution for this $\tilde{\phi}_{\rm free}(p)$, do this: \begin{align} \phi(x) = \int \frac{d^4 p}{(2\pi)^4} \tilde{j}(p) \tilde{\phi}_{\rm free}(p). \end{align} Where is the quantization, you may ask. It's in the solution of the free Green's function.
To explain the second part of your question, you need to include information about when your source was turned on. The $\phi(x)$ over time will also not be the one for the free field, since $\phi$ is a solution to the field equations, as we showed in 1. In other words, \begin{align} \phi(x) \neq \int \frac{d^3 p}{(2\pi)^3} a_{\vec{p}} e^{i p\cdot x} + a_{\vec{p}}^\dagger e^{-i p \cdot x} \end{align}.
Quick answer
My question is how does the presence of nonzero $J(x)$ results in a non-trivial spacetime dependent value of $\langle 0|\phi(x)|0\rangle$?
The equation $\phi(x)=\mathrm e^{-iPx}\phi(0)\mathrm e^{iPx}$ works both for $J=0$ and $J\neq 0$. Therefore, $$ \langle \phi(x)\rangle_J= {}_J\langle 0|\mathrm e^{-iPx}\phi(0)\mathrm e^{iPx}|0\rangle_J $$
What it is not true when $J\neq 0$ is that $P|0\rangle_J\stackrel{\text{no}}{=}0$ (because the source breaks the invariance), and therefore we cannot conclude that $$ \langle \phi(x)\rangle_J\stackrel{\text{no}}{=} {}_J\langle 0|\phi(0)|0\rangle_J $$
Therefore, if $J\neq 0$ the vev depends on position $x$.
To find the explicit dependence of $\langle \phi(x)\rangle_J$ with $x$, instead of using operators, it is easier to work with path integrals: $$ \langle \phi(x)\rangle_J=\frac{\delta}{\delta J(x)}\exp\left[-i\int \mathrm dy\,\mathrm dz\ J^*(y)\Delta(y-z)J(z)\right] $$ which I believe you can calculate yourself (note that the result is proportional to $J(x)$ and so the vev goes to zero as $J\to 0$, as expected).
The (somewhat) bigger picture
The first thing we have to do is to differentiate from internal sources and external ones:
An internal source is a term in the lagrangian that only includes dynamical fields, that is, fields that are part of the equations of motion. For example, you can have a KG theory, $$ \mathcal L\sim (\partial\phi)^2-m^2\phi^2+g\phi^3 $$ where the last term can be said to be an internal source (though the usual terminology is just interaction). This term is internal because it only depends on $\phi$, which is itself a dynamical field (determined from the EoM's). Another (more illustrative) example is the lagrangian for QED, $$ \mathcal L\sim \bar\psi(i\not\partial-m)\psi-F^2+eA_\mu \bar\psi\gamma^\mu\psi $$ Again, the last term is an internal source, because it only depends on dynamical fields, $\psi$ and $A$, which are determined from the EoM. I would like to stress that in general people don't say "internal source" but "interaction" instead.
An exernal source is a function in the lagrangian that is externally determined (fixed), that is, a function that is not dynamical (there is not an equation of motion for that function). Typical examples are the $J$'s that are used in path integrals, $$ \mathcal L\sim (\partial\phi)^2-m^2\phi^2+g\phi^3+\phi(x)J(x) $$ and fixed (background) functions in effective theories, such as, for example, the electromagnetic field in a low energy treatment of the Hydrogen atom: $$ \mathcal L\sim \bar\psi(i\not\partial-m)\psi+eA_\mu \bar\psi\gamma^\mu\psi $$ (here, $A_\mu$ is an external source, because there is not a kinetic term $F^2$ for it, and so the value of $A$ has to be written by hand, say, a Coulomb potential $A_0\sim e/r$.
Note that external sources break the translational invariance of the theory (because of the obvious reason: an external source has a fixed dependence on position, and so the "physics don't look the same everywhere"). Therefore, if there are external sources, $P_\mu|0\rangle\neq 0$ and vev's depend on position, as discussed in the first part of this answer.
On the other hand, internal sources don't break the translational invariance of the theory, because the sources themselves transform together with the fields. This might be easier to understand with an example. Consider first a theory with only internal sources: $$ S=\int\mathrm dx\ (\partial\phi(x))^2-m^2\phi(x)^2-g\phi(x)^3 $$ which, upon a translation $x\to x-a$ transforms into $$ S_a=\int\mathrm dx\ (\partial\phi(x-a))^2-m^2\phi(x-a)^2-g\phi(x-a)^3 $$ which is the same as before, $S_a=S$, because we integrate over all space and $\mathrm d(x-a)=\mathrm dx$.
On the other hand, consider a theory with an external source: $$ S=\int\mathrm dx\ (\partial\phi(x))^2-m^2\phi(x)^2-\phi(x)J(x) $$ which, upon a translation $x\to x-a$ transforms into $$ S_a=\int\mathrm dx\ (\partial\phi(x-a))^2-m^2\phi(x-a)^2-\phi(x-a)J(x) $$ which is not the same as before, because of the $J(x)$ term. The action is not the same as before, and so the translation changed the theory. At this point, you might want to read this answer of mine. In the notation of that post, the $(2)$ derivative of a lagrangian with external sources is non-zero.
To recapitulate,
If there are only internal sources, then the theory is translationally invariant, and so all the vev's are position independent (as can be easily shown using $P_\mu|0\rangle=0$ and $Q_\alpha(x)=\mathrm e^{-iPx}Q_\alpha(x)\mathrm e^{iPx}$, where $Q_\alpha(x)$ is any field). Most of the times we redefine every field $Q_\alpha(x)\to Q_\alpha(x)-\langle Q\rangle$ so that all the vev's are zero (this is relevant for renormalisation). In some cases (e.g., in the case of the Higgs field) a non-zero vev is physically relevant (but only makes sense because of the form of the lagrangian for the Higgs field, and wouldn't make sense for, say, a standard KG field). In any case, if the sources are internal then vev's are constant.
If there are only external sources, then the theory is free. Therefore, the vev's depend on position, but in the limit $J\to 0$ we must have $\langle\phi\rangle\to 0$, as it must be for a free theory.
If there are internal and external sources, the vev's are position-dependent and don't go to zero as the external sources go to zero (and therefore we must renormalise the fields).
Best Answer
TL;DR: The on-shellness of the Fourier decomposition (2.64) comes from the source-full/inhomogeneous Klein-Gordon equation (2.61).
P&S argue this in 2 steps:
First P&S consider the Fourier-transformed source-free/homogeneous Klein-Gordon equation, and argue that its solution $\phi_0$ is on-shell.
Next P&S introduce in eq. (2.63) a classical source $j$ via the retarded Greens function, which can also be brought on an on-shell form: $$\begin{align} D_R(x-y)~=~~~&\int\!\frac{d^4p}{(2\pi)^4}\frac{i}{(p^0+i\epsilon)^2-E_{\bf p}}e^{-ip\cdot(x-y)}\cr ~=~~~&-\int\!\frac{d^3p}{2E_{\bf p}(2\pi)^3} \frac{dp^0}{2\pi i}\left(\frac{1}{p^0+i\epsilon-E_{\bf p}}-\frac{1}{p^0+i\epsilon+E_{\bf p}}\right)e^{-ip\cdot(x-y)}\cr ~\stackrel{\text{Res.Thm.}}{=}&\theta(x^0-y^0)\int\!\frac{d^3p}{2E_{\bf p}(2\pi)^3} \left(\left. e^{-ip\cdot(x-y)}\right|_{p^0=E_{\bf p}} -\left. e^{-ip\cdot(x-y)}\right|_{p^0=-E_{\bf p}} \right)\cr ~=~~~&\theta(x^0-y^0)\int\!\frac{d^3p}{2E_{\bf p}(2\pi)^3} \left.\left( e^{-ip\cdot(x-y)} - e^{ip\cdot(x-y)} \right)\right|_{p^0=E_{\bf p}}, \end{align}\tag{2.58} $$ where the above $i\epsilon$ prescription is consistent with the figure on p.30 between eqs. (2.54-55).
References: