= Solution
Write $t=x^0-y^0$ and $\mathbf r=\mathbf x-\mathbf y$. By <time ordering> for a bosonic field,
$$
\Delta_F(x-y)=\theta(t)\int\frac{d^3p}{(2\pi)^3\,2E_p}e^{-iE_pt+i\mathbf p\cdot\mathbf r}
+\theta(-t)\int\frac{d^3p}{(2\pi)^3\,2E_p}e^{iE_pt-i\mathbf p\cdot\mathbf r}.
$$
Reverse $\mathbf p$ in the second integral. At fixed spatial momentum consider
$$
I(t)=\int_{-\infty}^{\infty}\frac{dp^0}{2\pi}\frac{i e^{-ip^0t}}{(p^0)^2-E_p^2+i0}.
$$
The <Feynman i-epsilon prescription> places the poles at $+E_p-i0$ and $-E_p+i0$. For $t>0$, close in the lower half-plane, clockwise. The positive-energy residue and the orientation give $e^{-iE_pt}/(2E_p)$. For $t<0$, close counterclockwise in the upper half-plane; the negative pole has derivative denominator $-2E_p$, and its residue gives $e^{iE_pt}/(2E_p)$. Thus $I(t)=e^{-iE_p|t|}/(2E_p)$, reproducing the time-ordered mode expression. We obtain
$$
\boxed{\Delta_F(x-y)=\int\frac{d^4p}{(2\pi)^4}\frac{i e^{-ip\cdot(x-y)}}{p^2-\mu^2+i0}.}
$$
This is a distributional boundary value, not an absolutely convergent Fourier integral. As an independent normalization check, $I'(0^+)-I'(0^-)=-i$, so $(\Box+\mu^2)\Delta_F=-i\delta^{(4)}$; the numerator is $i$ with the question's definition of the propagator.
Back to article page
