Where does the relation $\nabla^2(1/r)=-4\pi\delta^3({\bf r})$ between Laplacian and Dirac delta function come from? It is often quoted in physics textbooks for finding the electric potential using Green's function that 
$$\nabla ^2 \left(\frac{1}{r}\right)=-4\pi\delta^3({\bf r}),$$ 
or more generally 
$$\nabla ^2 \left(\frac{1}{|| \vec x - \vec x'||}\right)=-4\pi\delta^3(\vec x - \vec x'),$$
where $\delta^3$ is the 3-dimensional Dirac delta distribution.
However I don't understand how/where this comes from. Would anyone mind explaining?
 A: 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*}
A: 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})$
A: 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.
A: We can use the simplest method to display the results, as shown below : -
$$ \nabla ^2 \left(\frac{1}{r}\right) = \nabla \cdot \nabla \left( \frac 1 r \right) = \nabla \cdot \frac {-1 \mathbf {e_r}} {r^2} $$ 
Suppose there is a sphere centered on the origin, then the total flux on the surface of the sphere is : -
$$ \text {Total flux} = 4 \pi r^2 \frac {-1} {r^2} = -4 \pi $$
Suppose the volume of the sphere be $ \mathbf {v(r)}$, so by the definition, the divergence is : -
$$ \lim_{\text {volume} \to zero} \frac {\text {Total Flux}} {\text {Volume}} = \lim_{\text {v(r)} \to 0} \left(\frac {-4 \pi} {v(r)}\right)  $$
So obviously,
$$ \lim_{\text {r} \to 0} \left[ \nabla ^2 \left(\frac{1}{r}\right) \right]= \lim_{\text {r} \to 0} \left[ \nabla \cdot \nabla \left( \frac 1 r \right) \right]= \lim_{\text {r} \to 0} \left(\frac {-4 \pi} {v(r)} \right) = \text {infinite} $$
$$ \lim_{\text r\to 0} \int \nabla ^2 \left( \frac 1 r \right) dv(r) = \lim_{\text r\to 0}\int \frac {-4 \pi} {v(r)} dv(r) = -4\pi$$
Since the laplacian is zero everywhere except $r \to 0$, it is true and real that :-
$$ \nabla ^2 \left(\frac{1}{r}\right)=-4\pi\delta^3({\bf r}) $$
