Let $\psi$ be the distribution $\psi(r)=\frac1r$. Then, for $\phi \in C_C^\infty$, we have
$$\begin{align}
\langle \partial_{ij}\psi, \phi\rangle &=\langle \psi, \partial_{ij}\phi\rangle\\\\
&=\int_{\mathbb{R}^3}\frac1r \frac{\partial^2\phi(\vec r)}{\partial x_i\partial x_j}\,d^3\vec r\\\\
&=\underbrace{\int_{\mathbb{R}^3}\frac{\partial}{\partial x_i}\left(\frac1r \frac{\partial\phi(\vec r)}{\partial x_j}\right)\,d^3\vec r}_{=0\,\,\text{since}\,\,\phi\in C_C^\infty}-\int_{\mathbb{R}^3}\frac{\partial}{\partial x_i}\left(\frac1r \right)\left(\frac{\partial\phi(\vec r)}{\partial x_j}\right)\,d^3\vec r\\\\
&=-\lim_{\varepsilon\to 0^+}\int_{\mathbb{R}^3\setminus B(0,\varepsilon)}\frac{\partial}{\partial x_i}\left(\frac1r \right)\left(\frac{\partial\phi(\vec r)}{\partial x_j}\right)\,d^3 \vec r\tag1
\end{align}$$
where $B(\vec r_c,R)$ is a sphere of radius $R$ centered at $\vec r_C$. Now, integrating by parts again, we find that
$$\begin{align}
\langle \partial_{ij}\psi, \phi\rangle&=-\lim_{\varepsilon\to 0^+}\int_{\mathbb{R}^3\setminus B(0,\varepsilon)}\frac{\partial }{\partial x_j}\left(\phi(\vec r)\frac{\partial}{\partial x_i}\left(\frac1r \right)\right)\,d^3\vec r\\\\
&+\lim_{\varepsilon\to 0^+}\int_{\mathbb{R}^3\setminus B(0,\varepsilon)}\phi(\vec r)\frac{\partial^2 }{\partial x_i\partial x_j}\left(\frac1r\right)\,d^3\vec r\\\\
&=\lim_{\varepsilon\to0^+}\oint_{\partial B(0,\varepsilon)}(\hat n\cdot \hat x_j)\phi(\vec r)\frac{\partial }{\partial x_i}\left(\frac1r\right)\,dS\\\\
&+\lim_{\varepsilon\to 0^+}\int_{\mathbb{R}^3\setminus B(0,\varepsilon)}\phi(\vec r)\frac{\partial^2 }{\partial x_i\partial x_j}\left(\frac1r\right)\,d^3\vec r\\\\
&=-\lim_{\varepsilon\to0^+}\int_0^{2\pi}\int_0^\phi (\hat r\cdot \hat x_j)\phi(\vec r)\left.\left(\frac{x_i}{r^3}\right)\right|_{r=\varepsilon}\,\varepsilon^2\,\sin(\theta)\,d\theta\,d\phi\\\\
&+\lim_{\varepsilon\to 0^+}\int_{\mathbb{R}^3\setminus B(0,\varepsilon)}\phi(\vec r)\frac{\partial^2 }{\partial x_i\partial x_j}\left(\frac1r\right)\,d^3\vec r\\\\
&=-\frac{4\pi}{3}\delta_{ij}\phi(0)+\lim_{\varepsilon\to 0^+}\int_{\mathbb{R}^3\setminus B(0,\varepsilon)}\phi(\vec r)\frac{\partial^2 }{\partial x_i\partial x_j}\left(\frac1r\right)\,d^3\vec r\\\\
\end{align}$$
Therefore, in distribution we see that
$$\frac{\partial^2}{\partial x_i\partial x_j}\left(\frac1r\right)=-\frac{4\pi}{3}\delta_{ij}\delta(\vec r)+\text{PV}\left(\frac{3x_ix_j-\delta_{ij}r^2}{r^5}\right)$$
When $i=j$, we see that in distribution
$$\frac{\partial^2}{\partial x_i^2}\left(\frac1r\right)=-\frac{4\pi}{3}\delta(\vec r)+\text{PV}\left(\frac{3x_i^2-r^2}{r^5}\right)$$
whereupon summing over $i$ yields the familar result
$$\nabla^2 \frac1r =-4\pi \delta(\vec r)$$
$$\frac{\partial^2}{ \partial x^2} }\frac{1}{\sqrt{x^2+y^2+z^2}}} =\frac{3 x^2}{r^5}-\frac{1}{r^3}$$
there will be something wrong.
– Roland F Mar 11 '24 at 16:01