[proofplan]
We compare the two Busemann functions associated to the two ends of the line and use the standard Busemann-function part of the Cheeger-Gromoll argument: under $\operatorname{Ric}_g \geq 0$, each is weakly subharmonic and their sum is forced to vanish. Hence $b_+$ is both weakly subharmonic and weakly superharmonic, so it is weakly harmonic; elliptic regularity upgrades it to a smooth harmonic function. The global $1$-Lipschitz bound gives $|\nabla b|_g \leq 1$, while the identity $b(\gamma(s))=s$ gives equality along the line. The Bochner formula and the strong maximum principle then force $|\nabla b|_g^2$ to be constant, and the Bochner identity finally gives $\operatorname{Hess} b=0$.
[/proofplan]
[step:Introduce the opposite Busemann function and use the line comparison identity]
Define the backward Busemann function
\begin{align*}
b_-: M &\to \mathbb{R},\\
x &\mapsto \lim_{t \to \infty}\bigl(t - d_g(x,\gamma(-t))\bigr).
\end{align*}
The standard Busemann comparison theorem for a line in a complete manifold with nonnegative Ricci curvature states that $b_+$ and $b_-$ are locally Lipschitz, weakly subharmonic, and satisfy
\begin{align*}
b_+ + b_- = 0
\end{align*}
on $M$ (citing a result not yet in the wiki: Busemann comparison identity for a line under nonnegative Ricci curvature).
Since $b_-$ is weakly subharmonic and $b_+ = -b_-$, the function $b_+$ is weakly superharmonic. Since $b_+$ is also weakly subharmonic, $b=b_+$ is weakly harmonic:
\begin{align*}
\Delta b = 0
\end{align*}
in the distributional sense.
[/step]
[step:Upgrade the weakly harmonic Lipschitz function to a smooth harmonic function]
Because $b$ is locally Lipschitz, $b \in W_{\mathrm{loc}}^{1,\infty}(M) \subset W_{\mathrm{loc}}^{1,2}(M)$. Since $\Delta b=0$ weakly, interior elliptic regularity for the Laplace-Beltrami operator gives
\begin{align*}
b \in C^\infty(M)
\end{align*}
and the equality $\Delta b=0$ holds pointwise on $M$ (citing a result not yet in the wiki: elliptic regularity for weakly harmonic functions on smooth Riemannian manifolds).
[guided]
We have reached the point where the Busemann function is known only as a locally Lipschitz function. That regularity is enough to interpret its weak gradient and to formulate the weak equation $\Delta b=0$, but it is not enough yet to use the pointwise Bochner identity. The bridge is elliptic regularity for the Laplace-Beltrami operator.
Since $b$ is locally Lipschitz, for every relatively compact coordinate ball $U \subset M$ the coordinate representative of $b$ belongs to $W^{1,\infty}(U)$, hence also to $W^{1,2}(U)$. The weak harmonicity statement means that for every [test function](/page/Test%20Function) $\varphi \in C_c^\infty(U)$,
\begin{align*}
\int_U g(\nabla b,\nabla \varphi)_g \, d\operatorname{vol}_g = 0.
\end{align*}
In local coordinates this is a uniformly elliptic divergence-form equation with smooth coefficients, because the metric tensor $g$ is smooth and positive definite. Interior elliptic regularity therefore implies that $b$ is smooth on each such coordinate ball. Since these coordinate balls cover $M$, we obtain
\begin{align*}
b \in C^\infty(M).
\end{align*}
The same weak equation then agrees with the classical equation, so
\begin{align*}
\Delta b = 0
\end{align*}
pointwise on $M$.
[/guided]
[/step]
[step:Use the Lipschitz bound to obtain the global gradient estimate]
For each $t>0$, define
\begin{align*}
b_t: M &\to \mathbb{R},\\
x &\mapsto t - d_g(x,\gamma(t)).
\end{align*}
The [reverse triangle inequality](/theorems/2300) gives
\begin{align*}
|b_t(x)-b_t(y)| \leq d_g(x,y)
\end{align*}
for all $x,y \in M$. Passing to the limit $t \to \infty$ gives
\begin{align*}
|b(x)-b(y)| \leq d_g(x,y)
\end{align*}
for all $x,y \in M$, so $b$ is $1$-Lipschitz. Since $b$ is smooth, evaluating this Lipschitz inequality along smooth curves through any point gives
\begin{align*}
|\nabla b|_g \leq 1
\end{align*}
on $M$.
[/step]
[step:Show that the gradient has unit length along the line]
For every $s \in \mathbb{R}$ and every $t>s$,
\begin{align*}
d_g(\gamma(s),\gamma(t)) = t-s
\end{align*}
because $\gamma$ is a unit-speed line. Hence
\begin{align*}
b(\gamma(s))
&= \lim_{t \to \infty}\bigl(t-d_g(\gamma(s),\gamma(t))\bigr)\\
&= \lim_{t \to \infty}\bigl(t-(t-s)\bigr)\\
&= s.
\end{align*}
Differentiating the smooth identity $b(\gamma(s))=s$ with respect to $s$ gives
\begin{align*}
g_{\gamma(s)}\bigl(\nabla b|_{\gamma(s)},\dot{\gamma}(s)\bigr)=1.
\end{align*}
Since $|\dot{\gamma}(s)|_g=1$ and $|\nabla b|_g \leq 1$, the [Cauchy-Schwarz inequality](/theorems/432) forces
\begin{align*}
|\nabla b|_g(\gamma(s))=1
\end{align*}
for every $s \in \mathbb{R}$.
[/step]
[step:Apply Bochner and the strong maximum principle to make $|\nabla b|_g^2$ constant]
Define the smooth function
\begin{align*}
u: M &\to \mathbb{R},\\
x &\mapsto |\nabla b|_g^2(x).
\end{align*}
The Bochner identity for a smooth function on a Riemannian manifold gives
\begin{align*}
\frac{1}{2}\Delta u
=
|\operatorname{Hess} b|_g^2
+
g(\nabla b,\nabla \Delta b)_g
+
\operatorname{Ric}_g(\nabla b,\nabla b)
\end{align*}
(citing a result not yet in the wiki: Bochner formula for functions). Since $\Delta b=0$, the middle term vanishes. Since $\operatorname{Ric}_g \geq 0$, we obtain
\begin{align*}
\frac{1}{2}\Delta u
=
|\operatorname{Hess} b|_g^2
+
\operatorname{Ric}_g(\nabla b,\nabla b)
\geq 0.
\end{align*}
Thus $u$ is smooth and subharmonic. From the previous steps,
\begin{align*}
u \leq 1
\end{align*}
on $M$, and
\begin{align*}
u(\gamma(s))=1
\end{align*}
for every $s \in \mathbb{R}$. Therefore $u$ attains its global maximum at an interior point of the connected manifold $M$. By the strong maximum principle for smooth subharmonic functions (citing a result not yet in the wiki: strong maximum principle for subharmonic functions on connected Riemannian manifolds),
\begin{align*}
u \equiv 1
\end{align*}
on $M$.
[/step]
[step:Return to Bochner to force the Hessian to vanish]
Since $u=|\nabla b|_g^2 \equiv 1$, we have
\begin{align*}
\Delta u=0.
\end{align*}
Substituting this into the Bochner identity and using $\Delta b=0$ gives
\begin{align*}
0
=
|\operatorname{Hess} b|_g^2
+
\operatorname{Ric}_g(\nabla b,\nabla b).
\end{align*}
Both terms on the right are nonnegative: $|\operatorname{Hess} b|_g^2 \geq 0$ by definition of the tensor norm, and $\operatorname{Ric}_g(\nabla b,\nabla b)\geq 0$ by the Ricci curvature hypothesis. Therefore
\begin{align*}
|\operatorname{Hess} b|_g^2=0
\end{align*}
on $M$. Hence
\begin{align*}
\operatorname{Hess} b=0
\end{align*}
as a symmetric $2$-tensor on $M$. Together with the smoothness of $b$ and the identity $|\nabla b|_g^2=u\equiv 1$, this proves all asserted conclusions.
[/step]
