Divergence of vector fields and the trace of the Lie derivative of the metric

Let $(M,g)$ be a Riemannian manifold. Why does the following formula hold?

$$\operatorname{div} X= \tfrac12\,\mathrm{tr}_g\bigl(\mathcal{L}_Xg\bigr) $$

$\operatorname{div} X$ is defined by $ (\operatorname{div} X)\omega_g= d(i_X(\omega_g))$ where $i_X$ is interior multiplication by $X$.

Here is my tuppence worth. I hope it is what you are after and please excuse the perhaps unfamiliar notation - comments/questions/suggestions are welcome!

Let $M$ be a Riemannian manifold endowed with metric $g$ where

$$g \,=\, \delta_{ab}e^{a} \otimes e^{b}$$

with respect to a $g$-orthonormal coframe $\{e^{a}\}$. Denote the dual $g$-orthonormal frame $\{X_{b}\}$ (i.e. $e^{a}(X_{b})=\delta^{a}_{b}$). We introduce the operator

$$A_{X} \,\equiv\, \nabla_{X} - \mathcal{L}_{X}\quad \text{for any}\; X\in TM$$

in terms of the Levi-Civita connection $\nabla_{X}$ and Lie derivative $\mathcal{L}_{X}$. Note that since $\nabla_{X},\mathcal{L}_{X}$ commutes with contractions, so does $A_{X}$, and that for any function $f$ it should be clear that $A_{X}f=0$. In particular, for any $X,Y,Z\in TM$:

$$\begin{align}A_{X}[\,g(Y,Z)\,] \,=\, 0 &\,=\, (A_{X}g)(Y,Z) + g(A_{X}Y,Z) + g(Y,A_{X}Z) \\ &\,=\, -(\mathcal{L}_{X}g)(Y,Z) + g(A_{X}Y,Z) + g(Y,A_{X}Z), \end{align}$$

since $\nabla$ is metric compatible. Furthermore, since $\nabla$ is torsion-free, we can show $A_{X}Y=\nabla_{Y}X$. Putting all this together yields

$$(\mathcal{L}_{X}g)(Y,Z) \,=\, g(\nabla_{Y}X,Z) + g(Y,\nabla_{Z}X)$$

for any $X,Y,Z\in TM$. In particular, for $V\in TM$:

$$\begin{align} (\mathcal{L}_{V}g)(X_{k},X_{k}) \,=\, tr_{g}(\mathcal{L}_{V}g) &\,=\, g(\nabla_{X_{k}}V,X_{k}) + g(X_{k},\nabla_{X_{k}}V) \,=\, 2\,i_{X_{k}}\nabla_{X_{k}}\widetilde{V} \end{align}$$

where $\widetilde{X}=g(X,-)$ denotes the metric dual of any $X\in TM$ (we use the metric compatibility of $\nabla$ to get here). Thus we have

$$i_{X_{k}}\nabla_{X_{k}}\widetilde{V} \,=\, \frac{1}{2}tr_{g}(\mathcal{L}_{V}g)$$

Now for torsion-free $M$, we have

$$de^{a} + \omega^{a}_{\;b} \wedge e^{b} \,=\, 0$$

in terms of the connection 2-forms $\{\omega^{a}_{\;b}\}$ defined by:

$$\nabla_{X_{b}}e^{a} \,=\, -\omega^{a}_{\;c}(X_{b})e^{c}.$$

It is easy to show using these two relations that $d\equiv e^{a} \wedge \nabla_{X_{a}}$. Then using $\text{div}(V)=\star d\star \widetilde{V}$ in terms of the Hodge map $\star$:

$$\text{div}(V) \,=\, \star d\star \widetilde{V} \,=\, \star(e^{k} \wedge \nabla_{X_{k}}\star\widetilde{V}) \,=\,i_{X_{k}}\nabla_{X_{k}}\widetilde{V} \,=\, \frac{1}{2}tr_{g}(\mathcal{L}_{V}g) $$

where again you use the metric-compatibility of $\nabla$ to commute it with the Hodge map: $\nabla_{X}\,\star=\star\nabla_{X}$.

NB: If you are familar with it, you may notice that the object $i_{X_{k}}\nabla_{X_{k}}$ is, up to possibly a sign, related to the co-derivative $\delta$ on $M$. It should also be noted that most of this, again up to a possible sign, follows through for pseudo-Riemannian manifolds.

Fix a normal coordinate system $(x^1,\dots,x^n)$ centered at $p$ and write $X = X^i \partial_i$. Then (all calculations are evaluated only at $p$)

$$ \mathcal{L}_{\partial_i}{\partial_j} = \nabla_{\partial_i} \partial_j - \nabla_{\partial_j} \partial_i = [\partial_j, \partial_i ] = 0, \\ \mathcal{L}_X(\partial_j) = -\mathcal{L}_{\partial_j}(X^i \partial_i) = -\partial_j(X^i) \partial_i - X^i \mathcal{L}_{\partial_j}(\partial_i) = -\partial_j(X^i) \partial_i, \\ (\mathcal{L}_X g)(\partial_j, \partial_j) = Xg(\partial_j, \partial_j) - 2g(\mathcal{L}_X \partial_j, \partial_j) = 2g \left( \partial_j(X^i) \partial_i, \partial_j \right) = 2 \partial_j(X^i) \delta_i^j, \\ \frac{1}{2} \operatorname{tr}_g \left( \mathcal{L}_X g \right) = \frac{1}{2} \sum_{j=1}^n (\mathcal{L}_X g)(\partial_j, \partial_j) = \partial_j(X^j) = \operatorname{div}(X). $$

The latter formula for the divergence in a normal coordinate system can be verified similarly. Namely,

$$ d(i_X(\omega_g))(\partial_1, \dots, \partial_n) = \sum_{i=1}^n (-1)^{i-1} \partial_i ((i_X \omega_g)(\partial_1, \dots, \hat{\partial_i},\dots,\partial_n)) = \sum_{i=1}^n \partial_i \omega_g(\partial_1,\dots,\partial_{i-1},X,\partial_{i+1},\dots,\partial_n) = \partial_i (X^i) .$$