[proofplan] We proceed in three stages. First, in any coordinate chart $(U, \varphi)$ with coordinates $(x_1, \dots, x_n)$, we define $d$ by the coordinate formula $d(\sum_I \omega_I \, dx_I) = \sum_I d\omega_I \wedge dx_I$, where $d\omega_I$ is the differential of the function $\omega_I$, and verify that properties (1)–(3) hold locally. Second, we prove uniqueness: any operator $\tilde{d}$ satisfying (1)–(3) and agreeing with $d$ on functions must satisfy $\tilde{d}(dx_I) = 0$, and hence $\tilde{d}(f \, dx_I) = df \wedge dx_I$, forcing $\tilde{d} = d$ locally. Third, we verify coordinate independence by showing that uniqueness forces the coordinate formulas in overlapping charts to agree, so the local definitions glue to a well-defined global operator. Naturality (property 4) is checked by induction on degree. [/proofplan]

[step:Define $d$ in local coordinates and verify properties (1)–(3)]Let $(U, \varphi)$ be a coordinate chart on $M$ with local coordinates $(x_1, \dots, x_n)$. Every $k$-form $\omega \in \Omega^k(U)$ writes uniquely as \begin{align*} \omega = \sum_I \omega_I \, dx_I, \end{align*} where the sum is over strictly increasing multi-indices $I = (i_1 < \cdots < i_k)$ and $\omega_I \in C^\infty(U)$. Define \begin{align*} d_U\omega := \sum_I d\omega_I \wedge dx_I = \sum_I \sum_{m=1}^n \frac{\partial \omega_I}{\partial x_m} \, dx_m \wedge dx_I \in \Omega^{k+1}(U). \end{align*} **Property (1) — $\mathbb{R}$-linearity.** For $\omega, \eta \in \Omega^k(U)$ and $\lambda \in \mathbb{R}$, \begin{align*} d_U(\omega + \lambda\eta) = \sum_I d(\omega_I + \lambda\eta_I) \wedge dx_I = \sum_I (d\omega_I + \lambda \, d\eta_I) \wedge dx_I = d_U\omega + \lambda \, d_U\eta, \end{align*} where we used $\mathbb{R}$-linearity of $d$ on functions and of $\wedge$. **Property (2) — Leibniz rule.** For $\omega = f \, dx_I \in \Omega^k(U)$ and $\eta = g \, dx_J \in \Omega^l(U)$, \begin{align*} d_U(\omega \wedge \eta) = d_U(fg \, dx_I \wedge dx_J) = d(fg) \wedge dx_I \wedge dx_J. \end{align*} Since $d(fg) = f \, dg + g \, df$ (product rule on functions), this equals \begin{align*} (g \, df + f \, dg) \wedge dx_I \wedge dx_J = g \, df \wedge dx_I \wedge dx_J + f \, dg \wedge dx_I \wedge dx_J. \end{align*} Now $d_U\omega \wedge \eta = df \wedge dx_I \wedge g \, dx_J = g \, df \wedge dx_I \wedge dx_J$, and \begin{align*} (-1)^k \omega \wedge d_U\eta = (-1)^k f \, dx_I \wedge dg \wedge dx_J = f \, dg \wedge dx_I \wedge dx_J, \end{align*} where we moved $dg$ past $dx_I$ (a $k$-form), acquiring a sign $(-1)^k = (-1)^{|\omega|}$. Thus $d_U(\omega \wedge \eta) = d_U\omega \wedge \eta + (-1)^{|\omega|} \omega \wedge d_U\eta$, as required. **Property (3) — $d^2 = 0$.** For $\omega = \omega_I \, dx_I$, \begin{align*} d_U(d_U\omega) = d_U\!\left(\sum_m \frac{\partial \omega_I}{\partial x_m} \, dx_m \wedge dx_I\right) = \sum_{m,j} \frac{\partial^2 \omega_I}{\partial x_j \partial x_m} \, dx_j \wedge dx_m \wedge dx_I. \end{align*} The coefficient $\frac{\partial^2 \omega_I}{\partial x_j \partial x_m}$ is symmetric in $(j,m)$ by equality of mixed partial derivatives (since $\omega_I \in C^\infty$), while $dx_j \wedge dx_m$ is antisymmetric. Pairing each $(j,m)$ with $(m,j)$, the contributions cancel in pairs, giving $d_U^2\omega = 0$.[/step]

[guided]We want to define $d$ so that it satisfies properties (1)–(3) and extends the differential of functions. The coordinate formula $d(\sum_I \omega_I \, dx_I) = \sum_I d\omega_I \wedge dx_I$ is the canonical definition: it says "differentiate the coefficient functions, then wedge with the existing basis form." For property (2), the key is the product rule on coefficient functions: $d(fg) = f \, dg + g \, df$. When we pass $dg$ from the right of $dx_I$ to the left, we pick up a sign $(-1)^k$ because we are commuting a $1$-form past a $k$-form in the graded-commutative algebra $\Omega^*(M)$. For property (3), the cancellation $d^2 = 0$ is a direct consequence of two facts acting against each other: (i) mixed partials commute (symmetry of $\frac{\partial^2 \omega_I}{\partial x_j \partial x_m}$ in $j,m$), and (ii) the wedge product of two copies of the same $1$-form vanishes ($dx_j \wedge dx_m + dx_m \wedge dx_j = 0$). These two symmetry/antisymmetry conditions are exactly opposed, and their product sums to zero. This is the algebraic reason $d^2 = 0$ holds, and it is why $d$ is related to the boundary operator in algebraic topology.[/guided]

[step:Prove uniqueness of $d$ in local coordinates]Let $\tilde{d} : \Omega^*(U) \to \Omega^{*+1}(U)$ be any $\mathbb{R}$-linear operator satisfying properties (1)–(3) and agreeing with $d_U$ on $\Omega^0(U) = C^\infty(U)$. **Step 1: $\tilde{d}(dx_i) = 0$ for each coordinate function $x_i$.** Since $\tilde{d}$ agrees with $d_U$ on functions, $\tilde{d}(x_i) = dx_i$. Applying $\tilde{d}$ again and using $\tilde{d}^2 = 0$: \begin{align*} 0 = \tilde{d}(\tilde{d}(x_i)) = \tilde{d}(dx_i). \end{align*} **Step 2: $\tilde{d}(dx_I) = 0$ for every multi-index $I$.** By induction on $|I|$. The base case $|I| = 1$ is Step 1. For $|I| \ge 2$, write $dx_I = dx_{i_1} \wedge dx_{I'}$ where $I' = (i_2 < \cdots < i_l)$. Applying the Leibniz rule for $\tilde{d}$ and using the inductive hypothesis $\tilde{d}(dx_{i_1}) = 0$ and $\tilde{d}(dx_{I'}) = 0$: \begin{align*} \tilde{d}(dx_I) = \tilde{d}(dx_{i_1}) \wedge dx_{I'} + (-1)^1 dx_{i_1} \wedge \tilde{d}(dx_{I'}) = 0 \wedge dx_{I'} + (-1) dx_{i_1} \wedge 0 = 0. \end{align*} **Step 3: $\tilde{d}(f \, dx_I) = d_U(f \, dx_I)$ for all $f \in C^\infty(U)$.** Using the Leibniz rule for $\tilde{d}$ and Steps 1–2: \begin{align*} \tilde{d}(f \, dx_I) = \tilde{d}(f) \wedge dx_I + (-1)^0 f \wedge \tilde{d}(dx_I) = df \wedge dx_I + 0 = df \wedge dx_I = d_U(f \, dx_I). \end{align*} By $\mathbb{R}$-linearity, $\tilde{d}(\omega) = d_U(\omega)$ for all $\omega \in \Omega^k(U)$. Thus $d_U$ is the unique operator satisfying the stated properties in the chart $(U, \varphi)$.[/step]

[guided]The uniqueness argument is a clean cascade: $d^2 = 0$ pins down what happens to coordinate $1$-forms, and the Leibniz rule then pins down what happens to all forms. Why does $d^2 = 0$ kill $d(dx_i)$? Because $\tilde{d}$ agrees with $d$ on functions, so $\tilde{d}(x_i) = dx_i$. Applying $\tilde{d}$ again gives $\tilde{d}(dx_i) = \tilde{d}(\tilde{d}(x_i)) = \tilde{d}^2(x_i) = 0$ by property (3). No other input was needed — just $d^2 = 0$ and the agreement on functions. The inductive step for multi-indices then uses only the Leibniz rule. Once we know $\tilde{d}(dx_i) = 0$ for each $i$, the Leibniz rule forces $\tilde{d}(dx_I) = 0$ for each multi-index $I$, since a wedge product differentiates term by term and each term vanishes. Finally, $\tilde{d}(f \, dx_I) = df \wedge dx_I$ is forced by Leibniz applied to $f$ times the "constant" form $dx_I$: the second term $f \cdot \tilde{d}(dx_I) = 0$ drops out, and the first term $\tilde{d}(f) \wedge dx_I = df \wedge dx_I$ is the coordinate formula. So uniqueness follows from a two-step cascade: $d^2 = 0$ kills basis forms, then Leibniz reconstructs $\tilde{d}$ from $d$ on functions.[/guided]

[step:Verify coordinate independence so local definitions glue to a global $d$]Let $(U, \varphi)$ and $(V, \psi)$ be overlapping charts with coordinates $(x_1, \dots, x_n)$ and $(y_1, \dots, y_n)$ respectively. On $U \cap V$, each $y_j$ is a smooth function of $(x_1, \dots, x_n)$, and the two local operators $d_U$ and $d_V$ are defined by their respective coordinate formulas. We claim $d_U$ and $d_V$ agree on $\Omega^k(U \cap V)$. Both $d_U$ and $d_V$, when restricted to $U \cap V$, satisfy properties (1)–(3) and both agree with the differential of functions on $\Omega^0(U \cap V)$ (since the differential of a smooth function is coordinate-independent: if $f \in C^\infty(U \cap V)$ then $df = \sum_i \frac{\partial f}{\partial x_i} dx_i = \sum_j \frac{\partial f}{\partial y_j} dy_j$ by the chain rule). The uniqueness result of the previous step therefore forces \begin{align*} d_U\big|_{U \cap V} = d_V\big|_{U \cap V}. \end{align*} The local operators are thus consistent on all pairwise overlaps, and by partition of unity they assemble into a globally defined operator $d : \Omega^k(M) \to \Omega^{k+1}(M)$. In every chart the global $d$ coincides with the local coordinate formula, so properties (1)–(3) hold globally.[/step]

[guided]This step addresses a potential gap: we defined $d$ in each chart separately. For $d$ to be a well-defined global operator on $M$, we need the local definitions to agree on chart overlaps. The key observation is that the differential of a function is coordinate-independent. If $f \in C^\infty(U \cap V)$, then regardless of whether we use $x$-coordinates or $y$-coordinates, \begin{align*} df = \sum_{i=1}^n \frac{\partial f}{\partial x_i} \, dx_i = \sum_{j=1}^n \frac{\partial f}{\partial y_j} \, dy_j, \end{align*} as both expressions represent the same intrinsic object (the total derivative of $f$). This means both $d_U$ and $d_V$ agree on $\Omega^0(U \cap V)$. Since both $d_U|_{U \cap V}$ and $d_V|_{U \cap V}$ satisfy properties (1)–(3) and agree on functions, uniqueness (which we proved for any such operator in a single chart) forces them to be identical on $U \cap V$. Having verified consistency on all pairwise intersections, we may choose a partition of unity $\{\rho_\alpha\}$ subordinate to the cover $\{U_\alpha\}$ and define $d\omega = \sum_\alpha d_{U_\alpha}(\rho_\alpha \omega)$. Because any two local formulas agree on intersections, this sum is independent of the choice of partition of unity and defines a global operator.[/guided]

[step:Verify naturality: $d \circ f^* = f^* \circ d$ for smooth $f : M \to N$] Let $f : M \to N$ be a smooth map. We verify property (4) by induction on the degree $k$ of $\omega \in \Omega^k(N)$. **Base case $k = 0$:** For $\phi \in C^\infty(N) = \Omega^0(N)$, $f^*\phi = \phi \circ f \in C^\infty(M)$. For any vector field $X \in \Gamma(TM)$ and any $p \in M$, \begin{align*} (d(f^*\phi))_p(X_p) = X_p(\phi \circ f) = (df_p)^*(d\phi_{f(p)})(X_p) = (f^*(d\phi))_p(X_p), \end{align*} where $(df_p)^* : T^*_{f(p)}N \to T^*_p M$ is the pullback of the differential. Thus $d \circ f^* = f^* \circ d$ on $\Omega^0(N)$. **Inductive step:** Suppose naturality holds for all forms of degree less than $k$. Any $k$-form $\omega \in \Omega^k(N)$ is locally a sum of terms of the form $g \, dh_1 \wedge \cdots \wedge dh_k$ with $g, h_1, \dots, h_k \in C^\infty(N)$. By linearity it suffices to consider $\omega = g \, dh_1 \wedge \cdots \wedge dh_k$. Then \begin{align*} f^*\omega = (f^*g)(f^*dh_1) \wedge \cdots \wedge (f^*dh_k) = (g \circ f) \, d(h_1 \circ f) \wedge \cdots \wedge d(h_k \circ f), \end{align*} using the base case $f^*(dh_j) = d(f^*h_j) = d(h_j \circ f)$ in the last step. Applying $d$ to $f^*\omega$ and using the Leibniz rule: \begin{align*} d(f^*\omega) &= d(g \circ f) \wedge d(h_1 \circ f) \wedge \cdots \wedge d(h_k \circ f). \end{align*} On the other hand, $d\omega = dg \wedge dh_1 \wedge \cdots \wedge dh_k$, and \begin{align*} f^*(d\omega) = f^*(dg) \wedge f^*(dh_1) \wedge \cdots \wedge f^*(dh_k) = d(g \circ f) \wedge d(h_1 \circ f) \wedge \cdots \wedge d(h_k \circ f), \end{align*} again using the base case for each factor. The two expressions are equal, so $d(f^*\omega) = f^*(d\omega)$, completing the induction. [/step]