$\nabla \cdot \big(\frac{\hat{r}}{r^{2}}\big)$ and Dirac Delta Function
Solution 1:
Take $\phi \in C_c^\infty(\mathbb R^3)$ and consider the integral $$I := \int_{\mathbb R^3} \bigg( \nabla \cdot \frac{\mathbf x}{|\mathbf x|^3} \bigg) \phi(\mathbf x) \, d\mathbf x$$
We can rewrite this using $\nabla \cdot (\mathbf F \phi) = (\nabla \cdot \mathbf F) \phi + \mathbf F \cdot \nabla \phi$: $$I = \int_{\mathbb R^3} \nabla \cdot \bigg( \frac{\mathbf x}{|\mathbf x|^3} \phi(\mathbf x) \bigg) \, d\mathbf x - \int_{\mathbb R^3} \frac{\mathbf x}{|\mathbf x|^3} \cdot \nabla \phi(\mathbf x) \, d\mathbf x$$
Since $\phi(\mathbf x) = 0$ for big $|\mathbf x|$ the first term vanishes. Therefore, $$I = - \int_{\mathbb R^3} \frac{\mathbf x}{|\mathbf x|^3} \cdot \nabla \phi(\mathbf x) \, d\mathbf x$$
Now we change to polar coordinates, $\mathbf x = r \mathbf n$: $$I = - \int_0^\infty \int_{S^2} \frac{\mathbf n}{r^2} \cdot \nabla \phi(\mathbf x) \, r^2 d\Omega(\mathbf n) \, dr = - \int_0^\infty \int_{S^2} \mathbf n \cdot \nabla \phi(\mathbf x) \, d\Omega(\mathbf n) \, dr$$ where $d\Omega$ is the area measure on $S^2 = \{ \mathbf x \in \mathbb R^3 : |\mathbf x| = 1 \}$.
In polar coordinates, $\nabla \phi = \frac{\partial\phi}{\partial r} \mathbf n + (\text{terms orthogonal to $\mathbf n$})$ so $$I = - \int_0^\infty \int_{S^2} \frac{\partial\phi}{\partial r}(\mathbf x) \, d\Omega(\mathbf n) \, dr$$
Let us do the $r$ integral first: $$I = - \int_{S^2} \int_0^\infty \frac{\partial\phi}{\partial r}(\mathbf x) \, dr \, d\Omega(\mathbf n) = - \int_{S^2} \big[ \phi(r \mathbf n) \big]_{r=0}^{\infty} \, d\Omega(\mathbf n) = \int_{S^2} \phi(\mathbf 0) \, d\Omega(\mathbf n)$$
The left integral is trivial since we now have no dependency on $\mathbf n$ and the area of the unit sphere is $4\pi$: $$I = 4\pi \, \phi(\mathbf 0)$$ which can be written as $$I = 4\pi \int_{\mathbb R^3} \delta(\mathbf x) \phi(\mathbf x) \, d\mathbf x$$
Thus, $$\int_{\mathbb R^3} \bigg( \nabla \cdot \frac{\mathbf x}{|\mathbf x|^3} \bigg) \phi(\mathbf x) \, d\mathbf x = \int_{\mathbb R^3} 4\pi \, \delta(\mathbf x) \phi(\mathbf x) \, d\mathbf x$$
Since this is valid for all $\phi \in C_c^\infty(\mathbb R^3)$ we have $$\nabla \cdot \frac{\mathbf x}{|\mathbf x|^3} = 4\pi \, \delta(\mathbf x)$$
Edit
I got a question why $\int_{\mathbb R^3} \nabla \cdot \left( \frac{\mathbf x}{|\mathbf x|^3} \phi(\mathbf x) \right) \, d\mathbf x$ vanishes. It does so by divergence theorem. Since $\phi$ has compact support, for big enough $R$ we have $$ \int_{|\mathbf{x}|<R} \nabla \cdot \bigg( \frac{\mathbf x}{|\mathbf x|^3} \phi(\mathbf x) \bigg) \, d\mathbf x = \oint_{|\mathbf{x}|=R} \frac{\mathbf x}{|\mathbf x|^3} \phi(\mathbf x) \, \mathbf n \cdot dS(\mathbf x) = 0. $$ Therefore, $$ \int_{\mathbb R^3} \nabla \cdot \bigg( \frac{\mathbf x}{|\mathbf x|^3} \phi(\mathbf x) \bigg) \, d\mathbf x = \lim_{R\to\infty} \int_{|\mathbf{x}|<R} \nabla \cdot \bigg( \frac{\mathbf x}{|\mathbf x|^3} \phi(\mathbf x) \bigg) \, d\mathbf x = \oint_{|\mathbf{x}|=R} \frac{\mathbf x}{|\mathbf x|^3} \phi(\mathbf x) \, \mathbf n \cdot dS(\mathbf x) = 0. $$