Show that $\nabla\cdot\left(\dfrac{\mathbf{e}_r}{r^2}\right)=4\pi\delta(\mathbf{r})$ using the divergence theorem.

The surface integral can be evaluated in spherical coordinates:

1). Set

\begin{align} x&=r\sin\theta\cos\phi\;,\\ y&=r\sin\theta\sin\phi\;,\\ z&=r\cos\theta\;, \end{align}

2). Write

$\vec r(\theta,\phi)=r\sin\theta\cos\phi \vec i+r\sin\theta\sin\phi\vec j+r\cos\theta \vec k$

3). Then

$\vec r_{\theta}(\theta, \phi)=r\cos\theta\cos\phi \vec i+r\cos\theta\sin\phi\vec j-r\sin\theta \vec k$

and

$\vec r_{\phi}(\theta, \phi)=-r\sin\theta\sin\phi \vec i+r\sin\theta\cos\phi\vec j$

4). Find the magnitude of the cross product:

$\color{red}{\vert \vec r_{\phi}\times \vec r_{\theta}\vert =r^2\sin \theta.} $ $(\sin \theta \geq 0$ since $0\leq \theta <\pi)$,

so that

$d\sigma =r^2\sin \theta d\theta d\phi$.

5). Find the limits of integration:

$0\leq \phi <2\pi;\ 0\leq \theta <\pi$ because the surface is a sphere.

6). Substitute into the integral:

$\color{blue}{\iint_\limits{\large\text{surface enclosing}\, \tau}\dfrac{\mathbf{e}_r}{r^2}\cdot\mathbf{e}_r\,\mathrm{d}\sigma }=\int \int _{\sigma }\frac{1}{r^{2}}d\sigma =\int^{2\pi }_0 \int^{\pi }_0 \frac{1}{r^{2}}(r^{2}\sin \theta )d\theta d\phi=$

The vector $\vec r$ is a function of $\theta$ and $\phi$; which is given by $$\vec r=r\sin\theta\cos\phi \widehat i+r\sin\theta\sin\phi\widehat j+r\cos\theta \widehat k$$ and the unit vector $\widehat r$ is given by $$\color{teal}{\widehat r=\sin\theta\cos\phi \widehat i+\sin\theta\sin\phi\widehat j+\cos\theta \widehat k}\tag{2}$$

Therefore,

$$\frac{\partial r}{\partial \theta}=r\cos\theta\cos\phi \widehat{i}+r\cos\theta\sin\phi\widehat{j}-r\sin\theta \widehat{k}$$ and $$\frac{\partial r}{\partial \phi}=-r\sin\theta\sin\phi \widehat i+r\sin\theta\cos\phi\widehat j$$

So $$\begin{align}\frac{\partial r}{\partial \theta}\times\frac{\partial r}{\partial \phi}&=\begin{vmatrix} \widehat i & \widehat j & \widehat k \\ r\cos\theta\cos\phi & r\cos\theta\sin\phi & -r\sin\theta \\ -r\sin\theta\sin\phi & r\sin\theta\cos\phi & 0 \end{vmatrix}\\&=r^2\sin^2\theta\cos\phi\widehat i +r^2\sin^2\theta\sin\phi\widehat j+ r^2\cos\theta\sin\theta\cos^2\phi\widehat k+r^2\cos\theta\sin\theta\sin^2\phi\widehat k\\&=r^2\sin^2\theta\cos\phi\widehat i +r^2\sin^2\theta\sin\phi\widehat j+ r^2\cos\theta\sin\theta\left(\cos^2\phi+\sin^2\phi\right)\widehat k\\&=r^2{\left(\sin^2\theta\cos\phi\widehat i+\sin^2\theta\sin\phi\widehat j+\cos\theta\sin\theta\widehat k\right)}\\&=r^2\sin\theta\color{purple}{\underbrace{\left(\sin\theta\cos\phi\widehat i+\sin\theta\sin\phi\widehat j+\cos\theta\widehat k\right)}_{\color{teal}{\Large\text{From (2)}\,=\,\widehat r}}}\\&=r^2\sin\theta\,\widehat r\end{align}$$

Taking the magnitude gives $$\left|\frac{\partial r}{\partial \theta}\times\frac{\partial r}{\partial \phi}\right|=\sqrt{\left(r^2\sin\theta\,\widehat r\right)^2}={\sqrt{r^4\sin^2\theta \,\widehat r^2}}=\underbrace{r^2\sin\theta}_{\Large\text{Since}\,\color{teal}{\widehat r}^2=1}$$ Therefore $$\color{#F80}{\mathrm{d}\sigma=\left|\frac{\partial r}{\partial \theta}\times\frac{\partial r}{\partial \phi}\right|\,\mathrm{d}\theta\,\mathrm{d}\phi}\tag{1}$$ $$=r^2\sin\theta\,\mathrm{d}\theta\,\mathrm{d}\phi$$

As helpfully pointed out in a comment given by @Chilango we are integrating over a sphere of radius $r$. So the azimuthal angle $\phi$ has a limit range given by $0\leq \phi \lt 2\pi$ and the polar angle $\theta$ has a limit range given by $\ 0\leq \theta \lt\pi$.

Finally, substitution of $(1)$ and the limits into the $\color{blue}{\mathrm{blue}}$ integral gives the $\color{#180}{\mathrm{green}}$ integral as required.

A special thanks goes to @Chilango for his/her answer that gave me the ideas necessary to give this answer with full justification that $$\left|\frac{\partial r}{\partial \theta}\times\frac{\partial r}{\partial \phi}\right|=r^2\sin\theta$$

A note on what I think to be a misleading trend in some texts, where I have recently stumbled.

As a proper Riemann integral, $$\iiint_{\tau}\nabla\cdot\left(\frac{\mathbf{e}_r}{r^2}\right)d\tau$$ where $\frac{\mathbf{e}_r}{r^2}=\frac{\mathbf{x}-\mathbf{y}}{\|\mathbf{x}-\mathbf{y}\|^3}$ (it seems that $\mathbf{y}=\mathbf{0}$ in your case), is $0$ if $\mathbf{y}\notin\bar{\tau}$ and does not exist if $\mathbf{y}\in\bar{\tau}$, because the integrand of a Riemann integral, in the usual calculus definitions of it, has to be defined on all the domain.

As the limit $$\lim_{\varepsilon\to 0}\iiint_{\tau\setminus B(\mathbf{y},\varepsilon)}\nabla\cdot\left(\frac{\mathbf{x}-\mathbf{y}}{\|\mathbf{x}-\mathbf{y}\|^3}\right)dx_1dx_2dx_3$$ it is $0$ because $\lim_{\varepsilon\to 0}0=0$ and the same holds for the Lebesgue integral $$\int_{\tau}\nabla\cdot\left(\frac{\mathbf{x}-\mathbf{y}}{\|\mathbf{x}-\mathbf{y}\|^3}\right)d\mu_{\mathbf{x}}$$which is calculated as the preceding limit.

This shows that, under these definitions of the integral and the usual definition of the derivative, the divergence theorem, certainly valid if, instead of $\frac{\mathbf{e}_r}{r^2}$ you had a vector field $\mathbf{F}\in C^1(\mathring{A})$ with $\tau\subset\mathring{A}$ satisfying opportune assumptions, cannot be applied.

Since $\forall\mathbf{x}\in\mathbb{R}^3\setminus\{\mathbf{y}\}\quad\frac{\mathbf{x}-\mathbf{y}}{\|\mathbf{x}-\mathbf{y}\|^3}=-\nabla\left(\frac{1}{\|\mathbf{x}-\mathbf{y}\|}\right)$ and the divergence of the gradient is the Laplacian $\nabla\cdot\nabla=\nabla^2$, we see that $\nabla\cdot\left(\frac{\mathbf{x}-\mathbf{y}}{\|\mathbf{x}-\mathbf{y}\|^3}\right)=-\nabla^2\left(\frac{1}{\|\mathbf{x}-\mathbf{y}\|}\right)$. Then by reading the integral $$\int_{\tau}-\nabla^2\left(\frac{1}{r}\right)\varphi \,d\tau$$ where $\varphi\in C^2(\mathbb{R}^3)$ (typically $\varphi\in C^\infty(\mathbb{R}^3)$) is such that $\forall\mathbf{x}\notin\tau\quad\varphi(\mathbf{x})=0$, in the symbolic way of the Laplacian of the distribution defined by $-\frac{1}{r}$, it can be shown, as it is here, that $-\int_{\tau}\nabla^2\left(\frac{1}{r}\right)\varphi \,d\tau=4\pi\int_{\tau}\delta_{\mathbf{y}}\varphi \,d\tau:=4\pi\varphi(\mathbf{y})$ (where $\delta_{\mathbf{y}}(\mathbf{x}):=\delta(\mathbf{x}-\mathbf{y})$), but that is$$\lim_{\varepsilon\to 0}\iiint_{\tau\setminus B(\mathbf{y},\varepsilon)}\frac{-\nabla^2\varphi(\mathbf{x})}{\|\mathbf{x}-\mathbf{y}\|}dx_1dx_2dx_3=\int_{\tau}\frac{-\nabla^2\varphi(\mathbf{x})}{\|\mathbf{x}-\mathbf{y}\|}d\mu_{\mathbf{x}}$$if we use the usual Riemann (one th left) or Lebesgue (on the right) integrals, while $$\lim_{\varepsilon\to 0}\iiint_{\tau\setminus B(\mathbf{y},\varepsilon)}\nabla\cdot\left(\frac{\mathbf{x}-\mathbf{y}}{\|\mathbf{x}-\mathbf{y}\|^3}\right)\varphi(\mathbf{x})dx_1dx_2dx_3=\int_{\tau}\nabla\cdot\left(\frac{\mathbf{x}-\mathbf{y}}{\|\mathbf{x}-\mathbf{y}\|^3}\right)\varphi(\mathbf{x})d\mu_{\mathbf{x}}\equiv 0$$for all $\mathbf{y}$ and all functions $\varphi$.