[proofplan] We use the [Bochner formula for harmonic maps](/theorems/5688) to express $\Delta_g e(u)$ as a sum of three terms: the square norm $|\nabla du|^2$, a domain Ricci-curvature term, and a target curvature term with a minus sign. The Ricci term is nonnegative because $\operatorname{Ric}^M$ is a nonnegative symmetric [bilinear form](/page/Bilinear%20Form). The target curvature term is also nonnegative after the minus sign because nonpositive sectional curvature implies $h(R^N(a,b)b,a)\leq 0$ for all tangent vectors $a,b$ in $N$. Combining these sign estimates gives the inequality and hence subharmonicity. [/proofplan]

[step:Write the harmonic-map Bochner formula at a fixed point] Fix $p\in M$. Let $m=\dim M$, and choose a $g_p$-[orthonormal basis](/page/Orthonormal%20Basis) $(e_1,\dots,e_m)$ of $T_pM$. Let $\nabla du$ denote the covariant derivative of $du$, viewed as a section of $T^*M\otimes u^{-1}TN$, using the Levi-Civita connection of $g$ on $T^*M$ and the pullback Levi-Civita connection induced by $h$ on $u^{-1}TN$. Let $R^N$ denote the Riemann curvature tensor of $(N,h)$ with the convention that the sectional curvature of a $2$-plane spanned by linearly independent vectors $a,b\in T_qN$ is \begin{align*} K_N(a\wedge b)=\frac{h_q(R^N(a,b)b,a)}{|a|_h^2|b|_h^2-h_q(a,b)^2}. \end{align*} The hypotheses needed for the harmonic-map Bochner formula are exactly that $(M,g)$ and $(N,h)$ are smooth Riemannian manifolds and that $u:(M,g)\to(N,h)$ is a smooth harmonic map. These are part of the theorem statement, so the Bochner formula for harmonic maps gives, at $p$, \begin{align*} \Delta_g e(u)(p)=|\nabla du|_{g,h}^2(p)+\sum_{i=1}^{m} h_{u(p)}\bigl(du_p(\operatorname{Ric}^{M,\sharp}_p e_i),du_p(e_i)\bigr)-\sum_{i,j=1}^{m} h_{u(p)}\bigl(R^N(du_p(e_i),du_p(e_j))du_p(e_j),du_p(e_i)\bigr), \end{align*} where $\operatorname{Ric}^{M,\sharp}_p:T_pM\to T_pM$ is the $g_p$-self-adjoint [linear map](/page/Linear%20Map) determined by \begin{align*} g_p(\operatorname{Ric}^{M,\sharp}_p v,w)=\operatorname{Ric}^M_p(v,w) \end{align*} for all $v,w\in T_pM$. This is the standard Bochner identity for harmonic maps, cited here by name because the corresponding Androma theorem page is not yet available. [/step]

[step:Show that the Ricci-curvature contribution is nonnegative] Because $\operatorname{Ric}^M_p$ is a symmetric nonnegative bilinear form on the finite-dimensional [inner product space](/page/Inner%20Product%20Space) $(T_pM,g_p)$, there exists a $g_p$-orthonormal basis $(f_1,\dots,f_m)$ of $T_pM$ and numbers $\rho_1,\dots,\rho_m\in [0,\infty)$ such that \begin{align*} \operatorname{Ric}^{M,\sharp}_p f_i=\rho_i f_i \end{align*} for every $i\in\{1,\dots,m\}$. Rewriting the Ricci term in this basis gives \begin{align*} \sum_{i=1}^{m} h_{u(p)}\bigl(du_p(\operatorname{Ric}^{M,\sharp}_p f_i),du_p(f_i)\bigr)=\sum_{i=1}^{m} h_{u(p)}\bigl(du_p(\rho_i f_i),du_p(f_i)\bigr)=\sum_{i=1}^{m} \rho_i |du_p(f_i)|_h^2\geq 0. \end{align*} Thus the domain curvature contribution in the Bochner formula is nonnegative. [/step]

[step:Show that the target-curvature contribution is nonnegative after the Bochner sign]For $i,j\in\{1,\dots,m\}$, define $a_{ij}:=du_p(e_i)\in T_{u(p)}N$ and $b_{ij}:=du_p(e_j)\in T_{u(p)}N$. If $a_{ij}$ and $b_{ij}$ are linearly independent, the definition of sectional curvature gives \begin{align*} h_{u(p)}(R^N(a_{ij},b_{ij})b_{ij},a_{ij}) = K_N(a_{ij}\wedge b_{ij}) \bigl(|a_{ij}|_h^2|b_{ij}|_h^2-h_{u(p)}(a_{ij},b_{ij})^2\bigr). \end{align*} The Gram determinant \begin{align*} |a_{ij}|_h^2|b_{ij}|_h^2-h_{u(p)}(a_{ij},b_{ij})^2 \end{align*} is positive in the linearly independent case by the [Cauchy-Schwarz inequality](/theorems/432), and $K_N(a_{ij}\wedge b_{ij})\leq 0$ by hypothesis. Hence \begin{align*} h_{u(p)}(R^N(a_{ij},b_{ij})b_{ij},a_{ij})\leq 0. \end{align*} If $a_{ij}$ and $b_{ij}$ are linearly dependent, the alternating property of $R^N$ in its first two entries gives \begin{align*} R^N(a_{ij},b_{ij})=0, \end{align*} so the same inequality holds. Therefore \begin{align*} - \sum_{i,j=1}^{m} h_{u(p)} \bigl( R^N(du_p(e_i),du_p(e_j))du_p(e_j), du_p(e_i) \bigr) \geq 0. \end{align*}[/step]

[guided]The target term has the sign that is most important in the argument, so we check it one pair of indices at a time. For fixed $i,j\in\{1,\dots,m\}$, define the two tangent vectors at the point $u(p)\in N$ by $a_{ij}:=du_p(e_i)\in T_{u(p)}N$ and $b_{ij}:=du_p(e_j)\in T_{u(p)}N$. We need to determine the sign of \begin{align*} h_{u(p)}(R^N(a_{ij},b_{ij})b_{ij},a_{ij}). \end{align*} First suppose $a_{ij}$ and $b_{ij}$ are linearly independent. Then they span a $2$-plane in $T_{u(p)}N$, and the definition of sectional curvature gives \begin{align*} h_{u(p)}(R^N(a_{ij},b_{ij})b_{ij},a_{ij}) = K_N(a_{ij}\wedge b_{ij}) \bigl(|a_{ij}|_h^2|b_{ij}|_h^2-h_{u(p)}(a_{ij},b_{ij})^2\bigr). \end{align*} The second factor is the Gram determinant of the two vectors. By the Cauchy-Schwarz inequality, \begin{align*} h_{u(p)}(a_{ij},b_{ij})^2\leq |a_{ij}|_h^2|b_{ij}|_h^2, \end{align*} and because the vectors are linearly independent, equality does not occur. Thus \begin{align*} |a_{ij}|_h^2|b_{ij}|_h^2-h_{u(p)}(a_{ij},b_{ij})^2>0. \end{align*} The curvature hypothesis says $K_N(a_{ij}\wedge b_{ij})\leq 0$, so the product is nonpositive: \begin{align*} h_{u(p)}(R^N(a_{ij},b_{ij})b_{ij},a_{ij})\leq 0. \end{align*} Now suppose $a_{ij}$ and $b_{ij}$ are linearly dependent. If one of the two vectors is zero, then $R^N(a_{ij},b_{ij})=0$ by bilinearity. If neither is zero, there is a scalar $\lambda\in\mathbb{R}$ such that $b_{ij}=\lambda a_{ij}$, and the alternating property of the curvature tensor in its first two arguments gives \begin{align*} R^N(a_{ij},b_{ij}) = R^N(a_{ij},\lambda a_{ij}) = \lambda R^N(a_{ij},a_{ij}) = 0. \end{align*} Therefore the same inequality holds in the dependent case: \begin{align*} h_{u(p)}(R^N(a_{ij},b_{ij})b_{ij},a_{ij})\leq 0. \end{align*} Since this estimate holds for every pair $(i,j)$, summing preserves the inequality: \begin{align*} \sum_{i,j=1}^{m} h_{u(p)} \bigl( R^N(du_p(e_i),du_p(e_j))du_p(e_j), du_p(e_i) \bigr) \leq 0. \end{align*} Multiplying by the minus sign appearing in the Bochner formula reverses the sign and gives \begin{align*} - \sum_{i,j=1}^{m} h_{u(p)} \bigl( R^N(du_p(e_i),du_p(e_j))du_p(e_j), du_p(e_i) \bigr) \geq 0. \end{align*} This is exactly where the nonpositive curvature of the target is used.[/guided]

[step:Combine the sign estimates to obtain subharmonicity] Substituting the two nonnegativity estimates into the Bochner formula at $p$ gives \begin{align*} \Delta_g e(u)(p)=|\nabla du|_{g,h}^2(p)+\sum_{i=1}^{m} h_{u(p)}\bigl(du_p(\operatorname{Ric}^{M,\sharp}_p e_i),du_p(e_i)\bigr)-\sum_{i,j=1}^{m} h_{u(p)}\bigl(R^N(du_p(e_i),du_p(e_j))du_p(e_j),du_p(e_i)\bigr)\geq |\nabla du|_{g,h}^2(p). \end{align*} The square norm satisfies $|\nabla du|_{g,h}^2(p)\geq 0$, so \begin{align*} \Delta_g e(u)(p)\geq |\nabla du|_{g,h}^2(p)\geq 0. \end{align*} Since $p\in M$ was arbitrary, the inequality holds pointwise on $M$. By the definition of subharmonicity for smooth functions with respect to the positive Laplace-Beltrami operator, $\Delta_g e(u)\geq 0$ implies that $e(u)$ is subharmonic. [/step]