[proofplan]
We compute the first variation of the Dirichlet energy for an arbitrary smooth variation of $u$. The computation expresses the derivative of the energy as the negative $L^2$ pairing of the variation field with the tension field. If $\tau(u)=0$, this formula immediately gives stationarity. Conversely, if all first variations vanish, we use the exponential map of $N$ to realize the special variation field $V=\tau(u)$, forcing the integral of $|\tau(u)|_h^2$ to vanish and hence $\tau(u)=0$ pointwise.
[/proofplan]
[step:Compute the first variation of the Dirichlet energy]
Let $F:(-\varepsilon,\varepsilon)\times M\to N$ be a smooth variation of $u$, meaning $F(0,x)=u(x)$ for every $x\in M$. For each $t\in(-\varepsilon,\varepsilon)$, define the smooth map $u_t:M\to N$ by
\begin{align*}
u_t(x)=F(t,x).
\end{align*}
Define the variation field $V\in \Gamma(u^{-1}TN)$ by
\begin{align*}
V(x):=\left.\frac{\partial F}{\partial t}\right|_{(0,x)}\in T_{u(x)}N.
\end{align*}
We claim that
\begin{align*}
\left.\frac{d}{dt}\right|_{t=0} E(u_t)
=
-\int_M h_{u(x)}(\tau(u)(x),V(x))\,d\operatorname{vol}_g(x).
\end{align*}
Let $m=\dim M$. Since the identity is local on $M$, fix an [open set](/page/Open%20Set) $U\subset M$ and a smooth local $g$-orthonormal frame $e_1,\dots,e_m\in \Gamma(TU)$. On $U$, the energy density is
\begin{align*}
|du_t|_{g,h}^2=\sum_{i=1}^m h_{u_t(x)}(du_t(e_i),du_t(e_i)).
\end{align*}
Using the metric compatibility of the Levi-Civita connection of $h$, differentiating under the integral sign, and evaluating at $t=0$, we obtain
\begin{align*}
\left.\frac{d}{dt}\right|_{t=0}E(u_t)
&=
\int_M \sum_{i=1}^m h_{u(x)}\left(\nabla^N_{\partial_t}du_t(e_i)\big|_{t=0},du(e_i)\right)\,d\operatorname{vol}_g(x).
\end{align*}
The coordinate vector fields $\partial_t$ and $e_i$ commute on $(-\varepsilon,\varepsilon)\times U$, and the Levi-Civita connection of $h$ is torsion-free. Therefore
\begin{align*}
\nabla^N_{\partial_t}du_t(e_i)\big|_{t=0}
=
\nabla^{u^{-1}TN}_{e_i}V.
\end{align*}
Thus
\begin{align*}
\left.\frac{d}{dt}\right|_{t=0}E(u_t)
&=
\int_M \sum_{i=1}^m h_{u(x)}\left(\nabla^{u^{-1}TN}_{e_i}V,du(e_i)\right)\,d\operatorname{vol}_g(x).
\end{align*}
Define the smooth vector field $X\in \Gamma(TM)$ locally by
\begin{align*}
X:=\sum_{i=1}^m h_{u(x)}(V,du(e_i))\,e_i.
\end{align*}
At a chosen point of $U$, take the frame $e_1,\dots,e_m$ to be normal for $g$ at that point. Then $\nabla^M_{e_i}e_i=0$ at that point, and the divergence of $X$ at that point is
\begin{align*}
\operatorname{div}_g X
=
\sum_{i=1}^m e_i\left(h_{u(x)}(V,du(e_i))\right).
\end{align*}
By metric compatibility of the pullback connection,
\begin{align*}
\operatorname{div}_g X
=
\sum_{i=1}^m h_{u(x)}\left(\nabla^{u^{-1}TN}_{e_i}V,du(e_i)\right)
+
h_{u(x)}\left(V,\sum_{i=1}^m \nabla^{u^{-1}TN}_{e_i}du(e_i)\right)
\end{align*}
at that point. Since the frame is normal there, the second sum is exactly $h_{u(x)}(V,\tau(u))$ at that point. Because the point was arbitrary, the coordinate-invariant identity is
\begin{align*}
\sum_{i=1}^m h_{u(x)}\left(\nabla^{u^{-1}TN}_{e_i}V,du(e_i)\right)
=
\operatorname{div}_g X
-
h_{u(x)}(V,\tau(u)).
\end{align*}
Integrating this identity over $M$ with respect to $\operatorname{vol}_g$ gives
\begin{align*}
\left.\frac{d}{dt}\right|_{t=0}E(u_t)
&=
\int_M \operatorname{div}_g X\,d\operatorname{vol}_g(x)
-
\int_M h_{u(x)}(V,\tau(u))\,d\operatorname{vol}_g(x).
\end{align*}
Since $M$ is compact and has no boundary, the divergence term integrates to zero. Therefore
\begin{align*}
\left.\frac{d}{dt}\right|_{t=0}E(u_t)
=
-\int_M h_{u(x)}(\tau(u)(x),V(x))\,d\operatorname{vol}_g(x).
\end{align*}
[guided]
We begin with an arbitrary smooth variation $F:(-\varepsilon,\varepsilon)\times M\to N$ of $u$. For each time $t$, the slice $u_t:M\to N$ is the smooth map defined by
\begin{align*}
u_t(x)=F(t,x).
\end{align*}
The infinitesimal velocity at $t=0$ is the section
\begin{align*}
V(x):=\left.\frac{\partial F}{\partial t}\right|_{(0,x)}\in T_{u(x)}N.
\end{align*}
Thus $V\in \Gamma(u^{-1}TN)$.
The goal is to differentiate
\begin{align*}
E(u_t)=\frac{1}{2}\int_M |du_t|_{g,h}^2\,d\operatorname{vol}_g(x).
\end{align*}
The measure $\operatorname{vol}_g$ does not depend on $t$, so only the integrand is differentiated. Locally, choose a smooth $g$-orthonormal frame $e_1,\dots,e_m$ on an open set $U\subset M$, where $m=\dim M$. Then
\begin{align*}
|du_t|_{g,h}^2
=
\sum_{i=1}^m h_{u_t(x)}(du_t(e_i),du_t(e_i)).
\end{align*}
Differentiating this expression and using the metric compatibility of the Levi-Civita connection of $h$ gives
\begin{align*}
\left.\frac{d}{dt}\right|_{t=0}E(u_t)
&=
\int_M \sum_{i=1}^m h_{u(x)}\left(\nabla^N_{\partial_t}du_t(e_i)\big|_{t=0},du(e_i)\right)\,d\operatorname{vol}_g(x).
\end{align*}
The reason no factor of $2$ remains is that the energy includes the factor $\frac{1}{2}$.
Next we identify the differentiated term. The vector fields $\partial_t$ and $e_i$ commute on the product $(-\varepsilon,\varepsilon)\times U$, and the Levi-Civita connection of $h$ is torsion-free. Therefore differentiating first in $t$ and then in the spatial direction $e_i$ agrees with differentiating the variation field $V$ in the direction $e_i$:
\begin{align*}
\nabla^N_{\partial_t}du_t(e_i)\big|_{t=0}
=
\nabla^{u^{-1}TN}_{e_i}V.
\end{align*}
Substituting this into the energy derivative gives
\begin{align*}
\left.\frac{d}{dt}\right|_{t=0}E(u_t)
&=
\int_M \sum_{i=1}^m h_{u(x)}\left(\nabla^{u^{-1}TN}_{e_i}V,du(e_i)\right)\,d\operatorname{vol}_g(x).
\end{align*}
Now we integrate by parts on the compact manifold $M$. Define a local vector field $X\in \Gamma(TM)$ by
\begin{align*}
X:=\sum_{i=1}^m h_{u(x)}(V,du(e_i))\,e_i.
\end{align*}
At any chosen point of $U$, we may choose the orthonormal frame to be normal for $g$ at that point, so $\nabla^M_{e_i}e_i=0$ at that point. The divergence is then
\begin{align*}
\operatorname{div}_g X
=
\sum_{i=1}^m e_i\left(h_{u(x)}(V,du(e_i))\right).
\end{align*}
Using metric compatibility of the pullback connection,
\begin{align*}
\operatorname{div}_g X
=
\sum_{i=1}^m h_{u(x)}\left(\nabla^{u^{-1}TN}_{e_i}V,du(e_i)\right)
+
h_{u(x)}\left(V,\sum_{i=1}^m \nabla^{u^{-1}TN}_{e_i}du(e_i)\right).
\end{align*}
At the chosen point, the trace defining the tension field is
\begin{align*}
\tau(u)=\sum_{i=1}^m \nabla^{u^{-1}TN}_{e_i}du(e_i),
\end{align*}
because the frame is normal there. Hence
\begin{align*}
\sum_{i=1}^m h_{u(x)}\left(\nabla^{u^{-1}TN}_{e_i}V,du(e_i)\right)
=
\operatorname{div}_g X
-
h_{u(x)}(V,\tau(u)).
\end{align*}
This identity is tensorial, so it holds independently of the chosen normal frame.
Integrating over $M$ gives
\begin{align*}
\left.\frac{d}{dt}\right|_{t=0}E(u_t)
&=
\int_M \operatorname{div}_g X\,d\operatorname{vol}_g(x)
-
\int_M h_{u(x)}(V,\tau(u))\,d\operatorname{vol}_g(x).
\end{align*}
Because $M$ is compact and has no boundary, the integral of the divergence term is zero. Therefore the [first variation formula](/theorems/2728) is
\begin{align*}
\left.\frac{d}{dt}\right|_{t=0}E(u_t)
=
-\int_M h_{u(x)}(\tau(u)(x),V(x))\,d\operatorname{vol}_g(x).
\end{align*}
[/guided]
[/step]
[step:Deduce stationarity from the vanishing of the tension field]
Assume $\tau(u)=0$ as a section of $u^{-1}TN$. Let $F:(-\varepsilon,\varepsilon)\times M\to N$ be any smooth variation of $u$, and let $V\in \Gamma(u^{-1}TN)$ be its variation field. By the first variation formula,
\begin{align*}
\left.\frac{d}{dt}\right|_{t=0}E(u_t)
=
-\int_M h_{u(x)}(\tau(u)(x),V(x))\,d\operatorname{vol}_g(x)
=
0.
\end{align*}
Thus $u$ is a critical point of $E$ under every smooth variation.
[/step]
[step:Realize the tension field itself as a variation field]
Assume conversely that $u$ is a critical point of $E$ under every smooth variation. We show that the section $\tau(u)\in \Gamma(u^{-1}TN)$ is an admissible variation field.
Because $M$ is compact and $u$ is smooth, the image $u(M)\subset N$ is compact. The section $\tau(u):M\to u^{-1}TN$ is smooth, hence its image is compact in $TN$. The exponential map of $(N,h)$ is defined on an open neighbourhood of the zero section in $TN$, so there exists $\delta>0$ such that
\begin{align*}
t\,\tau(u)(x)\in \operatorname{Dom}(\exp)
\end{align*}
for every $x\in M$ and every $t\in(-\delta,\delta)$.
Define the smooth map $F:(-\delta,\delta)\times M\to N$ by
\begin{align*}
F(t,x)=\exp_{u(x)}(t\,\tau(u)(x)).
\end{align*}
This is a smooth variation of $u$, since $F(0,x)=\exp_{u(x)}(0)=u(x)$. Its variation field is
\begin{align*}
\left.\frac{\partial F}{\partial t}\right|_{(0,x)}
=
\tau(u)(x),
\end{align*}
because the differential of $\exp_{u(x)}$ at $0\in T_{u(x)}N$ is the identity map on $T_{u(x)}N$.
[/step]
[step:Test stationarity against the tension field and conclude pointwise vanishing]
Apply the first variation formula to the variation constructed above. Since $u$ is critical, the left-hand side vanishes, and since the variation field is $V=\tau(u)$, we obtain
\begin{align*}
0=\left.\frac{d}{dt}\right|_{t=0}E(u_t).
\end{align*}
The first variation formula with $V=\tau(u)$ gives
\begin{align*}
\left.\frac{d}{dt}\right|_{t=0}E(u_t)=-\int_M h_{u(x)}(\tau(u)(x),\tau(u)(x))\,d\operatorname{vol}_g(x).
\end{align*}
Since $h_{u(x)}(\tau(u)(x),\tau(u)(x))=|\tau(u)(x)|_h^2$, this is
\begin{align*}
0=-\int_M |\tau(u)(x)|_h^2\,d\operatorname{vol}_g(x).
\end{align*}
Therefore
\begin{align*}
\int_M |\tau(u)(x)|_h^2\,d\operatorname{vol}_g(x)=0.
\end{align*}
The function
\begin{align*}
x\mapsto |\tau(u)(x)|_h^2
\end{align*}
is continuous and non-negative on $M$. If it were positive at some point $x_0\in M$, continuity would give an open neighbourhood $U\subset M$ of $x_0$ on which it is positive, and the Riemannian volume measure satisfies $\operatorname{vol}_g(U)>0$. This would make the integral strictly positive, a contradiction. Hence $|\tau(u)(x)|_h^2=0$ for every $x\in M$, and therefore
\begin{align*}
\tau(u)=0.
\end{align*}
This proves the converse implication and completes the proof.
[/step]