Proof of Pullback by a Smooth Map Commutes with the Exterior Derivative

Step 0 (symmetry of mixed partials). Let $W\subseteq\mathbb{R}^m$ be open, $g:W\to\mathbb{R}$ a smooth map, and $i,j$ indices. For $a\in W$ and small $h,k>0$ consider

where $e_i,e_j$ are standard basis vectors. Applying the mean value theorem first in the $e_i$-variable and then in the $e_j$-variable gives $\Delta(h,k)=hk\,\partial_j\partial_i g(\xi)$ for some $\xi$ in the rectangle spanned by the four points; applying it in the other order gives $\Delta(h,k)=hk\,\partial_i\partial_j g(\eta)$ similarly. Letting $(h,k)\to(0,0)$ and using continuity of the second partial derivatives yields $\partial_j\partial_i g=\partial_i\partial_j g$ on $W$.

Step 1 (pullback of a decomposable form). A straightforward induction from Wedge Product of Differential Forms on Euclidean Space and Associativity of the Wedge Product of Differential Forms on Euclidean Space, using Sign of a Permutation, shows that for $1$-forms $\alpha_1,\dots,\alpha_k$ and vectors $v_1,\dots,v_k$,

with the determinant of the $k\times k$ matrix indicated. Write $F=(F_1,\dots,F_m)$ in components; each $F_i$ is a smooth $0$-form on $U$ and its exterior derivative is $dF_i=\sum_{l=1}^{n}\frac{\partial F_i}{\partial x_l}dx_l$. For a $C^1$ function $a:V\to\mathbb{R}$ and increasing indices $i_1<\cdots<i_k$, I claim

Indeed, by Pullback of a Differential Form by a C^1 Map the left side at $x$ on $(v_1,\dots,v_k)$ equals $a(F(x))$ times the value of $dy_{i_1}\wedge\cdots\wedge dy_{i_k}$ on the matrix-vector products $J_F(x)v_1,\dots,J_F(x)v_k$; by the determinant identity this is $a(F(x))\det\bigl[(J_F(x)v_s)_{i_r}\bigr]$, and $(J_F(x)v_s)_{i_r}=\sum_l \frac{\partial F_{i_r}}{\partial x_l}(x)(v_s)_l=(dF_{i_r})_x(v_s)$, which is the right side by the same determinant identity.

By Coordinate Expansion of Differential Forms on Euclidean Open Sets, $\omega=\sum_I a_I\,dy_{i_1}\wedge\cdots\wedge dy_{i_k}$ over increasing multi-indices $I=(i_1,\dots,i_k)$, and the pullback is additive in $\omega$ directly from Pullback of a Differential Form by a C^1 Map; hence

Evaluating on basis vectors $e_{j_1},\dots,e_{j_k}$ (increasing $J$) and using the uniqueness in Coordinate Expansion of Differential Forms on Euclidean Open Sets, the coefficient of $F^{*}\omega$ on $dx_{j_1}\wedge\cdots\wedge dx_{j_k}$ is

Each $a_I\circ F$ is $C^1$ by the chain rule (using C^1 Maps on Euclidean Open Sets are Differentiable), each $M_{I,J}$ is obtained from the smooth partial derivatives of $F$ by finitely many sums and products via Determinant of a Real Square Matrix, and products and sums of $C^1$ functions are $C^1$ by Products and Quotients of C^k Real-Valued Maps on Euclidean Open Sets Are C^k; hence $F^{*}\omega$ is a $C^1$ differential $k$-form.

Step 2 (the case $k=0$). For a $C^1$ function $b:V\to\mathbb{R}$, $F^{*}(db)=d(b\circ F)$: at $x$ on a vector $v$, the left side is $\sum_j \frac{\partial b}{\partial y_j}(F(x))(J_F(x)v)_j=\sum_{j,l}\frac{\partial b}{\partial y_j}(F(x))\frac{\partial F_j}{\partial x_l}(x)v_l$, which equals $d(b\circ F)_x(v)$ by the chain rule.

Step 3 (conclusion). Compute $d(F^{*}\omega)$ from Exterior Derivative of a C^1 Differential Form on a Euclidean Open Set:

Expand $\frac{\partial c_J}{\partial x_l}$ by the product rule for partial derivatives (the one-dimensional product rule applied to each coordinate slice): the terms in which the derivative falls on $a_I\circ F$ combine, by Step 2 and the two determinant identities of Step 1, into exactly the coordinate expansion of

which is the pullback of $d\omega$ by Step 1 applied to the expansion of $d\omega$. The remaining terms carry a derivative of an entry of $M_{I,J}$, i.e. a second partial $\frac{\partial^2 F_{i_r}}{\partial x_l\partial x_{j_s}}$; summed against $dx_l\wedge dx_{j_1}\wedge\cdots\wedge dx_{j_k}$ these assemble into expressions of the form $d\bigl(dF_{i}\bigr)\wedge(\text{remaining factors})$, whose coefficients are antisymmetric in the pair of differentiation indices by Permutation Rule for Wedge Products of Coordinate 1-Forms while the second partials are symmetric in that pair by Step 0; hence each such sum vanishes. Therefore $d(F^{*}\omega)=F^{*}(d\omega)$. $\blacksquare$