Solution (source code)

= Solution

Use units $\hbar=c=1$ and the <Minkowski metric> $g=\operatorname{diag}(1,-1,-1,-1)$. The <canonical momentum> of the <real scalar field> is $\pi=\partial\mathcal L/\partial\dot\phi=\dot\phi$. <Canonical quantization> replaces the classical fields by <operator-valued distributions> and imposes 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.
$$
These relations are imposed on an initial slice, with smeared fields understood. In the <Heisenberg picture>, simultaneous <unitary time evolution> of both operators preserves their <commutators>. Since the momentum is $\dot\phi$, \b[the requested equal-time commutator is] $[\phi(t,\mathbf x),\dot\phi(t,\mathbf y)]=i\delta^{(3)}(\mathbf x-\mathbf y)$.

The <Legendre transform in mechanics>, applied pointwise, gives the <canonical Hamiltonian density of a real scalar field>. Thus
$$
\boxed{H=\frac12\int d^3x\,[\pi^2+|\nabla\phi|^2+m^2\phi^2].}
$$
To obtain its dynamics rather than assume them, use the <Heisenberg equation of motion> $\dot O=i[H,O]$. The momentum-square term gives $[H,\phi]=-i\pi$. The field-square term gives $[H,\pi]=im^2\phi$, while differentiating the <Dirac delta function> in the gradient term and integrating by parts gives $[H,\pi]=-i\nabla^2\phi+im^2\phi$. Consequently
$$
\dot\phi=\pi,\qquad \dot\pi=\nabla^2\phi-m^2\phi,\qquad \boxed{(\Box+m^2)\phi=0.}
$$
Boundary terms vanish for suitably decaying smeared fields. This is the <Klein-Gordon equation> as an operator identity.

A spatial <Fourier transform> reduces the <Klein-Gordon equation> to $\ddot\phi_{\mathbf p}+E_p^2\phi_{\mathbf p}=0$, where $E_p=\sqrt{\mathbf p^2+m^2}$. Its two frequency branches and the <Hermitian adjoint> condition give
$$
\phi(x)=\int d\Pi_p\,[a(p)e^{-ip\cdot x}+a^\dagger(p)e^{ip\cdot x}],\qquad d\Pi_p=\frac{d^3p}{(2\pi)^3 2E_p},\quad p^0=E_p.
$$
The <invariant-normalized scalar mode extraction> can be written particularly simply as
$$
a(p)=\int d^3x\,e^{ip\cdot x}[E_p\phi(x)+i\pi(x)],\qquad a^\dagger(p)=\int d^3x\,e^{-ip\cdot x}[E_p\phi(x)-i\pi(x)].
$$
Substituting the <mode expansion of a free field> proves these inversion formulas: the opposite-frequency coefficient is proportional to $E_p-E_q$ and vanishes after spatial <Fourier inversion>. Alternatively they are the corresponding <Klein-Gordon inner products>. Taking their <commutators> using the equal-time fields gives $(E_p+E_q)(2\pi)^3\delta^{(3)}(\mathbf p-\mathbf q)$ in the mixed case. In the two equal-type cases the delta function enforces $\mathbf q=-\mathbf p$, multiplying $E_p-E_q=0$. Therefore
$$
\boxed{[a(p),a^\dagger(q)]=(2\pi)^3 2E_p\delta^{(3)}(\mathbf p-\mathbf q),\qquad [a(p),a(q)]=[a^\dagger(p),a^\dagger(q)]=0.}
$$
Conversely these mode <commutators> reproduce $[\phi,\dot\phi]=i\delta^{(3)}$: the two frequency contributions each supply half the spatial <Dirac delta function>.

For the <covariant oscillator Hamiltonian of a real scalar field>, insert the <mode expansion of a free field> into $H$. Spatial integration sets the two momenta equal in the mixed terms and opposite in the two-creator or two-annihilator terms. The latter have coefficient $-E_p^2+\mathbf p^2+m^2=0$. The mixed terms give
$$
H=\frac12\int d\Pi_p\,E_p[a^\dagger(p)a(p)+a(p)a^\dagger(p)]
=\int d\Pi_p\,E_p a^\dagger(p)a(p)+E_0,
\qquad E_0=\frac12\int d^3p\,E_p\delta^{(3)}(0).
$$
The constant $E_0$ is the divergent <zero-point energy>; a finite box and a momentum cutoff make its meaning explicit. <Normal ordering> sets the free <vacuum energy> to zero and leaves all <commutators> unchanged. Directly applying the mode <commutators> to the mixed expression yields
$$
\boxed{[H,a(p)]=-E_pa(p),\qquad [H,a^\dagger(p)]=E_pa^\dagger(p).}
$$
Choose the <Fock vacuum> by $a(p)|0\rangle=0$. Each <creation operator> adds a <boson> of energy $E_p$ and <four-momentum> $p$, while each <annihilation operator> removes one. Products of commuting <creation operators> give <bosonic statistics from commuting creation operators> and unrestricted <occupation numbers>. The <real scalar field> has one spin-zero species, identical to its <antiparticle>. Its <one-particle states> obey $\langle q|p\rangle=(2\pi)^3 2E_p\delta^{(3)}(\mathbf p-\mathbf q)$, precisely the <relativistic normalization of a one-particle state>.

For the <Feynman propagator>, put $z=x-y$ and contract the two fields in the <Fock vacuum>. Only $aa^\dagger$ survives, giving the <Wightman function>
$$
W(z)=\langle0|\phi(x)\phi(y)|0\rangle=\int d\Pi_p\,e^{-ip\cdot z}.
$$
The stipulated <time-ordered product> convention is $i\Delta_F(z)=\theta(z^0)W(z)+\theta(-z^0)W(-z)$. In a <contour integral> over $p^0$, the <Feynman i-epsilon prescription> places the positive-energy pole just below the real axis and the negative-energy pole just above it. For $z^0>0$ close below, clockwise, obtaining $-i e^{-iE_pz^0}/(2E_p)$; for $z^0<0$ close above, obtaining $-i e^{iE_pz^0}/(2E_p)$. Changing $\mathbf p\mapsto-\mathbf p$ in the second contribution therefore proves
$$
\boxed{\Delta_F(z)=\lim_{\epsilon\downarrow0}\int\frac{d^4p}{(2\pi)^4}\frac{e^{-ip\cdot z}}{p^2-m^2+i\epsilon}.}
$$
Acting with $\Box+m^2$ multiplies the integrand by $-(p^2-m^2)$. The remaining factor tends to $-1$ as a <distribution>, so <Fourier inversion> gives
$$
\boxed{(\Box+m^2)\Delta_F(z)=-\delta^{(4)}(z).}
$$
The sign also follows from the <derivative jump of a free scalar time-ordered two-point function>: the first time derivative of $i\Delta_F$ jumps by $\langle[\dot\phi,\phi]\rangle=-i\delta^{(3)}$, so that of $\Delta_F$ jumps by $-\delta^{(3)}$.