D'Alembertian of a delta-function of a space-time interval (i.e. on the light-cone)

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^3vec{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^3vec{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.

Using the theory of generalized functions, let us regularize the Dirac delta distribution$^1$ $$ delta(x^2)~=~lim_{varepsilonsearrow 0^+}delta_{varepsilon}(x^2), qquad x^2~:=~x^{mu}eta_{munu}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+ivarepsilon}, 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+ivarepsilon)^2},tag{3}$$ $$ partial_{mu}partial_{nu}delta_{varepsilon}(x^2) ~=~frac{1}{pi}{rm Im} left(frac{2eta_{munu}}{(x^2+ivarepsilon)^2}-frac{8x_{mu}x_{nu}}{(x^2+ivarepsilon)^3}right),tag{4}$$ $$begin{align} Box delta_{varepsilon}(x^2)~=~&frac{1}{pi}{rm Im} left(frac{8}{(x^2+ivarepsilon)^2}-frac{8x^2}{(x^2+ivarepsilon)^3}right) cr ~=~&frac{1}{pi}{rm Im}frac{8ivarepsilon}{(x^2+ivarepsilon)^3}cr ~longrightarrow~& -4pi 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)~=~-4pi delta^4(x).tag{6}$$

Sketched proof of the last line in eq. (5). Consider a test function $fin 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^2mathrm{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^2mathrm{d}r int_{mathbb{R}}!mathrm{d}t left(frac{2i}{t^2-a^2}right)^3cr ~=~~~& 4 f(0) {rm Im} int_{mathbb{R}_+}!r^2mathrm{d}r frac{-i}{a^3}int_{mathbb{R}}!mathrm{d}t left(frac{1}{t-a}-frac{1}{t+a}right)^3cr ~stackrel{text{residue thm.}}{=}&~ 4 f(0) {rm Im} int_{mathbb{R}_+}!r^2mathrm{d}r frac{-i}{a^3} 2pi i cr &qquad left(0+left.frac{3}{(t+a)^2}right|_{t=a}+left.frac{3}{(t-a)^2}right|_{t=-a}+0right)cr ~=~~~& 12pi f(0) {rm Im} int_{mathbb{R}_+}!r^2mathrm{d}r frac{1}{a^5} cr ~=~~~& 12pi f(0) {rm Im} int_{mathbb{R}_+}!frac{r^2mathrm{d}r}{(r^2+i)^{5/2}} cr ~=~~~& 12pi f(0) {rm Im} left[frac{r^3}{3i(r^2+i)^{3/2}}right]^{r=infty}_{r=0} cr ~=~~~& -4pi 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 $ivarepsilon$-prescription has an interpretation in terms of Wick rotation, cf. e.g. my Math.SE answer here.