Green's theorem and flux 
Given the vector field $\vec{F}(x,y) = (x^2+y^2)^{-1}\begin{bmatrix} x \\ y \end{bmatrix}$, calculate the flux of $\vec{F}$ across the circle $C$ of radius $a$ centered at the origin (with positive orientation).

It is my understanding that Green's theorem for flux and divergence says
$$
\int\limits_C \Phi_{\vec{F}} = \int\limits_C P\,dy - Q\,dx = \iint\limits_R\nabla\cdot\vec{F}\,dA
$$
if $\vec{F} = \begin{bmatrix} P & Q \end{bmatrix}$ (omitting other hypotheses of course). Note that $R$ is the region bounded by the curve $C$. For our $\vec{F}$, we have $\nabla\cdot\vec{F} = 0$. So shouldn't the flux of $\vec{F}$ through $C$ be zero? Everytime I try to compute it, I obtain zero. However, it appears that I am supposed to be getting $2\pi$. Can someone help me to understand where I am going wrong?
 A: 
Standard way:
Notice the interior of the area enclosed the curve $C$ has a singularity of the vector field, i.e., $\displaystyle \frac{\partial Q}{\partial x}$ or $\displaystyle \frac{\partial P}{\partial y}$ is not continuous. 
Standard procedure is to cut a hole with a small radius $r$ centered at this singularity , in this example it is $(0,0)$, such that outside this small disk of radius $r$ there is no singularity (See the picture above).
Denote $C'$ the boundary of this smaller disk $\Omega'$, rotating clockwisely. Then by $P = x/(x^2+y^2)$, $Q = y/(x^2+y^2)$:
\begin{align}
 &\oint_{C} \frac{xdy-ydx}{x^2+y^2} + \oint_{C'} \frac{xdy-ydx}{x^2+y^2} 
\\
=&
\int_{\Omega\backslash \Omega'}\left\{ \frac{\partial}{\partial x}\left(\frac{x}{x^2+y^2}\right)+\frac{\partial}{\partial y}\left(\frac{y}{x^2+y^2}\right)\right\} dx dy = 0.
\end{align}
Now parametrize the curve $C'$ using $t$, by letting $x = r\cos t$, $y = r\sin t$, so $dx = -r\sin t\, dt$, $dy = r \cos t \,dt$, $t$ from $2\pi$ to $0$ (clockwise).  The integral becomes:
\begin{align}
& \oint_{C} \frac{xdy-ydx}{x^2+y^2} = \oint_{C'}\frac{y }{x^2+y^2}dx-\frac{x}{x^2+y^2}dy
\\
 =& \int_{2\pi}^0 \left(\frac{r\sin t\, }{r^2\cos^2 t + r^2\sin^2 t} (-r\sin t)
- \frac{r\cos t\,}{r^2\cos^2 t + r^2\sin^2 t}r \cos t \right)dt 
=2\pi.
\end{align}

Cheating way: since Muphrid mentioned the delta function, we can evaluate using divergence theorem (a little cheated): Notice $F = \nabla \phi$, where $\phi = \big(\ln(x^2+y^2)\big)$.
$$
\mathrm{Flux} = \int_{\partial \Omega} F\cdot n \,ds = \frac{1}{2}\int_{\partial \Omega} \nabla \phi \cdot n \,ds = \frac{1}{2} \int_{\Omega} \Delta \phi\,dx = \frac{1}{2} \int_{\Omega} 4\pi \delta_0(x) = 2\pi.
$$  

Cheating way #2 : The domain has winding number being 1, so by definition:
$$
\frac{1}{2\pi} \oint_C \,\frac{x}{x^2+y^2}\,dy - \frac{y}{x^2+y^2}\,dx= 1.
$$
A: One way to get it is by using the middle expression for the line integral along the boundary $C$ parametrized by $x = a\,\cos\theta$ and $y = a\,\sin\theta$.  You can see how $a$ cancels out from the integral, as it should.  So we could even pick $a$ just barely greater than zero and the same value results.  This means that the non-vanishing contribution to the area integral in the third expression you give must come entirely from $\nabla\cdot F$ right at the origin.  Do you know about Dirac $\delta$ functions?  It vanishes everywhere except at a point but has non-vanishing integral.  Though it technically is not a function and is defined only through integral expressions, it's still common to write $\nabla\cdot F = 2\pi\,\delta(\vec{x})$ for the case of the $F$ you have given.  It's tricky, you might be tempted to think $\nabla\cdot F = 0$ and that's true everywhere except the origin, but it matters in the case of this integral.
A: Your vector field is exactly the Green's function for $\nabla$:  it is the unique vector field so that $\nabla \cdot F = 2\pi \delta$, where $\delta$ is the Dirac delta function.  Try to look at the limiting behavior at the origin; you should see that this diverges.
It is rather difficult to actually justify this calculation as a student, however; you're probably intended to do the line integral only.
