How to show the (Klein-Gordon) retarded propagator satisfies its equation of motion? The retarded propagator for a massless scalar field is
$$
G_R(t,\mathbf{x} ;t',\mathbf{x}' ) = \frac{ \Theta(t-t') \delta\big( - (t-t')^2 + |\mathbf{x} - \mathbf{x}'|^2 \big)}{2\pi} \tag{1}
$$
which is supposed to satisfy the equation
$$
( - \partial_t^2 + \nabla_{\mathbf{x}}^2 )G_R(t,\mathbf{x} ;t',\mathbf{x}' )  = - \delta^{(4)}(x-x')\ . \tag{2}
$$
How does one actual going about proving that the solution (1) actually satisfies (2)? I have tried inserting (1) into (2), and I get a lot of $\Theta$ and $\delta$ functions, as well as some derivatives $\delta'$ and $\delta''$, which I can not understand how this is supposed to result in the RHS of (2)?
EDIT: Being explicit, I find that (where $x^2 := -t^2 + |\mathbf{x}|^2$)
$$
\begin{align}
( - \partial_t^2 + \nabla_{\mathbf{x}}^2 )G_R(t,\mathbf{x} ;0,\mathbf{0} ) & = ( - \partial_t^2 + \nabla_{\mathbf{x}}^2 ) \big\{ \frac{\Theta(t) \delta(x^2 )}{2\pi} \big\} \\
& =  \frac{ \delta'(t) \delta(x^2) + 4 t \delta(t) \delta'(x^2) + 4 \Theta(t) \big[ 2 \delta'(x^2) + x^2 \delta''(x^2) \big] }{ 2\pi }
\end{align}
$$
This final expression does not seem to reduce to the claimed RHS.
 A: First, we need the expression for the propagator. Let's start with the expression for the Greens function
\begin{align}
\square G = \delta^4(x)\,.
\end{align}
We now perform a Fourier transform, so the right hand side of the equation becomes
\begin{align}
\int{\rm d}^4 x e^{{\rm i}(\omega t - \vec k\cdot\vec x)}\delta^4(x) = 1\,.
\end{align}
The left hand side of the equation is
\begin{align}
\square G(x) = (\partial_t^2 - \nabla^2)G(x)\to (-\omega^2 + (\vec k)^2)G(k)\,.
\end{align}
Thus,
\begin{align}
G(k) = \frac{-1}{\omega^2 - (\vec k)^2}\,.
\end{align}
The retarded propagator is obtained by the following prescription, $\omega\to \omega + {\rm i}\epsilon$. Thus, we can perform the inverse Fourier transform to obtain the position space representation of $G$:
\begin{align}
G_\text{ret}(r) = \int\frac{{\rm d}^4 k}{(2\pi)^4} \frac{e^{-{\rm i}(\omega t - \vec k\cdot\vec x)}}{(\omega + {\rm i}\epsilon)^2 - (\vec k)^2}\,.
\end{align}
Importantly, the poles of this function are in the lower half-plane. Therefore, the integral must close in the lower half plane in order to be non-zero (see the residue theorem for an explanation). When $t > 0$, the coefficient of $\omega$ is negative, so we close the integral in the lower half plane. Similarly, when $t < 0$, we close the integral in the upper half plane. The latter case does not contribute, which we express with a Heaviside function
\begin{align}
G_\text{ret} &= \theta(t)2\pi{\rm i}\int\frac{{\rm d}^3 k}{(2\pi)^4}\left( \frac{e^{-{\rm i}(-k t - \vec k\cdot\vec x)}}{2k} +  \frac{e^{-{\rm i}(k t - \vec k\cdot\vec x)}}{2k}\right)\,,\\
&=\theta(t)\int\frac{{\rm d}^3 k}{(2\pi)^3}e^{{\rm i} k x}\frac{\sin kt}{k}
\end{align}
We now take the angular integral
\begin{align}
G_\text{ret} &= \theta(t)\int \frac{4\pi k^2{\rm d} k}{(2\pi)^3}j_0(k r)\frac{\sin kt}{k}\,,\\
&=\frac{\theta(t)}{r}\int \frac{4\pi {\rm d} k}{(2\pi)^3}\sin(k r)\sin(k t)\,.  
\end{align}
We now observe that this is just an expression of the $\delta$-function
\begin{align}
G_\text{ret} &=\frac{\theta(t)}{4\pi r}\delta(t - r)\,.  
\end{align}
Now that we have the correct expression for the Greens function, let's apply the D'Alembertian to it and see that this is indeed the Greens function.
\begin{align}
\square\frac{\theta(t)}{4\pi r}\delta(t - r) &= \left(\partial_t^2 - \nabla^2\right)\frac{\theta(t)}{4\pi r}\delta(t - r)\,,\\
&=\frac{1}{4\pi r}\partial_t^2\theta(t)\delta(t - r) - \theta(t)\nabla^2\frac{1}{4\pi r}\delta(t - r)\,,\\
&=\frac{1}{4\pi r}\partial_t(\delta(t - r)\delta(t) + \theta(t)\delta'(t - r)) - \theta(t)\left(\delta(t - r)\nabla^2\frac{1}{4\pi r} + \frac{1}{4\pi r}\nabla^2\delta(t - r) + 2\partial_r\delta(t - r)\partial_r\frac{1}{4\pi r}\right)\,,\\
&=\frac{1}{4\pi r}(2\delta'(t - r)\delta(t)+\delta(t - r)\delta'(t) +\theta(t)\delta''(t - r)) - \theta(t)\left(-\delta(t - r)\delta(r) + \frac{1}{4\pi r}\delta''(t - r)- \frac{2}{4\pi r^2}\delta'(t - r) + 2\delta'(t - r)\frac{1}{4\pi r^2}\right)\,,\\
&=\frac{1}{4\pi r}(2\delta'(t - r)\delta(t)+\delta(t - r)\delta'(t) +\theta(t)\delta''(t - r)) - \theta(t)\left(-\delta(t - r)\delta(r) + \frac{1}{4\pi r}\delta''(t - r)\right)\,,\\
&=\frac{1}{4\pi r}(2\delta'(t - r)\delta(t)+\delta(t - r)\delta'(t)) - \theta(t)\left(-\delta(t - r)\delta(r)\right)\,,\\
&=\theta(t)\delta(t - r)\delta(r)\,.
\end{align}
Note, the identity that $\nabla^2 1/r = -4\pi\delta(r)$ is well explained elsewhere (I'll provide a link later).
