Why am I getting zero divergence of function $\vec{f} = \frac{1}{r^2} \hat{r}$, where $r$ is the distance from the origin and $\hat{r}$ is the unit vector in the radial direction. The divergence of this function over a sphere of radius $R$, which includes the origin.
$$\nabla \cdot f = \frac{1}{r^2} \frac{\partial}{\partial r}(r^2 f_r) = \frac{1}{r^2} \frac{\partial}{\partial r}(r^2 \frac{1}{r^2}) = 0$$
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)}$$
The calculation result is:
$$\nabla \cdot \frac{1}{r^2} \hat{r} = \frac{1}{r^2} \frac{\partial}{\partial r} \left[r^2 \left( \frac{1}{r^2} \right) \right] = \frac 0 {r^2} = \begin{cases}0&,r\ne 0\\\\ \text{undefined}&,r=0\end{cases}$$
It correctly shows that the divergence is zero everywhere except the origin.
However, unfortunately, it only says that the divergence is not defined at the origin and cannot provide more information, that is, $\nabla \cdot \frac{1}{r^2} \hat{r}$ is actually positive infinity at the origin.
Nevertheless, we can use the following simple mathematical way to obtain the situation at the origin point.
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 \left( \frac {1} {r^2} \right)= 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 {r} \to 0} \left(\frac {4 \pi} {v(r)}\right) $$
So obviously,
$$ \lim_{\text {r} \to 0} \left[ \nabla \cdot \left( \frac{1} {r^2} \hat r \right) \right] = \lim_{\text {r} \to 0} \left(\frac {4 \pi} {v(r)} \right) = \text {positive infinite} \qquad (1)$$
$$ \lim_{\text r\to 0} \int \nabla \cdot \left( \frac 1 {r^2} \hat r \right) dv(r) = \lim_{\text r\to 0}\int \frac {4 \pi} {v(r)} dv(r) = 4\pi \qquad (2)$$
Since the laplacian is zero everywhere except $r \to 0$, and according to equation (1) and (2), it is true and real that :-
$$ \nabla \cdot \left(\frac{1}{r^2} \hat r\right)=4\pi\delta^3({\bf r}) $$
$\hat{r} = (x,y,z)/\sqrt{x^2 + y^2 + z^2}$ $$ F(x,y,z) = \frac{1}{(x^2 + y^2 + z^2)^{3/2}} (x,y,z) \\ \frac{\partial}{\partial x} F = \frac{\partial}{\partial x} \frac{x}{(x^2 + y^2 + z^2)^{3/2}} = \frac{-2x^2 + y^2 + z^2}{(x^2 + y^2 + z^2)^{5/2}} \\ \frac{\partial}{\partial y} F = \frac{\partial}{\partial y} \frac{y}{(x^2 + y^2 + z^2)^{3/2}} = \frac{x^2 -2y^2 + z^2}{(x^2 + y^2 + z^2)^{5/2}} \\ \frac{\partial}{\partial z} F = \frac{\partial}{\partial z} \frac{z}{(x^2 + y^2 + z^2)^{3/2}} = \frac{x^2 + y^2 -2z^2}{(x^2 + y^2 + z^2)^{5/2}} \\ $$ Putting together : $$ \nabla\cdot F = \frac{\partial}{\partial x} F + \frac{\partial}{\partial y} F + \frac{\partial}{\partial z} F = \frac{0}{(x^2 + y^2 + z^2)^{5/2}} = 0 $$
