Derivatives of $ \frac{1}{r} $ and Dirac delta function 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?
 A: 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}
$$
A: 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}). $$
