I am trying understand how circulation density arises in Green's theorem and I'd like to know if my line of thinking is on the right track. Here it goes :).
Idea
We know that if we have a vector field $\vec F=P\hat x+Q\hat y\in\mathbb{C}^1$ on an open region $D$ containing the simple region $R$ whose boundary is $\Gamma$, then we can “cut” $\Gamma$ into $K$ positively-oriented rectangular paths such that $\Gamma=\bigcup_{k=1}^K\Gamma_k$.
Consequently,
$$\oint_\Gamma \vec F\cdot d\vec r=\sum_{k=1}^K \oint_{\Gamma_k}\vec F\cdot d\vec r $$
which says the the macroscopic circulation of $\vec F$ along the curve $\Gamma$ is equal to the sum of microscopic circulations over the region $R$ enclosed by $\Gamma$.
Moreover, we know that the circulation along one of these rectangles is
$$\oint_{\Gamma_k}\vec F\cdot d\vec r=\bigg(\frac{Q(x+\Delta x,\;y)-Q(x,y)}{\Delta x}-\frac{P(x,\;y+\Delta y)-P(x,y)}{\Delta y}\bigg)\;\Delta x\Delta y$$
(for proof, see Green's Theorem in the Plane: Circulation Density).
Hence, the circulation density, or curl, of $\vec F$ at any point $(x,y)$ in $R$ is given by
$$\text{curl}\vec F(x,y)=\lim_{\Delta A\rightarrow 0}\frac{1}{\Delta A}\oint_\Gamma \vec F \cdot d\vec r=Q_x-P_y,$$
which is a scalar function.
Now, if we let $\sigma (x,y) = Q_x-P_y$ and integrate over the $R$, we obtain
$$\sum_{i=1}^n\sum_{j=1}^m\sigma(x^*_{ij},y^*_{ij})\Delta x\Delta y $$
hence,
$$\lim_{n,m\rightarrow\infty}\sum_{i=1}^n\sum_{j=1}^m\sigma(x^*_{ij},y^*_{ij})\Delta x\Delta y =\iint_R \sigma(x,y) dxdy=\iint_R (Q_x-P_y )dA$$
whose RHS can be related back to the line integral $\oint_\Gamma\vec F\cdot d\vec r$ giving us what we know as Green's Theorem.