[proofplan]
We prove the comparison along each radial minimizing geodesic before the cut locus. On the smooth distance region, the Hessian of $r_p$ restricted to the orthogonal complement of the radial direction is the shape operator of geodesic spheres, and its trace is $\Delta r_p$. The Jacobi field equation gives a matrix Riccati equation for this shape operator; taking traces and using Cauchy-Schwarz plus the Ricci lower bound gives a scalar Riccati inequality. The model quantity $(n-1)\operatorname{ct}_k$ satisfies the corresponding equality with the same leading singular term near $0$, so a one-dimensional comparison argument yields the desired inequality.
[/proofplan]
[step:Reduce the pointwise inequality to a radial geodesic before the cut locus]
Fix $x \in \Omega_p$ with $r_p(x) \in I_k$, and set
\begin{align*}
r_0 := r_p(x).
\end{align*}
Since $x \in \Omega_p$, there is a unique unit-speed minimizing geodesic
\begin{align*}
\gamma: [0,r_0] &\to M \\
t &\mapsto \gamma(t)
\end{align*}
such that $\gamma(0)=p$, $\gamma(r_0)=x$, and $\gamma'(0)=\xi$ for some $\xi \in T_pM$ with $|\xi|_g=1$. Moreover $\gamma(t) \in \Omega_p$ for every $t \in (0,r_0]$, and
\begin{align*}
r_p(\gamma(t)) = t,
\qquad
\nabla r_p(\gamma(t)) = \gamma'(t)
\end{align*}
for every $t \in (0,r_0]$.
For $t \in (0,r_0]$, define the orthogonal hyperplane
\begin{align*}
E_t := \{V \in T_{\gamma(t)}M : g_{\gamma(t)}(V,\gamma'(t)) = 0\}.
\end{align*}
Let $\operatorname{Hess} r_p$ denote the symmetric covariant $2$-tensor on $\Omega_p$ defined by
\begin{align*}
(\operatorname{Hess} r_p)_y(V,W) := g_y(\nabla_V \nabla r_p,W)
\end{align*}
for $y\in\Omega_p$ and $V,W\in T_yM$. Define the shape operator of the geodesic sphere centered at $p$ along $\gamma$ by
\begin{align*}
S(t): E_t &\to E_t \\
V &\mapsto \nabla_V \nabla r_p.
\end{align*}
The map $S(t)$ is self-adjoint with respect to $g_{\gamma(t)}$. Define the radial Laplacian function
\begin{align*}
m:(0,r_0] &\to \mathbb{R} \\
t &\mapsto \operatorname{tr}_{E_t} S(t)=\Delta r_p(\gamma(t)).
\end{align*}
Thus it suffices to prove
\begin{align*}
m(r_0) \leq (n-1)\operatorname{ct}_k(r_0).
\end{align*}
[/step]
[step:Derive the trace Riccati inequality for the radial Laplacian]
Parallel transport along $\gamma$ identifies all $E_t$ with the fixed Euclidean space $E_0 := \xi^\perp \subset T_pM$. Under this identification, view $S(t)$ as a self-adjoint endomorphism
\begin{align*}
S(t): E_0 &\to E_0.
\end{align*}
Define the radial curvature endomorphism
\begin{align*}
R_\gamma(t): E_0 &\to E_0
\end{align*}
by the identity
\begin{align*}
g_{\gamma(t)}(R_\gamma(t)V,W)
=
g_{\gamma(t)}(R(V_t,\gamma'(t))\gamma'(t),W_t),
\end{align*}
where $V_t,W_t \in E_t$ are the parallel transports of $V,W \in E_0$ along $\gamma$ and $R$ is the Riemann curvature tensor with convention
\begin{align*}
R(X,Y)Z = \nabla_X\nabla_YZ-\nabla_Y\nabla_XZ-\nabla_{[X,Y]}Z.
\end{align*}
By the cut-locus characterization recorded on the [Cut Locus](/page/Cut%20Locus) page, applied to the complete Riemannian manifold $(M,g)$, the point $p$, and the unit initial vector $\xi$, the condition $x \in \Omega_p$ implies that the segment $\gamma|_{(0,r_0]}$ contains no cut point of $p$ and no conjugate point to $p$ along $\gamma$. Equivalently, the differential of the exponential map in directions orthogonal to $\xi$ is nonsingular for every $t \in (0,r_0]$. Define the Jacobi tensor
\begin{align*}
J:(0,r_0] &\to \mathcal{L}(E_0,E_0)
\end{align*}
by identifying, through parallel transport back to $E_0$, the Jacobi field $J_V(t)$ with initial conditions $J_V(0)=0$ and $\nabla_{\gamma'}J_V(0)=V$ for each $V\in E_0$. The nonsingularity of the orthogonal differential of $\exp_p$ implies that $J(t)$ is invertible for every $t\in(0,r_0]$. The shape operator is the logarithmic derivative
\begin{align*}
S(t)=J'(t)J(t)^{-1},
\end{align*}
where $J'(t)$ denotes covariant differentiation along $\gamma$ after the same parallel-transport identification. The Jacobi field equation for radial variations, as recorded on the [Jacobi Field](/page/Jacobi%20Field) page, is
\begin{align*}
J''(t)+R_\gamma(t)J(t)=0.
\end{align*}
Differentiating $S(t)=J'(t)J(t)^{-1}$ and using $(J^{-1})'=-J^{-1}J'J^{-1}$ gives
\begin{align*}
S'(t)
&=J''(t)J(t)^{-1}-J'(t)J(t)^{-1}J'(t)J(t)^{-1} \\
&=-R_\gamma(t)-S(t)^2.
\end{align*}
Thus the matrix Riccati equation for the shape operator of geodesic spheres is
\begin{align*}
S'(t) + S(t)^2 + R_\gamma(t) = 0.
\end{align*}
Taking the trace over the $(n-1)$-dimensional space $E_0$ gives
\begin{align*}
m'(t) + \operatorname{tr}_{E_0}(S(t)^2) + \operatorname{Ric}_{\gamma(t)}(\gamma'(t),\gamma'(t)) = 0.
\end{align*}
Since $S(t)$ is self-adjoint, the [Cauchy-Schwarz Inequality](/page/Cauchy-Schwarz%20Inequality) for its eigenvalues gives
\begin{align*}
\operatorname{tr}_{E_0}(S(t)^2) \geq \frac{m(t)^2}{n-1}.
\end{align*}
The geodesic has unit speed, so $|\gamma'(t)|_g=1$, and the Ricci lower bound gives
\begin{align*}
\operatorname{Ric}_{\gamma(t)}(\gamma'(t),\gamma'(t)) \geq (n-1)k.
\end{align*}
Combining these estimates yields the scalar inequality
\begin{align*}
m'(t) + \frac{m(t)^2}{n-1} + (n-1)k \leq 0
\end{align*}
for every $t \in (0,r_0]$.
[guided]
The geometric quantity we want to compare is $\Delta r_p$, and along the radial geodesic $\gamma$ this becomes a one-variable function
\begin{align*}
m(t) := \Delta r_p(\gamma(t)).
\end{align*}
The reason for introducing the shape operator is that $\Delta r_p$ is the trace of the Hessian of $r_p$. Here $\operatorname{Hess} r_p$ denotes the symmetric covariant $2$-tensor on $\Omega_p$ defined by
\begin{align*}
(\operatorname{Hess} r_p)_y(V,W) := g_y(\nabla_V \nabla r_p,W)
\end{align*}
for $y\in\Omega_p$ and $V,W\in T_yM$. Since $\nabla r_p=\gamma'$ along $\gamma$ and $|\nabla r_p|_g=1$, the radial component of the Hessian vanishes:
\begin{align*}
\operatorname{Hess} r_p(\gamma'(t),\gamma'(t)) = \frac{1}{2}\gamma'(t)\bigl(|\nabla r_p|_g^2\bigr)=0.
\end{align*}
Therefore the full trace of the Hessian of $r_p$ equals the trace over the orthogonal complement $E_t=\gamma'(t)^\perp$:
\begin{align*}
\Delta r_p(\gamma(t)) = \operatorname{tr}_{E_t}\bigl(V \mapsto \nabla_V\nabla r_p\bigr).
\end{align*}
To differentiate this operator in one fixed [vector space](/page/Vector%20Space), parallel transport $E_0=\xi^\perp$ along $\gamma$ and use it to identify $E_t$ with $E_0$. Under this identification the shape operator is a family of self-adjoint maps
\begin{align*}
S(t): E_0 &\to E_0.
\end{align*}
Because $x \in \Omega_p$, the segment from $p$ to $x$ stays before the cut locus. By the cut-locus characterization recorded on the [Cut Locus](/page/Cut%20Locus) page, applied to the complete Riemannian manifold $(M,g)$, the point $p$, and the unit initial vector $\xi$, for every $t \in (0,r_0]$ the radial geodesic is still minimizing and has no conjugate point to $p$ along $\gamma$; equivalently, the differential of the exponential map in directions orthogonal to $\xi$ is nonsingular. This nonsingularity is the hypothesis needed to pass from Jacobi fields to a logarithmic derivative: the Jacobi tensor along $\gamma$ is invertible on $E_0$, and the shape operator $S(t)$ is that logarithmic derivative. The Jacobi equation for radial variations says that this tensor satisfies a second-order equation involving the radial curvature endomorphism. Written in logarithmic derivative form, this is exactly the matrix Riccati equation
\begin{align*}
S'(t)+S(t)^2+R_\gamma(t)=0,
\end{align*}
where
\begin{align*}
g_{\gamma(t)}(R_\gamma(t)V,W)
=
g_{\gamma(t)}(R(V_t,\gamma'(t))\gamma'(t),W_t)
\end{align*}
for parallel fields $V_t,W_t$ along $\gamma$.
Now take traces. The trace of $S'(t)$ is $m'(t)$, the trace of $S(t)^2$ remains $\operatorname{tr}_{E_0}(S(t)^2)$, and the trace of $R_\gamma(t)$ over the hyperplane perpendicular to $\gamma'(t)$ is the Ricci curvature in the radial direction:
\begin{align*}
m'(t)+\operatorname{tr}_{E_0}(S(t)^2)
+\operatorname{Ric}_{\gamma(t)}(\gamma'(t),\gamma'(t))=0.
\end{align*}
This is the exact point where the hypothesis on Ricci curvature enters. Since $\gamma$ is unit speed, $|\gamma'(t)|_g=1$, and hence
\begin{align*}
\operatorname{Ric}_{\gamma(t)}(\gamma'(t),\gamma'(t))\geq (n-1)k.
\end{align*}
It remains to replace the full operator $S(t)$ by its trace $m(t)$. Because $S(t)$ is self-adjoint on the $(n-1)$-dimensional [inner product space](/page/Inner%20Product%20Space) $E_0$, choose an orthonormal eigenbasis with eigenvalues $\lambda_1(t),\dots,\lambda_{n-1}(t)$. Then
\begin{align*}
m(t)=\sum_{i=1}^{n-1}\lambda_i(t),
\qquad
\operatorname{tr}_{E_0}(S(t)^2)=\sum_{i=1}^{n-1}\lambda_i(t)^2.
\end{align*}
The [Cauchy-Schwarz Inequality](/page/Cauchy-Schwarz%20Inequality) applied to the vectors $(\lambda_1(t),\dots,\lambda_{n-1}(t))$ and $(1,\dots,1)$ gives
\begin{align*}
m(t)^2
=
\left(\sum_{i=1}^{n-1}\lambda_i(t)\right)^2
\leq
(n-1)\sum_{i=1}^{n-1}\lambda_i(t)^2.
\end{align*}
Thus
\begin{align*}
\operatorname{tr}_{E_0}(S(t)^2)\geq \frac{m(t)^2}{n-1}.
\end{align*}
Substituting this and the Ricci lower bound into the trace Riccati identity gives
\begin{align*}
m'(t)+\frac{m(t)^2}{n-1}+(n-1)k\leq 0.
\end{align*}
[/guided]
[/step]
[step:Identify the model solution and its initial singular behaviour]
Define
\begin{align*}
m_k: I_k &\to \mathbb{R} \\
t &\mapsto (n-1)\operatorname{ct}_k(t).
\end{align*}
Since $s_k$ satisfies
\begin{align*}
s_k''(t)+k s_k(t)=0,
\qquad
s_k(0)=0,
\qquad
s_k'(0)=1,
\end{align*}
and $\operatorname{ct}_k=s_k'/s_k$, differentiation gives
\begin{align*}
\operatorname{ct}_k'(t)+\operatorname{ct}_k(t)^2+k=0.
\end{align*}
Multiplying by $n-1$ gives
\begin{align*}
m_k'(t)+\frac{m_k(t)^2}{n-1}+(n-1)k=0.
\end{align*}
We use the following local meaning of big-$O$ notation: for functions defined near $0$, the assertion $f(t)=O(t^q)$ as $t\downarrow 0$ means that there exist constants $C>0$ and $\delta>0$ such that $|f(t)|\leq C t^q$ for every $t\in(0,\delta)$. We justify the singular expansion of $S(t)$ from the Jacobi tensor rather than using it as an unexplained normal-coordinate fact. For each $V\in E_0$, the Jacobi field $J_V$ satisfies $J_V(0)=0$, $\nabla_{\gamma'}J_V(0)=V$, and the Jacobi field equation recorded on the [Jacobi Field](/page/Jacobi%20Field) page. Let
\begin{align*}
I_{E_0}:E_0 &\to E_0 \\
V &\mapsto V
\end{align*}
denote the identity operator on $E_0$. Since the curvature endomorphism $R_\gamma(t)$ is smooth along $\gamma$, Taylor expansion of this linear ordinary differential equation gives, as endomorphisms of $E_0$,
\begin{align*}
J(t)=tI_{E_0}+O(t^3),
\qquad
J'(t)=I_{E_0}+O(t^2)
\end{align*}
as $t\downarrow 0$. Inverting the first expansion gives
\begin{align*}
J(t)^{-1}=\frac{1}{t}I_{E_0}+O(t).
\end{align*}
Using the logarithmic derivative formula $S(t)=J'(t)J(t)^{-1}$, we obtain
\begin{align*}
S(t)=\frac{1}{t}I_{E_0}+O(t)
\end{align*}
as $t \downarrow 0$. Taking the trace over the $(n-1)$-dimensional space $E_0$ gives
\begin{align*}
m(t)=\frac{n-1}{t}+O(t).
\end{align*}
Also $s_k(t)=t+O(t^3)$ and $s_k'(t)=1+O(t^2)$, so
\begin{align*}
m_k(t)=\frac{n-1}{t}+O(t).
\end{align*}
Thus, for
\begin{align*}
w:(0,r_0] &\to \mathbb{R} \\
t &\mapsto m(t)-m_k(t),
\end{align*}
we have
\begin{align*}
w(t)=O(t)
\end{align*}
as $t \downarrow 0$.
[/step]
[step:Compare the scalar Riccati inequalities from the common singular initial value]
Subtracting the equality for $m_k$ from the inequality for $m$ gives
\begin{align*}
w'(t)+a(t)w(t)\leq 0,
\end{align*}
where
\begin{align*}
a:(0,r_0] &\to \mathbb{R} \\
t &\mapsto \frac{m(t)+m_k(t)}{n-1}.
\end{align*}
Let $\mathcal{L}^1$ denote one-dimensional [Lebesgue measure](/page/Lebesgue%20Measure) on $\mathbb{R}$. For $0<\varepsilon<t\leq r_0$, define the integrating factor
\begin{align*}
A_{\varepsilon,t}:=\int_\varepsilon^t a(\tau)\,d\mathcal{L}^1(\tau).
\end{align*}
Multiplying the differential inequality by
\begin{align*}
\exp\left(\int_\varepsilon^t a(\tau)\,d\mathcal{L}^1(\tau)\right)
\end{align*}
and integrating from $\varepsilon$ to $t$ yields
\begin{align*}
w(t)\leq w(\varepsilon)\exp(-A_{\varepsilon,t}).
\end{align*}
The asymptotics of $m$ and $m_k$ give
\begin{align*}
a(\tau)=\frac{2}{\tau}+O(\tau)
\end{align*}
as $\tau\downarrow 0$. Hence, for each fixed $t\in(0,r_0]$,
\begin{align*}
\exp(-A_{\varepsilon,t}) = O(\varepsilon^2)
\end{align*}
as $\varepsilon\downarrow 0$. Since $w(\varepsilon)=O(\varepsilon)$, we obtain
\begin{align*}
w(\varepsilon)\exp(-A_{\varepsilon,t})\to 0.
\end{align*}
Letting $\varepsilon\downarrow 0$ gives
\begin{align*}
w(t)\leq 0.
\end{align*}
Therefore
\begin{align*}
m(t)\leq m_k(t)
\end{align*}
for every $t\in(0,r_0]$.
[guided]
The scalar comparison step turns the Riccati inequality into a first-order linear inequality for the difference between the actual radial Laplacian and the model radial Laplacian. Define
\begin{align*}
w:(0,r_0] &\to \mathbb{R} \\
t &\mapsto m(t)-m_k(t).
\end{align*}
The functions $m$ and $m_k$ satisfy
\begin{align*}
m'(t)+\frac{m(t)^2}{n-1}+(n-1)k\leq 0
\end{align*}
and
\begin{align*}
m_k'(t)+\frac{m_k(t)^2}{n-1}+(n-1)k=0.
\end{align*}
Subtracting the second identity from the [first inequality](/theorems/2897) gives
\begin{align*}
w'(t)+\frac{m(t)^2-m_k(t)^2}{n-1}\leq 0.
\end{align*}
Factor the difference of squares:
\begin{align*}
m(t)^2-m_k(t)^2=(m(t)+m_k(t))w(t).
\end{align*}
Thus, with
\begin{align*}
a:(0,r_0] &\to \mathbb{R} \\
t &\mapsto \frac{m(t)+m_k(t)}{n-1},
\end{align*}
we obtain the linear differential inequality
\begin{align*}
w'(t)+a(t)w(t)\leq 0.
\end{align*}
Fix $0<\varepsilon<t\leq r_0$. Since $a$ is continuous on the compact interval $[\varepsilon,t]$, the integrating factor is well-defined. Let
\begin{align*}
A_{\varepsilon,s}:=\int_\varepsilon^s a(\tau)\,d\mathcal{L}^1(\tau)
\end{align*}
for $s\in[\varepsilon,t]$. Multiplying by $e^{A_{\varepsilon,s}}$ gives
\begin{align*}
\frac{d}{ds}\left(e^{A_{\varepsilon,s}}w(s)\right)
=e^{A_{\varepsilon,s}}(w'(s)+a(s)w(s))\leq 0.
\end{align*}
Integrating this inequality over $[\varepsilon,t]$ with respect to $\mathcal{L}^1$ yields
\begin{align*}
e^{A_{\varepsilon,t}}w(t)\leq w(\varepsilon),
\end{align*}
and therefore
\begin{align*}
w(t)\leq w(\varepsilon)e^{-A_{\varepsilon,t}}.
\end{align*}
Now we use the common singular initial behaviour at $0$, and we restate its derivation here so that the guided comparison is self-contained. Let
\begin{align*}
I_{E_0}:E_0 &\to E_0 \\
V &\mapsto V
\end{align*}
denote the identity operator on $E_0$. The Jacobi tensor expansion along $\gamma$ gives
\begin{align*}
J(s)=sI_{E_0}+O(s^3),
\qquad
J'(s)=I_{E_0}+O(s^2),
\end{align*}
so $S(s)=J'(s)J(s)^{-1}=s^{-1}I_{E_0}+O(s)$ and hence
\begin{align*}
m(s)=\frac{n-1}{s}+O(s)
\end{align*}
as $s\downarrow 0$. The model function satisfies $s_k(s)=s+O(s^3)$ and $s_k'(s)=1+O(s^2)$, so
\begin{align*}
m_k(s)=\frac{n-1}{s}+O(s).
\end{align*}
Therefore
\begin{align*}
w(\varepsilon)=O(\varepsilon)
\end{align*}
as $\varepsilon\downarrow 0$, and, using the definition $a=(m+m_k)/(n-1)$,
\begin{align*}
a(\tau)=\frac{2}{\tau}+O(\tau)
\end{align*}
as $\tau\downarrow 0$. Hence, for fixed $t\in(0,r_0]$, there are constants $C>0$ and $\delta\in(0,t)$ such that for $0<\varepsilon<\delta$,
\begin{align*}
e^{-A_{\varepsilon,t}}\leq C\varepsilon^2.
\end{align*}
Combining the two estimates gives
\begin{align*}
w(\varepsilon)e^{-A_{\varepsilon,t}}\to 0
\end{align*}
as $\varepsilon\downarrow 0$. Passing to the limit in
\begin{align*}
w(t)\leq w(\varepsilon)e^{-A_{\varepsilon,t}}
\end{align*}
gives
\begin{align*}
w(t)\leq 0.
\end{align*}
Since $w(t)=m(t)-m_k(t)$, this proves
\begin{align*}
m(t)\leq m_k(t)
\end{align*}
for every $t\in(0,r_0]$.
[/guided]
[/step]
[step:Evaluate the comparison at the chosen point]
Taking $t=r_0$ in the scalar comparison gives
\begin{align*}
\Delta r_p(x)
=
m(r_0)
\leq
m_k(r_0)
=
(n-1)\operatorname{ct}_k(r_0).
\end{align*}
Since $r_0=r_p(x)$, this is
\begin{align*}
\Delta r_p(x)\leq (n-1)\operatorname{ct}_k(r_p(x)).
\end{align*}
The point $x\in\Omega_p$ with $r_p(x)\in I_k$ was arbitrary, so the asserted inequality holds throughout the stated smooth region.
[/step]