[proofplan]
Away from the cut locus of $p$, the distance function is smooth and the result is exactly the smooth [Laplacian comparison theorem](/theorems/5360). At a cut point $q$, we use Calabi's trick: move the base point a small distance $\delta$ along a minimizing geodesic from $p$ to $q$, and use the shifted function $\delta+d_g(p_\delta,\cdot)$ as a smooth upper support. Smooth Laplacian comparison applies to the shifted distance at $q$, and continuity of the model function $\operatorname{ct}_k$ lets us choose $\delta$ so that the right-hand side is evaluated at $r_p(q)$ up to the prescribed barrier tolerance.
[/proofplan]
[step:Apply smooth Laplacian comparison at regular distance points]
Let $q \in M \setminus \{p\}$ satisfy $r_p(q) \in I_k$. If $q$ is not in the cut locus of $p$, then there is an open neighbourhood $U \subset M \setminus \{p\}$ of $q$ on which
\begin{align*}
r_p|_U: U &\to \mathbb{R}
\end{align*}
is smooth. Since $(M,g)$ is complete and $\operatorname{Ric}_g \geq (n-1)kg$, the smooth Laplacian comparison theorem gives
\begin{align*}
\Delta r_p(q) \leq (n-1)\operatorname{ct}_k(r_p(q)).
\end{align*}
Here the cited comparison is the standard smooth distance comparison away from the cut locus, applied on the model interval $I_k$ where $\operatorname{sn}_k>0$ and $\operatorname{ct}_k$ is finite. Thus $r_p|_U$ itself is an admissible upper barrier, with zero tolerance error.
(given result not yet resolved in the wiki: Smooth Laplacian Comparison Theorem)
[/step]
[step:Construct a Calabi upper barrier at a cut point]
Assume now that $q$ lies in the cut locus of $p$. Set
\begin{align*}
L := r_p(q) = d_g(p,q).
\end{align*}
Since $(M,g)$ is complete, the Hopf-Rinow theorem gives a unit-speed minimizing geodesic
\begin{align*}
\gamma: [0,L] &\to M
\end{align*}
such that $\gamma(0)=p$ and $\gamma(L)=q$. Let $\delta$ be a number with $0<\delta<L$ and $L-\delta \in I_k$, and define
\begin{align*}
p_\delta &:= \gamma(\delta).
\end{align*}
Define the shifted distance support function
\begin{align*}
r_\delta: U_\delta &\to \mathbb{R}\\
x &\mapsto \delta+d_g(p_\delta,x),
\end{align*}
where $U_\delta$ is a neighbourhood of $q$ on which $d_g(p_\delta,\cdot)$ is smooth. Such a neighbourhood exists by Calabi's regularity lemma for the shifted base point along a minimizing geodesic.
For every $x \in M$, the triangle inequality gives
\begin{align*}
r_p(x)=d_g(p,x) \leq d_g(p,p_\delta)+d_g(p_\delta,x)=\delta+d_g(p_\delta,x)=r_\delta(x).
\end{align*}
At $x=q$, the segment $\gamma|_{[\delta,L]}$ is minimizing from $p_\delta$ to $q$, so
\begin{align*}
r_\delta(q)=\delta+d_g(p_\delta,q)=\delta+(L-\delta)=L=r_p(q).
\end{align*}
Thus $r_\delta$ touches $r_p$ from above at $q$.
(given result not yet resolved in the wiki: Hopf-Rinow Theorem)
(given result not yet resolved in the wiki: Calabi's Regularity Lemma for Distance Barriers)
[guided]
Because $q$ may lie in the cut locus of $p$, the original distance function $r_p$ may fail to be smooth at $q$, so it cannot necessarily be inserted directly into the Laplacian. The barrier method replaces $r_p$ near $q$ by a smooth function that sits above it and agrees with it at $q$.
Let
\begin{align*}
L := r_p(q)=d_g(p,q).
\end{align*}
Completeness of $(M,g)$ allows us to use the Hopf-Rinow theorem: there exists a unit-speed minimizing geodesic
\begin{align*}
\gamma: [0,L] &\to M
\end{align*}
with $\gamma(0)=p$ and $\gamma(L)=q$. Choose $0<\delta<L$ with $L-\delta \in I_k$, and define
\begin{align*}
p_\delta := \gamma(\delta).
\end{align*}
The point $p_\delta$ is obtained by moving the base point slightly forward along the chosen minimizing geodesic. Calabi's regularity lemma says that, for such a shifted base point, the distance function $d_g(p_\delta,\cdot)$ is smooth in some neighbourhood $U_\delta$ of $q$. Therefore the function
\begin{align*}
r_\delta: U_\delta &\to \mathbb{R}\\
x &\mapsto \delta+d_g(p_\delta,x)
\end{align*}
is a $C^2$ function near $q$.
We now verify that this smooth function is an upper support for $r_p$. For every $x \in M$, the triangle inequality gives
\begin{align*}
d_g(p,x) \leq d_g(p,p_\delta)+d_g(p_\delta,x).
\end{align*}
Since $\gamma$ is unit-speed and $p_\delta=\gamma(\delta)$, we have
\begin{align*}
d_g(p,p_\delta)=\delta.
\end{align*}
Hence
\begin{align*}
r_p(x)=d_g(p,x) \leq \delta+d_g(p_\delta,x)=r_\delta(x).
\end{align*}
At the point $q$, the subsegment $\gamma|_{[\delta,L]}$ is still minimizing, so
\begin{align*}
d_g(p_\delta,q)=L-\delta.
\end{align*}
Thus
\begin{align*}
r_\delta(q)=\delta+d_g(p_\delta,q)=\delta+(L-\delta)=L=r_p(q).
\end{align*}
So $r_\delta$ lies above $r_p$ and touches it exactly at $q$, which is precisely the geometric requirement for an upper barrier.
[/guided]
[/step]
[step:Estimate the Laplacian of the shifted barrier]
The function $r_\delta$ differs from $d_g(p_\delta,\cdot)$ by the constant $\delta$, so
\begin{align*}
\Delta r_\delta(q)=\Delta d_g(p_\delta,\cdot)(q).
\end{align*}
The radial segment from $p_\delta$ to $q$ has length
\begin{align*}
d_g(p_\delta,q)=L-\delta.
\end{align*}
Since the ambient Ricci lower bound $\operatorname{Ric}_g \geq (n-1)kg$ is unchanged by shifting the base point, and since $L-\delta \in I_k$, the smooth Laplacian comparison theorem applied to the smooth distance function $d_g(p_\delta,\cdot)$ at $q$ gives
\begin{align*}
\Delta r_\delta(q)
=
\Delta d_g(p_\delta,\cdot)(q)
\leq
(n-1)\operatorname{ct}_k(L-\delta).
\end{align*}
(given result not yet resolved in the wiki: Smooth Laplacian Comparison Theorem)
[/step]
[step:Choose the shifted base point to match the desired tolerance]
The number $L=r_p(q)$ belongs to $I_k$. Since $I_k$ is open and $\operatorname{ct}_k: I_k \to \mathbb{R}$ is smooth, it is continuous at $L$. Given $\varepsilon>0$, choose $0<\delta<L$ small enough that $L-\delta \in I_k$ and
\begin{align*}
(n-1)\operatorname{ct}_k(L-\delta)
\leq
(n-1)\operatorname{ct}_k(L)+\varepsilon.
\end{align*}
For this choice of $\delta$, the function $r_\delta \in C^2(U_\delta)$ satisfies
\begin{align*}
r_\delta(q)&=r_p(q),\\
r_\delta(x)&\geq r_p(x) \qquad \text{for every } x \in U_\delta,\\
\Delta r_\delta(q)
&\leq
(n-1)\operatorname{ct}_k(L-\delta)
\leq
(n-1)\operatorname{ct}_k(L)+\varepsilon\\
&=
(n-1)\operatorname{ct}_k(r_p(q))+\varepsilon.
\end{align*}
Thus $r_\delta$ is an admissible upper barrier at $q$ with tolerance $\varepsilon$.
[/step]
[step:Conclude the barrier comparison on the model interval]
For every $q \in M\setminus\{p\}$ with $r_p(q)\in I_k$, we have constructed, for each $\varepsilon>0$, a $C^2$ upper barrier whose Laplacian at $q$ is bounded above by
\begin{align*}
(n-1)\operatorname{ct}_k(r_p(q))+\varepsilon.
\end{align*}
Therefore
\begin{align*}
\Delta r_p \leq (n-1)\operatorname{ct}_k(r_p)
\end{align*}
on $M\setminus\{p\}$ in the upper barrier sense, at every point where $r_p$ lies in the interval $I_k$ on which the model expression is finite. This proves the theorem.
[/step]