If a scalar function $\phi$ satisfy the advection equation $$ \partial_t \phi + v_k \partial_k \phi = 0, $$ where $v_k$ is some vector field. Then the gradient $\partial_i \phi$ satisfy the following advection equation $$ \partial_t \partial_i \phi + v_k \partial_i \partial_k \phi + (\partial_i v_k) \partial_k \phi = 0. $$
With the help of Lie derivative we can rewrite the above equations as \begin{align} \partial_t \phi + \mathcal{L}_v \phi &= 0, \\ \partial_t \nabla \phi + \mathcal{L}_v \nabla \phi &= 0. \end{align}
What about the second derivative(Hessian) $\nabla^2 \phi$, what equation does it satisfy? I did the calculation in coordinates and I got $$ \partial_t \nabla^2 \phi + \mathcal{L}_v \nabla^2 \phi = - \nabla \phi \cdot \nabla^2 v, $$ just to clarify, we have $(\nabla \phi \cdot \nabla^2 v)_{ij} = (\partial_k \phi)(\partial_i \partial_j v_k)$.
Is there a nice way how to derive the equation for $\nabla^2 \phi$ without doing the calculation in coordinates? Does it scale to higher derivatives?