[proofplan]
The strategy is to test the proposed local formula against an arbitrary $1$-form $\alpha$ via the global $L^2$ inner product on $\Omega^1(M)$, and to use the formal-adjoint property $(\delta\beta, \alpha)_{L^2} = (\beta, d\alpha)_{L^2}$ defining $\delta$. We expand the inner product on the right using the local orthonormal coframe; combine the two halves of $d\alpha$ via the antisymmetrisation formula $d\alpha(X, Y) = (\nabla_X\alpha)(Y) - (\nabla_Y\alpha)(X)$ valid for the torsion-free Levi-Civita connection; and use the antisymmetry of $\beta$ in its two slots to combine the resulting two terms into $2\sum_{k} \beta(e_k, \cdot)\,(\nabla_{e_k}\alpha)$, paired against $\alpha$ via $\beta$. A "Leibniz + divergence theorem" manoeuvre — extracting a divergence and applying the [Divergence Theorem](/theorems/2754) on the closed manifold $M$ — moves the covariant derivative off $\alpha$ and onto $\beta$, producing the candidate formula. The fundamental lemma converts the resulting integral identity for all $\alpha$ into the pointwise statement $\delta\beta = -\sum_k i(e_k)\nabla_{e_k}\beta$.
[/proofplan]
[step:Set up notation and reduce to an integral identity for the codifferential]
Let $(M, g)$ be an oriented closed Riemannian manifold of dimension $n$, $\nabla$ the Levi-Civita connection, $\omega_g$ the Riemannian volume form, $\beta \in \Omega^2(M)$ a smooth $2$-form, and $\{e_k\}_{k=1}^n$ a local orthonormal frame field on an open set $U \subseteq M$ with dual coframe $\{\omega^k\}_{k=1}^n$ (so $\omega^k(e_l) = \delta^k_l$ and $g(e_k, e_l) = \delta_{kl}$ on $U$).
The codifferential $\delta : \Omega^2(M) \to \Omega^1(M)$ is the formal $L^2$-adjoint of the exterior derivative $d : \Omega^1(M) \to \Omega^2(M)$. By the [Co-differential is Formal Adjoint of $d$](/theorems/2742) identity applied to $\alpha \in \Omega^1(M)$ and $\beta \in \Omega^2(M)$ on the closed manifold $M$,
\begin{align*}
\int_M g(\delta\beta, \alpha)\,\omega_g = \int_M g(\beta, d\alpha)\,\omega_g \qquad \text{for every } \alpha \in \Omega^1(M).
\tag{A}
\end{align*}
Here $g(\cdot, \cdot)$ on $\Omega^1$ and $\Omega^2$ denotes the metric-induced inner product, given on $\Omega^1$ in an orthonormal coframe by $g(\xi, \eta) = \sum_k \xi(e_k)\eta(e_k)$, and on $\Omega^2$ by
\begin{align*}
g(\beta, \gamma) = \tfrac{1}{2}\sum_{k, l = 1}^n \beta(e_k, e_l)\,\gamma(e_k, e_l).
\end{align*}
(The factor $1/2$ accounts for the unordered pair convention; equivalently, summing over $k < l$ removes it.)
We will prove that
\begin{align*}
\int_M g(\delta\beta, \alpha)\,\omega_g = -\int_M \sum_{k=1}^n g\big( i(e_k)(\nabla_{e_k}\beta),\, \alpha \big)\,\omega_g \qquad \text{for every } \alpha \in \Omega^1(M).
\tag{T}
\end{align*}
By the fundamental lemma of the calculus of variations on the closed Riemannian manifold $M$, (T) — being valid for every smooth $\alpha$ — implies the pointwise identity
\begin{align*}
\delta\beta = -\sum_{k=1}^n i(e_k)(\nabla_{e_k}\beta) \qquad \text{on } U.
\end{align*}
Since the right-hand side, computed using a local orthonormal frame, is invariant under change of orthonormal frame (it is the metric trace of $\nabla\beta$ in the first slot, which is frame-independent), the formula extends globally and we may assume the frame is defined on all of $M$ for the purposes of the integral identity (using a partition of unity if necessary, or working in the open set $U$ paired with bump functions in $\alpha$).
The remaining task is to prove (T).
[/step]
[step:Expand $g(\beta, d\alpha)$ using $d\alpha(X, Y) = (\nabla_X\alpha)(Y) - (\nabla_Y\alpha)(X)$]
Fix $\alpha \in \Omega^1(M)$. The exterior derivative on a torsion-free connection satisfies the global identity
\begin{align*}
d\alpha(X, Y) = (\nabla_X \alpha)(Y) - (\nabla_Y \alpha)(X) \qquad \text{for all } X, Y \in \mathfrak{X}(M),
\end{align*}
which holds because $\nabla$ is torsion-free. Expand $g(\beta, d\alpha)$ in the orthonormal frame $\{e_k\}$ on $U$:
\begin{align*}
g(\beta, d\alpha) = \tfrac{1}{2}\sum_{k, l = 1}^n \beta(e_k, e_l)\,d\alpha(e_k, e_l) = \tfrac{1}{2}\sum_{k, l} \beta(e_k, e_l)\big[(\nabla_{e_k}\alpha)(e_l) - (\nabla_{e_l}\alpha)(e_k)\big].
\end{align*}
The two terms inside the bracket combine by antisymmetry of $\beta$. Reindex the second term by swapping $k \leftrightarrow l$:
\begin{align*}
\sum_{k, l} \beta(e_k, e_l)(\nabla_{e_l}\alpha)(e_k) = \sum_{k, l} \beta(e_l, e_k)(\nabla_{e_k}\alpha)(e_l) = -\sum_{k, l} \beta(e_k, e_l)(\nabla_{e_k}\alpha)(e_l),
\end{align*}
using $\beta(e_l, e_k) = -\beta(e_k, e_l)$. Therefore
\begin{align*}
g(\beta, d\alpha) = \tfrac{1}{2}\sum_{k, l}\beta(e_k, e_l)(\nabla_{e_k}\alpha)(e_l) - \tfrac{1}{2}\big(-\sum_{k, l}\beta(e_k, e_l)(\nabla_{e_k}\alpha)(e_l)\big) = \sum_{k, l}\beta(e_k, e_l)(\nabla_{e_k}\alpha)(e_l).
\end{align*}
Recognise the inner sum over $l$ as a contraction with the interior product. For fixed $k$, the $1$-form $i(e_k)\beta \in \Omega^1(M)$ is defined by $(i(e_k)\beta)(V) = \beta(e_k, V)$, and the metric inner product on $\Omega^1$ in the orthonormal frame is $g(\xi, \eta) = \sum_l \xi(e_l)\eta(e_l)$. Hence
\begin{align*}
\sum_{l=1}^n \beta(e_k, e_l)(\nabla_{e_k}\alpha)(e_l) = g\big(i(e_k)\beta,\; \nabla_{e_k}\alpha\big),
\end{align*}
and summing over $k$,
\begin{align*}
g(\beta, d\alpha) = \sum_{k=1}^n g\big( i(e_k)\beta,\; \nabla_{e_k}\alpha\big).
\tag{E}
\end{align*}
[/step]
[step:Move the covariant derivative off $\alpha$ via a divergence-on-a-divergence-free $1$-form trick]
For each $k = 1, \ldots, n$, the term $g(i(e_k)\beta,\, \nabla_{e_k}\alpha)$ on $U$ has the form "inner product of a $1$-form with the covariant derivative of $\alpha$". To integrate by parts, we apply the [Adjoint Formula for Covariant Derivative](/theorems/2758) — but here we work it out concretely for the $1$-form case.
Define the smooth vector field $Y_k \in \mathfrak{X}(M)$ to be the metric dual of the $1$-form $\alpha\cdot g(i(e_k)\beta, \cdot)$ — but more cleanly, we proceed by a direct Leibniz computation. Specifically: introduce the smooth function
\begin{align*}
f_k := g(i(e_k)\beta, \alpha) \in C^\infty(M).
\end{align*}
Differentiate $f_k$ in the direction $e_k$ using metric compatibility:
\begin{align*}
e_k(f_k) = g\big(\nabla_{e_k}(i(e_k)\beta),\, \alpha\big) + g\big( i(e_k)\beta,\, \nabla_{e_k}\alpha\big).
\end{align*}
Solving for the second term,
\begin{align*}
g\big(i(e_k)\beta,\, \nabla_{e_k}\alpha\big) = e_k(f_k) - g\big(\nabla_{e_k}(i(e_k)\beta),\, \alpha\big).
\end{align*}
The covariant derivative commutes with interior product up to a Leibniz rule:
\begin{align*}
\nabla_{e_k}(i(e_k)\beta) = i(\nabla_{e_k} e_k)\beta + i(e_k)(\nabla_{e_k}\beta).
\end{align*}
Substituting,
\begin{align*}
g\big(i(e_k)\beta, \nabla_{e_k}\alpha\big) = e_k(f_k) - g\big(i(\nabla_{e_k} e_k)\beta, \alpha\big) - g\big(i(e_k)(\nabla_{e_k}\beta), \alpha\big).
\tag{IBP}
\end{align*}
Sum (IBP) over $k = 1, \ldots, n$ and integrate against $\omega_g$ on $M$. We will show that the first two contributions on the right collapse into a single divergence term.
Define the vector field $Z := \sum_k f_k\, e_k \in \mathfrak{X}(M)$ (extended globally by partition of unity if the frame is local; for clarity, work first in $U$ and then patch). Compute its divergence using the [Divergence Leibniz Rule](/theorems/2752) on each summand:
\begin{align*}
\operatorname{div}(f_k\, e_k) = f_k\, \operatorname{div}(e_k) + \langle df_k, e_k\rangle = f_k\,\operatorname{div}(e_k) + e_k(f_k).
\end{align*}
Hence
\begin{align*}
\sum_{k=1}^n e_k(f_k) = \sum_{k=1}^n \operatorname{div}(f_k\, e_k) - \sum_{k=1}^n f_k\,\operatorname{div}(e_k) = \operatorname{div}(Z) - \sum_{k=1}^n f_k\,\operatorname{div}(e_k).
\end{align*}
Now $\operatorname{div}(e_k) = \sum_l g(\nabla_{e_l} e_k, e_l)$ in the orthonormal frame, and using metric compatibility together with $g(e_k, e_l) = \delta_{kl}$ — so $0 = e_l(g(e_k, e_l)) = g(\nabla_{e_l} e_k, e_l) + g(e_k, \nabla_{e_l} e_l)$ — one obtains $g(\nabla_{e_l} e_k, e_l) = -g(e_k, \nabla_{e_l} e_l)$. Summing,
\begin{align*}
\operatorname{div}(e_k) = -\sum_l g(e_k, \nabla_{e_l} e_l) = -g\big(e_k, \sum_l \nabla_{e_l} e_l\big).
\end{align*}
Likewise,
\begin{align*}
\sum_k f_k\,\operatorname{div}(e_k) = -\sum_k g(i(e_k)\beta, \alpha)\, g\big(e_k, \sum_l \nabla_{e_l} e_l\big).
\end{align*}
The frame-curvature terms — $\sum_k f_k\,\operatorname{div}(e_k)$ and $\sum_k g(i(\nabla_{e_k} e_k)\beta, \alpha)$ in (IBP) — cancel against each other. Indeed, using $i(\nabla_{e_k} e_k)\beta(V) = \beta(\nabla_{e_k} e_k, V)$ and the antisymmetry of $\beta$, one computes (with $W := \sum_k \nabla_{e_k} e_k$)
\begin{align*}
\sum_k g(i(\nabla_{e_k} e_k)\beta, \alpha) = \sum_k \sum_l \beta(\nabla_{e_k} e_k, e_l)\alpha(e_l) = \sum_l \beta(W, e_l)\alpha(e_l) = g(i(W)\beta, \alpha).
\end{align*}
And $\sum_k f_k\,\operatorname{div}(e_k) = -\sum_k g(i(e_k)\beta, \alpha)\cdot g(e_k, W) = -g(i(W)\beta, \alpha)$ by linearity of $i(\cdot)\beta$ in its argument. The two terms therefore add to zero. Combining,
\begin{align*}
\sum_{k=1}^n g\big(i(e_k)\beta, \nabla_{e_k}\alpha\big) = \operatorname{div}(Z) - \sum_{k=1}^n g\big(i(e_k)(\nabla_{e_k}\beta), \alpha\big).
\tag{IBP$'$}
\end{align*}
[/step]
[step:Integrate over $M$, apply the Divergence Theorem, and conclude]
Combine (E) and (IBP$'$):
\begin{align*}
g(\beta, d\alpha) = \operatorname{div}(Z) - \sum_{k=1}^n g\big(i(e_k)(\nabla_{e_k}\beta), \alpha\big).
\end{align*}
Integrating over $M$ against $\omega_g$, both sides are smooth on the compact $M$ and hence integrable:
\begin{align*}
\int_M g(\beta, d\alpha)\,\omega_g = \int_M \operatorname{div}(Z)\,\omega_g - \int_M \sum_{k=1}^n g\big(i(e_k)(\nabla_{e_k}\beta), \alpha\big)\,\omega_g.
\end{align*}
The vector field $Z = \sum_k f_k\, e_k$ is smooth and globally defined on $M$ (by globalising via a partition of unity, the resulting integral identity is independent of the choice). The hypotheses of the [Divergence Theorem](/theorems/2754) — $M$ oriented, compact, and without boundary — hold by hypothesis, so $\int_M \operatorname{div}(Z)\,\omega_g = 0$. Therefore
\begin{align*}
\int_M g(\beta, d\alpha)\,\omega_g = -\int_M \sum_{k=1}^n g\big(i(e_k)(\nabla_{e_k}\beta), \alpha\big)\,\omega_g.
\end{align*}
By (A), the left-hand side equals $\int_M g(\delta\beta, \alpha)\,\omega_g$. Hence
\begin{align*}
\int_M g(\delta\beta, \alpha)\,\omega_g = -\int_M \sum_{k=1}^n g\big(i(e_k)(\nabla_{e_k}\beta), \alpha\big)\,\omega_g \qquad \text{for every } \alpha \in \Omega^1(M).
\end{align*}
This is identity (T). The fundamental lemma of the calculus of variations — applied to the smooth $1$-form $\delta\beta + \sum_k i(e_k)(\nabla_{e_k}\beta) \in \Omega^1(M)$ — gives
\begin{align*}
\delta\beta = -\sum_{k=1}^n i(e_k)(\nabla_{e_k}\beta) \qquad \text{pointwise on } M.
\end{align*}
Evaluating both sides on a vector field $Y \in \mathfrak{X}(M)$ at a point of $U$,
\begin{align*}
(\delta\beta)(Y) = -\sum_{k=1}^n \big(i(e_k)(\nabla_{e_k}\beta)\big)(Y) = -\sum_{k=1}^n (\nabla_{e_k}\beta)(e_k, Y).
\end{align*}
This is the asserted local formula.
[/step]