[proofplan]
We approximate the Busemann function by the smooth-off-the-cut-locus functions $u_t(x)=t-d_g(x,\rho(t))$. Laplacian comparison under $\operatorname{Ric}_g \geq 0$ gives a distributional lower bound for $\Delta_g u_t$ whose error is bounded by
\begin{align*}
\frac{n-1}{d_g(x,\rho(t))}.
\end{align*}
On each compact set, the points $\rho(t)$ escape to infinity, so this error tends uniformly to $0$. Since $u_t \to b_\rho$ locally uniformly, we may pass the distributional inequality to the limit and obtain $\Delta_g b_\rho \geq 0$.
[/proofplan]
[step:Define the distance approximants and record their convergence]
For each $t>0$, define the distance function from $\rho(t)$ by
\begin{align*}
r_t: M &\to [0,\infty) \\
x &\mapsto d_g(x,\rho(t)),
\end{align*}
and define
\begin{align*}
u_t: M &\to \mathbb{R} \\
x &\mapsto t-r_t(x).
\end{align*}
By the definition of the Busemann function associated to the ray $\rho$, the functions $u_t$ converge pointwise to $b_\rho$ as $t\to\infty$.
Moreover, the convergence is locally uniform. Indeed, for $s,t \geq 0$ and $x \in M$, the triangle inequality gives
\begin{align*}
|r_t(x)-r_s(x)| \leq d_g(\rho(t),\rho(s)) = |t-s|.
\end{align*}
Thus the functions $u_t$ are $1$-Lipschitz in the spatial variable uniformly in $t$, since
\begin{align*}
|u_t(x)-u_t(y)| = |r_t(y)-r_t(x)| \leq d_g(x,y).
\end{align*}
The standard monotonicity argument for rays gives that $u_t(x)$ is nondecreasing in $t$: if $s \leq t$, then
\begin{align*}
r_t(x) \leq r_s(x)+d_g(\rho(s),\rho(t)) = r_s(x)+t-s,
\end{align*}
so
\begin{align*}
u_s(x)=s-r_s(x)\leq t-r_t(x)=u_t(x).
\end{align*}
We justify the local [uniform convergence](/page/Uniform%20Convergence) directly. Let $K_0\subset M$ be compact and fix $\varepsilon>0$. Since the functions $u_t$ are $1$-Lipschitz and the pointwise limit $b_\rho$ is also $1$-Lipschitz as a pointwise limit of $1$-Lipschitz functions, choose finitely many points $x_1,\dots,x_N\in K_0$ such that every $x\in K_0$ satisfies $d_g(x,x_j)<\varepsilon/3$ for some $j$. Since $u_t(x_j)\uparrow b_\rho(x_j)$ for each $j$, choose $T>0$ so that $0\leq b_\rho(x_j)-u_t(x_j)<\varepsilon/3$ for all $j$ and all $t\geq T$. Then, for $x\in K_0$ and a corresponding $x_j$,
\begin{align*}
0\leq b_\rho(x)-u_t(x)
&\leq |b_\rho(x)-b_\rho(x_j)|+|b_\rho(x_j)-u_t(x_j)|+|u_t(x_j)-u_t(x)| \\
&< \varepsilon.
\end{align*}
Thus $u_t\to b_\rho$ uniformly on $K_0$, and hence locally uniformly on $M$.
[guided]
The approximating functions are the finite-time Busemann functions. For each $t>0$, we first declare the maps
\begin{align*}
r_t: M &\to [0,\infty) \\
x &\mapsto d_g(x,\rho(t)),
\end{align*}
and
\begin{align*}
u_t: M &\to \mathbb{R} \\
x &\mapsto t-r_t(x).
\end{align*}
The Busemann function is defined as the limit of these maps:
\begin{align*}
b_\rho(x)=\lim_{t\to\infty}u_t(x).
\end{align*}
We need more than pointwise convergence because the final argument passes to the limit inside integrals over compact supports. The key compactness input is that all $u_t$ have the same Lipschitz constant. For any $x,y \in M$, the triangle inequality for the Riemannian distance gives
\begin{align*}
|r_t(x)-r_t(y)| \leq d_g(x,y).
\end{align*}
Since $u_t=t-r_t$, the constant $t$ cancels in differences, and hence
\begin{align*}
|u_t(x)-u_t(y)| \leq d_g(x,y).
\end{align*}
We also use the fact that the parameter family is monotone. If $s \leq t$, then the ray condition gives
\begin{align*}
d_g(\rho(s),\rho(t))=t-s.
\end{align*}
Applying the triangle inequality through $\rho(s)$,
\begin{align*}
r_t(x)=d_g(x,\rho(t)) \leq d_g(x,\rho(s))+d_g(\rho(s),\rho(t))
= r_s(x)+t-s.
\end{align*}
Rearranging gives
\begin{align*}
s-r_s(x)\leq t-r_t(x),
\end{align*}
which is exactly $u_s(x)\leq u_t(x)$.
Thus $(u_t)$ is a monotone family of functions with a uniform Lipschitz bound. We now prove the local uniform convergence rather than merely naming a compactness principle. Let $K_0\subset M$ be compact and let $\varepsilon>0$. Because $K_0$ is compact, choose finitely many points $x_1,\dots,x_N\in K_0$ such that every $x\in K_0$ satisfies $d_g(x,x_j)<\varepsilon/3$ for some $j$. The limit $b_\rho$ is also $1$-Lipschitz: for any $x,y\in M$,
\begin{align*}
|b_\rho(x)-b_\rho(y)|
=\lim_{t\to\infty}|u_t(x)-u_t(y)|
\leq d_g(x,y).
\end{align*}
For each fixed net point $x_j$, monotone pointwise convergence gives $u_t(x_j)\uparrow b_\rho(x_j)$. Since there are only finitely many net points, choose $T>0$ such that
\begin{align*}
0\leq b_\rho(x_j)-u_t(x_j)<\varepsilon/3
\end{align*}
for every $j\in\{1,\dots,N\}$ and every $t\geq T$. Now take any $x\in K_0$ and choose $x_j$ with $d_g(x,x_j)<\varepsilon/3$. Using the Lipschitz bounds for both $u_t$ and $b_\rho$,
\begin{align*}
0\leq b_\rho(x)-u_t(x)
&\leq |b_\rho(x)-b_\rho(x_j)|+|b_\rho(x_j)-u_t(x_j)|+|u_t(x_j)-u_t(x)| \\
&< \varepsilon.
\end{align*}
This proves uniform convergence on $K_0$. Since $K_0$ was an arbitrary compact subset of $M$, $u_t\to b_\rho$ locally uniformly.
[/guided]
[/step]
[step:Apply Laplacian comparison to the finite-time approximants]
Fix $t>0$. Let $\operatorname{Cut}(\rho(t))\subset M$ denote the cut locus of the point $\rho(t)$, namely the set of points where the minimizing geodesic from $\rho(t)$ ceases to be unique or ceases to minimize beyond the point. The distance function $r_t$ is smooth on $M\setminus(\{\rho(t)\}\cup \operatorname{Cut}(\rho(t)))$. Since $(M,g)$ is complete and $\operatorname{Ric}_g\geq 0$, the [Laplacian comparison theorem for distance functions](/page/Laplacian%20Comparison%20Theorem) applies to the distance from the point $\rho(t)$ and gives, in the barrier sense and hence in the distributional sense,
\begin{align*}
\Delta_g r_t \leq \frac{n-1}{r_t}
\end{align*}
on $M\setminus\{\rho(t)\}$.
Because $u_t=t-r_t$ and $t$ is constant in the spatial variable $x$, this gives
\begin{align*}
\Delta_g u_t = -\Delta_g r_t \geq -\frac{n-1}{r_t}
\end{align*}
in the distributional sense on $M\setminus\{\rho(t)\}$.
Equivalently, for every nonnegative $\varphi\in C_c^\infty(M\setminus\{\rho(t)\})$,
\begin{align*}
\int_M u_t(x)\,\Delta_g\varphi(x)\,d\operatorname{vol}_g(x)
\geq
-\int_M \frac{n-1}{r_t(x)}\varphi(x)\,d\operatorname{vol}_g(x).
\end{align*}
[guided]
Fix $t>0$. The point $\rho(t)$ is the pole of the distance function, and we declare
\begin{align*}
r_t: M &\to [0,\infty) \\
x &\mapsto d_g(x,\rho(t)),
\end{align*}
so $r_t$ is the Riemannian distance from $\rho(t)$. Let $\operatorname{Cut}(\rho(t))\subset M$ denote the [cut locus](/page/Cut%20Locus) of $\rho(t)$. On the [open set](/page/Open%20Set) $M\setminus(\{\rho(t)\}\cup \operatorname{Cut}(\rho(t)))$, the distance function is smooth.
We apply the [Laplacian comparison theorem for distance functions](/page/Laplacian%20Comparison%20Theorem). Its relevant hypotheses are completeness of $(M,g)$ and the lower Ricci bound $\operatorname{Ric}_g\geq 0$, both of which are hypotheses of the theorem. Applied to the distance from $\rho(t)$, it yields the barrier inequality
\begin{align*}
\Delta_g r_t \leq \frac{n-1}{r_t}
\end{align*}
on $M\setminus\{\rho(t)\}$. The barrier formulation implies the same inequality in the distributional sense, so it may be tested against nonnegative smooth compactly supported functions away from the pole.
Now define
\begin{align*}
u_t: M &\to \mathbb{R} \\
x &\mapsto t-r_t(x).
\end{align*}
The scalar $t$ is constant as a function of the spatial variable $x$, so $\Delta_g t=0$ distributionally. Therefore
\begin{align*}
\Delta_g u_t=\Delta_g(t-r_t)=-\Delta_g r_t\geq -\frac{n-1}{r_t}
\end{align*}
in the distributional sense on $M\setminus\{\rho(t)\}$. Equivalently, for every nonnegative $\varphi\in C_c^\infty(M\setminus\{\rho(t)\})$,
\begin{align*}
\int_M u_t(x)\,\Delta_g\varphi(x)\,d\operatorname{vol}_g(x)
\geq
-\int_M \frac{n-1}{r_t(x)}\varphi(x)\,d\operatorname{vol}_g(x).
\end{align*}
This is the finite-time subharmonicity estimate; the negative right-hand side is the error that must vanish as the pole $\rho(t)$ goes to infinity.
[/guided]
[/step]
[step:Move the pole away from the compact support]
Let $\varphi\in C_c^\infty(M)$ be nonnegative, and define its compact support by
\begin{align*}
K:=\operatorname{supp}\varphi \subset M.
\end{align*}
Choose a point $p_0\in M$ and define
\begin{align*}
A_K:=\sup_{x\in K} d_g(x,p_0).
\end{align*}
This number is finite because $K$ is compact and $x\mapsto d_g(x,p_0)$ is continuous.
For every $x\in K$, the [reverse triangle inequality](/theorems/2300) gives
\begin{align*}
r_t(x)=d_g(x,\rho(t))
\geq d_g(p_0,\rho(t))-d_g(x,p_0).
\end{align*}
Since $\rho$ is a unit-speed ray,
\begin{align*}
d_g(p_0,\rho(t))\geq d_g(\rho(0),\rho(t))-d_g(p_0,\rho(0))
= t-d_g(p_0,\rho(0)).
\end{align*}
Therefore, for every $x\in K$,
\begin{align*}
r_t(x)\geq t-d_g(p_0,\rho(0))-A_K.
\end{align*}
In particular, for all sufficiently large $t$, the point $\rho(t)$ is not in $K$, so the distributional inequality from the previous step applies to $\varphi$.
For such $t$, set
\begin{align*}
E_K(t):=t-d_g(p_0,\rho(0))-A_K.
\end{align*}
Then $E_K(t)\to\infty$ as $t\to\infty$, and on $K$ we have $r_t(x)\geq E_K(t)>0$. Hence
\begin{align*}
0\leq
\int_M \frac{n-1}{r_t(x)}\varphi(x)\,d\operatorname{vol}_g(x)
\leq
\frac{n-1}{E_K(t)}
\int_M \varphi(x)\,d\operatorname{vol}_g(x).
\end{align*}
The right-hand side tends to $0$ as $t\to\infty$.
[guided]
Let $\varphi\in C_c^\infty(M)$ be nonnegative, and define the compact set
\begin{align*}
K:=\operatorname{supp}\varphi\subset M.
\end{align*}
The distributional inequality from the previous step can only be tested against functions whose support avoids the pole $\rho(t)$. Thus we must prove that, for large $t$, the point $\rho(t)$ lies outside $K$, and we also need a quantitative lower bound for $r_t$ on $K$.
Choose a point $p_0\in M$ and define
\begin{align*}
A_K:=\sup_{x\in K}d_g(x,p_0).
\end{align*}
The number $A_K$ is finite because $K$ is compact and the map $x\mapsto d_g(x,p_0)$ is continuous. For any $x\in K$, the reverse triangle inequality gives
\begin{align*}
r_t(x)=d_g(x,\rho(t))
\geq d_g(p_0,\rho(t))-d_g(x,p_0).
\end{align*}
We next estimate the first term from below. Since $\rho$ is a unit-speed ray,
\begin{align*}
d_g(\rho(0),\rho(t))=t.
\end{align*}
Applying the reverse triangle inequality again,
\begin{align*}
d_g(p_0,\rho(t))
\geq d_g(\rho(0),\rho(t))-d_g(p_0,\rho(0))
= t-d_g(p_0,\rho(0)).
\end{align*}
Combining the two estimates gives, for every $x\in K$,
\begin{align*}
r_t(x)\geq t-d_g(p_0,\rho(0))-A_K.
\end{align*}
Define the positive lower-bound function for large $t$ by
\begin{align*}
E_K(t):=t-d_g(p_0,\rho(0))-A_K.
\end{align*}
Then $E_K(t)\to\infty$ as $t\to\infty$, and for all sufficiently large $t$ we have $E_K(t)>0$. Hence $\rho(t)\notin K$ and $r_t(x)\geq E_K(t)$ on $K$.
For those large values of $t$, the nonnegative [test function](/page/Test%20Function) $\varphi$ is admissible in the distributional inequality from the previous step. Since $\varphi$ is supported in $K$ and $r_t(x)\geq E_K(t)$ on $K$,
\begin{align*}
0\leq
\int_M \frac{n-1}{r_t(x)}\varphi(x)\,d\operatorname{vol}_g(x)
\leq
\frac{n-1}{E_K(t)}
\int_M \varphi(x)\,d\operatorname{vol}_g(x).
\end{align*}
The integral of $\varphi$ is finite because $\varphi$ is smooth with compact support. Since $E_K(t)\to\infty$, the right-hand side tends to $0$. This proves that the Laplacian comparison error vanishes on every fixed compact support.
[/guided]
[/step]
[step:Pass the distributional inequality to the Busemann limit]
For all sufficiently large $t$, combining the previous two steps gives
\begin{align*}
\int_M u_t(x)\,\Delta_g\varphi(x)\,d\operatorname{vol}_g(x)
\geq
-\int_M \frac{n-1}{r_t(x)}\varphi(x)\,d\operatorname{vol}_g(x).
\end{align*}
The right-hand side tends to $0$ as $t\to\infty$.
Since $\operatorname{supp}\varphi=K$ is compact, $\Delta_g\varphi$ is smooth with compact support contained in $K$. We justify the limit by a uniform estimate on $K$:
\begin{align*}
\left|
\int_M (u_t(x)-b_\rho(x))\,\Delta_g\varphi(x)\,d\operatorname{vol}_g(x)
\right|
&\leq
\sup_{x\in K}|u_t(x)-b_\rho(x)|
\int_K |\Delta_g\varphi(x)|\,d\operatorname{vol}_g(x).
\end{align*}
The integral on the right is finite because $\Delta_g\varphi$ is smooth and compactly supported, and the supremum tends to $0$ by locally uniform convergence. Therefore
\begin{align*}
\lim_{t\to\infty}
\int_M u_t(x)\,\Delta_g\varphi(x)\,d\operatorname{vol}_g(x)
=
\int_M b_\rho(x)\,\Delta_g\varphi(x)\,d\operatorname{vol}_g(x).
\end{align*}
Taking the limit inferior in the preceding distributional inequality yields
\begin{align*}
\int_M b_\rho(x)\,\Delta_g\varphi(x)\,d\operatorname{vol}_g(x)\geq 0.
\end{align*}
Since $\varphi\in C_c^\infty(M)$ was arbitrary and nonnegative, this is precisely $\Delta_g b_\rho\geq 0$ in the [distributional sense](/page/Distribution). Therefore $b_\rho$ is [subharmonic](/page/Subharmonic%20Function).
[guided]
Let $\varphi\in C_c^\infty(M)$ be nonnegative, and let
\begin{align*}
K:=\operatorname{supp}\varphi\subset M.
\end{align*}
For all sufficiently large $t$, the pole $\rho(t)$ is outside $K$, so the finite-time distributional inequality applies to $\varphi$:
\begin{align*}
\int_M u_t(x)\,\Delta_g\varphi(x)\,d\operatorname{vol}_g(x)
\geq
-\int_M \frac{n-1}{r_t(x)}\varphi(x)\,d\operatorname{vol}_g(x).
\end{align*}
The previous step proved that the right-hand side tends to $0$ as $t\to\infty$.
It remains to justify the limit on the left. Since $\Delta_g\varphi$ is smooth and supported in the compact set $K$, we estimate the difference of the two distributional pairings directly:
\begin{align*}
\left|
\int_M (u_t(x)-b_\rho(x))\,\Delta_g\varphi(x)\,d\operatorname{vol}_g(x)
\right|
&\leq
\sup_{x\in K}|u_t(x)-b_\rho(x)|
\int_K |\Delta_g\varphi(x)|\,d\operatorname{vol}_g(x).
\end{align*}
The integral $\int_K |\Delta_g\varphi(x)|\,d\operatorname{vol}_g(x)$ is finite because $\Delta_g\varphi$ is continuous and $K$ is compact. The supremum tends to $0$ because $u_t\to b_\rho$ locally uniformly. Hence
\begin{align*}
\lim_{t\to\infty}
\int_M u_t(x)\,\Delta_g\varphi(x)\,d\operatorname{vol}_g(x)
=
\int_M b_\rho(x)\,\Delta_g\varphi(x)\,d\operatorname{vol}_g(x).
\end{align*}
Taking the limit inferior in the finite-time inequality gives
\begin{align*}
\int_M b_\rho(x)\,\Delta_g\varphi(x)\,d\operatorname{vol}_g(x)\geq 0.
\end{align*}
Because this holds for every nonnegative $\varphi\in C_c^\infty(M)$, it is exactly the definition of $\Delta_g b_\rho\geq 0$ in the [distributional sense](/page/Distribution). Therefore the [Busemann function](/page/Busemann%20Function) $b_\rho$ is [subharmonic](/page/Subharmonic%20Function).
[/guided]
[/step]