D'Alembertian of a delta-function of a space-time interval (i.e. on the light-cone) How one differentiates a delta-function of a space-time interval? Namely,
$$[\partial_t^2 - \partial_x^2 - \partial_y^2 - \partial_z^2] \, \delta(t^2-x^2-y^2-z^2) \, .$$
Somewhere I saw that the result was $$4\pi\delta^{4}(x).$$ However, have no idea on how to obtain it.
 A: Using the theory of generalized functions, let us regularize the Dirac delta distribution$^1$
$$ \delta(x^2)~=~\lim_{\varepsilon\searrow 0^+}\delta_{\varepsilon}(x^2), \qquad x^2~:=~x^{\mu}\eta_{\mu\nu}x^{\nu}, \tag{1}$$
via a smooth $C^{\infty}$-function$^2$
$$ \delta_{\varepsilon}(x^2)
~=~ \frac{1}{\pi}\frac{\varepsilon}{(x^2)^2+\varepsilon^2} 
~=~-\frac{1}{\pi}{\rm Im} \frac{1}{x^2+i\varepsilon}, \qquad \varepsilon~>~0. \tag{2}$$
Eqs. (1)-(2) are a well-known representation of the Dirac delta distribution.
It is now mathematically well-defined to apply the d'Alembert operator on the smooth function (2). We calculate:
$$  \partial_{\mu}\delta_{\varepsilon}(x^2)~=~\frac{1}{\pi}{\rm Im} \frac{2x_{\mu}}{(x^2+i\varepsilon)^2},\tag{3}$$
$$  \partial_{\mu}\partial_{\nu}\delta_{\varepsilon}(x^2)
~=~\frac{1}{\pi}{\rm Im} \left(\frac{2\eta_{\mu\nu}}{(x^2+i\varepsilon)^2}-\frac{8x_{\mu}x_{\nu}}{(x^2+i\varepsilon)^3}\right),\tag{4}$$
$$\begin{align}  \Box \delta_{\varepsilon}(x^2)~=~&\frac{1}{\pi}{\rm Im} \left(\frac{8}{(x^2+i\varepsilon)^2}-\frac{8x^2}{(x^2+i\varepsilon)^3}\right) \cr
~=~&\frac{1}{\pi}{\rm Im}\frac{8i\varepsilon}{(x^2+i\varepsilon)^3}\cr
~\longrightarrow~& -4\pi  \delta^4(x) 
\qquad \mathrm{for} \qquad \varepsilon ~\searrow~ 0^+.
\end{align}\tag{5}$$
Eq. (5) implies OP's sought-for proposition.

Proposition.
$$ \Box \delta (x^2)~=~-4\pi  \delta^4(x).\tag{6}$$

Sketched proof of the last line in eq. (5). Consider a test function $f\in C^{\infty}_c(\mathbb{R}^4)$, i.e., an infinitely often differentiable function $f$ with compact support. Then
$$\begin{align}
\int_{\mathbb{R}^4}\! \mathrm{d}^4x & ~ f(x)~ \Box \delta_{\varepsilon}(x^2) \cr
~\stackrel{x=\sqrt{\varepsilon} y}{=}~& 
\int_{\mathbb{R}^4}\!\frac{\mathrm{d}^4y}{\pi} ~ f(\sqrt{\varepsilon} y)~{\rm Im}\frac{8i}{(y^2+i)^3}\cr  
~\longrightarrow~&  f(0) \int_{\mathbb{R}^4}\!\frac{\mathrm{d}^4y}{\pi} ~ {\rm Im}\frac{8i}{(y^2+i)^3}\cr
~\stackrel{y=(t,\vec{r})}{=}~&~ 4 f(0) {\rm Im} \int_{\mathbb{R}_+}\!r^2\mathrm{d}r \int_{\mathbb{R}}\!\mathrm{d}t \frac{8i}{(r^2-t^2+i)^3}\cr
~\stackrel{a:=\sqrt{r^2+i}}{=}&~ 4 f(0) {\rm Im} \int_{\mathbb{R}_+}\!r^2\mathrm{d}r \int_{\mathbb{R}}\!\mathrm{d}t \left(\frac{2i}{t^2-a^2}\right)^3\cr
~=~~~& 4 f(0) {\rm Im} \int_{\mathbb{R}_+}\!r^2\mathrm{d}r \frac{-i}{a^3}\int_{\mathbb{R}}\!\mathrm{d}t \left(\frac{1}{t-a}-\frac{1}{t+a}\right)^3\cr
~\stackrel{\text{residue thm.}}{=}&~ 4 f(0) {\rm Im} \int_{\mathbb{R}_+}\!r^2\mathrm{d}r \frac{-i}{a^3} 2\pi i \cr
&\qquad \left(0+\left.\frac{3}{(t+a)^2}\right|_{t=a}+\left.\frac{3}{(t-a)^2}\right|_{t=-a}+0\right)\cr
~=~~~& 12\pi f(0) {\rm Im} \int_{\mathbb{R}_+}\!r^2\mathrm{d}r \frac{1}{a^5} \cr
~=~~~& 12\pi f(0) {\rm Im} \int_{\mathbb{R}_+}\!\frac{r^2\mathrm{d}r}{(r^2+i)^{5/2}} \cr
~=~~~& 12\pi f(0) {\rm Im} \left[\frac{r^3}{3i(r^2+i)^{3/2}}\right]^{r=\infty}_{r=0} \cr
~=~~~& -4\pi f(0)   \qquad \mathrm{for} \qquad \varepsilon ~\searrow~ 0^+,
\end{align}\tag{7}$$
because of, e.g., Lebesgue's dominated convergence theorem. Here we have implicitly assumed that $a:=\sqrt{r^2+i}$ is the square root branch in the upper complex plane. $\Box$
--
$^1$ In this answer we work in units where the speed-of-light $c=1$ is one, and we use the Minkowski sign convention $(−,+,+,+)$.
$^2$ The $i\varepsilon$-prescription has an interpretation in terms of Wick rotation, cf. e.g. my Math.SE answer here.
A: The action of any distributional "function" such as a delta function is really only determined when you integrate it with a test function;  and two distributions are equal if and only if their actions agree on all test functions.  So let's see what this distribution does when we integrate it with a smooth test function $f(t, \vec{r})$.
We start with the integral
$$
I[f] = \iiiint d^4 x \, \left[ f(t, \vec{r}) \Box  \delta(\Delta s^2) \right]
$$
where $\Box$ stands for the wave operator under the metric convention $(+ - - \, -)$.  Integrating by parts twice and assuming that $f$ has compact support (as all good test functions should), this becomes
$$
I[f] = \iiiint dt d^3\vec{r} \, \left[ (\Box f)  \delta(t^2 - x^2 - y^2 - z^2) \right]
$$
We can then use the rules concerning composition of the delta-function with a function to rewrite this as
$$
I[f] = \iiiint dt d^3\vec{r} \, (\Box f) \left[ \frac{\delta(t - |\vec{r}|)}{2 |t|} + \frac{\delta(t + |\vec{r}|)}{2 |t|} \right].
$$
Finally, we can integrate over $t$ to obtain
$$
I[f] = \iiint d^3 \vec{r} \, \frac{ \Box f(|\vec{r}|, \vec{r}) + \Box f(-|\vec{r}|, \vec{r})}{2 |\vec{r}|}.
$$
The notation here for $\Box f(\pm|\vec{r}|, \vec{r})$ is awkward, but it means "$\Box f$ evaluated at $t = \pm|\vec{r}|$."
This action will differ from the action of the distribution $4 \pi \delta^4(x^\mu)$.  In particular, we have
$$
J[f] = \iiiint d^4 x \, f(t,\vec{r}) \left[ 4 \pi \delta^4(x^\mu) \right] = 4 \pi f(0, \vec{0}).
$$
$I$ and $J$ are only equal if their results agree on all test functions $f$.  But if we choose a test function that vanishes at the origin and for which $\Box f$ is positive on some portion of the light cone, then $I$ and $J$ will yield a positive result and a zero result respectively.  Thus, the distributions are not equal.
It is possible that the actions of these two distributions will be equal on test functions $f$ that have special properties (for example, if $f$ is the Green's function for the wave operator or something like that.)  This may be all that is necessary in a particular context.  But the two distributions are not equal in a strict sense.
