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.



Unexpected regularity of the distance from a $C^2$ submanifold
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).


