$\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.