Let $(M,g)$ be a Riemannian manifold and $f\in C^{\infty}(M)$. Why does the equality $$\nabla^j\nabla_i\nabla_jf-\nabla_i\nabla^j\nabla_jf=R_{ij}\nabla^jf$$ hold?
I tried writing $$\nabla^j\nabla_i-\nabla_i\nabla^j=g^{jk}(\nabla_k\nabla_i-\nabla_i\nabla_k)=g^{jk}R(e_k,e_i)$$ assuming $[e_i,e_j]=0$, but now I'm stuck.
For any vector field $Z$, we have that $$R(X,Y)Z = \nabla_X\nabla_YZ - \nabla_Y\nabla_XZ - \nabla_{[X,Y]}Z \tag{1}$$may be rewritten as $$R(X,Y)Z = \nabla_X\nabla_YZ - \nabla_{\nabla_XY}Z - (\nabla_Y\nabla_XZ - \nabla_{\nabla_YX}Z),\tag{2}$$i.e., $$R(X,Y)Z = [\nabla_X(\nabla Z)](Y) - [\nabla_Y(\nabla Z)](X),\tag{3}$$which then reads $$R_{ijk}^{\phantom{ijk}\ell}Z^k = \nabla_i\nabla_jZ^\ell - \nabla_j\nabla_i Z^\ell.\tag{4}$$Now make $i=\ell$ to obtain $$R_{jk}Z^k = \nabla_i\nabla_jZ^i - \nabla_j\nabla_iZ^i,\tag{5}$$which, in a coordinate-free way, means $${\rm Ric}(\cdot, Z) = {\rm div}(\nabla Z) - {\rm d}({\rm div}\,Z).\tag{6}$$For $Z = \nabla f$, this reduces to $${\rm Ric}(\cdot,\nabla f) = {\rm div}\,{\rm Hess}\,f - {\rm d}(\triangle f),\tag{7}$$where ${\rm Hess}\,f$ and $\triangle f$ are the Hessian and Laplacian of $f$, respectively.
N.B.: if your sign convention for the curvature is the opposite of the one I used in $(1)$, then signs get reversed on $(2)$-$(4)$, but not on $(5)$, which is valid for any torsionfree connection $\nabla$ on $M$.