Reasoning behind the divergence: $\overrightarrow{\nabla} \cdot (\frac{1}{r} \overrightarrow{e_r}) = \frac{1}{r}\delta(r)$ Short summary:
I want to calculate the divergence of $\overrightarrow{A}=\overrightarrow{(\frac{x}{x^2+y^2},\frac{y}{x^2+y^2})}^T=\frac{1}{r} \overrightarrow{e_r} $. So I calculate $\nabla \overrightarrow{A}= \frac{1}{r}\frac{\partial}{\partial r}(r\cdot\frac{1}{r})+\frac{1}{r}\frac{\partial v_\phi}{\partial \phi}=\frac{1}{r}(\frac{r-r}{r^2})$ by using the quotient rule. Then i see that the derivative is always zero except for r=0 where its infinity, does this mean its the dirac delta function? And why do I apply the product rule and dont say $\frac{\partial}{\partial r}(r\cdot\frac{1}{r})=\frac{\partial}{\partial r}1=0$?
Longer background not absolutley necessary for the question:
If I calculate the divergence in cartesian coordinates I get: $\nabla \overrightarrow{A}= \frac{y^2-x^2+x^2-y^2}{(x^2+y^2)^2}$ wich means i cant do this $\frac{\partial}{\partial r}(r\cdot\frac{1}{r})=\frac{\partial}{\partial r}1=0$, but I still dont get the mathematical reason behind it.
The exercise wich made me think about the question was the following exercise about the 2d divergence theorem:
If I calculate the flux of A trough the boundary of the square with the corners $(\pm1,\pm'1)$ I get 2π. Normaly Id need to make a second boundary around r=0 wich would give me -2π cause the normal vector points in the other direction. Therefore the flux of A trough the boundary would be 0, just like the integral of the divergence of A over the area of the square without the pole point. Then the divergence theorem would apply. My exercise wants me to calculate the integral of the divergence of A over the area of the square with the pole point. So i need to know if it is the dirac delta function, wich would make sense since: $\int_0^{2\pi}\int_0^{r(\phi)}\frac{1}{r}\delta(r)r \ dr \ d\phi=2\pi$. But this would mean that the divergence theorem still applies, wich it shouldnt because there is a pole point.
Thanks in advance.
 A: Let $\vec r=\hat x x+\hat yy$ and $r=\sqrt{x^2+y^2}$. In terms of classical real analysis, we have for $\vec A=\frac{\vec r}{r^2}$, $r\ne 0$
$$\nabla \cdot \vec A=0$$
That is to say, that the divergence of $\vec A$ is zero in the domain of $\vec A$.

METHODOLOGY $1$:
As a distribution, we adopt the approach used in THIS ANSWER for the 3-dimensional case, to show that $\nabla \cdot \vec A=2\pi \delta(\vec r)$.  
Proceeding, we regularize the function $\left(\frac{\vec 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} \tag 1$$
Taking the divergence of $(1)$ reveals that 
$$\nabla \cdot \vec \psi(\vec r; a)=\frac{2a^2}{(r^2+a^2)^2}$$
Now, for any sufficiently smooth test function $\phi$, we have that
$$\begin{align}
\lim_{a \to 0}\int_S \nabla \cdot \vec \psi(\vec r; a)\phi(\vec r)\,dS&=\lim_{a \to 0}\int_S \frac{2a^2}{(r^2+a^2)^2}\phi(\vec r)\,dS\\\\
&=0
\end{align}$$
if $S$ does not include the origin.  

Now, suppose that $S$ does include the origin.  Then, we have
$$\begin{align}
\lim_{a \to 0}\int_S \nabla \cdot \vec \psi(\vec r; a)\phi(\vec r)\,dS&=\lim_{a\to 0}\int_{S-S_{\delta}} \frac{2a^2}{(r^2+a^2)^2}\phi(\vec r)\,dS+\lim_{a\to 0}\int_{S_{\delta}} \frac{2a^2}{(r^2+a^2)^2}\phi(\vec r)\,dS\\\\
&=\lim_{a\to 0}\int_{S_{\delta}} \frac{2a^2}{(r^2+a^2)^2}\phi(\vec r)\,dS
\end{align}$$
where $S_{\delta}$ is the circular 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/(2\pi)$ whenever $0<|\vec r|< \delta$.  Then, we have
$$\begin{align}
\lim_{a \to 0}\left|\int_S \nabla \cdot \vec \psi(\vec r; a)(\phi(\vec r)-\phi(0))\,dS\right|&\le \lim_{a\to 0} \int_{S_{\delta}} \left|\phi(\vec r)-\phi(0)\right|\frac{2a^2}{(r^2+a^2)^2}\,dS\\\\
&\le \left(\frac{\epsilon}{2\pi}\,2\pi\right)\,\lim_{a \to 0}\int_{0}^{\infty}\frac{2a^2}{(r^2+a^2)^2}\,r\,dr\\\\
&\le \epsilon
\end{align}$$
Thus, we have for any test function $\phi$, 
$$\begin{align}
\lim_{a \to 0}\int_S \nabla \cdot \vec \psi(\vec r; a)\phi(\vec r)\,dS&=2\pi \phi(0)
\end{align}$$
and it is in this sense (i.e., as a distribution) that 
$$\bbox[5px,border:2px solid #C0A000]{\lim_{a\to 0} \nabla \cdot \vec \psi(\vec r;a)=2\pi \delta(\vec r)}$$

METHODOLOGY $2$:
Alternatively, we write in distribution 
$$\begin{align}
\langle \nabla \cdot \vec A, \phi \rangle&=-\langle \vec A, \nabla \phi \rangle\\\\
&=-\int_0^{2\pi}\int_0^\infty \left(\frac{\vec r}{r^2}\right)\cdot \nabla \phi(\vec r)\,r\,dr\,d\phi\\\\
&=-2\pi \int_0^\infty \frac{\partial \phi(\vec r)}{\partial r}\,dr\\\\
&=2\pi \phi(0)
\end{align}$$
as expected!
