How to prove $$\textbf{n}d^2\sigma=\frac{\partial r}{\partial u} \times \frac{\partial r}{\partial v}dudv$$ where $\sigma$ is an open surface that has a smooth parameterization $\textbf{r}(u,v)$ and $\textbf{n}$ is the unit vector normal to the oriented surface, using the notion of a Jacobian?
I understand that the cross product gives a vector normal to the surface, but I am not sure how to justify this rigorously. On the internet, most sources use this fact as if it is obvious so I am finding it difficult to find a formal argument that justifies this claim. Perhaps someone has an idea how to prove this?
It is difficult to prove the formula $${\bf n}\,{\rm d}\sigma={\bf r}_u\times{\bf r}_v\ {\rm d}(u,v)\tag{1}$$ rigorously without going into "geometric measure theory". In the first place: Do you have a clear idea what the "vectorial surface element" $\ {\bf n}\,{\rm d}\sigma\ $ is as a mathematical object? In the following I shall try to make $(1)$ plausible and prove that it behaves properly under parameter transformations.
A linear map $(u,v)\mapsto L(u,v)=u{\bf a}+v{\bf b}\in{\mathbb R}^3$ maps the unit square of the $(u,v)$-plane onto a parallelogram of area $|{\bf a}\times{\bf b}|$. In fact such a map multiplies the area of any large or small domain in the $(u,v)$-plane by this same factor $|{\bf a}\times{\bf b}|$. It follows that a general map $(u,v)\mapsto{\bf r}(u,v)\in{\mathbb R}^3$ maps a fine grid of mesh sizes $\Delta u$, $\Delta v$ in the $(u,v)$-plane to a mesh in space consisting of tiny parallelograms of area $|{\bf r}_u\times{\bf r}_v|\Delta u\Delta v$. This idea is encoded in the formula $${\rm d}\sigma=|{\bf r}_u\times{\bf r}_v|\ {\rm d}(u,v)$$ which leads to the well known formula for computing areas of surfaces in space. Given that at the same time the unit normal vector of the surface at the point ${\bf r}(u,v)$ is $${\bf n}(u,v)={{\bf r}_u\times{\bf r}_v\over|{\bf r}_u\times{\bf r}_v|}\ ,$$ we can say that the "vectorial surface element" encoding information not only about area, but also about the spacial orientation of this element, is given by $(1)$.
If $u$ and $v$ are replaced by new coordinates $x$ and $y$ via $$\psi:\qquad u=u(x,y),\quad v=v(x,y)$$ then we obtain a new parametrization $$(x,y)\mapsto{\bf r}\bigl(u(x,y),v(x,y)\bigr)\ ,$$ again denoted by ${\bf r}$. One has $${\bf r}_x={\bf r}_u u_x+{\bf r}_v v_x,\quad {\bf r}_y={\bf r}_u u_y+{\bf r}_v v_y$$ and therefore $${\bf r}_x\times{\bf r}_y=(u_xv_y-u_yv_x){\bf r}_u\times{\bf r}_v=J_\psi(x,y) \>{\bf r}_u\times{\bf r}_v\>\ .$$ When we are integrating such things over corresponding domains this amounts to $${\bf r}_x\times{\bf r}_y\>{\rm d}(x,y)={\bf r}_u\times{\bf r}_v\>{\rm d}(u,v)\ ,$$ as it should be.