[proofplan]
The proof is a pointwise computation in a normal orthonormal frame on $M$. We differentiate the identity $e(u)=\frac12\sum_i |du(e_i)|^2$ twice, using metric compatibility to produce the positive term $|\nabla du|^2$ and a remaining rough-Laplacian contraction. We then commute the covariant derivatives of $du$; the domain curvature contracts to the Ricci term, the target curvature gives the negative sectional-curvature term, and the derivative of the tension field vanishes because $u$ is harmonic. Since every term is tensorial, the pointwise normal-frame calculation proves the formula globally.
[/proofplan]
[step:Declare the induced connection and reduce to a normal-frame computation]
Let $E:=u^*TN$ be the [pullback vector bundle](/page/Pullback%20Bundle) over $M$. The connection induced by the Levi-Civita connection $\nabla^N$ on $E$ is denoted by $\nabla^E$, and the induced tensor-product connection on $T^*M\otimes E$ is denoted by $\nabla$. Thus $du$ is a smooth section
\begin{align*}
du \in \Gamma(T^*M\otimes E),
\end{align*}
and its covariant derivative is the section
\begin{align*}
\nabla du \in \Gamma(T^*M\otimes T^*M\otimes E)
\end{align*}
defined by
\begin{align*}
(\nabla du)(X,Y)
:=
\nabla^E_X(du(Y))-du(\nabla^M_XY)
\end{align*}
for smooth vector fields $X,Y\in\mathfrak X(M)$.
The tension field of $u$ is the section
\begin{align*}
\tau(u)\in\Gamma(E),
\qquad
\tau(u):=\sum_{i=1}^m(\nabla du)(e_i,e_i),
\end{align*}
where $e_1,\dots,e_m$ is any local $g$-orthonormal frame. The harmonicity assumption means
\begin{align*}
\tau(u)=0.
\end{align*}
Fix a point $p\in M$. Choose a local $g$-orthonormal frame $e_1,\dots,e_m$ near $p$ such that
\begin{align*}
\nabla^M_{e_i}e_j(p)=0
\end{align*}
for all $1\leq i,j\leq m$. Because $\nabla^M$ is torsion-free, this also gives
\begin{align*}
[e_i,e_j](p)=0.
\end{align*}
The desired formula is tensorial in the frame and in the point $p$, so it suffices to prove it at $p$ in this frame.
[/step]
[step:Differentiate the energy density twice]
For $1\leq j\leq m$, define the local section $A_j\in \Gamma(E)$ by
\begin{align*}
A_j:=du(e_j).
\end{align*}
For $1\leq i,j\leq m$, define the local section $B_{ij}\in\Gamma(E)$ by
\begin{align*}
B_{ij}:=(\nabla du)(e_i,e_j).
\end{align*}
By the definition of the induced connection on $T^*M\otimes E$,
\begin{align*}
B_{ij}=\nabla^E_{e_i}A_j-du(\nabla^M_{e_i}e_j).
\end{align*}
At the point $p$, the normal-frame condition gives
\begin{align*}
B_{ij}(p)=\nabla^E_{e_i}A_j(p).
\end{align*}
For $1\leq i,j,k\leq m$, define the smooth real-valued connection coefficient $\Gamma_{ij}^k$ by
\begin{align*}
\nabla^M_{e_i}e_j=\sum_{k=1}^m \Gamma_{ij}^k e_k.
\end{align*}
Metric compatibility of $\nabla^M$ and the $g$-orthonormality of the frame give
\begin{align*}
0=e_i g(e_j,e_k)=\Gamma_{ij}^k+\Gamma_{ik}^j,
\end{align*}
so $\Gamma_{ij}^k=-\Gamma_{ik}^j$. Since
\begin{align*}
e(u)=\frac12\sum_{j=1}^m h(A_j,A_j),
\end{align*}
metric compatibility of $\nabla^E$ gives, on the neighbourhood where the frame is defined,
\begin{align*}
e_i e(u)=\sum_{j=1}^m h(\nabla^E_{e_i}A_j,A_j).
\end{align*}
Using $\nabla^E_{e_i}A_j=B_{ij}+du(\nabla^M_{e_i}e_j)$, this becomes
\begin{align*}
e_i e(u)=\sum_{j=1}^m h(B_{ij},A_j)+\sum_{j,k=1}^m \Gamma_{ij}^k h(A_k,A_j).
\end{align*}
The last sum vanishes identically for each $i$, because $\Gamma_{ij}^k$ is skew-symmetric in $j,k$ while $h(A_k,A_j)$ is symmetric in $j,k$. Hence
\begin{align*}
e_i e(u)=\sum_{j=1}^m h(B_{ij},A_j)
\end{align*}
as an identity of smooth functions on the frame neighbourhood. Since $\Delta_g f(p)=\sum_{i=1}^m e_i e_i f(p)$ for every smooth real-valued function $f:M\to\mathbb R$ in this normal frame, differentiating this identity and using metric compatibility of $\nabla^E$ gives
\begin{align*}
\Delta_g e(u)(p)=\sum_{i,j=1}^m h(\nabla^E_{e_i}B_{ij},A_j)(p)+\sum_{i,j=1}^m h(B_{ij},\nabla^E_{e_i}A_j)(p).
\end{align*}
At $p$, the normal-frame condition gives $\nabla^E_{e_i}A_j(p)=B_{ij}(p)$, so
\begin{align*}
\Delta_g e(u)(p)
=
\sum_{i,j=1}^m h(\nabla^E_{e_i}B_{ij},A_j)(p)
+
\sum_{i,j=1}^m h(B_{ij},B_{ij})(p).
\end{align*}
The second sum is the pointwise tensor norm $|\nabla du|^2(p)$ in the orthonormal frame. Thus
\begin{align*}
\Delta_g e(u)(p)
=
\sum_{j=1}^m h(C_j,A_j)(p)+|\nabla du|^2(p),
\end{align*}
where the vector $C_j\in E_p$ is defined by
\begin{align*}
C_j:=\sum_{i=1}^m\nabla^E_{e_i}B_{ij}(p).
\end{align*}
[guided]
We now compute the Laplacian of the scalar function $e(u)$ without differentiating the normal-frame identity $B_{ij}(p)=\nabla^E_{e_i}A_j(p)$ as if it held to first order. The reason for introducing $B_{ij}$ is precisely to keep the calculation covariant: $B_{ij}$ is the component of the tensor $\nabla du$, while $\nabla^E_{e_i}A_j$ alone would miss derivatives of the moving frame.
For each $j$, define $A_j:=du(e_j)\in\Gamma(E)$, and for each pair $i,j$, define $B_{ij}:=(\nabla du)(e_i,e_j)\in\Gamma(E)$. By the definition of the induced connection on $T^*M\otimes E$,
\begin{align*}
B_{ij}=\nabla^E_{e_i}(du(e_j))-du(\nabla^M_{e_i}e_j).
\end{align*}
At the fixed point $p$, the normal-frame condition gives $\nabla^M_{e_i}e_j(p)=0$, so
\begin{align*}
B_{ij}(p)=\nabla^E_{e_i}A_j(p).
\end{align*}
This equality is used only at the point $p$, not differentiated as an identity of sections.
The energy density is
\begin{align*}
e(u)=\frac12|du|^2=\frac12\sum_{j=1}^m h(A_j,A_j).
\end{align*}
The subtle point is that $B_{ij}(p)=\nabla^E_{e_i}A_j(p)$ is only a pointwise equality. Away from $p$ there is an extra moving-frame term. To account for it, define smooth real-valued functions $\Gamma_{ij}^k$ by
\begin{align*}
\nabla^M_{e_i}e_j=\sum_{k=1}^m \Gamma_{ij}^k e_k.
\end{align*}
Because the frame is $g$-orthonormal and $\nabla^M$ is metric-compatible,
\begin{align*}
0=e_i g(e_j,e_k)=g(\nabla^M_{e_i}e_j,e_k)+g(e_j,\nabla^M_{e_i}e_k)=\Gamma_{ij}^k+\Gamma_{ik}^j.
\end{align*}
Thus $\Gamma_{ij}^k=-\Gamma_{ik}^j$: for fixed $i$, the connection coefficients are skew-symmetric in the frame indices $j,k$.
Metric compatibility of $\nabla^E$ gives
\begin{align*}
e_i e(u)=\sum_{j=1}^m h(\nabla^E_{e_i}A_j,A_j).
\end{align*}
Using the identity $\nabla^E_{e_i}A_j=B_{ij}+du(\nabla^M_{e_i}e_j)$, we obtain
\begin{align*}
e_i e(u)=\sum_{j=1}^m h(B_{ij}+du(\nabla^M_{e_i}e_j),A_j).
\end{align*}
Substituting $\nabla^M_{e_i}e_j=\sum_{k=1}^m\Gamma_{ij}^k e_k$ gives
\begin{align*}
e_i e(u)=\sum_{j=1}^m h(B_{ij},A_j)+\sum_{j,k=1}^m \Gamma_{ij}^k h(A_k,A_j).
\end{align*}
Why does the second sum disappear? The coefficient $\Gamma_{ij}^k$ changes sign when $j$ and $k$ are interchanged, while $h(A_k,A_j)=h(A_j,A_k)$ is symmetric. Pairing the $(j,k)$ and $(k,j)$ terms therefore cancels the whole sum, including the diagonal terms because $\Gamma_{ij}^j=0$. Hence, as an identity of functions on the frame neighbourhood,
\begin{align*}
e_i e(u)=\sum_{j=1}^m h(B_{ij},A_j).
\end{align*}
This is the identity we may safely differentiate.
At $p$, the scalar Laplacian in a normal orthonormal frame is
\begin{align*}
\Delta_g f(p)=\sum_{i=1}^m e_i e_i f(p)
\end{align*}
for every smooth real-valued function $f:M\to\mathbb R$. Differentiating the preceding identity gives
\begin{align*}
\Delta_g e(u)(p)=\sum_{i,j=1}^m e_i h(B_{ij},A_j)(p).
\end{align*}
Using metric compatibility of $\nabla^E$ on each summand, this becomes
\begin{align*}
\Delta_g e(u)(p)=\sum_{i,j=1}^m h(\nabla^E_{e_i}B_{ij},A_j)(p)+\sum_{i,j=1}^m h(B_{ij},\nabla^E_{e_i}A_j)(p).
\end{align*}
Only now do we use the normal-frame equality at the point $p$: since $\nabla^M_{e_i}e_j(p)=0$, we have $\nabla^E_{e_i}A_j(p)=B_{ij}(p)$. Therefore
\begin{align*}
\Delta_g e(u)(p)=\sum_{j=1}^m h\left(\sum_{i=1}^m\nabla^E_{e_i}B_{ij},A_j\right)(p)+\sum_{i,j=1}^m h(B_{ij},B_{ij})(p).
\end{align*}
The last sum is $|\nabla du|^2(p)$ by the definition of the pointwise tensor norm in the orthonormal frame. Therefore
\begin{align*}
\Delta_g e(u)(p)=\sum_{j=1}^m h(C_j,A_j)(p)+|\nabla du|^2(p),
\end{align*}
where $C_j:=\sum_{i=1}^m\nabla^E_{e_i}B_{ij}(p)\in E_p$. The remaining task is to identify $C_j$ by commuting the covariant derivatives of the tensor $du$.
[/guided]
[/step]
[step:Commute the covariant derivatives of $du$]
At $p$, write
\begin{align*}
C_j
=
\sum_{i=1}^m\nabla^E_{e_i}B_{ij}(p).
\end{align*}
We first justify the symmetry of $\nabla du$. For smooth vector fields $X,Y\in\mathfrak X(M)$, the torsion-free property of $\nabla^M$ gives $[X,Y]=\nabla^M_XY-\nabla^M_YX$, while the torsion-free property of $\nabla^N$ gives
\begin{align*}
\nabla^E_X(du(Y))-\nabla^E_Y(du(X))=du([X,Y]).
\end{align*}
Substituting the definition of $\nabla du$ therefore gives
\begin{align*}
(\nabla du)(X,Y)=(\nabla du)(Y,X).
\end{align*}
Thus $B_{ij}=B_{ji}$, and hence
\begin{align*}
C_j
=
\sum_{i=1}^m\nabla^E_{e_i}B_{ji}(p).
\end{align*}
Now compare $\nabla^E_{e_i}B_{ji}$ with $\nabla^E_{e_j}B_{ii}$. For a section $\alpha\in\Gamma(T^*M\otimes E)$, the curvature of the [tensor-product connection](/page/Connection) on $T^*M\otimes E$ is
\begin{align*}
((R^{(T^*M)\otimes E}(X,Y)\alpha)(Z))=R^E(X,Y)(\alpha(Z))-\alpha(R^M(X,Y)Z).
\end{align*}
Indeed, the first term is the curvature action on the $E$ factor, and the second term is the dual curvature action on the covector factor. Applying this identity to $\alpha=du$, $X=e_i$, $Y=e_j$, and $Z=e_i$, we evaluate the left-hand side at $p$. The terms coming from differentiating the tensor slot vanish there because $\nabla^M_{e_i}e_k(p)=0$ for every $k$, and the Lie-bracket term in the curvature commutator vanishes because $[e_i,e_j](p)=0$. Hence
\begin{align*}
\nabla^E_{e_i}B_{ji}-\nabla^E_{e_j}B_{ii}=R^N(A_i,A_j)A_i-du(R^M(e_i,e_j)e_i).
\end{align*}
Therefore
\begin{align*}
C_j=\sum_{i=1}^m\nabla^E_{e_j}B_{ii}(p)+\sum_{i=1}^m R^N(A_i,A_j)A_i(p)-\sum_{i=1}^m du(R^M(e_i,e_j)e_i)(p).
\end{align*}
Since $\nabla^M_{e_j}e_i(p)=0$, the first sum is $\nabla^E_{e_j}\tau(u)(p)$. Therefore
\begin{align*}
C_j=\nabla^E_{e_j}\tau(u)(p)+\sum_{i=1}^m R^N(A_i,A_j)A_i(p)-du\left(\sum_{i=1}^m R^M(e_i,e_j)e_i\right)(p).
\end{align*}
Since $u$ is harmonic, $\tau(u)=0$, so $\nabla^E_{e_j}\tau(u)(p)=0$. With the stated curvature convention,
\begin{align*}
\sum_{i=1}^m R^M(e_i,e_j)e_i=-\operatorname{Ric}^M(e_j).
\end{align*}
Thus
\begin{align*}
C_j
=
du(\operatorname{Ric}^M(e_j))(p)
+
\sum_{i=1}^m R^N(A_i,A_j)A_i(p).
\end{align*}
[/step]
[step:Convert the curvature contraction to the stated sign]
Substituting the expression for $C_j$ into the second-derivative identity gives
\begin{align*}
\Delta_g e(u)(p)=|\nabla du|^2(p)+\sum_{j=1}^m h(du(\operatorname{Ric}^M(e_j)),A_j)(p)+\sum_{i,j=1}^m h(R^N(A_i,A_j)A_i,A_j)(p).
\end{align*}
Using metric compatibility of the [Levi-Civita connection](/page/Levi-Civita%20Connection) on $N$, the curvature endomorphism $R^N(X,Y):T_{u(p)}N\to T_{u(p)}N$ is skew-adjoint with respect to $h$, so
\begin{align*}
h(R^N(X,Y)Z,W)=-h(R^N(X,Y)W,Z).
\end{align*}
We apply this identity with $X=A_i$, $Y=A_j$, $Z=A_i$, and $W=A_j$. This gives
\begin{align*}
h(R^N(A_i,A_j)A_i,A_j)
=
-h(R^N(A_i,A_j)A_j,A_i).
\end{align*}
Since $A_i=du(e_i)$ and $A_j=du(e_j)$, this becomes
\begin{align*}
\Delta_g e(u)(p)=|\nabla du|^2(p)+\sum_{j=1}^m h(du(\operatorname{Ric}^M(e_j)),du(e_j))(p)-\sum_{i,j=1}^m h(R^N(du(e_i),du(e_j))du(e_j),du(e_i))(p).
\end{align*}
This is the desired identity at $p$.
[/step]
[step:Use tensoriality to obtain the global formula]
The point $p\in M$ was arbitrary. Each term in the identity is independent of the chosen normal frame: $\Delta_g e(u)$ is a scalar function, $|\nabla du|^2$ is the pointwise norm of a tensor, the Ricci term is the trace of the tensor pairing between $du\circ\operatorname{Ric}^M$ and $du$, and the target-curvature term is the full orthonormal-frame contraction of the curvature tensor of $N$ against $du$. Hence the identity proved at $p$ holds at every point of $M$ and in every local $g$-orthonormal frame. This proves the Bochner formula for harmonic maps.
[/step]