What is the gradient of length of gradient of function $f$? Suppose $f\in C^3(M)$ and $\nabla f$ denotes the gradient of $f$ w.r.t Riemannian metric $g$. Then what is the equivalent expression of the following?
$$\nabla \langle \nabla f, \nabla f\rangle=?$$
Background problem: By definition $\mathrm{Hess} f(X, Y) = \langle \nabla _X(\nabla f), Y\rangle$; $\triangle f= \mathrm{tr}(\mathrm{Hess} f)$. Now I want to calculate $\triangle|\nabla f|^2$. By above definitions we have:
$$\triangle|\nabla f|^2=\mathrm{tr}(\mathrm{Hess} \langle\nabla f,\nabla f\rangle)=\langle \nabla _{X_i}(\nabla \langle \nabla f, \nabla f\rangle), X_i\rangle.$$
But I saw somewhere that 
$$\triangle|\nabla f|^2=\sum_iX_iX_i\langle\nabla f,\nabla f\rangle.$$
Why these two (my calculation and last one) are equal?
 A: I'm using tensor notation.
$$
\nabla\langle\nabla f,\nabla f\rangle\equiv \nabla^a(\nabla_b f\:\nabla^bf)=(\nabla^a\nabla_bf)\nabla^bf+\nabla_bf(\nabla^a\nabla^bf)=2(\nabla^a\nabla^bf)\nabla_bf
$$
$$
\nabla_XZ \equiv X^a\nabla_aZ^b
$$
$$
\mathbf{Hess}\langle\nabla f,\nabla f\rangle(X,Y) \equiv
X^a\{\nabla_a \nabla_c(\nabla_b f\:\nabla^b f)\}Y^c =
X^a(\nabla_a \nabla_c h)Y^c
$$
To compute trace of (2,0)-tensor, we have to rise one index and then contract them. To do this we need dual basis to $X_i$, and when 
$\langle X_i,X_j\rangle=\delta_{ij}$ we get it by lowering index of $(X_i)^a$.
So $\mathbf{Tr\:T}=\sum_i{T_a}^b(X_i)^a(X_i)_b=\sum_iT_{ab}(X_i)^a(X_i)^b$
Action of a vector field on a function is $Xf\equiv X^a\nabla_a f$, so
$$
X\:Y h\equiv X^c\nabla_c(Y^a\nabla_ah)=(X^c\nabla_cY^a)\nabla_ah+X^c Y^a \nabla_c \nabla_a h
$$
If $X_i$ satisfy geodesic equation then $\nabla_{X_i}X_i\equiv (X_i)^c\nabla_c(X_i)^a=0$ and we are done, because
$$
\mathbf{Tr}(\mathbf{Hess}\:h)=\sum_i{X_i X_i}\:h 
$$
when $X_i$ are orthonormal, geodesic coordinate vector fields. 
A: I think I got the answer. Thanks for posting the question I have faced the same problem in the last few days. Here is a generalized solution. If I made any error kindly let me know.
We know that if $g$ is the Riemannian metric then $\langle\nabla f,\nabla h\rangle =g(\nabla f,\nabla h)=g_{pq}\nabla^p f\nabla^q h=\nabla_q f\nabla^q h$ (called lowering the index [similarly raising the index is another term]) thus
$$\begin{eqnarray}
\nabla\langle\nabla f,\nabla h\rangle &\equiv& \nabla^a(\nabla_b f\nabla^b h)\\ &=& \nabla^a(\nabla_b f)\nabla^b h+\nabla^a(\nabla^b h)\nabla_b f\\ &=& \nabla^a(\nabla_b f)g^{kb}\nabla_k h+\nabla^a(\nabla^b h)\nabla_b f\\ &=& \nabla^a(g^{kb}\nabla_b f)\nabla_k h+\nabla^a(\nabla^b h)\nabla_b f\text{ (Since $\nabla g=0$)}\\ &=& \nabla^a(\nabla^k f)\nabla_k h+\nabla^a(\nabla^b h)\nabla_b f
\end{eqnarray},$$
Changing the dummy index $k$ to $b$ we obtain.
$$\begin{eqnarray}&=&\nabla^a(\nabla^b f)\nabla_b h+\nabla^a(\nabla^b h)\nabla_b f\\ &\equiv& \text{Hess }f\ \nabla h+\text{Hess }h\ \nabla f\end{eqnarray}$$
In particular for $f=h$ we have the answer by user $87091403130$ which is $\nabla \langle\nabla f,\nabla f\rangle=2\text{Hess }f\ \nabla f$
You've asked another question about calculating $\Delta |\nabla f|^2$, I think we can use the Bochner formula to calculate over Euclidean space, Riemannian manifold and other suitable spaces. One can try to see the proof of it for better understanding.
