[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]