Laplacians and Dirac delta functions

The gradient of $\frac1r$ (noting that $r=\sqrt{x^2+y^2+z^2}$) is

$$ \nabla \frac1r = -\frac{\mathbf{r}}{r^3} $$ when $r\neq 0$, where $\mathbf{r}=x\mathbf{i}+y\mathbf{j}+z\mathbf{k}$. Now, the divergence of this is

$$ \nabla\cdot \left(-\frac{\mathbf{r}}{r^3}\right) = 0 $$ when $r\neq 0$. Therefore, for all points for which $r\neq 0$,

$$ \nabla^2\frac1r = 0 $$ However, if we integrate this function over a sphere, $S$, of radius $a$, then, applying Gauss's Theorem, we get

$$ \iiint_S \nabla^2\frac1rdV = \iint_{\Delta S} -\frac{\mathbf{r}}{r^3}.d\mathbf{S} $$ where $\Delta S$ is the surface of the sphere, and is outward-facing. Now, $d\mathbf{S}=\mathbf{\hat r}dA$, where $dA=r^2\sin\theta d\phi d\theta$. Therefore, we may write our surface integral as $$\begin{align} \iint_{\Delta S} -\frac{\mathbf{r}}{r^3}.d\mathbf{S}&=-\int_0^\pi\int_0^{2\pi}\frac{r}{r^3}r^2\sin\theta d\phi d\theta\\ &=-\int_0^\pi\sin\theta d\theta\int_0^{2\pi}d\phi\\ &= -2\cdot 2\pi = -4\pi \end{align}$$ Therefore, the value of the laplacian is zero everywhere except zero, and the integral over any volume containing the origin is equal to $-4\pi$. Therefore, the laplacian is equal to $-4\pi \delta(\mathbf{r})$.

EDIT: The general case is then obtained by replacing $r=|\mathbf{r}|$ with $s=|\mathbf{r}-\mathbf{r_0}|$, in which case the function shifts to $-4\pi \delta(\mathbf{r}-\mathbf{r_0})$

I'm new around here (so suggestions about posting are welcome!) and want to give my contribution to this question, even though a bit old. I feel I need to because using the divergence theorem in this context is not quite rigorous. Strictly speaking $1/r$ is not even differentiable at the origin. So here's a proof using limits of distributions.

Let $\mathbf{x}\in\mathbb{R}^{3}$ and $r=|\mathbf{x}|=\sqrt{x^2+y^2+z^2}$. It is evident from direct calculation that $\nabla^{2}\left(\frac{1}{r}\right)=0$ everywhere except in $\mathbf{x}=0$, where it is in fact not defined. Thus, the integral of $\nabla^{2}\left(\frac{1}{r}\right)$ over any volume non containing the origin is zero.

So let $\eta>0$ and $r_{\eta}=\sqrt{x^2+y^2+z^2+\eta^2}$. Obviously $\lim\limits_{\eta\rightarrow 0}r_{\eta}=r$. Direct calculation brings \begin{equation} \nabla^{2}\left(\frac{1}{r_{\eta}}\right) = \frac{-3\eta^{2}}{r_{\eta}^5} \end{equation} Now let us consider the distribution represented by $\nabla^{2}\left(\frac{1}{r_{\eta}}\right)$ and let $\rho$ be a test function (for example in the Schwartz space). I use Dirac's bra-ket notation to express the action of a distribution over a test function. Let $S^{2}$ be the unit sphere. Thus we calculate \begin{align} \lim\limits_{\eta\rightarrow 0}\left.\left\langle \nabla^{2}\left(\frac{1}{r_{\eta}}\right)\right|\rho\right\rangle &= \lim\limits_{\eta\rightarrow 0}\iiint\limits_{\mathbb{R}^3}\mathrm{d}^{3}x\, \nabla^{2}\left(\frac{1}{r_{\eta}}\right)\rho(\mathbf{x})\\ &=\lim\limits_{\eta\rightarrow 0}\left\{\iiint\limits_{\mathbb{R}^3\setminus S^{2}}\mathrm{d}^{3}x\, \frac{-3\eta^{2}}{r_{\eta}^5}\rho(\mathbf{x})+ \iiint\limits_{S^{2}}\mathrm{d}^{3}x\, \frac{-3\eta^{2}}{r_{\eta}^5}\rho(\mathbf{x})\right\}\\ &=\lim\limits_{\eta\rightarrow 0}\iiint\limits_{S^{2}}\mathrm{d}^{3}x\, \frac{-3\eta^{2}}{r_{\eta}^5}\rho(\mathbf{x}) \end{align} Where the limit of the first of the integrals in the curly braces is zero (easy to show, referring to the laplacian of $1/r$, no need for $\eta$ in sets not containing the origin). Now Taylor expand $\rho$ at $\mathbf{x}=0$ and integrate using spherical polar coordinates: \begin{equation} \lim\limits_{\eta\rightarrow 0}\int_{0}^{\pi}\mathrm{d}\theta\,\sin\theta\int_{0}^{2\pi}\mathrm{d}\varphi\,\int_{0}^{1}\mathrm{d}t\,\frac{-3\eta^{2}t^{2}}{(t^{2}+\eta^{2})^{5/2}}(\rho(0)+O(t^{2})) \end{equation} Integrating you will get that all the terms contained in $O(t^2)$ vanish as $\eta\rightarrow 0$, while the term with $\rho(0)$ remains. In fact you get \begin{align} \lim\limits_{\eta\rightarrow 0}\frac{-4\pi\rho(0)}{\sqrt{1+\eta^{2}}}&=-4\pi\rho(0)\\ &=\langle -4\pi\delta_{0}^{(3)}|\rho\rangle \end{align}

From this argument one defines the limit distribution \begin{equation} \nabla^{2}\left(\frac{1}{r}\right):=\lim\limits_{\eta\rightarrow 0}\nabla^{2}\left(\frac{1}{r_{\eta}}\right)=-4\pi\delta_{0}^{(3)} \end{equation} The generalization to $r=|\mathbf{x}-\mathbf{x}_{0}|$ is obvious.

Here is a back-alley derivation using Fourier transform properties. Take the Fourier transform of $\frac{1}{4 \pi r}$ to get $\frac{1}{k^2}$. Therefore, $$ \frac{1}{4\pi r} = \int \frac{d^3k}{(2\pi)^3} \frac{e^{-ik\cdot r}}{k^2} $$ Now, apply $-\nabla^2$ to get $$ -\nabla^2 \frac{1}{4\pi r} = \int \frac{d^3k}{(2\pi)^3} e^{-ik\cdot r} $$ This is the delta function, but we can explicitly put in our test function: \begin{align*} \int d^3r \left(-\nabla^2 \frac{1}{4 \pi r}\right) f(r) &= \int \frac{d^3k}{(2\pi)^3} \int d^3r e^{-ik\cdot r} f(r) \\ &= \int \frac{d^3k}{(2\pi)^3} \tilde{f}(k) \\ & = f(r=0). \end{align*}