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?
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.