Let $M$ be an oriented Riemannian manifold. Let $\nu$ be a differential form of degree $n=\dim M$ defined in the following way: \begin{align} \nu(v_1,\ldots,v_n)(p)&=\pm\sqrt{\det(\langle v_i,v_j \rangle)} \\ &= \text{orient. vol. } \{v_1,\ldots,v_n\},\quad p\in M, \end{align} where $v_1,\ldots,v_n \in T_p(M)$ are linearly independent, and the oriented volume is affected by the sign $+$ or $-$ depending on whether or not the basis $\{v_1,\ldots,v_n\}$ belongs to the orientation of $M$; $\nu$ is called the volume element of $M$. For a vector field $X \in \mathcal X(M)$ define the interior product $i(X)\nu$ of $X$ with $\nu$ as the $(n-1)$-form: $$ i(X)\nu(Y_2,\ldots,Y_n)=\nu(X,Y_2,\ldots,Y_n), \quad Y_2,\ldots,Y_n \in \mathcal X(M). $$ Prove that $d(i(X)\nu)=\text{div} X\nu$.
Hint: Let $p \in M$ and let $E_i$ be a geodesic frame at $p$. Write $X$ as a sum $X=\sum_i f_i E_i$ and let $\omega_i$ be differential forms of degree one defined on a neighborhood of $p$ by $\omega_i(E_j)=\delta_{ij}$. Show that $\omega_i \wedge \cdots \wedge \omega_n$ is a volume form of $\nu$ on $M$. Next put $\theta_i=\omega_1 \wedge \cdots \wedge \hat \omega_i \wedge \cdots \wedge \omega_n$, where $\hat \omega_i$ signifies that the factor $\hat \omega_i$ is not present. $\color{red}{\textrm{Prove that $i(X) \nu = \sum_i (-1)^{i+1} f_i \theta_i$.}}$ It then follows that \begin{align} d(i(X)\nu)&=\sum_i (-1)^{i+1} df_i \wedge \theta_i + \sum_i (-1)^{i+1} f_i \wedge d\theta_i \\ &= (\sum_i E_i(f_i))\nu+\sum_i (-1)^{i+1} f_i \wedge d\theta_i. \end{align} But $d\theta_i=0$ at $p$, since \begin{align} \color{blue}{\textrm d \omega_k(E_i,E_j)} &= \color{blue}{\textrm E_i \omega_k(E_j)-E_j \omega_k(E_i)-\omega_k([E_i,E_j])} \\ &= \omega_k(\nabla_{E_i} E_j - \nabla_{E_j} E_i). \end{align} Therefore $$ d(i(X)\nu)(p)=(\sum_i E_i(f_i)(p))\nu=\text{div} X(p)\nu, $$ and since $p$ is arbitrary, this completes the proof.
I cannot figure out the part---which I highlighted in $\color{red}{\textrm{red}}$ from the hint---where we have to establish $i(X)\nu = \sum_i (-1)^{i+1} f_i \theta_i$. Because I was not sure what I can do with the given $(n-1)$-form $$ i(X)\nu(Y_2,\ldots,Y_n) = \nu(X,Y_2,\ldots,Y_n). $$
Edit: It turns out that I also did not quite understand what I highlighted above in $\color{blue}{\textrm{blue}}$: $$ d \omega_k(E_i,E_j) = E_i \omega_k(E_j)-E_j \omega_k(E_i)-\omega_k([E_i,E_j]). $$ Why is this equality true? I am not sure how was one able to take the exterior derivative of $d\omega_k(E_i,E_j)$. (Though, the equality does keep reminding me of the curvature tensor introduced in Chapter 4 of Do Carmo...)
Other than these, I understand the details of all other steps given in the hint.
Let $\theta_i=\omega_i\wedge\dots\hat\omega_i\wedge\dots\omega_n$, where $\hat\omega_i$ denotes that $\omega_i$ is removed. It suffices to check the desired result on vectors $v_2,\dots,v_n$: \begin{align*} i_X(\nu)(v_2,\dots,v_n) & = \nu(X,v_2,\dots,v_n) \\ & = \sum_if_i\nu(E_i,v_2,\dots,v_n) \\ & = \sum_if_i\omega_1\wedge\dots\wedge\omega_n(E_i,v_2,\dots,v_n) \\ & = \sum_i(-1)^{i+1}f_i\theta_i(v_2,\dots,v_n), \end{align*} where $(-1)^{i+1}$ arises from the number of swaps needed to shift $E_i$ to the $i^{\mathrm{th}}$ coordinate, which is the only permutation with a nonzero evaluation since $\omega_i(E_j)=\delta^i_j$.
The equality in blue follows from the following result: Let $\alpha$ be a smooth 1-form and let $X,Y\in\mathfrak{X}(M)$ be smooth vector field. Then $d\alpha(X,Y)=X\alpha(Y)-Y\alpha(X)-\alpha([X,Y])$.
By definition we have $$d\alpha(X,Y) = (i_X(d\alpha))(Y),$$ and recall Cartan's formula: $$\mathcal{L}_X\alpha = d(i_X\alpha)+i_X(d\alpha),$$ where $\mathcal{L}_X$ is the Lie derivative with respect to $X$. Now we see that \begin{align*} d\alpha(X,Y) & = (i_X(d\alpha))(Y) \\ & = (\mathcal{L}_X\alpha)(Y) - d(i_X\alpha)(Y) \\ & = (\mathcal{L}_X\alpha)(Y) - Y(\alpha(X)) \\ & = \left(\mathcal{L}_X(\alpha(Y)) - \alpha(\mathcal{L}_XY)\right) - Y(\alpha(X)) \\ & = X(\alpha(Y)) - \alpha([X,Y]) - Y(\alpha(X)). \end{align*} Alternatively, if you haven't worked too much with differential forms (in particular with Lie derivatives and Cartan's formula) you should also be able to prove this result by noting that every smooth 1-form $\alpha$ can locally be written as $\alpha = g~df$ where $f$ and $g$ are smooth functions and then showing that both sides of the desired equality are indeed equal.