Divergence of $\vec{F} = \frac{\hat{\mathrm{r}}}{r^{2}}$

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

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.

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.