[proofplan] The divergence on a Riemannian manifold is defined as the trace of the covariant derivative: $\operatorname{div} X = \operatorname{tr}(\nabla X)$, where $\nabla X \in \Gamma(T^*M \otimes TM)$ is regarded as an endomorphism of $TM$ via the canonical isomorphism $T^*M \otimes TM \cong \operatorname{End}(TM)$. The Leibniz rule for $\nabla$ on the product $fX$ produces $\nabla(fX) = df \otimes X + f\,\nabla X$. Taking the pointwise trace and identifying $\operatorname{tr}(df \otimes X) = \langle df, X \rangle$ via the same canonical isomorphism gives the result. The argument is purely tensorial and pointwise. [/proofplan]

[step:Set up the divergence as the trace of the covariant derivative] Let $(M, g)$ be a Riemannian manifold, $\nabla$ the Levi-Civita connection, $f \in C^\infty(M)$, and $X \in \mathfrak{X}(M)$. The covariant derivative of $X$ is the $(1,1)$-tensor field \begin{align*} \nabla X : \mathfrak{X}(M) &\to \mathfrak{X}(M), \\ Y &\mapsto \nabla_Y X, \end{align*} which we regard as a section of $T^*M \otimes TM$, equivalently as a section of $\operatorname{End}(TM)$ via the trace-contraction pairing $T^*M \otimes TM \cong \operatorname{End}(TM)$. The divergence is defined pointwise by \begin{align*} \operatorname{div} X := \operatorname{tr}(\nabla X) \in C^\infty(M), \end{align*} so $(\operatorname{div} X)(p) = \operatorname{tr}\big( v \mapsto \nabla_v X \big)$ at each $p \in M$, where the trace is taken on the linear map $T_p M \to T_p M$, $v \mapsto (\nabla_v X)(p)$. [/step]

[step:Apply the Leibniz rule for $\nabla$ to the product $fX$] The Levi-Civita connection $\nabla$ is a derivation: for any $f \in C^\infty(M)$, $X \in \mathfrak{X}(M)$, and $Y \in \mathfrak{X}(M)$, \begin{align*} \nabla_Y(fX) = (Yf)\,X + f\,\nabla_Y X. \end{align*} Recognising $Yf = df(Y)$, this rewrites at the level of $(1,1)$-tensors as \begin{align*} \nabla(fX) = df \otimes X + f\,\nabla X, \tag{L} \end{align*} where $df \otimes X$ is the section of $T^*M \otimes TM$ acting on $Y \in \mathfrak{X}(M)$ by $(df \otimes X)(Y) = df(Y) \, X$. [/step]

[step:Take the pointwise trace and identify each term]Apply the trace, viewed as the canonical contraction $T^*M \otimes TM \to C^\infty(M)$, $\theta \otimes Z \mapsto \theta(Z)$, to both sides of (L). The trace is $C^\infty(M)$-linear, so \begin{align*} \operatorname{div}(fX) = \operatorname{tr}\big(\nabla(fX)\big) = \operatorname{tr}(df \otimes X) + f\,\operatorname{tr}(\nabla X). \end{align*} The second term is $f\,\operatorname{div} X$ by definition of $\operatorname{div}$. For the first term, fix $p \in M$ and a basis $\{e_i\}_{i=1}^n$ of $T_p M$ with dual basis $\{e^i\}$ of $T_p^* M$. Writing $df_p = \sum_i a_i\, e^i$ and $X(p) = \sum_j b_j\, e_j$, the endomorphism $df_p \otimes X(p) : T_p M \to T_p M$ acts by \begin{align*} v \mapsto df_p(v) \, X(p) = \big(\sum_i a_i\, e^i(v)\big) \big(\sum_j b_j\, e_j\big), \end{align*} whose matrix in the basis $\{e_i\}$ has $(i,j)$-entry $b_i a_j$. Its trace is \begin{align*} \operatorname{tr}(df_p \otimes X(p)) = \sum_i b_i a_i = df_p(X(p)) = \langle df, X \rangle(p), \end{align*} where $\langle df, X \rangle := df(X)$ is the natural pairing of a $1$-form with a vector field. Since $p$ was arbitrary, \begin{align*} \operatorname{tr}(df \otimes X) = \langle df, X \rangle. \end{align*} Combining, \begin{align*} \operatorname{div}(fX) = \langle df, X \rangle + f\,\operatorname{div} X = f\,\operatorname{div} X + \langle df, X \rangle. \end{align*} This is the asserted identity.[/step]

[guided]We have the tensor-level Leibniz identity from the previous step: \begin{align*} \nabla(fX) = df \otimes X + f\,\nabla X. \tag{L} \end{align*} Our goal is to take the trace of both sides and recognise each term as a recognisable object — the divergence on the left, and the function-derivative pairing $\langle df, X \rangle$ on the right. The trace we use is the canonical contraction \begin{align*} \operatorname{tr} : T^*M \otimes TM &\to C^\infty(M), \\ \theta \otimes Z &\mapsto \theta(Z), \end{align*} extended $C^\infty(M)$-linearly. Why is this the "right" trace? Because under the trace-contraction isomorphism $T^*M \otimes TM \cong \operatorname{End}(TM)$, this contraction is exactly the trace of the corresponding endomorphism. Applying it to (L) and using $C^\infty(M)$-linearity (the function $f$ pulls out of the trace), \begin{align*} \operatorname{div}(fX) = \operatorname{tr}\big(\nabla(fX)\big) = \operatorname{tr}(df \otimes X) + f\,\operatorname{tr}(\nabla X). \end{align*} The second term is recognisable immediately: by the very definition $\operatorname{div} X := \operatorname{tr}(\nabla X)$, we have $f\,\operatorname{tr}(\nabla X) = f\,\operatorname{div} X$. The first term requires more thought. What is $\operatorname{tr}(df \otimes X)$? Since trace is defined pointwise, fix $p \in M$ and compute in coordinates. Choose a basis $\{e_i\}_{i=1}^n$ of $T_p M$ with dual basis $\{e^i\}$ of $T_p^* M$, and expand \begin{align*} df_p = \sum_i a_i\, e^i, \qquad X(p) = \sum_j b_j\, e_j. \end{align*} The endomorphism $df_p \otimes X(p) : T_p M \to T_p M$ acts on $v \in T_p M$ by \begin{align*} v \mapsto df_p(v) \, X(p) = \Big(\sum_i a_i\, e^i(v)\Big)\Big(\sum_j b_j\, e_j\Big). \end{align*} What is its matrix? Reading off the coefficient of $e_i$ in the image of $e_j$, the $(i,j)$-entry is $b_i a_j$. The trace is the sum of diagonal entries: \begin{align*} \operatorname{tr}(df_p \otimes X(p)) = \sum_i b_i a_i = df_p(X(p)) = \langle df, X \rangle(p), \end{align*} where we recognise $\sum_i a_i b_i$ as the natural pairing $df_p(X(p))$, written $\langle df, X \rangle(p)$. This is exactly the abstract statement that the trace-contraction pairing $T^*M \otimes TM \cong \operatorname{End}(TM)$ sends the simple tensor $\theta \otimes Z$ to $\theta(Z)$ — here $\theta = df_p$ and $Z = X(p)$. Since $p \in M$ was arbitrary, \begin{align*} \operatorname{tr}(df \otimes X) = \langle df, X \rangle. \end{align*} Combining the two terms, \begin{align*} \operatorname{div}(fX) = \langle df, X \rangle + f\,\operatorname{div} X = f\,\operatorname{div} X + \langle df, X \rangle, \end{align*} which is the asserted identity. A remark on what was and was not used. The argument is purely tensorial and pointwise: nowhere did we invoke metric compatibility ($\nabla g = 0$) or torsion-freeness of $\nabla$. The only property of $\nabla$ used was that it is a derivation on $\mathfrak{X}(M)$ over $C^\infty(M)$ — this is what produced (L) — together with the algebraic identification of the trace of a simple tensor as the pairing. The function-derivative term $\langle df, X \rangle$ here plays the role of the elementary Euclidean expression $X(f)$ from vector calculus.[/guided]