Divergence of $\vec{f} = \frac{1}{r^2} \hat{r}$

The divergence of $\vec F(\vec r)=\frac{\vec r}{r^3}$ is not $4\pi$. By direct computation, we have

$$\nabla \cdot \vec F(\vec r)=\begin{cases}0&,r\ne 0\\\\\text{undefined}&,r=0\end{cases}$$

However, we can assign meaning to $\nabla \cdot \vec F(\vec r)$ at the origin and write

$$\bbox[5px,border:2px solid #C0A000]{\nabla \cdot \vec F(\vec r)\sim 4\pi \delta(\vec r)}$$

where $\delta(\vec r)$ is the Dirac Delta distribution. To do this, we introduce the regularization function

$$\psi(\vec r;a)=\frac{\vec r}{(r^2+a^2)^{3/2}}\tag 1$$

Note that for $r\ne 0$, $\lim_{a\to 0}\psi(\vec r;a)=\vec F(\vec r)$. Taking the divergence of $\psi(\vec r;a)$ yields

$$\nabla \cdot \psi(\vec r;a)=\frac{3a^2}{(r^2+a^2)^{5/2}}$$

Now, in THIS ANSWER, I showed that for any smooth test function $\phi(\vec r)$

$$\begin{align} \lim_{a \to 0}\int_V \nabla \cdot \vec \psi(\vec r; a)\phi(\vec r)\,dV&=\begin{cases}4\pi \phi(0)&,0\in V\\\\0&,0\notin V\end{cases} \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)}$$