Unexpected regularity of the distance from a $C^2$ submanifold

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.

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).