Let $(M,g)$ be a Riemannian manifold. Let $\nabla_v$ be the covariant derivative in the $v$ direction for all $v\in T_xM$, and denote with $\nabla^k h$ the $(k,0)$-tensor field defined in local coordinates inductively by $$ \nabla^0h=dh,\quad(\nabla^kh)_{i_1,\dots,i_k}=(\nabla_{\partial_{i_1}}h)_{i_2,\dots,i_k}. $$ for any smooth function $h$.
My question is: is there a nice way to express the difference $\nabla\nabla_udh-\nabla_u\nabla dh$?
To avoid confusion, I am considering the expression given by $$ \nabla(\nabla_udh)(X,Y)-\nabla_u(\nabla dh)(X,Y)=\nabla_X(\underbrace{\nabla_udh}_{(1,0) -tensor\,field})(Y)-\nabla_u(\underbrace{\nabla dh}_{(2,0)-tensor\,field})(X,Y). $$ This looks somehow similar to the Riemannian curvature tensor applied to forms. I have tried to develop the difference, but I cannot see anything familiar. More generally (but maybe I am asking too much), is there a nice way to write $$ \nabla^k\nabla_udh-\nabla_u\nabla^kdh=? $$
Write $\nabla_u dh = c^1_1 ( u\otimes \nabla dh)$, where $c^1_1$ is the contraction, then
\begin{align} \nabla (\nabla_u dh ) &= \nabla(c^1_1 ( u\otimes \nabla dh)) \\ &=c^1_1 \nabla (u\otimes \nabla dh) \\ &= c^1_1( \nabla u \otimes \nabla dh + u \otimes \nabla \nabla dh) \end{align}
In particular, it means for all $X, Y$ and using Ricci identity,
\begin{align} \nabla (\nabla_u dh ) (X, Y) &= (\nabla_{\nabla_X u} dh) (Y)+ \nabla_X \nabla_u dh (Y)\\ &= (\nabla_{\nabla_X u} dh) (Y)+ \nabla_u \nabla_X dh (Y) + R(u, X)dh (Y) \end{align}
thus
$$\big( \nabla (\nabla_u dh ) - \nabla_u \nabla dh \big)(X, Y) = (\nabla_{\nabla_X u} dh) (Y)+ R(u, X)dh (Y).$$
so as expected the curvature terms come out. Also we have $\nabla u$. In general, when calculating $$ \nabla^k \nabla_u dh- \nabla _u \nabla^k dh,$$ you have to differentiate $u$ $k$-times and use Ricci identity $k$-times. I guess there won't be a nice formula.