Unexpected regularity of the distance from a $C^2$ submanifold Let $\Gamma$ be a $C^2$ compact submanifold of $\mathbb{R}^n$. Consider the distance function $\delta$ from $\Gamma$. It is well known that, for sufficiently small $\varepsilon>0$, $\delta$ is $C^2$ on $\{ 0<\delta < \varepsilon\}$, and that it satisfies the eikonal equation
$$ \| \nabla \delta \| = 1, \qquad \text{with} \qquad  \delta|_{\Gamma} = 0. $$
Now recall Bochner's formula, valid for a smooth function $u \in C^\infty(M)$ on a general Riemannian manifold. It reads:
$$ \frac{1}{2} \Delta\left( \| \nabla u\|^2 \right) = \nabla u \cdot \nabla \Delta u + \| \mathrm{Hess} (u) \|^2_{\mathrm{HS}} + \mathrm{Ric}(\nabla u, \nabla u). \qquad (\star)$$
Here $\| \cdot \|_{\mathrm{HS}}$ is the Hilbert-Schmidt norm and $\Delta$ is the Laplace-Beltrami operator. When we specify $(\star)$ to $M = \mathbb{R}^n$ and $u$ is a smooth solution of the eikonal equation, we have the following:
$$ \nabla u \cdot \nabla \Delta u = - \| \mathrm{Hess} (u) \|^2_{\mathrm{HS}}. \qquad (\star\star)$$
Observe that the r.h.s. requires only $C^2$ regularity of $u$, while the l.h.s., a priori, requires a further derivative ($\nabla u \cdot \nabla \Delta u$ is indeed the directional derivative of $\Delta u$ in the direction $\nabla u$).
Q: Is it true, even if $\delta$ is a priori only a $C^2$ solution of the eikonal equation, that $\Delta \delta$ admits a directional derivative in the direction of $\nabla \delta$? Or, in other words, is it true that $(\star\star)$ holds for $\delta$?
The same question indeed can be posed in the Riemannian setting.
P.S: A direct computation with the distance from the $C^2$ (but not $C^3$) surface in $\mathbb{R}^2$ given by $y= x^{5/2}$ seems to support this claim.
This could be a standard fact about the regularity of solutions of the eikonal equation, but I have not been able to find any reference on this precise point.
 A: If $u$ is smooth, the left-hand-side in $(\star\star)$ could be rewritten as $(\nabla u)(\Delta u)$. Assuming that $u$ is $C^2$ this expression $(\nabla u)(\Delta u)$ is always defined, while your original expression $\nabla u \cdot \nabla \Delta u$ might be undefined.
Indeed assume $\gamma$ is a unit-speed geodesic such that $u\circ\gamma(t)\equiv t +\mathrm{const}$, or equivalently $\dot\gamma(t)=\nabla_{\gamma(t)}u$.
Then $H(t)=\mathrm{Hess}_{\gamma(t)}u$ satisfies Riccati equation
$$H'+H^2=0$$
Taking the trace your equation follows.
Hope it helps.
A: Since this question seems to have attracted some interest, I will post my own answer (which is a proof of the statement in the comments by Anton).
If $u : M \to \mathbb{R}$ is a $C^2$ solution of the Eikonal equation
$$ \| \nabla u\| = 1, $$
then $\mathrm{Hess}(u)$, which a priori is only continuous, is smooth along integral lines of $\nabla u$ (and analytic, in the Euclidean case).
The proof boils down to show that the Riccati equation for $\mathrm{Hess}(u)$ holds true in the distributional sense (any "simpler" proof is welcome). Assume $M = \mathbb{R}^n$. The eikonal equation is (sum over repeated indices is assumed)
$$(\partial_i u) (\partial_i u) = 1.$$
Its first derivative gives 
$$(\partial_{ij} u) (\partial_i u) = 0.$$
We can indeed differentiate again, but not use Leibnitz since $\partial_{ij} u$ might not be differentiable. Nevertheless, as distributions we obtain that
$$\partial_k (\partial_i \partial_j u) (\partial_i u)= - (\partial_{k\ell} u)( \partial_{\ell j} u) .$$ 
(To be precise one must consider $\partial_{i}\partial_j u \in C_c^0(\mathbb{R}^n)^*$ so that its distributional derivative is in $C_c^1(\mathbb{R}^n)^*$, and so the multiplication for $C^1(\mathbb{R}^n)$ functions, such as $\partial_i u$, is well defined). We can interchange the derivatives w.r.t. $k$ and $i$ in the distributional sense (indeed we can do it for test functions in $C^\infty_c(\mathbb{R}^n)$, and then by density also for test functions in $C^1_c(\mathbb{R}^n)$):
$$\partial_i (\partial_k \partial_j u) (\partial_i u)  = - (\partial_{k\ell} u)( \partial_{\ell j} u). $$ 
Since the r.h.s. is continuous, the identity holds also in the strong sense. In particular, the Hessian $H_{kj} := \partial_{k}\partial_j u$ admits derivative in the direction of $v = \nabla u$ and it satisfies the Riccati equation
$$ D_v H = -H^2 .$$
One can compute explicit solutions of this along the integral lines of $v=\nabla u$, and they are analytic.
Riemannian case. The argument generalizes in a straightforward way to the Riemannian setting, with the non-trivial occurrence of the curvature due to the interchange of covariant derivatives. Denoting with $S$ the symmetric $(1,1)$ tensor associated with $\mathrm{Hess}(u)$ (the shape operator of the level sets of $u$):
$$ \nabla_{v} S + S^2 + R = 0,$$
where $R$ is the curvature operator $g(R X,Y) = \mathrm{Rm}(x,\nabla u,\nabla u,Y)$. In particular, $S$ is smooth along integral lines of $\nabla u$. Taking the trace, we get:
$$(\nabla u) ( \Delta u) + \|\mathrm{Hess}(u)\|_{\mathrm{HS}}^2 + \mathrm{Ric}(\nabla u,\nabla u) = 0,$$
which is the particular case of Bochner's formula for $C^2$ solutions of the Eikonal equation.
Comments. These equations are indeed classical, the only non-trivial part (at least to me), was the unexpected extra regularity of $\mathrm{Hess}(u)$ along the integral curve $\gamma(t)$ of $\nabla u$ (which is indeed a geodesic).
Remark. This means that the mean curvature $\Delta u$ of a $C^2$ surface $\{u = 0\}$ is actually smooth in the normal direction (and analytic if the ambient space is Euclidean).
