Solution (source code)

= Solution

Use natural units $\hbar=c=1$ and the <Minkowski metric> $g=\operatorname{diag}(1,-1,-1,-1)$. The <canonical momentum> is $\pi=\partial\mathcal L/\partial\dot\phi=\dot\phi$, and the <Legendre transform in mechanics> gives
$$
H=\frac12\int d^3x\,[\pi^2+(\nabla\phi)^2+m^2\phi^2].
$$
For <canonical quantization of a real scalar field>, impose the equal-time <canonical commutation relations>
$$
[\phi(t,\mathbf x),\pi(t,\mathbf y)]=i\delta^3(\mathbf x-\mathbf y),\qquad
[\phi(t,\mathbf x),\phi(t,\mathbf y)]=[\pi(t,\mathbf x),\pi(t,\mathbf y)]=0.
$$
The <Heisenberg picture> equation $\dot O=i[H,O]$ then gives $\dot\phi=\pi$. For the other equation, the momentum term commutes with $\pi$; commuting the gradient term past $\pi$ and integrating the derivative of the <Dirac delta function> by parts gives $\dot\pi=\nabla^2\phi-m^2\phi$. Consequently
$$
\boxed{(\Box+m^2)\phi=0,\qquad\Box=\partial_t^2-\nabla^2.}
$$
The <Fourier transform> reduces this <Klein-Gordon equation> to oscillators of frequency $E_{\mathbf p}=\sqrt{\mathbf p^2+m^2}$. Hermiticity of the <real scalar field> pairs their positive- and negative-frequency coefficients, so its general mode expansion can be normalized as
$$
\phi(x)=\int\frac{d^3p}{(2\pi)^3\sqrt{2E_{\mathbf p}}}
[a(\mathbf p)e^{-ip\cdot x}+a^\dagger(\mathbf p)e^{ip\cdot x}],\qquad p^0=E_{\mathbf p}.
$$
To verify its normalization rather than guess the oscillator algebra, use <scalar field oscillator inversion> at $t=0$:
$$
a(\mathbf p)=\int d^3x\,e^{-i\mathbf p\cdot\mathbf x}
\left[\sqrt{E_{\mathbf p}/2}\,\phi(0,\mathbf x)+\frac{i}{\sqrt{2E_{\mathbf p}}}\pi(0,\mathbf x)\right].
$$
Its adjoint has the conjugate phase and the opposite sign in front of $i\pi$. The two mixed field-momentum <commutators> give
$$
[a(\mathbf p),a^\dagger(\mathbf q)]
=\frac{E_{\mathbf p}+E_{\mathbf q}}{2\sqrt{E_{\mathbf p}E_{\mathbf q}}}
\int d^3x\,e^{i(\mathbf q-\mathbf p)\cdot\mathbf x}
=\boxed{(2\pi)^3\delta^3(\mathbf p-\mathbf q)}.
$$
The <Fourier representation of the Dirac delta function> sets the energies equal in the last step. Similarly, $[a(\mathbf p),a(\mathbf q)]$ has a factor $E_{\mathbf q}-E_{\mathbf p}$ multiplying $\delta^3(\mathbf p+\mathbf q)$ and vanishes; the two-creator commutator vanishes as well.

Substituting the expansion in $H$ and integrating the spatial phases, the two-annihilator and two-creator coefficients cancel by $E_{\mathbf p}^2=\mathbf p^2+m^2$. The <free real scalar Hamiltonian in oscillator variables> is
$$
H=\frac12\int\frac{d^3p}{(2\pi)^3}E_{\mathbf p}
[a^\dagger(\mathbf p)a(\mathbf p)+a(\mathbf p)a^\dagger(\mathbf p)].
$$
It differs from $:H:=\int d^3p\,E_{\mathbf p}a^\dagger a/(2\pi)^3$ by its constant <vacuum energy>; <normal ordering> removes that constant. Using the oscillator commutators, for either Hamiltonian,
$$
\boxed{[H,a(\mathbf p)]=-E_{\mathbf p}a(\mathbf p),\qquad
[H,a^\dagger(\mathbf p)]=E_{\mathbf p}a^\dagger(\mathbf p).}
$$
Choose the <Fock vacuum> annihilated by every $a(\mathbf p)$. Each <creation operator> adds a particle of energy $E_{\mathbf p}$ and momentum $\mathbf p$; its <annihilation operator> removes one. Their mutual commutativity gives <bosonic Fock space> and symmetric many-particle states. Thus the scalar mass-shell frequency has a particle interpretation, rather than merely describing a classical oscillation.

For $t=x^0-y^0$ and $\mathbf r=\mathbf x-\mathbf y$, vacuum contraction and <time ordering> give the <Feynman propagator>
$$
D_F(t,\mathbf r)=\int\frac{d^3p}{(2\pi)^3}\frac{e^{i\mathbf p\cdot\mathbf r}}{2E_{\mathbf p}}
[\theta(t)e^{-iE_{\mathbf p}t}+\theta(-t)e^{iE_{\mathbf p}t}].
$$
In the <scalar Feynman propagator pole prescription>, the integral of $i e^{-ip^0t}/((p^0)^2-E_{\mathbf p}^2+i0)$ has poles at $E_{\mathbf p}-i0$ and $-E_{\mathbf p}+i0$. Close clockwise below for $t>0$ and counterclockwise above for $t<0$. The <residue theorem> gives exactly the corresponding term above, proving
$$
\boxed{D_F(x-y)=\int\frac{d^4p}{(2\pi)^4}\frac{i\,e^{-ip\cdot(x-y)}}{p^2-m^2+i0}.}
$$
As a direct check of the distributional source, $g_E(t)=e^{-iE|t|}/(2E)$ is continuous, but $g'_E(0^+)-g'_E(0^-)=-i$. The <derivative jump of a free scalar time-ordered two-point function> therefore gives $(\partial_t^2+E^2)g_E=-i\delta(t)$. Spatial Fourier inversion yields
$$
\boxed{(\Box+m^2)D_F(x-y)=-i\delta^4(x-y).}
$$