[proofplan]
The Cartan-Hadamard theorem gives a global smoothly endpoint-dependent geodesic between any two points of $N$, so the pointwise geodesics define a smooth geodesic homotopy. The energy along such a homotopy is convex because the second variation integrand is a sum of a square and a nonnegative curvature contribution. Since the endpoint maps are harmonic and $M$ has no boundary, the first variation gives zero derivative at both endpoints; convexity then forces the energy to be constant. In negative curvature, equality in the second variation forces all spatial derivatives of the homotopy to be collinear with the variation field, so a nonzero variation field would put the whole image in one geodesic; the stated exception therefore implies the variation field vanishes.
[/proofplan]
[step:Construct the smooth pointwise geodesic homotopy]
By the [Cartan-Hadamard theorem](/page/Cartan-Hadamard%20Theorem), applied to the complete simply connected Riemannian manifold $(N,h)$ with $K_N\leq 0$, the exponential map
\begin{align*}
\exp_q: T_qN \to N
\end{align*}
is a diffeomorphism for every $q\in N$. Hence for every ordered pair $(q_0,q_1)\in N\times N$ there is a unique constant-speed geodesic
\begin{align*}
\Gamma(q_0,q_1,\cdot):[0,1]\to N
\end{align*}
with $\Gamma(q_0,q_1,0)=q_0$ and $\Gamma(q_0,q_1,1)=q_1$. The endpoint geodesic map
\begin{align*}
\Gamma:N\times N\times[0,1]\to N
\end{align*}
is given explicitly by
\begin{align*}
\Gamma(q_0,q_1,t)=\exp_{q_0}\left(t\,\exp_{q_0}^{-1}(q_1)\right).
\end{align*}
Cartan-Hadamard gives that $\exp_{q_0}$ is a global diffeomorphism for every $q_0\in N$, so the inverse endpoint vector $\exp_{q_0}^{-1}(q_1)\in T_{q_0}N$ depends smoothly on $(q_0,q_1)$. Therefore $\Gamma$ is smooth.
Define the homotopy map
\begin{align*}
H:M\times[0,1]\to N
\end{align*}
by $H(p,t)=\Gamma(u_0(p),u_1(p),t)$. Since $u_0$, $u_1$, and $\Gamma$ are smooth, $H$ is smooth. For each $p\in M$, define
\begin{align*}
H_p:[0,1]\to N
\end{align*}
by $H_p(t)=H(p,t)$. This curve is the unique constant-speed geodesic from $u_0(p)$ to $u_1(p)$. Thus $H$ is the required [geodesic homotopy](/page/Geodesic%20Homotopy).
[/step]
[step:Use first variation at the harmonic endpoints]
For $t\in[0,1]$, define
\begin{align*}
H_t:M\to N
\end{align*}
by $H_t(p)=H(p,t)$, and define the variation field
\begin{align*}
V_t\in \Gamma(H_t^*TN), \qquad V_t(p):=\partial_t H(p,t).
\end{align*}
Let
\begin{align*}
F:[0,1]\to \mathbb{R}
\end{align*}
be the [energy](/page/Energy%20of%20a%20Map) function defined by
\begin{align*}
F(t)=E(H_t):=\frac{1}{2}\int_M |dH_t|_{g,h}^2\, d\operatorname{vol}_g(p).
\end{align*}
The [first variation formula](/theorems/2728) for the energy of maps from a compact manifold without boundary gives, for $0<t<1$,
\begin{align*}
F'(t)=-\int_M h(\tau(H_t),V_t)\, d\operatorname{vol}_g(p),
\end{align*}
and the same formula gives the right derivative $F'_+(0)$ at $0$ and the left derivative $F'_-(1)$ at $1$. Here $\tau(H_t)\in\Gamma(H_t^*TN)$ denotes the [tension field](/page/Tension%20Field). Since $u_0=H_0$ and $u_1=H_1$ are [harmonic maps](/page/Harmonic%20Map), $\tau(H_0)=0$ and $\tau(H_1)=0$. Therefore
\begin{align*}
F'_+(0)=0,\qquad F'_-(1)=0.
\end{align*}
[guided]
We now use exactly where the harmonicity of the endpoint maps enters. Define the energy function
\begin{align*}
F:[0,1]\to \mathbb{R}
\end{align*}
by
\begin{align*}
F(t)=\frac{1}{2}\int_M |dH_t|_{g,h}^2\, d\operatorname{vol}_g(p).
\end{align*}
The velocity of the variation is the section
\begin{align*}
V_t\in\Gamma(H_t^*TN), \qquad V_t(p)=\partial_tH(p,t).
\end{align*}
The first variation formula for harmonic-map energy on a compact manifold without boundary says, for $0<t<1$,
\begin{align*}
F'(t)=-\int_M h(\tau(H_t),V_t)\, d\operatorname{vol}_g(p),
\end{align*}
because there is no boundary term when $\partial M=\varnothing$. Since $F$ is defined on the closed interval $[0,1]$, the endpoint derivatives are one-sided: the same first variation computation gives the right derivative $F'_+(0)$ and the left derivative $F'_-(1)$. Here $\tau(H_t)$ is the [tension field](/page/Tension%20Field) of $H_t$. The endpoint hypotheses say precisely that $H_0=u_0$ and $H_1=u_1$ are [harmonic maps](/page/Harmonic%20Map), meaning
\begin{align*}
\tau(H_0)=0,\qquad \tau(H_1)=0.
\end{align*}
Substituting these equalities into the endpoint first variation formulas gives
\begin{align*}
F'_+(0)=0,\qquad F'_-(1)=0.
\end{align*}
This is the mechanism that turns convexity into constancy: convexity alone would only give monotonicity information, but zero one-sided derivative at both endpoints pins the whole convex function flat.
[/guided]
[/step]
[step:Apply convexity of energy along a geodesic homotopy]
Because $H$ is smooth, $M$ is compact without boundary, and each curve $t\mapsto H(p,t)$ is a geodesic, the [second variation formula](/theorems/2729) for energy has no boundary term and no acceleration term. Thus, for every $t\in(0,1)$,
\begin{align*}
F''(t)=\int_M \sum_{i=1}^m
\left(
|\nabla_{e_i}^{H_t}V_t|_h^2
-
h\left(R^N(V_t,dH_t(e_i))dH_t(e_i),V_t\right)
\right)\, d\operatorname{vol}_g(p),
\end{align*}
where $m=\dim M$, $(e_1,\dots,e_m)$ is any local $g$-orthonormal frame on $M$, $\nabla^{H_t}$ is the pullback connection on $H_t^*TN$, and $R^N$ is the curvature tensor of $(N,h)$.
Since $K_N\leq 0$, for every pair of tangent vectors $\xi,\eta\in T_qN$ at any point $q\in N$,
\begin{align*}
-h(R^N(\xi,\eta)\eta,\xi)\geq 0.
\end{align*}
Thus every summand in the integrand defining $F''(t)$ is nonnegative, and hence
\begin{align*}
F''(t)\geq 0
\end{align*}
for all $t\in(0,1)$. Therefore $F$ is convex on $[0,1]$.
Since $F$ is convex, its one-sided derivatives are monotone nondecreasing in the sense that, for $0<s<t<1$,
\begin{align*}
F'_+(0)\leq F'(s)\leq F'(t)\leq F'_-(1).
\end{align*}
Together with $F'_+(0)=F'_-(1)=0$, this implies $F'(t)=0$ for all $t\in(0,1)$, and therefore $F$ is constant on $[0,1]$.
[/step]
[step:Extract the equality conditions from the second variation]
Since $F$ is constant, $F''(t)=0$ for every $t\in(0,1)$. The integrand in the second variation formula is a sum of nonnegative continuous functions. Therefore, for every $p\in M$, every $t\in(0,1)$, and every local $g$-orthonormal frame vector $e_i$ at $p$,
\begin{align*}
|\nabla_{e_i}^{H_t}V_t|_h^2=0
\end{align*}
and
\begin{align*}
-h\left(R^N(V_t,dH_t(e_i))dH_t(e_i),V_t\right)=0.
\end{align*}
Hence
\begin{align*}
\nabla_{e_i}^{H_t}V_t=0.
\end{align*}
If $K_N<0$, the curvature equality also implies that $V_t(p)$ and $dH_t(e_i)$ are linearly dependent in $T_{H(p,t)}N$ for every $i$, because strict negative sectional curvature makes
\begin{align*}
-h(R^N(\xi,\eta)\eta,\xi)>0
\end{align*}
whenever $\xi$ and $\eta$ span a two-dimensional plane.
[/step]
[step:Show that a nonzero variation field forces the image into one geodesic]
Assume $K_N<0$. Define the global variation field along $H$ by
\begin{align*}
V\in \Gamma(H^*TN), \qquad V(p,t):=\partial_tH(p,t)=V_t(p).
\end{align*}
Suppose there exists $(p_0,t_0)\in M\times(0,1)$ such that $V(p_0,t_0)\neq 0$. Since
\begin{align*}
\nabla_{e_i}^{H_t}V_t=0
\end{align*}
for all spatial directions $e_i$, the function
\begin{align*}
p\mapsto |V_t(p)|_h^2
\end{align*}
has zero derivative in every direction on the connected manifold $M$. Hence $|V_t|_h$ is constant on $M$ for each fixed $t$. Since each curve $t\mapsto H(p,t)$ is a constant-speed geodesic, $\nabla_{\partial_t}^{H}V=0$, so $|V(p,t)|_h$ is also independent of $t$. Thus $V(p,t)\neq 0$ for every $(p,t)\in M\times[0,1]$.
By the equality condition from negative curvature, for every tangent vector $X\in T_pM$ there exists a scalar $a_X(p,t)\in\mathbb{R}$ such that
\begin{align*}
dH_{(p,t)}(X,0)=a_X(p,t)V(p,t).
\end{align*}
Also
\begin{align*}
dH_{(p,t)}(0,\partial_t)=V(p,t).
\end{align*}
Thus every differential of $H$ is tangent to the one-dimensional subspace spanned by $V(p,t)$.
We now prove that this one-dimensional tangency forces the image to lie in one geodesic. Let
\begin{align*}
\alpha:[0,1]\to M\times[0,1]
\end{align*}
be any smooth path, and define
\begin{align*}
\beta:[0,1]\to N
\end{align*}
by $\beta(s)=H(\alpha(s))$. Since $dH_{\alpha(s)}(\alpha'(s))$ is collinear with $V(\alpha(s))$, there is a smooth scalar function
\begin{align*}
b:[0,1]\to\mathbb{R}
\end{align*}
such that
\begin{align*}
\beta'(s)=b(s)V(\alpha(s)).
\end{align*}
The field $V$ is parallel along $\alpha$ because $\nabla_{\partial_t}^{H}V=0$ and $\nabla_X^{H_t}V_t=0$ for every $X\in T_pM$. Therefore the vector field $W$ along $\beta$ defined by $W(s)=V(\alpha(s))$ is parallel. Let $\mathcal{L}^1$ denote one-dimensional [Lebesgue measure](/page/Lebesgue%20Measure) on $[0,1]$, and define
\begin{align*}
r:[0,1]\to\mathbb{R}, \qquad r(s):=\int_0^s b(\sigma)\,d\mathcal{L}^1(\sigma).
\end{align*}
Parallel transport along $\beta$ carries $W(s)$ back to $W(0)$, so the velocity direction of $\beta$ remains fixed after this transport. Let $\gamma:\mathbb{R}\to N$ be the geodesic with initial data $\gamma(0)=\beta(0)$ and $\gamma'(0)=W(0)$. Define
\begin{align*}
\widetilde{\beta}:[0,1]\to N, \qquad \widetilde{\beta}(s):=\gamma(r(s)).
\end{align*}
Along $\widetilde{\beta}$, define $\widetilde{W}(s)$ to be the parallel translate of $W(0)$ along $\gamma$ to the point $\gamma(r(s))$. Then
\begin{align*}
\widetilde{\beta}'(s)=r'(s)\widetilde{W}(s)=b(s)\widetilde{W}(s),
\end{align*}
and $\widetilde{W}$ is parallel along $\widetilde{\beta}$. Thus the pair $(\widetilde{\beta},\widetilde{W})$ solves the same first-order system
\begin{align*}
\zeta'(s)=b(s)Y(s),\qquad \nabla_{\zeta'}Y=0
\end{align*}
with initial data $\zeta(0)=\beta(0)$ and $Y(0)=W(0)$ as the pair $(\beta,W)$. By uniqueness for this coupled curve-and-parallel-field system, $\beta(s)=\widetilde{\beta}(s)=\gamma(r(s))$ for all $s\in[0,1]$. Hence the image of $\beta$ is contained in the image of the single geodesic $\gamma$.
Because $M$ is connected and $[0,1]$ is connected, the product $M\times[0,1]$ is path-connected. Applying the preceding paragraph to paths from a fixed base point in $M\times[0,1]$ to arbitrary points shows that every point of $H(M\times[0,1])$ lies on the same geodesic image in $N$.
[/step]
[step:Conclude uniqueness outside the geodesic-image exception]
Assume now that $K_N<0$ and that $H(M\times[0,1])$ is not contained in the image of a single geodesic. The previous step shows that the alternative $V_t(p)\neq 0$ somewhere would force the whole image of $H$ to be contained in one geodesic. Hence
\begin{align*}
V_t(p)=0
\end{align*}
for every $(p,t)\in M\times[0,1]$.
Therefore each curve
\begin{align*}
H_p:[0,1]\to N
\end{align*}
defined by $H_p(t)=H(p,t)$ has zero velocity, and so it is constant. For every $p\in M$,
\begin{align*}
u_0(p)=H(p,0)=H(p,1)=u_1(p).
\end{align*}
Thus $u_0=u_1$, completing the proof.
[/step]