[proofplan] We use geodesic polar coordinates centered at $p$ on the pre-cut domain. The exponential map is a diffeomorphism from the set of pairs $(t,v)$ with $v \in S_pM$ and $0<t<c(v)$ onto $\Omega_p$, and in these coordinates the distance from $p$ is exactly the radial coordinate $t$. Smoothness follows because the radial coordinate is smooth. The gradient formula then follows from the Gauss lemma: radial coordinate vectors are unit length and orthogonal to angular coordinate vectors. [/proofplan] [step:Pass to geodesic polar coordinates on the pre-cut domain] Let \begin{align*} D_p=\{(t,v)\in (0,\infty)\times S_pM: 0<t<c(v)\}. \end{align*} Define the polar exponential map \begin{align*} E_p: D_p &\to M \\ (t,v) &\mapsto \exp_p(tv). \end{align*} The standard structure theorem for the cut locus states that $D_p$ is open in $(0,\infty)\times S_pM$ and that \begin{align*} E_p: D_p \to \Omega_p \end{align*} is a smooth diffeomorphism. This result is obtained by combining the [inverse function theorem](/page/Inverse%20Function%20Theorem) for smooth manifolds with the fact that, before the cut time, the radial geodesic has no conjugate obstruction and is the unique minimizing geodesic from $p$ to its endpoint. We record the citation status explicitly: (citing a result not yet in the wiki: Smooth Polar Coordinate Theorem Away from the Cut Locus). Thus every $x \in \Omega_p$ has a unique representation \begin{align*} x=E_p(t,v)=\exp_p(tv) \end{align*} with $v \in S_pM$ and $0<t<c(v)$. [/step] [step:Identify the distance function with the radial coordinate] Let \begin{align*} \pi_1: D_p &\to (0,\infty) \\ (t,v) &\mapsto t \end{align*} be the first-coordinate projection. Since $E_p: D_p \to \Omega_p$ is a diffeomorphism, it is enough to compute $r_p \circ E_p$. Fix $(t,v)\in D_p$. By the definition of the cut time, the geodesic segment \begin{align*} \gamma_v|_{[0,t]}: [0,t] \to M \end{align*} is minimizing from $p$ to $\gamma_v(t)=E_p(t,v)$. Since $\gamma_v$ has unit speed, \begin{align*} d_g(p,E_p(t,v)) = d_g(p,\gamma_v(t)) = t. \end{align*} Therefore \begin{align*} (r_p\circ E_p)(t,v)=t=\pi_1(t,v). \end{align*} Hence \begin{align*} r_p|_{\Omega_p}=\pi_1\circ E_p^{-1}. \end{align*} Because $\pi_1$ and $E_p^{-1}$ are smooth maps, $r_p$ is smooth on $\Omega_p$. [/step] [step:Compute the radial and angular coordinate vectors] Fix $v\in S_pM$, fix $t$ with $0<t<c(v)$, and set \begin{align*} x=E_p(t,v)=\exp_p(tv). \end{align*} Let \begin{align*} \partial_t \in T_xM \end{align*} denote the radial coordinate vector at $x$, defined by \begin{align*} \partial_t = d(E_p)_{(t,v)}(1,0). \end{align*} By differentiating the curve \begin{align*} s \mapsto E_p(s,v)=\gamma_v(s) \end{align*} at $s=t$, we obtain \begin{align*} \partial_t=\dot{\gamma}_v(t). \end{align*} Since $\gamma_v$ is a unit-speed geodesic, \begin{align*} |\partial_t|_g = |\dot{\gamma}_v(t)|_g = 1. \end{align*} For each $\xi \in T_vS_pM$, define the corresponding angular coordinate vector \begin{align*} A_\xi = d(E_p)_{(t,v)}(0,\xi)\in T_xM. \end{align*} The Gauss lemma for the exponential map gives \begin{align*} g_x(\partial_t,A_\xi)=0 \end{align*} for every $\xi \in T_vS_pM$. We record the citation status explicitly: (citing a result not yet in the wiki: Gauss Lemma for the Exponential Map). [/step] [step:Use the Gauss lemma to identify the gradient with the radial vector] We prove that $\nabla r_p(x)=\partial_t$. Since $E_p$ is a diffeomorphism, every tangent vector $w \in T_xM$ has a unique decomposition \begin{align*} w=a\,\partial_t + A_\xi \end{align*} for some $a\in \mathbb{R}$ and some $\xi \in T_vS_pM$. Because $r_p\circ E_p=\pi_1$, the differential of $r_p$ at $x$ satisfies \begin{align*} dr_p|_x(w) &=d\pi_1|_{(t,v)}(a,\xi) \\ &=a. \end{align*} On the other hand, using $|\partial_t|_g=1$ and the Gauss lemma orthogonality, \begin{align*} g_x(\partial_t,w) &=g_x(\partial_t,a\,\partial_t+A_\xi) \\ &=a\,g_x(\partial_t,\partial_t)+g_x(\partial_t,A_\xi) \\ &=a. \end{align*} Thus \begin{align*} dr_p|_x(w)=g_x(\partial_t,w) \end{align*} for every $w\in T_xM$. By the defining property of the Riemannian gradient, this implies \begin{align*} \nabla r_p(x)=\partial_t=\dot{\gamma}_v(t). \end{align*} [guided] We now identify the gradient by testing it against an arbitrary tangent vector. The definition of the Riemannian gradient says that $\nabla r_p(x)$ is the unique vector in $T_xM$ satisfying \begin{align*} dr_p|_x(w)=g_x(\nabla r_p(x),w) \end{align*} for every $w\in T_xM$. Therefore it is enough to prove that the radial vector $\partial_t$ has this property. Since $E_p: D_p\to \Omega_p$ is a diffeomorphism, its differential \begin{align*} d(E_p)_{(t,v)}: T_{(t,v)}D_p \to T_xM \end{align*} is a linear isomorphism. The tangent space of $D_p\subset (0,\infty)\times S_pM$ at $(t,v)$ is naturally \begin{align*} T_{(t,v)}D_p \cong \mathbb{R}\times T_vS_pM. \end{align*} Hence every $w\in T_xM$ can be written uniquely in the form \begin{align*} w=d(E_p)_{(t,v)}(a,\xi)=a\,\partial_t + A_\xi, \end{align*} where $a\in\mathbb{R}$ and $\xi\in T_vS_pM$. The distance function is especially simple in these coordinates: from the previous step, \begin{align*} r_p\circ E_p=\pi_1, \end{align*} where $\pi_1(t,v)=t$. Applying the chain rule to $w=d(E_p)_{(t,v)}(a,\xi)$ gives \begin{align*} dr_p|_x(w) &=d(r_p\circ E_p)|_{(t,v)}(a,\xi) \\ &=d\pi_1|_{(t,v)}(a,\xi) \\ &=a. \end{align*} Now we compute the [inner product](/page/Inner%20Product) of $w$ with the radial vector. The radial vector is \begin{align*} \partial_t=d(E_p)_{(t,v)}(1,0)=\dot{\gamma}_v(t), \end{align*} so it has unit length because $\gamma_v$ is parametrized by arclength: \begin{align*} g_x(\partial_t,\partial_t)=1. \end{align*} The angular part is \begin{align*} A_\xi=d(E_p)_{(t,v)}(0,\xi). \end{align*} The Gauss lemma states that radial and angular coordinate vectors are orthogonal in geodesic polar coordinates: \begin{align*} g_x(\partial_t,A_\xi)=0. \end{align*} Therefore \begin{align*} g_x(\partial_t,w) &=g_x(\partial_t,a\,\partial_t+A_\xi) \\ &=a\,g_x(\partial_t,\partial_t)+g_x(\partial_t,A_\xi) \\ &=a. \end{align*} We have proved that \begin{align*} dr_p|_x(w)=g_x(\partial_t,w) \end{align*} for every $w\in T_xM$. By uniqueness in the definition of the Riemannian gradient, \begin{align*} \nabla r_p(x)=\partial_t. \end{align*} Finally, since $\partial_t=\dot{\gamma}_v(t)$, we obtain \begin{align*} \nabla r_p(x)=\dot{\gamma}_v(t). \end{align*} [/guided] [/step] [step:Conclude the smoothness and gradient formula] The previous steps show that $r_p|_{\Omega_p}=\pi_1\circ E_p^{-1}$, so $r_p$ is smooth on $\Omega_p$. For every $v\in S_pM$ and every $t$ with $0<t<c(v)$, the point $x=\exp_p(tv)$ lies in $\Omega_p$, and the gradient computation gives \begin{align*} \nabla r_p(x)=\dot{\gamma}_v(t). \end{align*} This is exactly the asserted formula. [/step]