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?

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}). $$
