[proofplan]
We first use the triangle inequality along the line $\gamma$ to prove the global upper bound $b_+ + b_- \leq 0$. We then compute the two Busemann functions on the line itself, obtaining equality there. Under the hypothesis $\operatorname{Ric}_g \geq 0$, each Busemann function is distributionally subharmonic, so their sum is a Lipschitz distributionally subharmonic function with a global maximum. The strong maximum principle for Lipschitz distributionally subharmonic functions forces this sum to be constant on the connected manifold, and the value on the line identifies the constant as $0$.
[/proofplan]
[step:Use the triangle inequality to bound the sum from above]
Define the function $u: M \to \mathbb{R}$ by
\begin{align*}
u(x) := b_+(x) + b_-(x).
\end{align*}
Fix $x \in M$ and $t > 0$. Since $\gamma$ is a line,
\begin{align*}
d_g(\gamma(t),\gamma(-t)) = 2t.
\end{align*}
The triangle inequality for the Riemannian distance $d_g$ gives
\begin{align*}
2t
= d_g(\gamma(t),\gamma(-t))
\leq d_g(\gamma(t),x) + d_g(x,\gamma(-t)).
\end{align*}
Rearranging,
\begin{align*}
\bigl(t - d_g(x,\gamma(t))\bigr)
+
\bigl(t - d_g(x,\gamma(-t))\bigr)
\leq 0.
\end{align*}
Taking the limit as $t \to \infty$ in the two defining Busemann limits gives
\begin{align*}
u(x) = b_+(x) + b_-(x) \leq 0.
\end{align*}
Since $x \in M$ was arbitrary, $u \leq 0$ on $M$.
[/step]
[step:Compute the sum on the line]
Fix $s \in \mathbb{R}$. For every $t > |s|$, the line property gives
\begin{align*}
d_g(\gamma(s),\gamma(t)) &= t-s,\\
d_g(\gamma(s),\gamma(-t)) &= t+s.
\end{align*}
Therefore
\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*}
and
\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*}
Thus
\begin{align*}
u(\gamma(s)) = b_+(\gamma(s)) + b_-(\gamma(s)) = s - s = 0.
\end{align*}
In particular, $u$ attains its global maximum value $0$ at every point of $\gamma(\mathbb{R})$.
[/step]
[step:Show that the sum is distributionally subharmonic]
Each Busemann function $b_+,b_-:M \to \mathbb{R}$ is locally Lipschitz, because each function $x \mapsto t - d_g(x,\gamma(\pm t))$ is $1$-Lipschitz and the Busemann functions are pointwise limits of such functions. Since $\operatorname{Ric}_g \geq 0$, the Busemann subharmonicity theorem gives that $b_+$ and $b_-$ are distributionally subharmonic on $M$ (citing a result not yet in the wiki: Busemann functions are distributionally subharmonic under nonnegative Ricci curvature).
Distributional subharmonicity is closed under addition: if $\Delta_g b_+ \geq 0$ and $\Delta_g b_- \geq 0$ in the sense of distributions, then
\begin{align*}
\Delta_g u = \Delta_g b_+ + \Delta_g b_- \geq 0
\end{align*}
in the sense of distributions. Hence $u$ is a locally Lipschitz distributionally subharmonic function on $M$.
[guided]
The maximum principle will apply to the function
\begin{align*}
u: M \to \mathbb{R}, \qquad x \mapsto b_+(x)+b_-(x).
\end{align*}
For that reason, we must verify that $u$ has the correct regularity and subharmonicity.
First, $b_+$ and $b_-$ are locally Lipschitz. Indeed, for each fixed $t>0$, the maps
\begin{align*}
x &\mapsto t-d_g(x,\gamma(t)),\\
x &\mapsto t-d_g(x,\gamma(-t))
\end{align*}
are $1$-Lipschitz, because the distance function is $1$-Lipschitz in each variable. Passing to the pointwise Busemann limits preserves the same Lipschitz bound, so both $b_+$ and $b_-$ are $1$-Lipschitz on $M$. Therefore their sum $u=b_++b_-$ is locally Lipschitz.
Next, we use the curvature hypothesis. Since $(M,g)$ is complete and satisfies $\operatorname{Ric}_g \geq 0$, the Busemann subharmonicity theorem states that Busemann functions associated to rays are subharmonic in the distributional sense (citing a result not yet in the wiki: Busemann functions are distributionally subharmonic under nonnegative Ricci curvature). Applying this theorem to the ray $t \mapsto \gamma(t)$ gives
\begin{align*}
\Delta_g b_+ \geq 0
\end{align*}
distributionally, and applying it to the ray $t \mapsto \gamma(-t)$ gives
\begin{align*}
\Delta_g b_- \geq 0
\end{align*}
distributionally.
Because the distributional Laplacian is linear, the distributional Laplacian of $u$ is
\begin{align*}
\Delta_g u = \Delta_g b_+ + \Delta_g b_-.
\end{align*}
The sum of two nonnegative distributions is nonnegative, so
\begin{align*}
\Delta_g u \geq 0
\end{align*}
in the sense of distributions. Thus $u$ is a locally Lipschitz distributionally subharmonic function on $M$.
[/guided]
[/step]
[step:Apply the strong maximum principle to force constancy]
From the first step, $u \leq 0$ on $M$. From the second step, $u(\gamma(0))=0$. Thus $u$ attains a global maximum at the point $\gamma(0) \in M$.
Since $u$ is locally Lipschitz and distributionally subharmonic, the strong maximum principle for Lipschitz distributionally subharmonic functions applies (citing a result not yet in the wiki: strong maximum principle for distributionally subharmonic functions). It implies that if such a function attains a local maximum on a connected domain, then it is constant on that domain. Because $M$ is connected, $u$ is constant on all of $M$.
The constant is determined by evaluating at $\gamma(0)$:
\begin{align*}
u(x) = u(\gamma(0)) = 0
\end{align*}
for every $x \in M$. Therefore
\begin{align*}
b_+(x)+b_-(x)=0
\end{align*}
for every $x \in M$, as required.
[/step]
