[proofplan] The proof is the metric part of Calabi's trick. The triangle inequality gives a global upper support function for the nonsmooth distance function $r_p$. Equality at $q$ follows because the chosen minimizing geodesic from $p$ to $q$ passes through $p_\delta$ and remains minimizing on every subinterval. The only smoothness input is the defining property of the cut locus: away from $p_\delta$ and its cut locus, the distance from $p_\delta$ is smooth in a neighbourhood of the point. [/proofplan]

[step:Use the triangle inequality to place $r_\delta$ above $r_p$]Fix $\delta \in (0,R)$. Since $p_\delta = \gamma(\delta)$ and $\gamma$ is unit speed, \begin{align*} d(p,p_\delta) = d(\gamma(0),\gamma(\delta)) = \delta. \end{align*} For an arbitrary point $x \in M$, the triangle inequality for the Riemannian distance $d: M \times M \to \mathbb{R}$ gives \begin{align*} r_p(x) = d(p,x) \le d(p,p_\delta) + d(p_\delta,x) = \delta + d(p_\delta,x) = r_\delta(x). \end{align*} Thus $r_\delta \ge r_p$ on all of $M$, hence in particular on every neighbourhood of $q$.[/step]

[guided]Fix $\delta \in (0,R)$. The point $p_\delta$ lies on the chosen geodesic from $p$ to $q$: \begin{align*} p_\delta = \gamma(\delta). \end{align*} Because $\gamma$ is parametrized by arclength and is minimizing on $[0,R]$, the distance from $p=\gamma(0)$ to $p_\delta=\gamma(\delta)$ is exactly the length of the subsegment: \begin{align*} d(p,p_\delta)=d(\gamma(0),\gamma(\delta))=\delta. \end{align*} Now take any point $x \in M$. The function $r_\delta$ is designed to compare the broken path from $p$ to $p_\delta$ and then from $p_\delta$ to $x$ with the shortest path from $p$ to $x$. The triangle inequality gives \begin{align*} d(p,x) \le d(p,p_\delta)+d(p_\delta,x). \end{align*} Substituting $d(p,p_\delta)=\delta$ yields \begin{align*} r_p(x) = d(p,x) \le \delta+d(p_\delta,x) = r_\delta(x). \end{align*} This proves the stronger global comparison $r_\delta \ge r_p$ on $M$, not merely a local comparison near $q$.[/guided]

[step:Evaluate the comparison function at $q$ along the minimizing geodesic] Since $\gamma|_{[\delta,R]}$ is a unit-speed minimizing geodesic from $p_\delta$ to $q$, we have \begin{align*} d(p_\delta,q)=R-\delta. \end{align*} Therefore \begin{align*} r_\delta(q) = \delta+d(p_\delta,q) = \delta+(R-\delta) = R = r_p(q). \end{align*} Hence $r_\delta$ touches $r_p$ from above at $q$. [/step]

[step:Use absence from the cut locus to obtain a smooth local barrier] Assume that $q$ is not in the cut locus of $p_\delta$. Since $\delta \in (0,R)$, we have $p_\delta \ne q$. By the standard local smoothness property of the Riemannian distance away from the base point and its cut locus, there exists an open neighbourhood $U_\delta \subset M$ of $q$ such that the map \begin{align*} x &\mapsto d(p_\delta,x) \end{align*} is smooth on $U_\delta$. Since $r_\delta$ is obtained from this map by adding the constant $\delta$, the restriction $r_\delta|_{U_\delta}$ is smooth. From the first step, \begin{align*} r_\delta(x) \ge r_p(x) \quad \text{for all } x \in U_\delta, \end{align*} and from the second step, \begin{align*} r_\delta(q)=r_p(q). \end{align*} Thus $r_\delta|_{U_\delta}$ is a smooth upper barrier for $r_p$ at $q$. [/step]