[proofplan]
We compare the two top-degree forms $(\operatorname{div} X)\,\omega_g$ and $d(i(X)\omega_g)$ by evaluating both on a frame at an arbitrary point $p \in M$. The Leibniz rule for $\nabla$ applied to the contraction $i(X)$, combined with metric compatibility $\nabla\omega_g = 0$, reduces $\nabla_Y(i(X)\omega_g)$ to $i(\nabla_Y X)\omega_g$. Choosing a normal orthonormal frame at $p$ converts the exterior derivative $d$ into an antisymmetrisation of $\nabla$ at $p$, and the trace structure of the resulting expression collapses it to $\operatorname{tr}(\nabla X)\,\omega_g = (\operatorname{div} X)\,\omega_g$. Both sides are tensorial in $X$, so the pointwise identity at every $p$ gives the global identity.
[/proofplan]
[step:Reduce to a pointwise computation in a normal orthonormal frame]
Let $(M, g)$ be an oriented Riemannian manifold, $\nabla$ the Levi-Civita connection, $\omega_g$ the Riemannian volume form, $X \in \mathfrak{X}(M)$, and $i(X) : \Omega^k(M) \to \Omega^{k-1}(M)$ the interior product (contraction with $X$ in the first slot). Both sides of the identity
\begin{align*}
(\operatorname{div} X)\,\omega_g = d(i(X)\,\omega_g)
\tag{$\star$}
\end{align*}
are top-degree forms on $M$, and both depend pointwise on $X$ and its first covariant derivative. Hence ($\star$) is an equality of tensorial expressions, and it suffices to verify it at an arbitrary point $p \in M$.
Fix $p \in M$ and choose a **normal orthonormal frame** $\{e_1, \ldots, e_n\}$ at $p$: a local orthonormal frame on a neighbourhood of $p$ with $\nabla_Y e_i|_p = 0$ for all $Y \in T_p M$ and all $i = 1, \ldots, n$. Such a frame exists; one can construct it by parallel-transporting an orthonormal basis of $T_p M$ along radial geodesics in a normal coordinate chart centred at $p$. Let $\{\omega^1, \ldots, \omega^n\}$ be the dual coframe, so $\omega^i(e_j) = \delta_{ij}$. In this frame at $p$,
\begin{align*}
\omega_g|_p = \omega^1 \wedge \cdots \wedge \omega^n|_p.
\end{align*}
[/step]
[step:Compute $\nabla_Y(i(X)\omega_g)$ via the Leibniz rule and metric compatibility]
The covariant derivative $\nabla$ extends to all tensor bundles as a derivation that commutes with contractions. Applied to the contraction $i(X)\psi$ of a vector field $X$ with a $k$-form $\psi$, this gives the Leibniz identity
\begin{align*}
\nabla_Y(i(X)\psi) = i(\nabla_Y X)\psi + i(X)\nabla_Y\psi \qquad \text{for all } Y \in \mathfrak{X}(M).
\tag{L}
\end{align*}
Specialise (L) to $\psi := \omega_g$. The Levi-Civita connection is metric-compatible, and the Riemannian volume form is determined by $g$ together with the orientation; consequently $\nabla\omega_g = 0$, i.e., $\nabla_Y \omega_g = 0$ for all $Y \in \mathfrak{X}(M)$. The second term in (L) drops, leaving
\begin{align*}
\nabla_Y(i(X)\omega_g) = i(\nabla_Y X)\,\omega_g.
\tag{L$'$}
\end{align*}
[/step]
[step:Express $d$ as an antisymmetrisation of $\nabla$ at $p$]
For any $(n-1)$-form $\eta \in \Omega^{n-1}(M)$, the exterior derivative on a torsion-free connection satisfies the global identity
\begin{align*}
d\eta(Y_0, Y_1, \ldots, Y_{n-1}) = \sum_{j=0}^{n-1} (-1)^j (\nabla_{Y_j}\eta)(Y_0, \ldots, \widehat{Y_j}, \ldots, Y_{n-1}),
\end{align*}
where $\widehat{\cdot}$ denotes omission. This is the standard formula expressing $d$ as the antisymmetrisation of $\nabla$ (which holds because $\nabla$ is torsion-free).
Apply this to $\eta := i(X)\omega_g$ at $p$ with $(Y_0, \ldots, Y_{n-1}) := (e_1, \ldots, e_n)$:
\begin{align*}
d(i(X)\omega_g)(e_1, \ldots, e_n)|_p = \sum_{j=1}^{n} (-1)^{j-1} (\nabla_{e_j}(i(X)\omega_g))(e_1, \ldots, \widehat{e_j}, \ldots, e_n)|_p.
\end{align*}
By (L$'$) at $p$, $\nabla_{e_j}(i(X)\omega_g)|_p = i(\nabla_{e_j}X)\omega_g|_p$, so
\begin{align*}
d(i(X)\omega_g)(e_1, \ldots, e_n)|_p = \sum_{j=1}^{n} (-1)^{j-1} \big(i(\nabla_{e_j}X)\omega_g\big)(e_1, \ldots, \widehat{e_j}, \ldots, e_n)|_p.
\tag{D}
\end{align*}
[/step]
[step:Evaluate the contraction sum and identify the trace]
Write $\nabla_{e_j} X|_p = \sum_{k=1}^n A_{kj}\, e_k|_p$ in the frame at $p$, so that the matrix of $\nabla X|_p : T_p M \to T_p M$ in the basis $\{e_k\}$ is $(A_{kj})$. By definition of $\operatorname{div}$,
\begin{align*}
(\operatorname{div} X)(p) = \operatorname{tr}(\nabla X|_p) = \sum_{k=1}^n A_{kk}.
\end{align*}
Compute one term of the sum (D). For fixed $j$,
\begin{align*}
\big(i(\nabla_{e_j}X)\omega_g\big)(e_1, \ldots, \widehat{e_j}, \ldots, e_n)|_p &= \omega_g(\nabla_{e_j}X, e_1, \ldots, \widehat{e_j}, \ldots, e_n)|_p \\
&= \sum_{k=1}^n A_{kj}\, \omega_g(e_k, e_1, \ldots, \widehat{e_j}, \ldots, e_n)|_p.
\end{align*}
The form $\omega_g|_p = \omega^1 \wedge \cdots \wedge \omega^n|_p$ evaluated on the tuple $(e_k, e_1, \ldots, \widehat{e_j}, \ldots, e_n)$ vanishes unless $\{k, 1, \ldots, \widehat{j}, \ldots, n\} = \{1, \ldots, n\}$, i.e., unless $k = j$. When $k = j$, the tuple $(e_j, e_1, \ldots, \widehat{e_j}, \ldots, e_n)$ is the cyclic shift bringing $e_j$ to the front of $(e_1, \ldots, e_n)$, achieved by $j-1$ transpositions, so
\begin{align*}
\omega_g(e_j, e_1, \ldots, \widehat{e_j}, \ldots, e_n)|_p = (-1)^{j-1} \omega_g(e_1, \ldots, e_n)|_p = (-1)^{j-1}.
\end{align*}
Therefore
\begin{align*}
\big(i(\nabla_{e_j}X)\omega_g\big)(e_1, \ldots, \widehat{e_j}, \ldots, e_n)|_p = (-1)^{j-1} A_{jj}.
\end{align*}
Substituting into (D),
\begin{align*}
d(i(X)\omega_g)(e_1, \ldots, e_n)|_p = \sum_{j=1}^n (-1)^{j-1} \cdot (-1)^{j-1} A_{jj} = \sum_{j=1}^n A_{jj} = (\operatorname{div} X)(p).
\end{align*}
On the other hand,
\begin{align*}
\big((\operatorname{div} X)\,\omega_g\big)(e_1, \ldots, e_n)|_p = (\operatorname{div} X)(p) \cdot \omega_g(e_1, \ldots, e_n)|_p = (\operatorname{div} X)(p).
\end{align*}
The two top-degree forms agree on the frame $(e_1, \ldots, e_n)$ at $p$. Since both are alternating $n$-forms on the $n$-dimensional space $T_p M$, agreement on one ordered basis forces agreement as forms at $p$:
\begin{align*}
\big(d(i(X)\omega_g)\big)|_p = \big((\operatorname{div} X)\,\omega_g\big)|_p.
\end{align*}
The point $p \in M$ was arbitrary, so this is the asserted global identity ($\star$).
[guided]
We have arrived at formula (D),
\begin{align*}
d(i(X)\omega_g)(e_1, \ldots, e_n)|_p = \sum_{j=1}^{n} (-1)^{j-1} \big(i(\nabla_{e_j}X)\omega_g\big)(e_1, \ldots, \widehat{e_j}, \ldots, e_n)|_p,
\end{align*}
and we want to evaluate this and recognise the right-hand side as $(\operatorname{div} X)(p)$. To do so we must extract numerical content from each $j$-th summand. The natural strategy is to expand $\nabla_{e_j} X$ in the orthonormal frame, plug into the volume form, and watch the alternating signs collide with the signs coming from the cyclic structure of $\omega_g$.
**Expand $\nabla_{e_j} X$ in the frame.** The vector $\nabla_{e_j} X|_p \in T_p M$ has coordinates in the basis $\{e_k|_p\}$:
\begin{align*}
\nabla_{e_j} X|_p = \sum_{k=1}^n A_{kj}\, e_k|_p,
\end{align*}
where the matrix $(A_{kj})$ is precisely the matrix of the endomorphism $\nabla X|_p : T_p M \to T_p M$ in the basis $\{e_k\}$. Why does this help? Because the divergence is by definition the trace of this endomorphism:
\begin{align*}
(\operatorname{div} X)(p) = \operatorname{tr}(\nabla X|_p) = \sum_{k=1}^n A_{kk}.
\end{align*}
So our goal is to show that the right-hand side of (D) collapses to $\sum_k A_{kk}$ — only diagonal entries should survive.
**Compute one summand of (D).** Fix $j$. Using the definition $i(V)\omega_g(\cdot) = \omega_g(V, \cdot)$ and the linearity of $\omega_g$ in its first argument,
\begin{align*}
\big(i(\nabla_{e_j}X)\omega_g\big)(e_1, \ldots, \widehat{e_j}, \ldots, e_n)|_p
&= \omega_g(\nabla_{e_j}X, e_1, \ldots, \widehat{e_j}, \ldots, e_n)|_p \\
&= \sum_{k=1}^n A_{kj}\, \omega_g(e_k, e_1, \ldots, \widehat{e_j}, \ldots, e_n)|_p.
\end{align*}
Now we ask: which $k$ contribute? At $p$ the form $\omega_g|_p = \omega^1 \wedge \cdots \wedge \omega^n|_p$ is alternating, so it vanishes on any tuple with a repeated entry. The tuple $(e_k, e_1, \ldots, \widehat{e_j}, \ldots, e_n)$ already contains every $e_i$ with $i \neq j$. If $k \neq j$, then $e_k$ duplicates one of these, and $\omega_g$ vanishes. The only surviving term is $k = j$ — which is exactly the diagonal entry $A_{jj}$. This is the mechanism by which off-diagonal entries are killed.
**Compute the sign for $k = j$.** When $k = j$, we evaluate $\omega_g$ on $(e_j, e_1, \ldots, \widehat{e_j}, \ldots, e_n)$. To compare this with the standard ordering $(e_1, \ldots, e_n)$ on which $\omega_g|_p$ takes the value $+1$, count the transpositions needed. The tuple is obtained from $(e_1, \ldots, e_n)$ by extracting $e_j$ from its position and pulling it to the front — that is, by $j-1$ adjacent transpositions. Since $\omega_g$ is alternating, each transposition changes the sign:
\begin{align*}
\omega_g(e_j, e_1, \ldots, \widehat{e_j}, \ldots, e_n)|_p = (-1)^{j-1} \omega_g(e_1, \ldots, e_n)|_p = (-1)^{j-1}.
\end{align*}
Combining,
\begin{align*}
\big(i(\nabla_{e_j}X)\omega_g\big)(e_1, \ldots, \widehat{e_j}, \ldots, e_n)|_p = (-1)^{j-1} A_{jj}.
\end{align*}
**Watch the signs cancel.** Substitute back into (D). The summand carries two sign factors of $(-1)^{j-1}$: one from the antisymmetrisation in the $d$-formula, and one from the cyclic permutation we just computed. They square to $+1$:
\begin{align*}
d(i(X)\omega_g)(e_1, \ldots, e_n)|_p
= \sum_{j=1}^n (-1)^{j-1} \cdot (-1)^{j-1} A_{jj}
= \sum_{j=1}^n A_{jj}
= (\operatorname{div} X)(p).
\end{align*}
This is the heart of the proof: the alternating signs from $d$ and the alternating signs from $\omega_g$ exactly compensate, leaving the bare trace.
**Compare with the left-hand side.** The form $(\operatorname{div} X)\,\omega_g$ is, by definition, the function $\operatorname{div} X$ multiplied by the volume form. Evaluating on the frame at $p$,
\begin{align*}
\big((\operatorname{div} X)\,\omega_g\big)(e_1, \ldots, e_n)|_p
= (\operatorname{div} X)(p) \cdot \omega_g(e_1, \ldots, e_n)|_p
= (\operatorname{div} X)(p),
\end{align*}
using $\omega_g(e_1, \ldots, e_n)|_p = 1$ since $\{e_i\}$ is a positively oriented orthonormal basis at $p$.
**Conclude pointwise then globally.** We have shown
\begin{align*}
d(i(X)\omega_g)(e_1, \ldots, e_n)|_p = (\operatorname{div} X)(p) = \big((\operatorname{div} X)\,\omega_g\big)(e_1, \ldots, e_n)|_p.
\end{align*}
Why does agreement on a single ordered basis suffice? Both sides are alternating $n$-forms on the $n$-dimensional vector space $T_p M$, and the space of alternating $n$-forms on an $n$-dimensional space is one-dimensional. An alternating $n$-form is determined by its value on any single ordered basis, so equality on $(e_1, \ldots, e_n)$ forces
\begin{align*}
\big(d(i(X)\omega_g)\big)|_p = \big((\operatorname{div} X)\,\omega_g\big)|_p.
\end{align*}
Finally, the point $p \in M$ was arbitrary. As discussed in the first step, both sides of ($\star$) are tensorial in $X$ and depend only on $X$ and $\nabla X$ pointwise, so a pointwise identity at every $p$ is exactly the global identity ($\star$). The proof is complete.
Two retrospective remarks on the strategy. First, the choice of a normal orthonormal frame served a single purpose: it made $\nabla_{e_j} e_k|_p = 0$ so that the global Koszul-type formula for $d$ reduced at $p$ to the naive antisymmetrisation of $\nabla$. Without the normal-frame condition, (D) would have carried extra terms of the form $\eta(\ldots, \nabla_{e_j} e_k, \ldots)$ that complicate the cyclic-shift bookkeeping. Second, the entire computation is a "Cramer-style" extraction of a trace: the contraction $i(X)$ followed by $d$ produces, via the alternating structure of $\omega_g$, exactly the trace of $\nabla X$ — the same trace that defines $\operatorname{div} X$. This is why the divergence and $d \circ i(X)$ on the volume form are the same object.
[/guided]
[/step]
