This is a cross-post from MSE (no answer there).
Let $M,N$ be oriented $d$-dimensional Riemannian manifolds, $M$ compact*, and let $f:M \to N$ be smooth.
Consider the Dirichlet energy functional: $E_{M,N}(f)=\int_M \|df\|^2 \operatorname{Vol}_g$. ($\operatorname{Vol}_g$ is the Riemannian volume form of $M$).
Fix another $d$-dimensional Riemannian manifold $W$. For every isometric immersion $\phi:N \to W$,
$$ E_{M,N}(f)=E_{M,W}(\phi \circ f) \, \text{ for any } \, f:M \to N. \tag{1}$$
Let us call a smooth map $\phi$ which sarisfies $(1)$ a symmetry of the Dirichlet's integral.
Question: Is every symmetry an isometric immersion?
Does anything change if we restrict the symmetries to be invertible maps? (or diffeomorphisms)?
Has this notion of symmetry been studied in the context of mappings between manifolds? (In more "classical Euclidean" settings I found similar notions by the name of "variational symmetry groups").
*We can relax the compactness assumption, but then we need to restrict the discussion to maps which are constant outside a compact domain.
Comments:
$(1)$ The "pointwise analog" is trivial:
If $B \in M_d$, and for every $A \in M_d$, $ \|BA\|=\|A\|$ (where $\| \|$ is the Euclidean norm), then $B$ is an orthogonal matrix.
The challenge is that the notion of a Dirichlet's-symmetry is a global one, while the concept of isometric immersion is local.
I think the rough idea should be to choose "test maps" $f$ which are very "localized" (are constant outside small balls). However, this does not seem trivial, since the differential would have to pass from a given linear map to zero, so its norm would vary. (Think of a bump function which goes quickly from one to zero, you can't make the integral of the derivative small).
$(2)$ The answer could a-priori depend on the manifolds. Even the case where $M$ is an Euclidean ball in $\mathbb{R}^d$, $N=W=\mathbb{R}^d$, the answer does not seem to be trivial (see the previous comment).