Are distance functions necessarily nonsmooth on the cut locus? Let $(M,g)$ be a complete Riemannian manifold and fix $p\in M$. Consider the distance function $r(x):=d(p,x)$. It is well-known that $r$ is smooth outside $\operatorname{cut}(p)\cup\{p\}$ where $\operatorname{cut}(p)$ is the cut locus of $p$. My question is:

Is $r$ necessarily nonsmooth on every point of $\operatorname{cut}(p)$?

It is well-known that $x\in\operatorname{cut}(p)$ if and only if either (a) there are two distinct unit-speed minimizing geodesics $\gamma_1,\gamma_2:[0,\ell]\to M$ joining $p$ and $x$, or (b) $x$ is a critical value of $\exp_p$. In Peter Petersen's Riemannian Geometry, the author gave a remark on this: In case (a), $\nabla r$ could be either $\gamma_1'(\ell)$ or $\gamma_2'(\ell)$ and hence does not exist; in case (b), $\operatorname{Hess}r$ is undefined since it must tend to $-\infty$ along certain
fields.
I know that the part about (a) is intuitive, but is there any way to make the argument rigorous? O the other hand, I don't see why $\operatorname{Hess}r$ must blow up.
 A: Answer: Yes, they are non-smooth.
For a): Let $U$ be the set on which $r$ is differentiable. Since $r$ is 1-Lipschitz, we have $\Vert \nabla r \Vert \le 1$. Anyway, what I want to show is that for any shortest geodesic $\gamma$ with $\gamma(0) = p$, we have $$(\nabla r)_{\gamma(t)} = \gamma'(t).$$
For this, let $v \in T_{\gamma(t)}M$ be arbitrary and $\tilde{\gamma}$ the geodesic with $\tilde{\gamma}(0) = \gamma(t)$ and $\tilde{\gamma}'(0) = v$. Then we can compute $$\langle(\nabla r)_{\gamma(t)},v\rangle = (dr)_{\gamma(t)} \cdot v = \frac{d}{dt}_{\vert t=0} r(\tilde{\gamma}(t)) = \frac{d}{dt}_{\vert t=0} d(p,\tilde{\gamma}(t)) = \langle \gamma'(t), \tilde{\gamma}'(0)\rangle,$$
where the last equality follows from the first variation formula. By the uniqueness of the gradient we gain our claim.
Also, here's another way of computing $(\nabla r)_{\gamma(t)}$ without using the first variation formula:
$$\langle \nabla r, \gamma' \rangle = \frac{d}{dt} r(\gamma(t)) = \frac{d}{dt} t = 1$$
but also by Cauchy Schwarz
$$\langle \nabla r, \gamma' \rangle \le 1 \cdot 1 = 1$$
and hence we have $\nabla r = \gamma'$.
Notice that we computed $\frac{d}{dt} r(\gamma(t))$ with the limit from below, assuming it was differentiable. So, if you are in case a) and assumed $r$ was differentiable in $x$, then you would get $\nabla r(x) = \gamma_1'(x)$ but also $\nabla r(x) = \gamma_2'(x)$, which is a contradiction.
For b): If $\gamma:[0,L] \to M^n$ is a geodesic and $x = \gamma(L)$ its first conjugated point to $p = \gamma(0)$, then the Weingarten map $A(t) = \nabla_\cdot N$ (where $N = \nabla r$ is a normed normal field along the distance spheres $S_t(p)$) has a pole in $t = L$. This is because $A(t) \cdot J(t) = J'(t)$ for $0<t<L$ and $J$ any Jacobi field along $\gamma$ with $J(0)=0$ and $J'(0) \neq 0$. But since the two points are conjugated, there exists such a $J$ with also $J(L)=0$ and $J'(L) \neq 0$ (otherwise $J\equiv 0$). Thus $$\lim\limits_{t \to L} A(t) \cdot J(t) = \lim\limits_{t \to L} J'(t) = J'(L) \neq 0,$$ but $$\lim\limits_{t \to L} J(t) = 0,$$ so $A(t)$ must blow up for $t \rightarrow L$.
Since the Hessian of $r$ is (tangentially to the distance spheres) given by the Weingarten map, the claim follows, since if $r$ were smooth at $x = \gamma(T)$ then $A(t)$ was continuous which is impossible since $\lim\limits_{t \to T} A(t)$ blows up. Thus, $r$ can't be contiuously differentiable at $x$.
A: I'd like to offer another answer on why $\mathrm{Hess}r$ blows up. Actuarally the following arguement has already been contained on Petersen's book.
Suppose $\gamma:[0,l]\to M$ be a unique minimizing geodesic connecting $\gamma(0)=p$ and $\gamma(l)$ with $\gamma(l)$ being the first conjugate point to $p$. Suppose $J$ is a nontrivial Jacobi field along $\gamma$ such that $J(0)=0=J(l)$.
If $r$ is smooth at $\gamma(0)$, then as (a) implies
$$\nabla r=\gamma'(l),$$
\begin{equation}
\begin{aligned}
\mathrm{Hess}r(J,J)&=JJr-\nabla_JJr=J\langle\gamma'(l),J\rangle-\langle\gamma'(l),\nabla_JJ\rangle\\
&=\langle \nabla_J\gamma'(l),J\rangle=\langle \nabla_{\gamma'(l)}J,J\rangle\\
&=\frac{d}{dt}\Big{|}_{t=l}\langle J,J\rangle.
\end{aligned}
\end{equation}
But also note that $|J(t)|^2\to 0$ as $t\to l$ which implies $\log|J|^2\to -\infty$ as $t\to l$. By elementary calculus there exists a sequence $\{t_k\}$ such that $$\lim\limits_{k\to\infty}\frac{2\mathrm{Hess}r(J(t_k),J(t_k))}{|J(t_k)|^2}=\lim\limits_{k\to\infty}\frac{d}{dt}\Big{|}_{t=t_k}\log |J|^2= -\infty.$$
Define a smooth vector field along $\gamma(t)$ by $$X(t)=\frac{J(t)}{|J(t)|}.$$ The norm of $X(t)$ equals $1$.
In local coordinates around $\gamma(l)$, we can assume $X(t)=\sum\limits_{i=1}^nX^i(t)\frac{\partial}{\partial x_i}\Big{|}_{\gamma(t)}$. Since the norm of $X$ is bounded, we can subtract a subsequence $\{t_{k_l}\} $of $\{t_k\}$ such that $X^i(t_{k_l})$ converges to some number $X^i$. Now $X(t_{k_l})$ converges to some $X=X^i\frac{\partial}{\partial x_i}\Big{|}_{\gamma(l)}$
However, by what we previously argued $$\lim\limits_{l\to\infty}\mathrm{Hess}r(X(t_{k_l}),X(t_{k_l}))=-\infty,$$ thus $\mathrm{Hess}r(X,X)$ is not well defined.
