[proofplan]
The proof uses the Bochner identity for harmonic maps. Since the domain has nonnegative Ricci curvature and the target has nonpositive sectional curvature, the curvature terms in the Bochner formula have the favorable sign, giving $\Delta_g e(u) \geq |\nabla du|_{g,h}^2$. Integrating over the compact boundaryless manifold kills the Laplacian term by [Stokes' theorem](/theorems/1530), forcing $|\nabla du|_{g,h}^2$ to vanish everywhere. Once $du$ is parallel, the energy density is parallel as a scalar function and hence constant.
[/proofplan]
[step:Apply the harmonic map Bochner identity with the curvature signs]
Let $\nabla^M$ and $\nabla^N$ denote the Levi-Civita connections of $(M,g)$ and $(N,h)$, respectively. Let
$E := u^{-1}TN$
denote the pullback tangent bundle over $M$, equipped with the pullback connection induced by $\nabla^N$. The differential $du$ is a smooth section of $T^*M \otimes E$, and its covariant derivative $\nabla du$ is a smooth section of $T^*M \otimes T^*M \otimes E$. Let $\Delta_g: C^\infty(M) \to C^\infty(M)$ denote the Laplace-Beltrami operator associated to $g$, with the sign convention $\Delta_g f := \operatorname{div}_g(\nabla f)$ for every $f \in C^\infty(M)$, and let $R^N$ denote the Riemann curvature tensor of $\nabla^N$.
By the [Eells-Sampson Bochner formula for harmonic maps](/page/Eells-Sampson%20Bochner%20Formula) applied to the smooth harmonic map $u$, define the energy density as the smooth map $e(u): M \to \mathbb{R}$ whose value at $p \in M$ is
\begin{align*}
e(u)(p) := \frac{1}{2}|du_p|_{g,h}^2.
\end{align*}
The stated Riemannian metrics provide the Levi-Civita and pullback connections required in the formula. We use the curvature convention
\begin{align*}
K_N(A,B) := \frac{h(R^N(A,B)B,A)}{|A|_h^2|B|_h^2 - h(A,B)^2}
\end{align*}
for linearly independent vectors $A,B \in T_qN$ at a point $q \in N$. With $m := \dim M$ and any local $g$-orthonormal frame $(e_1,\dots,e_m)$ on $M$, the identity gives
\begin{align*}
\Delta_g e(u)
=
|\nabla du|_{g,h}^2
+
\sum_{i=1}^{m} h\bigl(du(\operatorname{Ric}^M(e_i)),du(e_i)\bigr)
-
\sum_{i,j=1}^{m} h\bigl(R^N(du(e_i),du(e_j))du(e_j),du(e_i)\bigr).
\end{align*}
Since $\operatorname{Ric}^M \geq 0$, the Ricci endomorphism is self-adjoint and positive semidefinite. At each point $p \in M$, choose a $g$-orthonormal eigenbasis $(E_1,\dots,E_m)$ of $T_pM$ with eigenvalues $\lambda_1,\dots,\lambda_m \geq 0$. Then the Ricci contribution at $p$ is
\begin{align*}
\sum_{a=1}^{m} h\bigl(du(\operatorname{Ric}^M(E_a)),du(E_a)\bigr)
=
\sum_{a=1}^{m} \lambda_a |du(E_a)|_h^2
\geq 0.
\end{align*}
Since $K_N \leq 0$, the target curvature term appears with a nonnegative sign. If $du(e_i)$ and $du(e_j)$ are linearly independent, the definition of sectional curvature gives
\begin{align*}
-
h\bigl(R^N(du(e_i),du(e_j))du(e_j),du(e_i)\bigr)
\geq 0.
\end{align*}
If $du(e_i)$ and $du(e_j)$ are linearly dependent, then skew-symmetry and multilinearity of $R^N$ in its first two arguments give
\begin{align*}
h\bigl(R^N(du(e_i),du(e_j))du(e_j),du(e_i)\bigr)=0.
\end{align*}
Thus the same inequality holds for every pair $i,j$. Therefore
\begin{align*}
\Delta_g e(u) \geq |\nabla du|_{g,h}^2 \geq 0.
\end{align*}
[guided]
Let $\nabla^M$ and $\nabla^N$ denote the Levi-Civita connections of $(M,g)$ and $(N,h)$, respectively. The energy density is the smooth map $e(u): M \to \mathbb{R}$ whose value at $p \in M$ is
\begin{align*}
e(u)(p) := \frac{1}{2}|du_p|_{g,h}^2.
\end{align*}
The object whose vanishing we want is not $du$ itself, but its covariant derivative. More precisely, if $E := u^{-1}TN$ is the pullback tangent bundle equipped with the pullback connection induced by $\nabla^N$, then $du$ is a smooth section of $T^*M \otimes E$, and $\nabla du$ is a smooth section of $T^*M \otimes T^*M \otimes E$. Let $\Delta_g: C^\infty(M) \to C^\infty(M)$ denote the Laplace-Beltrami operator associated to $g$, with the sign convention $\Delta_g f := \operatorname{div}_g(\nabla f)$ for every $f \in C^\infty(M)$, and let $R^N$ denote the Riemann curvature tensor of $\nabla^N$.
We now use the [Eells-Sampson Bochner formula for harmonic maps](/page/Eells-Sampson%20Bochner%20Formula). The hypotheses needed for that identity are that $u$ is smooth and harmonic and that the Riemannian metrics provide the Levi-Civita and pullback connections appearing above; smoothness and harmonicity are part of the theorem statement, while the connections are supplied by the Riemannian structures on $(M,g)$ and $(N,h)$. We use the curvature convention
\begin{align*}
K_N(A,B) := \frac{h(R^N(A,B)B,A)}{|A|_h^2|B|_h^2 - h(A,B)^2}
\end{align*}
for linearly independent vectors $A,B \in T_qN$ at a point $q \in N$. For a local $g$-orthonormal frame $(e_1,\dots,e_m)$ on $M$, where $m := \dim M$, the identity gives
\begin{align*}
\Delta_g e(u)
=
|\nabla du|_{g,h}^2
+
\sum_{i=1}^{m} h\bigl(du(\operatorname{Ric}^M(e_i)),du(e_i)\bigr)
-
\sum_{i,j=1}^{m} h\bigl(R^N(du(e_i),du(e_j))du(e_j),du(e_i)\bigr).
\end{align*}
The point of the curvature assumptions is exactly to control the two curvature terms. The condition $\operatorname{Ric}^M \geq 0$ says that the Ricci endomorphism is self-adjoint and positive semidefinite. To see the sign directly, fix $p \in M$ and choose a $g$-orthonormal eigenbasis $(E_1,\dots,E_m)$ of $T_pM$ with eigenvalues $\lambda_1,\dots,\lambda_m \geq 0$. In that basis the Ricci contribution is
\begin{align*}
\sum_{a=1}^{m} h\bigl(du(\operatorname{Ric}^M(E_a)),du(E_a)\bigr)
=
\sum_{a=1}^{m} \lambda_a |du(E_a)|_h^2
\geq 0.
\end{align*}
The condition $K_N \leq 0$ says that the sectional curvature of every $2$-plane in $T_{u(p)}N$ is nonpositive. If $du(e_i)$ and $du(e_j)$ are linearly independent, this gives
\begin{align*}
-
h\bigl(R^N(du(e_i),du(e_j))du(e_j),du(e_i)\bigr)
\geq 0.
\end{align*}
If they are linearly dependent, then skew-symmetry and multilinearity of $R^N$ in its first two arguments force
\begin{align*}
h\bigl(R^N(du(e_i),du(e_j))du(e_j),du(e_i)\bigr)=0,
\end{align*}
so the same nonnegativity conclusion still holds. Thus all curvature contributions have the favorable sign, and the identity implies the pointwise estimate
\begin{align*}
\Delta_g e(u) \geq |\nabla du|_{g,h}^2 \geq 0.
\end{align*}
[/guided]
[/step]
[step:Integrate the Bochner inequality over the compact boundaryless domain]
Let $d\operatorname{vol}_g$ denote the Riemannian volume measure on $M$. Since $M$ is compact and has no boundary, [Stokes' theorem](/page/Stokes%20Theorem) applied to the vector field $\nabla e(u)$ gives
\begin{align*}
\int_M \Delta_g e(u)\,d\operatorname{vol}_g = 0.
\end{align*}
Integrating the pointwise inequality from the previous step with respect to $d\operatorname{vol}_g$ gives
\begin{align*}
0
=
\int_M \Delta_g e(u)\,d\operatorname{vol}_g
\geq
\int_M |\nabla du|_{g,h}^2\,d\operatorname{vol}_g
\geq 0.
\end{align*}
Therefore
\begin{align*}
\int_M |\nabla du|_{g,h}^2\,d\operatorname{vol}_g = 0.
\end{align*}
[/step]
[step:Conclude that $du$ is parallel everywhere]
The function $|\nabla du|_{g,h}^2: M \to \mathbb{R}$ defined by
\begin{align*}
|\nabla du|_{g,h}^2(p) := |\nabla du(p)|_{g,h}^2
\end{align*}
is continuous and nonnegative because $u$ is smooth. The Riemannian volume measure $d\operatorname{vol}_g$ has full support on $M$: every nonempty open subset of $M$ has positive Riemannian volume. Since the integral of this continuous nonnegative function over $M$ with respect to $d\operatorname{vol}_g$ is zero, it vanishes identically:
\begin{align*}
|\nabla du|_{g,h}^2 = 0
\end{align*}
on $M$. Hence
\begin{align*}
\nabla du = 0
\end{align*}
as a section of $T^*M \otimes T^*M \otimes u^{-1}TN$.
[/step]
[step:Differentiate the energy density using the parallelness of $du$]
Let $X \in \mathfrak{X}(M)$ be an arbitrary smooth vector field. Using metric compatibility of the induced connection on $T^*M \otimes u^{-1}TN$, we compute
\begin{align*}
X(e(u))
=
X\left(\frac{1}{2}|du|_{g,h}^2\right)
=
\bigl(\nabla_X du,du\bigr)_{g,h}.
\end{align*}
Since $\nabla du = 0$, the right-hand side is zero for every smooth vector field $X$. Thus $\nabla e(u)=0$ as a smooth real-valued function on $M$. Since $M$ is connected by the revised theorem statement, every smooth real-valued function with zero gradient is constant on $M$.
This proves both conclusions: the differential $du$ is parallel, and the energy density $e(u)$ is constant.
[/step]