Derivatives of $ \frac{1}{r} $ and Dirac delta function

Question

I am trying to understand the formula

\begin{equation} \nabla^2\left(\frac{1}{|{\bf r}-{\bf r}'|}\right) = - 4 \pi \delta(\bf{r}-\bf{r}'), \qquad\qquad {\rm (I)} \end{equation}

where ${\bf r}=(x,y,z)$. This is something heavily used in electrostatics and the steps to 'show' this is often the following:

The first derivative reads \begin{equation} \nabla \frac{1}{| {\bf r} - {\bf r}' |} = - \frac{ {\bf r} - {\bf r}'}{| {\bf r} - {\bf r}'|^3} \end{equation} And taking the second derivative gives zero, except for the singularity at ${\bf r} = {\bf r'}$. Then from the divergence theorem we have \begin{equation} \int dV \, \nabla^2 \frac{1}{| {\bf r} - {\bf r'}|} = \int dS \,\,{\bf n} \cdot ( \nabla\frac{1}{|{\bf r} - {\bf r}'|}) = -4 \pi \end{equation} where the integration is performed over a sphere centered at ${\bf r}'$.

Q1: Is there a more direct proof for equation (I)?

Then my main question is about the separate second-order differentials. For instance, we can obtain, by direct computation

\begin{equation} \partial_x^2 \, \frac{1}{| {\bf r} - {\bf r}'|} = \frac{ 3 (x-x')^2 }{| {\bf r} - {\bf r'} |^5} - \frac{1}{| {\bf r} - {\bf r}'|^3} \end{equation}

Q2: Should there be a $\delta$ function on the r.h.s of this equation?

Answer 1

A1. If you are not familiar with distribution theory, we might consider an alternative approach using the idea of approximate Dirac delta function. Indeed, define

$$ f_{\epsilon}(\mathbf{x}) = \frac{1}{\sqrt{\|\mathbf{x}\|^2+\epsilon^2}}=\frac{1}{\sqrt{x^2+y^2+z^2+\epsilon^2}}. $$

Then its Laplacian is

$$ \Delta f_{\epsilon}(\mathbf{x}) = -\frac{3\epsilon^2}{(x^2+y^2+z^2+\epsilon^2)^{5/2}}. $$

So, if $\varphi$ is any compactly supported smooth function on $\mathbb{R}^3$, then

\begin{align*} \int_{\mathbb{R}^3} \varphi(\mathbf{x}) \Delta f_{\epsilon}(\mathbf{x}) \, \mathrm{d}\mathbf{x} &= - \int_{\mathbb{R}^3} \varphi(\mathbf{x}) \frac{3\epsilon^2}{(x^2+y^2+z^2+\epsilon^2)^{5/2}} \, \mathrm{d}\mathbf{x} \\ &= - \int_{0}^{\infty} \int_{\mathbb{S}^2} \varphi(r\omega) \frac{3\epsilon^2 r^2}{(r^2+\epsilon^2)^{5/2}}\, \sigma(\mathrm{d}\omega)\mathrm{d}r \tag{$\mathbf{x}=r\omega$} \\ &= - \int_{0}^{\infty} \int_{\mathbb{S}^2} \varphi(\epsilon s \omega) \frac{3s^2}{(s^2+1)^{5/2}}\, \sigma(\mathrm{d}\omega)\mathrm{d}s, \tag{$r=\epsilon s$} \end{align*}

where $\mathbb{S}^2$ is the unit sphere centered at the origin and $\sigma$ is the surface measure of $\mathbb{S}^2$. (If this sounds a bit abstract, just think of the spherical coordinates change!) Now letting $\epsilon \to 0^+$, the dominated convergence theorem tells that switching the order of limit and integration is valid in this case, hence the integral converges to

\begin{align*} \lim_{\epsilon \to 0^+} \int_{\mathbb{R}^3} \varphi(\mathbf{x}) \Delta f_{\epsilon}(\mathbf{x}) \, \mathrm{d}\mathbf{x} = - \int_{0}^{\infty} \int_{\mathbb{S}^2} \varphi(0) \frac{3s^2}{(s^2+1)^{5/2}}\, \sigma(\mathrm{d}\omega)\mathrm{d}s = - 4\pi \varphi(0). \end{align*}

Here, we utilized $\int_{\mathbb{S}^2} \sigma(\mathrm{d}\omega) = 4\pi$ and $\int_{0}^{\infty} \frac{3s^2}{(s^2+1)^{5/2}} \, \mathrm{d}s = 1$.

---

A2. Still using the above setting, we have

\begin{align*} \partial^2_x f_{\epsilon}(\mathbf{x}) = \frac{2x^2-y^2-z^2-\epsilon^2}{(\|\mathbf{x}\|+\epsilon^2)^{5/2}} = \frac{2x^2-y^2-z^2}{(\|\mathbf{x}\|^2+\epsilon^2)^{5/2}} + \frac{1}{3}\Delta f_{\epsilon}(\mathbf{x}) \end{align*}

So it suffices to analyze the contribution of the first term in the last line. To this end, note that if $B_r$ denotes the ball of radius $r$ centered at the origin, then

$$ \int_{B_r} \frac{2x^2-y^2-z^2}{(\|\mathbf{x}\|^2+\epsilon^2)^{5/2}} \, \mathrm{d}\mathbf{x} = 0 $$

by the symmetry, and so, we may write

\begin{align*} &\int_{\mathbb{R}^3} \varphi(\mathbf{x}) \frac{2x^2-y^2-z^2}{(\|\mathbf{x}\|^2+\epsilon^2)^{5/2}} \, \mathrm{d}\mathbf{x} \\ &= \int_{\mathbb{R}^3} \left( \varphi(\mathbf{x}) - \varphi(0)\mathbf{1}_{B_r}(\mathbf{x}) \right) \frac{2x^2-y^2-z^2}{(\|\mathbf{x}\|^2+\epsilon^2)^{5/2}} \, \mathrm{d}\mathbf{x} \end{align*}

Introducing the regularizing term $- \varphi(0)\mathbf{1}_{B_r}(\mathbf{x})$ makes the integrand decay fast enough, i.e.,

$$ \left( \varphi(\mathbf{x}) - \varphi(0)\mathbf{1}_{B_r}(\mathbf{x}) \right) (2x^2-y^2-z^2) = \mathcal{O}(\|\mathbf{x}\|^3) $$

as $\|\mathbf{x}\| \to 0$, and so, we can utilize the dominated convergence theorem to conclude that

\begin{align*} &\lim_{\epsilon \to 0^+} \int_{\mathbb{R}^3} \varphi(\mathbf{x}) \frac{2x^2-y^2-z^2}{(\|\mathbf{x}\|^2+\epsilon^2)^{5/2}} \, \mathrm{d}\mathbf{x} \\ &= \int_{\mathbb{R}^3} \left( \varphi(\mathbf{x}) - \varphi(0)\mathbf{1}_{B_r}(\mathbf{x}) \right) \frac{2x^2-y^2-z^2}{\|\mathbf{x}\|^5} \, \mathrm{d}\mathbf{x}. \end{align*}

This defines a distribution on $\mathbb{R}^3$ which we may write

$$ \operatorname{p.v.}\left(\frac{2x^2-y^2-z^2}{\|\mathbf{x}\|^5}\right) $$

by analogy with the Cauchy principal value in the one-dimensional setting. In conclusion, we get

$$ \partial_x^2 \frac{1}{\|\mathbf{x}\|} = \operatorname{p.v.}\left(\frac{2x^2-y^2-z^2}{\|\mathbf{x}\|^5}\right) - \frac{4\pi}{3}\delta(\mathbf{x}). $$

Answer 2

So, first, in the classical sense, this function is not derivable at $x=0$. For $x≠0$, the classical calculus tells you that $$ \Delta (\tfrac{1}{|x|}) = 0 $$

To understand why the Dirac delta appears, one has to get a new notion of derivatives allowing to retrieve information about what is happening at $x=0$. This is called the theory of distribution. In the sense of distributions, $$ \Delta (\tfrac{1}{|x|}) = -4π\, \delta_0 $$ A way to prove it is to come back to the definition of derivative in the sense of distributions $$ \begin{align*} \langle \Delta (\tfrac{1}{|x|}),\varphi\rangle &= \int_{\mathbb{R}^3} \tfrac{\Delta \varphi(x)}{|x|}\,\mathrm{d}x = -\int_{\mathbb{R}^3} \tfrac{x\cdot\nabla \varphi(x)}{|x|^3}\,\mathrm{d}x \\ &= -\int_{|x|<1} \tfrac{x\cdot\nabla (\varphi(x)-\varphi(0))}{|x|^3}\,\mathrm{d}x -\int_{|x|>1} \tfrac{x\cdot\nabla \varphi(x)}{|x|^3}\,\mathrm{d}x \end{align*} $$ On the two last integrals, we can integrate by parts one more time, and find $-4π\, \varphi(0)$. Similar computations gives you the Hessian $\nabla^2$ (warning, here I am using the mathematical notation, $\nabla^2 = \nabla\nabla$ and $\Delta = \nabla\cdot\nabla$) $$ \nabla^2 (\tfrac{1}{|x|}) = \mathrm{pv.}\left(\frac{3\,x\otimes x - |x|^2\,\mathrm{Id}}{|x|^5}\right) - \frac{4π}{3}\, \delta_0 \,\mathrm{Id} $$