How to interpret the meaning of $\mathbf n d\sigma$ in terms of differential forms? We often write stokes's formula in $\mathbb R^n$ as
$$
\int_\Omega \nabla\cdot \mathbf f d\mu=\int_{\partial \Omega} \mathbf f \cdot \mathbf n d\sigma.
$$
My questions is: what does $\mathbf n d\sigma$ mean? it is written as if it is a vector, but $d \sigma$ is essentially a $(n-1)$-covector.  Multiplying a covector by $\mathbf n$ does not make sense.
Interpreting $\mathbf f \cdot \mathbf n d\sigma$ as the product of a function and a $(n-1)$-covector does not work as well, because when I calculate it in two dimensions, it does not work. Let's say that the normal vector is $(\sin \theta, -\cos \theta)$, then $\mathbf f \cdot \mathbf n = f_1 \sin \theta -f_2 \cos \theta$, and $d \sigma = \cos \theta dx_1 + \sin \theta dx_2$. Multiplying $\mathbf f \cdot \mathbf n$ and $d \sigma$ does NOT lead to the expected expression $f_1 dx_2 -f_2 dx_1$ or something similar.
So how could I translate $\mathbf n d\sigma$ to the language of differential form?
 A: The other answer gives you a nice way to formulate the divergence theorem abstractly on a Riemannian manifold. The issue you're facing is that your expression for $d\sigma$ is not quite correct; you've written it using the forms $dx$ and $dy$, but these are technically speaking differential forms on $\Bbb{R}^2$, while $d\sigma$ is supposed to be a differential form on a $2-1=1$-dimensional embedded submanifold $N$ of $\Bbb{R}^2$.
So, technically, if you let $\iota:N \to \Bbb{R}^2$ be the inclusion mapping, then you should be using $\iota^*(dx) = d(\iota^*x)$ and $\iota^*(dy) = d(\iota^*y)$ when writing down $d\sigma$, NOT $dx$ and $dy$. I'll use $\xi:= \iota^*x$ and $\eta:= \iota^*y$, note that these are functions on $N$, which shouldn't be confused with the coordinate functions $x$ and $y$, which are defined on all of $\Bbb{R}^2$. Now, let $\mathbf{f}$ be the vector field and $\mathbf{n}$ the unit normal to $N$.
Then, $d\sigma = n_1\,d\eta - n_2\,d\xi$ is the line element to $N$ (in abstract language, if $dV$ is the volume form on the big manifold $M$, then the volume form $d\sigma$ on the submanifold is $\iota^*(\mathbf{n} \mathbin\lrcorner dV)$; i.e you first take the interior product of the normal $\mathbf{n}$ with the volume form $dV$, and then pull-back the whole thing to $N$). Next, $\langle \mathbf{f},\mathbf{n}\rangle = f_1n_1 + f_2 n_2$. If we multiply these together and use the fact $(n_1)^2 + (n_2)^2 = 1$, then a few lines of algebra will show that
\begin{align}
\langle \mathbf{f},\mathbf{n}\rangle\, d\sigma &= (f_1 n_1 + f_2 n_2)\cdot (n_1\,d\eta - n_2\,d\xi) \\
&= (f_1 d \eta - f_2 d\xi) + (f_2 n_1 - f_1 n_2)\underbrace{(n_1\, d\xi + n_2\, d\eta)}_{=0}
\end{align}
That last term is zero, because if $\mathbf{n} = n_1\frac{\partial}{\partial x} + n_2 \frac{\partial}{\partial y}$ is normal to $N$, then $\tau = -n_2\frac{\partial}{\partial x} + n_1 \frac{\partial}{\partial y}$ is tangent to $N$. So, if you evaluate the underlined $1$-form on any vector which is tangent to $N$, the result is zero.
I'm sure that when you wrote things out using $dx$ and $dy$, you got something similar, but were just unsure of how to "get rid of" the ugly term. Well, the point is that you have to look at the differential forms when pulled back to $N$, because if you look at $n_1dx + n_2 dy$, and you take a point $p\in N$, then you get a covector $\mu_p = n_1(p) dx_p + n_2(p) dy_p: T_pM=T_p\Bbb{R}^2 \to \Bbb{R}$. Unfortunately this map $\mu_p$ isn't necessarily the zero function. It is only when you restrict $\mu_p$ to the smaller subspace  $T_pN$ (which is what pulling-back using $\iota$ does), then you get the zero map.

By the way if you refer to the other answer, the place where you're running into trouble in your example is with the formula $\iota_X(dV) = \langle X, \mathbf{n}\rangle \, dS$ (this again isn't technically right, on the LHS, we have a form on $M$, so to get a true equality, it needs to be pulled back to $N$). If you take a book on differential geometry and see how this formula is proven, we first decompose $X = X^{\perp} + X^{\parallel}$ into a part which is perpendicular to $N$ and a part which is parallel to $N$, then it will turn out that the pull-back of the parallel part's interior product with $dV$ evaluates to zero (similar to how it was zero above).
A: Domains in $\mathbb{R}^n$ are naturally Riemannian manifolds, so I may answer this question in terms of differential forms on Riemannian manifolds. Let $(M,g)$ be a Riemannian manifold of dimension $n$,for simplicity, let assume $M$ is compact. In this case, your formula simply reads
$$\int_M \mathrm{div}(X)dV = \int_{\partial M} \left \langle X, \mathbf{n} \right \rangle dS$$
where $\mathbf{n}$ is outward normal vector field, that is, if $\partial_1,...,\partial_{n-1}$ is a basis of the tangent space $T_x(\partial M)$ then $\mathbf{n},\partial_1,...,\partial_{n-1}$ is a basic of $T_x(M)$; $dS$ is the  volume form of its boundary $\partial M$ and $X$ is just a vector field. The divergence operator $\mathrm{div}$ is
$$\mathrm{div}(X) = \sum_{i=1}^{\mathrm{dim}(M)}\frac{\partial X_i}{\partial x_i}$$ in terms of a normal coordinate $(U,x_i)$. By definition of the volume form
$$\iota_X(dV) = \left \langle X, \mathbf{n} \right \rangle dS$$
where $\iota_X$ is contraction operator. That is, $\iota_X(\omega)(Y_1,...Y_{r-1}) = \omega(X,Y_1,...,Y_{r-1})$ for every $r$-form $\omega$. Therefore, in terms of differential forms, the div-formula is
$$\int_M \mathrm{div}(X)dV = \int_{\partial_M}\iota_X(dV).$$
Actually, the above formula is called divergence theorem. The key ideas in its proof are Stoke's theorem and Cartan magic formula. In case of closed manifolds, I could quickly give a proof here. What we need to prove is
$$\int_M \mathrm{div}(X)dV = 0$$
since the boudary-term vanishes, hence we would find a form $\eta$ such that $d\eta = \mathrm{div}(X)dV$. I claim that $\eta = \iota(X)dV$ is such a form. Using a normal coordinate $(U,x_i)$ we find out that
$$\begin{cases} dV = dx_1 \wedge ... \wedge dx_n \\ \iota(X)dV = \sum_{i=1}^n (-1)^{i-1}X_idx_1 \wedge ... \hat{dx_i} \wedge ... \wedge dx_n & \mathrm{with} \ X = \sum_{i=1}^n X_i \frac{\partial}{\partial x_i}. \end{cases}$$
The rest of the proof follows easily by direct computations.
