I am looking to derive an identity I found which commutes the Hessian and the Laplacian of a scalar function $f \in C^4(M)$ on a Riemannian manifold $(M,g).$
$$\Delta \nabla_i \nabla_j f = \nabla_i \nabla_j \Delta f + (R_{jp}g_{ik} + R_{ip} g_{jk} - 2R_{kipj}) \nabla_k \nabla_p f + (\nabla_i R_{jp} + \nabla_j R_{pi} - \nabla_p R_{ij}) \nabla^p f. $$
One can note you have the Laplacian of the Hessian on the LHS of the equality, and the Hessian of the Laplacian on the RHS + some extra terms. This formula
It can be found in page 28 of these Lecture Notes: https://www.math.uci.edu/~jviaclov/courses/865_F07.html
The problem with understanding the derivation he does there is the notation, as I am more familiar with index notation, so that the Laplacian of a scalar function for example is expressed as $\Delta f = g^{ik}u_{,ik}$ and the Hessian is $u_{,ij}.$
First, $\Delta \nabla_i \nabla_j f$ is the Laplacian of the Hessian, so that in index notation it would be $g^{kl}(f_{,ij})_{,kl},$ whilst he writes $g^{kl}\nabla_k \nabla_l \nabla_i \nabla_j f.$ From this I understand that he is doing $g^{kl}\nabla_k \nabla_l (\nabla_i \nabla_j f)$ and he operates from the left. But then he does $g^{kl}\nabla_k \nabla_l \nabla_i \nabla_j f = g^{kl} \nabla_k (\nabla_l \nabla_i \nabla_j f),$ which would mean $g^{kl}(f_{,lij})_{,k}$ but this does not seem to work well in index notation and I can no longer follow through his steps.
The following question helps me but I would need two covariant derivatives instead of one: Commutator of laplacian and covariant derivative of a tensor