The divergence on a Riemannian manifold is defined as $$\operatorname{div}X=\operatorname{tr}(Y\mapsto \nabla_YX)$$ for $X,Y\in \mathfrak{X}(M)$.
I found online that this can be apparently defined also as $$d(\iota_XdV_g)= (\operatorname{div}X)dV_g.$$ Here $dV_g$ is the volume form with the metric $g$.
I have a hard time understanding why these two definitions are equivalent. It is also furthermore said that $\mathcal{L}_X(dV_g)=(\operatorname{div}X)dV_g$ which only adds to my confusion. Could anyone elaborate on why these are all equivalent definitions?
Let $\{E_1,\ldots,E_n\}$ be an oriented orthonormal frame, and$\{\theta^1,\ldots,\theta^n\}$ be its dual coframe, that is $\theta^j = g(\cdot,E_j)$. $\newcommand{\d}{\mathrm{d}} \newcommand{\L}{\mathscr{L}}$ Then locally, $$ \d V_g = \theta^1\wedge \cdots \wedge \theta^n. $$ Let $X$ be a vector field. Since $\d(\d V_g) = 0$, by Cartan's magic formula $$ \d(\iota_X\d V_g) = \d(\iota_X\d Vg) + \iota_X\d(\d V_g) = \L_X(\d V_g). $$ The Lie derivative is a derivation of the exterior algebra. Hence $$ \L_X(\d V_g) = \sum_{j=1}^n \theta^1\wedge\cdots\wedge \L_X\theta^j\wedge\cdots\wedge\theta^n. $$ For $j\in\{1,\ldots,n\}$, and $Y$ any vector field, Leibniz rule together with the torsion-free condition gives \begin{align} (\L_X\theta^j)(Y) &= \L_X(\theta^j(Y)) - \theta^j(\L_XY) \\ &= X(g(Y,E_j)) - g([X,Y],E_j)\\ &= g(\nabla_XY-[X,Y],E_j) + g(Y,\nabla_XE_j)\\ &= \theta^j(\nabla_YX) + \alpha^j(Y), \end{align} with $\alpha^j(Y) = g(Y,\nabla_XE_j)$. Notice that $\nabla_XE_j \perp E_j$ by differentiating $\|E_j\|^2=1$ with respect to $X$. It follows that $\alpha^j \in \mathrm{span}\{\theta^k \mid k\neq j\}$. Therefore, $\theta^1\wedge\cdots\wedge \theta^{j-1}\wedge \alpha^j\wedge \theta^{j+1}\wedge\cdots\wedge\theta^n = 0$, and it follows that $$ \L_X(\d V_g) = \sum_{j=1}^n \theta^1\wedge\cdots\wedge (\theta^j\circ \nabla X) \wedge \cdots \wedge \theta^n. $$ Finally, since $\Lambda^n(T^*M)$ has rank $1$, there exists a function $f$ such that $$ \L_X(\d V_g) = f\, \d V_g. $$ Applying to $\{E_1,\ldots,E_n\}$, we obtain \begin{align} f &= \left(\sum_{j=1}^n \theta^1\wedge\cdots\wedge (\theta^j\circ \nabla X)\wedge\cdots\wedge\theta^n\right)(E_1,\ldots,E_n) \\ &= \sum_{j=1}^n \theta^j(\nabla_{E_j}X)\\ &= \sum_{j=1}^n g(\nabla_{E_j}X,E_j) \\ &= \mathrm{trace}(\nabla X). \end{align}
This shows that both definitions of the divergence, namely $\mathrm{trace}(\nabla X)$ or the function $f$ above, are equivalent.
Here is a different but somewhat similar computation for $f$. We have \begin{align} f &= (\L_X(\d V_g))(E_1,\ldots,E_n) \\ &= \L_X(\d V_g(E_1,\ldots,E_n)) - \sum_{j=1}^n \d V_g(E_1,\ldots,[X,E_j],\ldots,E_n) \\ &= 0 - \sum_{j=1}^n \det_{\{E_1,\ldots,E_n\}}(E_1,\ldots,[X,E_j],\ldots,E_n) \\ &= -\sum_{j=1}^n g([X,E_j],E_j) \\ &= \sum_{j=1}^n g(\nabla_{E_j}X,E_j) - g(\nabla_X E_j,E_j)\\ &= \mathrm{trace}(\nabla X), \end{align} where we have used successively the Leibniz rule, the fact that $\d V_g(E_1,\ldots,E_n) = 1$ is constant, the fact that the volume form is the determinant in any orthonormal basis, the torsion-free condition, and finally the fact that $\nabla_XE_j\perp E_j$.