Divergence of $\vec{F} = \frac{\hat{\mathrm{r}}}{r^{2}}$ Consider the vector field
$$\vec{\mathrm{F}} = \frac{\hat{\mathrm{r}}}{r^{2}},$$
then the divergence of this field is:
$$\vec{\nabla}\cdot\left(\frac{\hat{\mathrm{r}}}{r^{2}}\right) = 4\pi\delta^{3}(\mathrm{\vec{r}})$$
What is the proof of this relation?
 A: What does "$4\pi\delta^3(r)$" mean? It means if you integrate any compactly supported test function $f$ against it, the result will be $4\pi f(0)$. 
So what you really want to show is that for any compactly supported test function $f$, you have
$$ \int \bigg(\nabla\cdot \frac{\hat{r}}{r^2}\bigg)f(r)\ dr = 4\pi f(0). $$
Integrate by parts and take a limit around the singularity at $0$ to get the answer.
A: Can't really give "the" proof if I don't know what you're allowed to assume.
But, to suit my personal taste I would start with the well known identity
$ \int d^3x \mathbf{\nabla\cdot F} = \oint d^2x \mathbf{\hat{n}}\cdot\mathbf{F}$  Insert $\mathbf{F}=\mathbf{\hat{r}}/r^2$ and take it from there.
A: A common way to show that $\nabla \cdot \left(\frac{\hat r}{r^2}\right)=4\pi \delta (\vec r)$ is to regularize the function $\left(\frac{\hat r}{r^2}\right)$ in terms of a parameter, say $a$.  To that end, let $\vec \psi$ be the regularized function given by 
$$\vec \psi(\vec r;a)=\frac{\vec r}{(r^2+a^2)^{3/2}} \tag 1$$
Taking the divergence of $(1)$ reveals that 
$$\nabla \cdot \vec \psi(\vec r; a)=\frac{3a^2}{(r^2+a^2)^{5/2}}$$
Now, for any sufficiently smooth test function $\phi$, we have that
$$\begin{align}
\lim_{a \to 0}\int_V \nabla \cdot \vec \psi(\vec r; a)\phi(\vec r)dV&=\lim_{a \to 0}\int_V \frac{3a^2}{(r^2+a^2)^{5/2}}\phi(\vec r)dV\\\\
&=0
\end{align}$$
if $V$ does not include the origin.  
Now, suppose that $V$ does include the origin.  Then, we have
$$\begin{align}
\lim_{a \to 0}\int_V \nabla \cdot \vec \psi(\vec r; a)\phi(\vec r)dV&=\lim_{a\to 0}\int_{V-V_{\delta}} \frac{3a^2}{(r^2+a^2)^{5/2}}\phi(\vec r)dV+\lim_{a\to 0}\int_{V_{\delta}} \frac{3a^2}{(r^2+a^2)^{5/2}}\phi(\vec r)dV\\\\
&=\lim_{a\to 0}\int_{V_{\delta}} \frac{3a^2}{(r^2+a^2)^{5/2}}\phi(\vec r)dV
\end{align}$$
where $V_{\delta}$ is a spherical region centered at $\vec r=0$ with radius $\delta$.  For any $\epsilon>0$, take $\delta>0$ such that $|\phi(\vec r)-\phi(0)|\le \epsilon/(4\pi)$ whenever $0<|\vec r|< \delta$.  Then, we have
$$\begin{align}
\lim_{a \to 0}\left|\int_V \nabla \cdot \vec \psi(\vec r; a)(\phi(\vec r)-\phi(0))\,dV\right|&\le \lim_{a\to 0} \int_{V_{\delta}} \left|\phi(\vec r)-\phi(0)\right|\frac{3a^2}{(r^2+a^2)^{5/2}}dV\\\\
&\le \left(\frac{\epsilon}{4\pi}\,4\pi\right)\,\lim_{a \to 0}\int_{0}^{\infty}\frac{3a^2}{(r^2+a^2)^{5/2}}r^2\,dr\\\\
&\le \epsilon
\end{align}$$
Thus, we have for any test function $\phi$, 
$$\begin{align}
\lim_{a \to 0}\int_V \nabla \cdot \vec \psi(\vec r; a)\phi(\vec r)\,dV&=4\pi \phi(0)
\end{align}$$
and it is in this sense that 
$$\bbox[5px,border:2px solid #C0A000]{\lim_{a\to 0} \nabla \cdot \vec \psi(\vec r;a)=4\pi \delta(\vec r)}$$
A: The below is a more detailed version of bob.sacamento's answer. I will use boldface for vector quantities.
Suppose we integrate over a sphere of radius $R$ centered at the origin. The surface
integral is
$$
\begin{aligned}
\oint \mathbf{F} \cdot d \mathbf{a} &=\int\left(\frac{1}{R^{2}} \hat{\mathbf{r}}\right) \cdot\left(R^{2} \sin \theta d \theta d \phi \hat{\mathbf{r}}\right) \\
&=\left(\int_{0}^{\pi} \sin \theta d \theta\right)\left(\int_{0}^{2 \pi} d \phi\right)=4 \pi
\end{aligned}
$$
Using the divergence theorem:
$$
\begin{aligned}
\oint \mathbf{F} \cdot d \mathbf{a} &=\int \nabla\cdot{\bf F}\, d{\mathbf{\tau}}
\end{aligned}
$$
where ${\mathbf{\tau}}$ is the volume differential. The above two equations imply that
$$
\begin{aligned}
\int \nabla\cdot{\bf F}\, d{\mathbf{\tau}} &=4\pi\, 
\end{aligned}
$$
Therefore,
$$
\begin{aligned}
\nabla\cdot\left(\frac{1}{r^{2}} \hat{\mathbf{r}}\right) &=4\pi\delta^3({\bf r}) \, .
\end{aligned}
$$
A: IMHO I think that all the other answers of this question are excessively complex. A very simple explanation, just based in the generalized Stokes theorem, is as follow: suppose that $R$ is an open, bounded and convex region in $\mathbb{R}^3$ and let $\mathbf{F}:=\frac{\mathbf{\hat r}}{r^2}$, then it follows from Stokes theorem that if $\mathbf{0}\notin R$ then $\operatorname{div}\mathbf{F}$ is well-defined and equal to zero at every point of $R$, then from Stokes theorem we have that
$$
\int_{\partial R}\mathbf{F}\cdot d\mathbf{S}=\int_{R}\operatorname{div} \mathbf{F}\, dV=\int_{R}0\, dV=0\tag1
$$
Now, to handle the case when $\mathbf{0}\in R$ we can exploit (1) and the simple result
$$
\int_{\partial \mathbb{B}(\mathbf 0,\epsilon )}\mathbf{F} \cdot d\mathbf{S}=4\pi\tag2
$$
what follows easily using spherical coordinates in $\mathbb{B}(\mathbf 0,\epsilon )$. Then as $R$ is open by assumption then there exists some $\epsilon >0$ such that $\mathbb{B}(\mathbf{0},\epsilon )\subset R$, therefore
$$
0=\int_{R\setminus \mathbb{B}(\mathbf{0},\epsilon )}\operatorname{div} \mathbf{F}\,dV=\int_{\partial (R\setminus \mathbb{B}(\mathbf 0,\epsilon ))}\mathbf{F}\cdot d\mathbf{S}=\int_{\partial R}\mathbf{F}\cdot d\mathbf{S}-\int_{\partial \mathbb{B}(\mathbf 0,\epsilon )}\mathbf{F}\cdot  d\mathbf{S}\\
\therefore\quad \int_{\partial R}\mathbf{F}\cdot d\mathbf{S}=4\pi\tag3
$$
Then, to extend the Stokes theorem to the case where $\mathbf{0}\in R$ we can set $\operatorname{div}\left(\frac{\mathbf{\hat r}}{r^2}\right):=4\pi \delta ^3$ and with this definition we have that
$$
\int_{\partial R}\mathbf{F}\cdot d\mathbf{S}=4\pi=\int_{R}\operatorname{div}\mathbf{F}\,dV\tag4
$$
∎
