$\newcommand{\Ric}{\text{Ric}}$ $\newcommand{\div}{\text{div}}$ $\newcommand{\tr}{\text{trace}}$ $\newcommand{\g}{\mathfrak{g}}$
Let $M$ be a smooth closed oriented Riemannian manifold.
I am looking for a reference (or a sketch of proof) for the following equality:
For every vector field $V \in \Gamma(TM)$, $$ \int_M-\Ric(V,V)+ |\nabla V|^2 =\int_M \frac{1}{2}|L_Vg|^2-\big(\div(V)\big)^2,$$
where $ \nabla$ is the Levi-Civita connection, and all integrations are against the Riemannian volume form.
Motivation:
The LHS is $ \int_M \langle JV , V \rangle$ where $J:\Gamma(TM) \to \Gamma(TM)$ is the Jacobi operator, i.e.
$$ J(V)=J_{\text{id}_M}(V)=-\tr_{\g} R^{M}\big(V,(\cdot)\big)(\cdot)+\nabla^*{\nabla}V. $$
So, the question is to prove another formula for the hessian of the Dirichlet energy at the identity, namely:
$$ H(E)_{\text{id}_M}(V,V)=\int_M \langle JV , V \rangle=\int_M \frac{1}{2}|L_Vg|^2-\big(\div(V)\big)^2.$$
This formula is useful in showing various results regarding the stability of the identity map (for instance w.r.t to conformal variations).
I found a source claiming that there is a reference for this in "K. Yano and S. Bochner, Curvature and Betti numbers" (1953) - pg 57, but I could not find it there.
Another source cited "On harmonic and Killing vector fields, Yano" (formula 4.4) but again, I did not find it there.
Perhaps I am missing something, or having trouble with the notation (both references are from the 50's).
$\newcommand{\M}{\mathcal{M}}$ $\newcommand{\Ric}{\text{Ric}}$ $\newcommand{\div}{\text{div}}$ $\newcommand{\tr}{\text{trace}}$ $\newcommand{\g}{\mathfrak{g}}$
This answer is based on Anthony Carapetis's comment.
We shall use the following identity (See equation $(3)$ here):
$$\int_{\M} -\Ric(V,V)= \int_{\M} -(\div V)^2 +\tr ( \nabla V \circ \nabla V) . \tag{1}$$
Using this, the required identity $$\int_M-\Ric(V,V)+ |\nabla V|^2 =\int_M \frac{1}{2}|L_Vg|^2-\big(\div V \big)^2,$$ reduces to
$$\int_{\M} \tr ( \nabla V \circ \nabla V)+ |\nabla V|^2 =\int_M \frac{1}{2}|L_Vg|^2. \tag{2}$$
Now, using the fact $$ L_Vg(X,Y)=\langle \nabla_X V,Y \rangle+\langle \nabla_Y V,X \rangle,$$
and choosing an orthonormal basis $e_i$ for $T\M$, we obtain
$$ \frac{1}{2}|L_Vg|^2=\frac{1}{2}\sum_{ij} \big(L_Vg(e_i,e_j)\big)^2=\frac{1}{2}\sum_{ij} \big(\langle \nabla_{e_i} V,e_j \rangle+\langle \nabla_{e_j} V,e_i \rangle \big)^2$$
$$ =\sum_{ij} \langle \nabla_{e_i} V,e_j \rangle^2+ \sum_{ij} \langle \nabla_{e_i} V,e_j \rangle \langle \nabla_{e_j} V,e_i \rangle =\sum_{i} |\nabla_{e_i} V|^2+ \sum_{ij} \langle \nabla_{e_i} V,e_j \rangle \langle \nabla_{e_j} V,e_i \rangle$$
$$ = |\nabla V|^2+ \sum_{ij} \langle \nabla_{e_i} V,e_j \rangle \langle \nabla_{e_j} V,e_i \rangle.$$
Using the last equation, the required result $(2)$ reduces to
$$ \int_{\M} \tr ( \nabla V \circ \nabla V)=\int_{\M} \sum_{ij} \langle \nabla_{e_i} V,e_j \rangle \langle \nabla_{e_j} V,e_i \rangle,$$
But now the integrands are equal; This follows from a pointwise identity. Given a linear map $T:V \to V$, $$\tr(T^2)=\sum_{ij} \langle T(e_i) ,e_j \rangle \langle T(e_j) ,e_i \rangle,$$ where $e_i$ is an orthonormal basis. (Apply this for $T=\nabla V$).
Indeed, write $ A_{ij}:=\langle T(e_i) ,e_j \rangle$. $(A_{ij})$ is the representing matrix of $T$ w.r.t the basis $e_i$. Then
$$\sum_{ij} \langle T(e_i) ,e_j \rangle \langle T(e_j) ,e_i \rangle=\sum_{ij} A_{ij} A_{ji}=\langle A,A^T \rangle=\tr((A^T)^TA)=\tr(A^2).$$
Since $A^2$ is the representing matrix of $T$, we deduce
$$\sum_{ij} \langle T(e_i) ,e_j \rangle \langle T(e_j) ,e_i \rangle=\tr(T^2)$$ as required.