[proofplan]
We compute the Hessian of
\begin{align*}
\frac{1}{2}r^2
\end{align*}
by varying the endpoint $x$ and taking the geodesics from $S$ to the moving endpoint. The associated variational field is a Jacobi field satisfying the natural moving-boundary condition at $S$, and the second variation of energy gives the index form; total geodesy of $S$ removes the boundary term. The relation between the $[0,1]$ and arclength parametrisations is then a direct change of variables. Finally, Fermi coordinates along the totally geodesic submanifold show that, on each relatively compact tubular subneighbourhood, $r^2$ agrees with the squared normal coordinate up to second-order metric errors with uniformly bounded constants, yielding the stated local lower bound.
[/proofplan]
[step:Represent nearby endpoints by normal geodesics from $S$]
Because $U$ is a tubular neighbourhood on which $\pi:U\to S$ is smooth and uniquely defined, every $y\in U$ has a unique representation
\begin{align*}
y=\exp_{\pi(y)}(\rho(y)\, n(y)),
\end{align*}
where $\rho(y)=r(y)$, $n(y)\in \nu_{\pi(y)}(S)$, and $|n(y)|_g=1$ whenever $y\notin S$.
Fix $x\in U\setminus S$ and $X\in T_xM$. Choose the geodesic
\begin{align*}
c:(-\varepsilon,\varepsilon)&\to U \\
\tau&\mapsto \exp_x(\tau X)
\end{align*}
with $\varepsilon>0$ small enough that $c(( -\varepsilon,\varepsilon))\subset U$. Then $c(0)=x$, $\dot c(0)=X$, and $\nabla_{\dot c}\dot c(0)=0$.
Define the maps
\begin{align*}
s:(-\varepsilon,\varepsilon)&\to S, & \tau&\mapsto \pi(c(\tau)),\\
\ell:(-\varepsilon,\varepsilon)&\to (0,\infty), & \tau&\mapsto r(c(\tau)),\\
n:(-\varepsilon,\varepsilon)&\to \nu(S), & \tau&\mapsto \ell(\tau)^{-1}\exp_{s(\tau)}^{-1}(c(\tau)).
\end{align*}
After shrinking $\varepsilon>0$, these maps are smooth, with $n(\tau)\in \nu_{s(\tau)}(S)$ and $|n(\tau)|_g=1$.
Define the geodesic variation
\begin{align*}
F:(-\varepsilon,\varepsilon)\times[0,1]&\to M \\
(\tau,u)&\mapsto \exp_{s(\tau)}(u\ell(\tau)n(\tau)).
\end{align*}
For each $\tau$, the curve $u\mapsto F(\tau,u)$ is the unique normal geodesic from $S$ to $c(\tau)$, parametrised on $[0,1]$. Its variational field along $\alpha(u)=F(0,u)$ is
\begin{align*}
J:[0,1]&\to TM\\
u&\mapsto \partial_\tau F(0,u).
\end{align*}
Since $F(\tau,1)=c(\tau)$, we have $J(1)=X$. Since $F(\tau,0)=s(\tau)\in S$, we have $J(0)=\dot s(0)\in T_sS$.
[guided]
The nearest-point projection gives a controlled way to vary $x$: rather than choosing arbitrary curves from $S$ to $c(\tau)$, we choose the distinguished shortest normal geodesics. We also choose the endpoint curve itself to be the geodesic $c(\tau)=\exp_x(\tau X)$, so $\nabla_{\dot c}\dot c(0)=0$; this is the condition that removes the acceleration term in the Hessian computation below. Because $U$ lies inside the tubular neighbourhood where this representation is unique, each point $c(\tau)$ can be written as
\begin{align*}
c(\tau)=\exp_{s(\tau)}(\ell(\tau)n(\tau)),
\end{align*}
where $s:(-\varepsilon,\varepsilon)\to S$ is $s(\tau)=\pi(c(\tau))$, $\ell:(-\varepsilon,\varepsilon)\to(0,\infty)$ is $\ell(\tau)=r(c(\tau))$, and $n:(-\varepsilon,\varepsilon)\to\nu(S)$ is the unit normal field satisfying $n(\tau)\in \nu_{s(\tau)}(S)$.
Now define
\begin{align*}
F:(-\varepsilon,\varepsilon)\times[0,1]&\to M \\
(\tau,u)&\mapsto \exp_{s(\tau)}(u\ell(\tau)n(\tau)).
\end{align*}
For fixed $\tau$, the map $u\mapsto F(\tau,u)$ is a geodesic because it is obtained from the exponential map along a fixed normal vector. It starts at $s(\tau)\in S$ and ends at $c(\tau)$:
\begin{align*}
F(\tau,0)=s(\tau),\qquad F(\tau,1)=c(\tau).
\end{align*}
The variational field along the central geodesic $\alpha(u)=F(0,u)$ is
\begin{align*}
J(u)=\partial_\tau F(0,u).
\end{align*}
Differentiating the endpoint identity $F(\tau,1)=c(\tau)$ at $\tau=0$ gives
\begin{align*}
J(1)=\partial_\tau F(0,1)=\dot c(0)=X.
\end{align*}
Differentiating $F(\tau,0)=s(\tau)\in S$ gives
\begin{align*}
J(0)=\dot s(0)\in T_sS.
\end{align*}
Thus the variation produces precisely a Jacobi field with the required endpoint value and with initial value tangent to $S$.
[/guided]
[/step]
[step:Verify the moving-boundary transversality condition]
For every $\tau$, the initial velocity
\begin{align*}
\partial_uF(\tau,0)=\ell(\tau)n(\tau)
\end{align*}
lies in $\nu_{s(\tau)}(S)$. Hence, for every smooth tangent field $Y(\tau)\in T_{s(\tau)}S$ along $s(\tau)$,
\begin{align*}
\langle \partial_uF(\tau,0),Y(\tau)\rangle_g=0.
\end{align*}
Differentiate at $\tau=0$. Using the symmetry of covariant derivatives in a smooth two-parameter variation,
\begin{align*}
\nabla_{\partial_\tau}\partial_uF=\nabla_{\partial_u}\partial_\tau F,
\end{align*}
we obtain
\begin{align*}
0
&=
\frac{d}{d\tau}\bigg|_{\tau=0}
\langle \partial_uF(\tau,0),Y(\tau)\rangle_g \\
&=
\langle \nabla_{\dot\alpha}J(0),Y(0)\rangle_g
+
\langle \dot\alpha(0),\nabla_{\dot s}Y\rangle_g .
\end{align*}
Since $S$ is totally geodesic, $\nabla_{\dot s}Y\in T_sS$ whenever $Y$ is tangent to $S$, while $\dot\alpha(0)=\ell n\in \nu_s(S)$. Thus the second [inner product](/page/Inner%20Product) vanishes. Since $Y(0)\in T_sS$ is arbitrary,
\begin{align*}
\bigl(\nabla_{\dot\alpha}J(0)\bigr)^\top=0.
\end{align*}
[/step]
[step:Apply the second variation of energy to obtain the index form]
Let $\mathcal{L}^1$ denote one-dimensional [Lebesgue measure](/page/Lebesgue%20Measure) on $[0,1]$. Define the energy map
\begin{align*}
E:(-\varepsilon,\varepsilon)&\to \mathbb{R} \\
\tau&\mapsto \frac{1}{2}\int_0^1 |\partial_uF(\tau,u)|_g^2\,d\mathcal{L}^1(u).
\end{align*}
Because $u\mapsto F(\tau,u)$ is the constant-speed minimizing normal geodesic from $S$ to $c(\tau)$, its speed is $\ell(\tau)=r(c(\tau))$, and hence
\begin{align*}
E(\tau)=\frac{1}{2}r(c(\tau))^2.
\end{align*}
Therefore, using $\nabla_{\dot c}\dot c(0)=0$ for the chosen geodesic $c$,
\begin{align*}
E''(0)
&=\frac{d^2}{d\tau^2}\bigg|_{\tau=0}\left(\frac{1}{2}r(c(\tau))^2\right)\\
&=\operatorname{Hess}_x\left(\frac{1}{2}r^2\right)(X,X)
+d\left(\frac{1}{2}r^2\right)_x\bigl(\nabla_{\dot c}\dot c(0)\bigr)\\
&=\operatorname{Hess}_x\left(\frac{1}{2}r^2\right)(X,X)
=\frac{1}{2}\operatorname{Hess}_x(r^2)(X,X).
\end{align*}
The [second variation formula](/theorems/2729) for energy with the initial endpoint constrained to lie in $S$ and the terminal endpoint equal to $c(\tau)$ gives
\begin{align*}
E''(0)
=
\int_0^1
\left(
|\nabla_{\dot\alpha}J|_g^2
-
\langle R(J,\dot\alpha)\dot\alpha,J\rangle_g
\right)
\,d\mathcal{L}^1(u)
-
\langle A_{\dot\alpha(0)}J(0),J(0)\rangle_g
+
\langle \nabla_{\dot c}\dot c(0),\dot\alpha(1)\rangle_g,
\end{align*}
where
\begin{align*}
A_{\dot\alpha(0)}:T_sS&\to T_sS
\end{align*}
is the shape operator of $S$ in the normal direction $\dot\alpha(0)$. Let
\begin{align*}
\mathrm{II}:T_sS\times T_sS&\to \nu_s(S)
\end{align*}
denote the second fundamental form of $S$ at $s$. We use the sign convention
\begin{align*}
\langle A_\eta V,W\rangle_g=\langle \mathrm{II}(V,W),\eta\rangle_g
\end{align*}
for $\eta\in\nu_s(S)$ and $V,W\in T_sS$. This is the standard [constrained-endpoint second variation formula](/page/Second%20Variation) for a smooth geodesic variation; its hypotheses hold because $F$ is smooth, $\alpha$ is a geodesic, $J(0)\in T_sS$, and the transversality condition $(\nabla_{\dot\alpha}J(0))^\top=0$ was verified above. The terminal boundary term vanishes because $c$ is a geodesic. Since $S$ is totally geodesic, its second fundamental form vanishes, so $A_{\dot\alpha(0)}=0$. Hence
\begin{align*}
\frac{1}{2}\operatorname{Hess}_x(r^2)(X,X)
=
\int_0^1
\left(
|\nabla_{\dot\alpha}J|_g^2
-
\langle R(J,\dot\alpha)\dot\alpha,J\rangle_g
\right)
\,d\mathcal{L}^1(u)
=
I_\alpha(J,J).
\end{align*}
[guided]
The function
\begin{align*}
\frac{1}{2}r^2
\end{align*}
is naturally an energy. For each $\tau$, the curve $u\mapsto F(\tau,u)$ has constant speed $\ell(\tau)=r(c(\tau))$, so
\begin{align*}
E(\tau)
&=
\frac{1}{2}\int_0^1 |\partial_uF(\tau,u)|_g^2\,d\mathcal{L}^1(u)\\
&=
\frac{1}{2}\int_0^1 \ell(\tau)^2\,d\mathcal{L}^1(u)\\
&=
\frac{1}{2}r(c(\tau))^2.
\end{align*}
Thus differentiating twice at $\tau=0$ computes the Hessian because $c$ is a geodesic. For a general curve one would have an additional first-derivative term, namely
\begin{align*}
\frac{d^2}{d\tau^2}\bigg|_{\tau=0}\left(\frac{1}{2}r(c(\tau))^2\right)
=
\operatorname{Hess}_x\left(\frac{1}{2}r^2\right)(X,X)
+d\left(\frac{1}{2}r^2\right)_x\bigl(\nabla_{\dot c}\dot c(0)\bigr).
\end{align*}
Our chosen curve has $\nabla_{\dot c}\dot c(0)=0$, so
\begin{align*}
E''(0)=\operatorname{Hess}_x\left(\frac{1}{2}r^2\right)(X,X)
=
\frac{1}{2}\operatorname{Hess}_x(r^2)(X,X).
\end{align*}
We now use the standard [constrained-endpoint second variation formula](/page/Second%20Variation) for energy for geodesic variations whose initial point moves inside the submanifold $S$ and whose terminal point is the prescribed endpoint curve $c(\tau)$. We check its hypotheses in this setting. The map $F$ is smooth by smoothness of the tubular parametrisation, the central curve $\alpha=F(0,\cdot)$ is a geodesic, the variational field is $J=\partial_\tau F(0,\cdot)$, the initial value satisfies $J(0)\in T_sS$, and the preceding step proved the transversality condition $(\nabla_{\dot\alpha}J(0))^\top=0$. The initial velocity $\dot\alpha(0)$ is normal to $S$, and the terminal acceleration is $\nabla_{\dot c}\dot c(0)$. The formula gives
\begin{align*}
E''(0)
=
\int_0^1
\left(
|\nabla_{\dot\alpha}J|_g^2
-
\langle R(J,\dot\alpha)\dot\alpha,J\rangle_g
\right)
\,d\mathcal{L}^1(u)
-
\langle A_{\dot\alpha(0)}J(0),J(0)\rangle_g
+
\langle \nabla_{\dot c}\dot c(0),\dot\alpha(1)\rangle_g.
\end{align*}
Here
\begin{align*}
A_{\dot\alpha(0)}:T_sS\to T_sS
\end{align*}
is the shape operator of $S$ in the normal direction $\dot\alpha(0)$. The term involving $A_{\dot\alpha(0)}$ is the correction from allowing the initial point to slide along $S$, while the final inner product is the correction from the moving terminal endpoint.
Both boundary corrections vanish in the present situation. The terminal correction vanishes because $c$ is a geodesic:
\begin{align*}
\langle \nabla_{\dot c}\dot c(0),\dot\alpha(1)\rangle_g=0.
\end{align*}
The hypothesis that $S$ is totally geodesic is also used. A totally geodesic submanifold has vanishing second fundamental form, so every shape operator of $S$ is the zero map. Therefore
\begin{align*}
\langle A_{\dot\alpha(0)}J(0),J(0)\rangle_g=0.
\end{align*}
Consequently,
\begin{align*}
\frac{1}{2}\operatorname{Hess}_x(r^2)(X,X)
=
\int_0^1
\left(
|\nabla_{\dot\alpha}J|_g^2
-
\langle R(J,\dot\alpha)\dot\alpha,J\rangle_g
\right)
\,d\mathcal{L}^1(u),
\end{align*}
which is the claimed index form identity.
[/guided]
[/step]
[step:Convert the index form to arclength parametrisation]
Let
\begin{align*}
\beta:[0,\ell]&\to M\\
t&\mapsto \exp_s(tn)
\end{align*}
be the arclength parametrisation of the same geodesic, and define
\begin{align*}
\widetilde J:[0,\ell]&\to TM\\
t&\mapsto J(t/\ell).
\end{align*}
Then $\alpha(u)=\beta(\ell u)$ and $\dot\alpha(u)=\ell\dot\beta(\ell u)$. Also,
\begin{align*}
\nabla_{\dot\alpha}J(u)=\ell\,\nabla_{\dot\beta}\widetilde J(\ell u).
\end{align*}
Substituting $t=\ell u$, so that $d\mathcal{L}^1(t)=\ell\,d\mathcal{L}^1(u)$, gives
\begin{align*}
I_\alpha(J,J)
&=
\int_0^1
\left(
\ell^2|\nabla_{\dot\beta}\widetilde J(\ell u)|_g^2
-
\ell^2\langle R(\widetilde J(\ell u),\dot\beta(\ell u))\dot\beta(\ell u),\widetilde J(\ell u)\rangle_g
\right)
\,d\mathcal{L}^1(u)\\
&=
\ell
\int_0^\ell
\left(
|\nabla_{\dot\beta}\widetilde J(t)|_g^2
-
\langle R(\widetilde J(t),\dot\beta(t))\dot\beta(t),\widetilde J(t)\rangle_g
\right)
\,d\mathcal{L}^1(t)\\
&=
\ell\, I_\beta(\widetilde J,\widetilde J).
\end{align*}
[/step]
[step:Compare the Hessian with the product model in Fermi coordinates]
Let $U'\Subset U$ be a relatively compact tubular subneighbourhood. Choose a relatively compact coordinate patch $V\subset S$ and an orthonormal local frame
\begin{align*}
e_1,\dots,e_k
\end{align*}
for the normal bundle $\nu(S)$ over $V$. Define Fermi coordinates
\begin{align*}
\Phi: V\times B_{\mathbb{R}^k}(0,\delta)&\to M\\
(s,z)&\mapsto \exp_s\left(\sum_{a=1}^k z_a e_a(s)\right).
\end{align*}
After shrinking $\delta>0$, $\Phi$ is a diffeomorphism onto its image and
\begin{align*}
r(\Phi(s,z))^2=|z|^2.
\end{align*}
In these coordinates, the normal-fibre component of a tangent vector corresponds to its $z$-component.
We use the [Fermi-coordinate expansion](/page/Fermi%20Coordinates) near a totally geodesic submanifold. Along $z=0$, the metric coefficients agree with the product metric, the tangential first normal derivatives vanish because the second fundamental form of $S$ is zero, and the Christoffel symbols of the product model have no normal component. More precisely, writing $g_{AB}(s,z)$ for the coordinate coefficients of $g$ and $g^{\mathrm{prod}}_{AB}(s,z)$ for the product coefficients determined by $g|_S$ and the chosen orthonormal normal frame, there are constants $K_V>0$ and $\delta_V>0$ such that, for $|z|<\delta_V$,
\begin{align*}
|g_{AB}(s,z)-g^{\mathrm{prod}}_{AB}(s,z)|\leq K_V |z|^2,
\qquad
|\partial_C g_{AB}(s,z)-\partial_C g^{\mathrm{prod}}_{AB}(s,z)|\leq K_V |z|
\end{align*}
for all coordinate indices $A,B,C$. In addition, the normal Christoffel contraction satisfies
\begin{align*}
\left|\sum_{c=1}^k \Gamma_{AB}^{c}(s,z)z_c\right|\leq K_V |z|^2
\end{align*}
for all $A,B$. For tangential-tangential indices this follows from total geodesy, which gives $\Gamma_{ij}^{c}(s,0)=0$; for mixed and normal-normal indices it follows from the Fermi-coordinate normalisation along each normal fibre, which gives $\Gamma_{ia}^{c}(s,0)=0$ and $\Gamma_{ab}^{c}(s,z)=0$ along the fibre. The constant $K_V$ is finite because the patch is relatively compact and these estimates depend only on bounded curvature and bounded first covariant derivatives of the relevant coordinate coefficients on the patch.
Let $Y\in T_{\Phi(s,z)}M$, and write its coordinate representation as
\begin{align*}
Y=\sum_A Y^A\partial_A.
\end{align*}
For the function $f:V\times B_{\mathbb{R}^k}(0,\delta_V)\to\mathbb{R}$ given by $f(s,z)=|z|^2$, the coordinate Hessian formula gives
\begin{align*}
\operatorname{Hess}_g(f)(Y,Y)
&=
\sum_{A,B}Y^A Y^B\left(\partial_A\partial_B f-\sum_C\Gamma_{AB}^{C}\partial_C f\right)\\
&=
2\sum_{a=1}^k (Y^a)^2
-
2\sum_{A,B}Y^A Y^B\sum_{c=1}^k\Gamma_{AB}^{c}(s,z)z_c.
\end{align*}
The Christoffel contraction estimate gives
\begin{align*}
\left|
2\sum_{A,B}Y^A Y^B\sum_{c=1}^k\Gamma_{AB}^{c}(s,z)z_c
\right|
\leq C_{1,V}|z|^2|Y|_g^2,
\end{align*}
where $C_{1,V}>0$ also includes the uniform equivalence between coordinate component norms and the Riemannian norm on the relatively compact coordinate patch. The metric comparison with the product model gives another constant $C_{2,V}>0$ such that
\begin{align*}
\left|2\sum_{a=1}^k (Y^a)^2-2|Y^\perp|_g^2\right|
\leq C_{2,V}|z|^2|Y|_g^2.
\end{align*}
Defining
\begin{align*}
C_V:=C_{1,V}+C_{2,V},
\end{align*}
we obtain throughout $V\times B_{\mathbb{R}^k}(0,\delta_V)$
\begin{align*}
\left|
\operatorname{Hess}_g(|z|^2)(Y,Y)-2|Y^\perp|_g^2
\right|
\leq C_V |z|^2 |Y|_g^2.
\end{align*}
Since $r(\Phi(s,z))=|z|$, this gives
\begin{align*}
\operatorname{Hess}_g(r^2)(Y,Y)
\geq
2|Y^\perp|_g^2-C_V r(\Phi(s,z))^2|Y|_g^2.
\end{align*}
The set $\overline{U'}$ is compact. Hence its projection to $S$ is compact, and finitely many of the relatively compact Fermi-coordinate patches above cover $\overline{U'}$. Taking
\begin{align*}
C_{U'}:=\max_V C_V
\end{align*}
over this finite cover, with each $C_V$ already incorporating the coordinate-to-Riemannian norm comparison on its patch, gives
\begin{align*}
\operatorname{Hess}_x(r^2)(X,X)
\geq
2|X^\perp|_g^2-C_{U'}r(x)^2|X|_g^2
\end{align*}
for every $x\in U'$ and $X\in T_xM$.
In the exact product model $S\times\mathbb{R}^k$, the metric is the product metric and $r(s,z)^2=|z|^2$. Therefore the Hessian has no mixed or tangential contribution:
\begin{align*}
\operatorname{Hess}(r^2)(X,X)=2|X^\perp|^2.
\end{align*}
This proves both the local lower bound on relatively compact tubular subneighbourhoods and the equality in the product case.
[/step]
