I am attempting to prove that divergence of a smooth vector field $X$ over a $n$-dimensional Riemannian manifold $(M,g)$ is invariant under change of coordinates (or invariant of the choice of frame).
I am stuck in the proof and here is my attempt.
I am starting from the formula given in [p. 33, 1]. So, we have for $$X=X^i \frac{\partial}{\partial x^i}$$ the divergence is $$ \text{div}\left( X^i \frac{\partial}{\partial x^i}\right) = \frac{1}{\sqrt{\det g_{ij}} }\frac{\partial}{\partial x^i}\big( X^i \sqrt{\det g_{ij}}\big). $$
As a first step I am trying to prove that the RHS of that formula is the same for two different orthonormal frames.
So, I assume an orthonormal frame $\{E_i,E_2,...,E_n\}$ with dual frame $\{e^1,e^2,...,e^n\}$. First of all in this frame we have $$X^i = e^i(X)$$ and due to orthonormality we obtain $$ \det g_{ij} = 1.$$ I shall use $\text{div}(X)_{Ee}$ to emphasize the dependence on the frame. So we obtain $$\text{div}(X)_{Ee} = E_i(e^i(X))$$
Now I consider another orthonormal frame $\{F_i,F_2,...,F_n\}$ with dual frame $\{f^1,f^2,...,f^n\}$. For this frame we have $$\text{div}(X)_{Ff} = F_i(f^i(X))$$
Now, to prove invariance under the two frames, I need to prove $$\text{div}(X)_{Ff} = \text{div}(X)_{Ee}$$
I expand the second frame in terms of the first as $$\begin{eqnarray} F_i &=& F_i^l E_l \\ f^i &=& f^i_m e^m\end{eqnarray}$$
Then we have $$\begin{eqnarray} \text{div}(X)_{Ff} &=& F_i(f^i(X)) \nonumber \\ &=& F_i^lE_l(f^i_m e^m(X)) \nonumber \\ &=& F_i^lE_l(f^i_m) e^m(X) + F_i^lf^i_m E_l(e^m(X)) \\ &=& F_i^lE_l(f^i_m) e^m(X) + \delta^l_m E_l(e^m(X)) \\ &=& F_i^lE_l(f^i_m) e^m(X) + \text{div}(X)_{Ee} \end{eqnarray}$$
Where the third step is from Leibniz's rule and the fourth step is because the matrices $[f_m^i]$ and $[F_i^l]$ are inverses of each other.
So essentially I need to prove $$F_i^lE_l(f^i_m) e^m(X) = 0$$ and I am not able to do this, can someone help?
[1] Lee, John M. Introduction to Riemannian manifolds. Vol. 176. Springer, 2018.
P.S. I did go through Divergence of a smooth vector field and that answer uses the covariant derivative. I am attempting to avoid using it.
I don't understand what you're trying to do here. The formula you're using ($\text{div}(X)_{Ee} = E_i(e^i(X))$) doesn't appear in my book, and it's not correct, except in a flat Riemannian manifold with a coordinate frame.
Page 32 of my book gives a definition of $\text{div}\, X$ that's manifestly coordinate-invariant: it's the unique function that satisfies $$ d(X\mathrel{\lrcorner} dV_g) = (\text{div}\, X) dV_g.$$ Then on page 33, there's a coordinate formula that works in any local coordinates: $$ \text{div}\left( X^i \frac{\partial}{\partial x^i}\right) = \frac{1}{\sqrt{\det g} }\frac{\partial}{\partial x^i}\big( X^i \sqrt{\det g}\big). $$ This reduces to the familiar vector calculus formula in the special case of standard coordinates on Euclidean space.
It's an interesting exercise to prove that the general coordinate formula follows from the coordinate-free definition. You can also try to prove directly that the coordinate formula transforms correctly when you change coordinates, but that's a lot messier and not necessary.