Tweetable way to see Riemannian isometries are harmonic? $\newcommand{\al}{\alpha}$
$\newcommand{\euc}{\mathcal{e}}$
$\newcommand{\Cof}{\operatorname{Cof}}$
$\newcommand{\Det}{\operatorname{Det}}$
Smooth Riemannian isometries are harmonic. Can one conclude this  just from "staring" at the Dirichlet energy functional? (Ideally without even computing the Euler-Lagrange equations).
(As is the case of constant mappings which are global minimizers)
I am interested in finding a very simple explanation, which involves a minimal amount of computation. 
To be honest, I suspect such a simple argument cannot exist (I describe below why), but maybe I am overly pessimistic.
The difficulties:
(1) It's false for isometric embeddings, i.e the result only holds for mappings between equidimensional manifolds. So any argument must be involved enough to address this distinction.
(2) It does not hold in a weak sense. ("Weak isometries are not weakly harmonic"). If we assume $f$ is only weakly differentiable and $df \in O(n)$ almost everywhere, then $f$ need not be weakly harmonic.
Indeed, Gromov constructed arcwise isometries between equidimensional non-flat manifolds and flat ones. Such maps are $1$-Lipschitz, hence weakly differentiable and classically differentiable a.e, and their differentials are isometries a.e. However, they cannot be weakly harmonic, since continuous weakly harmonic maps are smooth, which contradicts the mismatch in curvatures. 
Gromov's pathological maps cannot be orientation-preserving or orientation-reversing on any open subset of the domain. My colleagues and I showed that any weakly differentiable map $f \in W^{1,\infty}(M,N)$ satisfying $df \in SO(n)$ a.e is weakly harmonic. 
The two points discussed above hint that it's unresonable to expect a "pointwise identity" that will show "isometry $\Rightarrow$ harmonicity".

I know essentally of 3 different ways to prove isometries are harmonic (all of them requires too much computation, in my opinion):
(I) Prove every totally geodesic map is harmonic ($f:M \to N$ is totally geodesic if it maps geodesics to geodesics). 
The equidimensionality  enters the game when showing our isometry  preserve geodesics. However, the claim "geodesic $\Rightarrow$ harmonic" requires a computation in order to pass from "preservation of geodesics" to $\nabla df =0$ which implies harmonicity since $\tau(f)=\operatorname{trace}_g(\nabla df)$).
(II) Prove every minimal isometric immersion is harmonic. Again, this requires a computation (and the equidimensionality also plays its role).
(III) Showing that for every smooth map $f:M \to N$, $\delta (\Cof df)=0$ where
$\Cof df:TM \to f^*(TN)$ is the cofactor map of $df$ defined by $$ \Cof df= (-1)^{d-1} \star_{f^*TN}^{d-1} (\wedge^{d-1} df) \star_{TM}^1. $$
($\Cof df$ represents the action of $df$ on $d-1$ dimensional parallelepipeds).
My colleagues and I showed that $\delta (\Cof df)=0$ for arbitrary (smooth) maps $M \to N$ (for details see Proposition 3.4 here). Again, this requires computations (If someone has a more heuristic argument, I would be glad to know).
Since $df$ is an isometry, $df=\Det df \cdot \Cof df$ (See here for heuristics). The smoothnes of $f$ implies $\Det df$ is constant, so 
$$ \delta ( df)=\Det df \cdot \delta (\Cof df)=0 $$
This point of view enlightens the importance of not switching orientation in order to be harmonic.

"Comment:"
In the case where the target space is Euclidean there is a relatively simple proof:
Suppose $f:(M,g) \to (\mathbb{R}^d,e)$ is a smooth isometry.
Harmonicity of $f$ amounts to the component-wise harmonicity $\Delta f^\al= 0$. 
This follows from the pullback property,
$$
\star_{f^*\euc} df^\al = \star_{f^*\euc}(f^*\,dx^\al) = \pm  f^*(\star_\euc dx^\al), \tag{1}
$$
Where the sign is determined according to whether $f$ preserves or reverses orientation. Since $f$ is smooth this sign is constant (if $M$ is connected).
This implies that
$$
d\star_{f^*\euc} df^\al = d\big(\pm  f^*(\star_\euc dx^\al)\big) \stackrel{(a)}{=} \pm d\big(  f^*(\star_\euc dx^\al)\big) = \pm f^*(d\star_\euc dx^\al) \stackrel{(b)}{=} 0, \tag{2}
$$
Where in equality $(a)$ we used the fact the sign is constant, and equality $(b)$ follows since $$d\star_\euc dx^\al=d \big((-1)^{\al+1}dx^1 \wedge \dots \wedge \widehat{dx^\al} \wedge \dots \wedge dx^d\big)=0.$$
Since $f^*e=g$, equation $(2)$ becomes $ d\star_{g} df^\al = 0 $ which implies $ \delta (df^\al)=-\star_{g} d\star_{g} df^\al = 0$, so $f^\al$ is harmonic.
Perhaps there is a way to generalize this proof to the general case, using some adapted coordinates.
 A: Let $u: M\rightarrow N$ be a diffeomorphism. By staring at the Dirichlet energy formula and knowing that integration by parts works just as well here as for the classical case for functions, you know that a map is harmonic if and only if $\Delta u = 0$.
However, if $u$ is isometric, then you can see in your head that
$$
0 = \nabla_k(\partial_iu\cdot\partial_ju) = \nabla^2u_{ik}\cdot\partial_ju  + \partial_iu\cdot\nabla^2_{kj}u.
$$
But you also know, given the symmetries of the indices that the standard argument proving the uniqueness of the Levi-Civita connection (also known as Cartan's lemma) implies that
$$
\partial_ku\cdot\nabla^2_{ij}u = 0.
$$
Specifically,
\begin{align*}
  \partial_ku\cdot\nabla^2_{ij}u
& = \frac{1}{2}(\nabla_i(\partial_ku\cdot\partial_ju)
  + \nabla_j(\partial_ku\cdot\partial_iu)
  - \nabla_k(\partial_iu\cdot\partial_ju))\\
&  = \frac{1}{2}(\nabla_ig_{kj} + \nabla_jg_{ki} - \nabla_kg_{ij})\\
&= 0.
\end{align*}
Since $u$ is a diffeomophism, $\partial_1 u, \cdots, \partial_n u$ span $T_*N$ and therefore, $\nabla^2u = 0$ and $u$ is harmonic.
If $\dim M < \dim N$, then the argument above shows that $\nabla^2u$ is normal to $T_*M$. In fact, it is the second fundamental form (viewed as a normal-vector-valued symmetric tensor).
A: Not exactly 'tweetable', but perhaps the identity (1) may help, if all you want to do is avoid the Euler-Lagrange equations.  For simplicity, assume that $M^n$ is oriented.  (One can write the identity (1) below as an identity on densities, so the orientability hypothesis is not essential, but I'll leave that detail for the interested.)  
If $g$ is a metric on $M$, let $\mathrm{vol}_g>0$ denote its associated volume form (positively oriented).  Then for any vector field $X$ on $M$, one has the $n$-form identity
$$
\tfrac12\,\mathrm{tr}_g\bigl(\mathcal{L}_Xg\bigr)\,\mathrm{vol}_g 
= \mathrm{d}\bigl(\iota(X)\,\mathrm{vol}_g\bigr),\tag1
$$
where $\mathcal{L}_X$ denotes Lie derivative with respect to $X$ and $\iota(X)\,\mathrm{vol}_g$ denotes the $(n{-}1)$-form that is the interior product of $X$ with $\mathrm{vol}_g$.
Since, for any vector field $X$ on $M$ (say, with compact support), one has, by definition
$$
\left.\frac{d\ }{dt}\right|_{t=0}\bigl(\mathrm{exp}_{tX}^*g\bigr) 
= \mathcal{L}_Xg\,,
$$
one has, for the $1$-parameter family of diffeomorphisms $\mathrm{exp}_{tX}:M\to M$,
$$
\left.\frac{d\ }{dt}\right|_{t=0}\bigl(e(\mathrm{exp}_{tX})\bigr)
=\left.\frac{d\ }{dt}\right|_{t=0}\bigl(\tfrac12\mathrm{tr}_g(\mathrm{exp}_{tX}^*g)\,\mathrm{vol}_g\bigr) 
= \tfrac12\mathrm{tr}_g(\mathcal{L}_Xg)\,\mathrm{vol}_g 
= \mathrm{d}\bigl(\iota(X)\,\mathrm{vol}_g\bigr).
$$
Integrating this over $M$,
$$
\left.\frac{d\ }{dt}\right|_{t=0}E(\mathrm{exp}_{tX})
=  \int_M \left(\left.\frac{d\ }{dt}\right|_{t=0}e(\mathrm{exp}_{tX})\right)
= \int_M \mathrm{d}\bigl(\iota(X)\,\mathrm{vol}_g\bigr) = 0,
$$
by Stokes' Theorem.  Thus, the identity $(= \mathrm{exp}_{0X})$
is $E$-critical for all compactly supported variations, as desired.
