The concept of harmonicity of a map between Riemannian manifolds depends only on the Levi-Civita connection of the target manifold, as can be seen from the Euler-Lagrange equations (details below).
This is not trivial from the definition*, since the Dirichlet integral is explicitly dependent upon the metric of the taregt, not just its induced connection.
Is there a way to "see in advance", without deriving the Euler-Lagrange equations, that the harmonicity will "turn out" to depend only upon the connection?
Another way to view this phenomena, is to think of it as an additional symmetry which is "created" when passing from the Dirichlrt energy functional to its associated E-L equations; While the original functional was not invariant under connection-preserving transformations, its set of critical points is invariant under such transformations.
Is there a way then to see "in advance" the set of solutions to the E-L equations will have "more symmetries" than the functional? (inclusion in one way is obvious)
Suggested approach:
Suppose we somehow knew that the Euler-Lagrange equations, evaluated at a map $\phi:M \to N$ are of the form $\alpha_{M,N}(\phi)=0$, where in general $\alpha_{M,N}(\psi)$ takes values in $\Gamma(\psi^*TN)$. (Perhaps there is a way to deduce that somehow but let's leave it for now).
This raises the following question:
Is this assumption, together with the known symmetry of isometry-invariance forces our operator to be the "right" one?
More precisely, consider pairs $(M,g),(N,h)$ of smooth Riemannian manifolds.
Are there any non-zero operators $\alpha_{M,N}$ which takes $\phi \in C^{\infty}(M,N)$ to $\Gamma(\phi^*TN)$ and are invariant under "scaled isometries", other than $\alpha_{M,N}(\phi)=\operatorname{trace}_g(\nabla d\phi)=\delta(d\phi)$ (up to scaling)?
($\nabla$ is the tensor product connection of $\nabla^{T^*M},\phi^*(\nabla^{TN})$, where $\nabla^{T^*M}$ is the dual connection to the L-C connection on $TM$, $\nabla^{TN}$ is the L-C connection on $TN$).
By "invariant", I mean that if $\phi:N \to \tilde N$ is a scaled isometry then
$$ d \tilde \phi \big( \alpha_{M,N} (\phi) \big)=\alpha_{M,\tilde N} (\tilde \phi \circ \phi)$$
where I consider $ d \tilde \phi$ to as a map $\phi^*TN \to (\tilde \phi \circ \phi)^*T\tilde N$,
and I require similar invariance under compositions with scaled isometries from the left.
Perhaps I should also require that $\alpha_{M,N}$ will be dependent only upon $d\phi$ and its derivative.
Comment: Requiring only invariance under isometries (not including all scaled isometries), is not restrictive enough, e.g take
$$ \beta_{M,N} (\phi)=\| d\phi\|\operatorname{trace}_g(\nabla d\phi)=\| d\phi\|\delta(d\phi)$$
Then $\beta$ commutes with isometries,
$$ \beta_{M,\tilde N} (\tilde \phi \circ \phi)=\| d \tilde \phi \circ d\phi\|\delta(d \tilde \phi \circ d\phi)=\| d\phi\|d \tilde \phi\big(\delta(d\phi)\big)$$
$$ = d \tilde \phi\big(\| d\phi\|\delta(d\phi)\big)=d \tilde \phi \big( \beta_{M,N} (\phi) \big)$$
but not with scaled isometries: If $\tilde \phi$ is a scaled isometry with conformal factor $\lambda \in \mathbb{R}$, then
$$ \beta_{M,\tilde N} (\tilde \phi \circ \phi)=\| d \tilde \phi \circ d\phi\|\delta(d \tilde \phi \circ d\phi)=\lambda\| d\phi\|d \tilde \phi\big(\delta(d\phi)\big)$$
$$ = \lambda d \tilde \phi\big(\| d\phi\|\delta(d\phi)\big)=\lambda d \tilde \phi \big( \beta_{M,N} (\phi) \big)$$
If the answer to the above question is positive, then this gives a nice explanation for the additional symmetry "created". It is the only possible outcome...
*I am using the definition of harmonic maps as critical points of the Dirichlet functional.
Proposition:
Let $h,\tilde h$ be Riemannian metrics on a manifold $N$, with the same Levi-Civita connection.
Then $f:(M,g) \to (N,h)$ is harmonic $\iff$ $f:(M,g) \to (N,\tilde h)$ is harmonic.
Proof:
Let $\nabla=\nabla_h=\nabla_{\tilde{h}}$ be the L-C connection of $h,\tilde h$ on $TN$. Denote by $\nabla^{f^*TN}$ the pullback connection by $f$. The corresponding well-known E-L equations imply
$f:(M,g) \to (N,h)$ is harmonic $\iff \delta_h(df)=0$, where
$\delta_h:\Omega^1(M,f^*TN) \to \Omega^0(M,f^*TN)=\Gamma(f^*TN)$, is the adjoint of $\nabla^{f^*TN} $, taken w.r.t metrics $g,f^*h$ on $TM,f^*TN$.
Similarly,
$f:(M,g) \to (N,\tilde h)$ is harmonic $\iff \delta_{\tilde h}(df)=0$
Since both $f^*h,f^*{\tilde h}$ are metric w.r.t $\nabla^{f^*TN}$, and since the codifferential of $d_{\nabla^{f^*TN}}$ is the same for every $\nabla^{f^*TN}$-compatible metric on $f^*TN$, we get that $\delta_{\tilde h}=\delta_h$.
