[proofplan] The proof is a direct comparison of definitions. On the open normal domain, the distance function is smooth and its gradient $\nu=\nabla r$ is the outward unit normal to the geodesic spheres. For tangent vectors to a geodesic sphere, the Hessian of $r$ is exactly the metric pairing of $\nabla_X\nu$ with the second tangent vector, which is the definition of the shape operator used here. The radial vanishing follows by differentiating the identity $|\nabla r|_g^2=1$ and translating the result into the Hessian. [/proofplan] [step:Identify the Hessian in tangent directions with the derivative of the unit normal] Fix $x \in \Omega_p$ and $X,Y \in T_xS(p,r(x))$. Let $\widetilde X$ and $\widetilde Y$ be smooth local vector fields along $S(p,r(x))$ such that $\widetilde X_x=X$ and $\widetilde Y_x=Y$. Since $\nu=\nabla r$ on $\Omega_p$, the definition of the Hessian of a smooth function gives \begin{align*} \operatorname{Hess} r_x(X,Y) &= g_x\bigl(\nabla_{\widetilde X}\nabla r,\widetilde Y\bigr) \\ &= g_x\bigl(\nabla_{\widetilde X}\nu,\widetilde Y\bigr). \end{align*} By the definition of the shape operator $A_x$ with respect to the outward normal $\nu$, \begin{align*} A_xX=\nabla_{\widetilde X}\nu\big|_x. \end{align*} Therefore \begin{align*} \operatorname{Hess} r_x(X,Y)=g_x(A_xX,Y). \end{align*} [guided] Fix a point $x \in \Omega_p$ and tangent vectors $X,Y \in T_xS(p,r(x))$. To evaluate the Hessian at tangent vectors, choose smooth local tangent extensions $\widetilde X$ and $\widetilde Y$ along the geodesic sphere $S(p,r(x))$ such that $\widetilde X_x=X$ and $\widetilde Y_x=Y$. The key observation is that the outward unit normal field to the geodesic spheres is the gradient of the distance function: \begin{align*} \nu=\nabla r. \end{align*} The Hessian of a smooth function $r$ is defined by covariantly differentiating its gradient and pairing with the metric. Hence \begin{align*} \operatorname{Hess} r_x(X,Y) &= g_x\bigl(\nabla_{\widetilde X}\nabla r,\widetilde Y\bigr) \\ &= g_x\bigl(\nabla_{\widetilde X}\nu,\widetilde Y\bigr). \end{align*} Now compare this with the definition of the shape operator. With the sign convention in the statement, the shape operator of the geodesic sphere at $x$ is \begin{align*} A_xX=\nabla_{\widetilde X}\nu\big|_x. \end{align*} Substituting this definition into the preceding identity gives \begin{align*} \operatorname{Hess} r_x(X,Y)=g_x(A_xX,Y). \end{align*} Thus the Hessian of the distance function restricted to tangent directions is the second fundamental form of the geodesic sphere, written through the shape operator. [/guided] [/step] [step:Differentiate the unit gradient identity to obtain radial vanishing] Fix $x \in \Omega_p$ and $Z \in T_xM$. Let $\widetilde Z$ be a smooth local vector field on $\Omega_p$ with $\widetilde Z_x=Z$. Since $|\nabla r|_g^2=1$ on $\Omega_p$, differentiating this identity in the direction $\widetilde Z$ gives \begin{align*} 0 &= \widetilde Z\bigl(g(\nabla r,\nabla r)\bigr)\big|_x. \end{align*} Using metric compatibility of the Levi-Civita connection, \begin{align*} \widetilde Z\bigl(g(\nabla r,\nabla r)\bigr) &= g(\nabla_{\widetilde Z}\nabla r,\nabla r)+g(\nabla r,\nabla_{\widetilde Z}\nabla r) \\ &= 2g(\nabla_{\widetilde Z}\nabla r,\nabla r). \end{align*} Evaluating at $x$ yields \begin{align*} 0=2g_x(\nabla_{\widetilde Z}\nabla r,\nabla r(x)) =2\operatorname{Hess} r_x(Z,\nabla r(x)). \end{align*} Therefore \begin{align*} \operatorname{Hess} r_x(Z,\nabla r(x))=0. \end{align*} Since the Levi-Civita connection is torsion-free, the Hessian of a smooth function is symmetric, so \begin{align*} \operatorname{Hess} r_x(\nabla r(x),Z)=0. \end{align*} This proves the radial vanishing identity. [/step]