$\newcommand{\Cof}{\text{Cof}}$ $\newcommand{\R}{\mathbb{R}}$ Let $M,N$ be smooth oriented manifolds, M compact, $\partial N=\emptyset$. Let $\omega \in \Omega^m(N)$ be closed. ($m=\dim M$). I don't assume that $\dim M=\dim N$.
Consider the functional $E:C^{\infty}(M,N) \to \mathbb{R}$, defined by $E(f)=\int_M f^*\omega$. Can we obtain an "explicit" form of its Euler-Lagrange equation?
I imagine something of the form $\Phi_{\omega}(f)=0$, where $\Phi_{\omega}$ is some second order differential operator (which depends on $\omega$).
One idea I have is to introduce more structure, e.g. endow $M$ and $N$ with Riemannian metrics $g_M,g_N$. Now we might be able to find a "metric-dependent" formulation of the Euler-Lagrange equation. (Even though our functional does not depend on the Riemannian structures, of course).
More specifically, let $V \in \Gamma(f^*TN)$ be compactly supported in the interior $M^o$, and let $f_t$ be a variation of $f$ satisfying $\left. \frac{\partial}{\partial t} \right|_{t=0} f_t=V$. My goal is to express
$$ \left. \frac{d}{dt} \right|_{t=0} E(f_t)=\int_M \langle V,\Phi_{\omega,g_M,g_N}(f)\rangle Vol_M,$$ where $\Phi_{\omega,g_M,g_N}:C^{\infty}(M,N) \to \Gamma(f^*TN)$ is some (second order) differential operator, which depends on $\omega$ and on the metrics $g_M,g_N$. Actually, since only the metric on $N$ appears in $\langle V,\Phi_{\omega,g_M,g_N}(f)\rangle$, $\Phi_{\omega,g_M,g_N}$ depends "in essence" only on $g_N$; however, $g_M$ may "get in" in a way which "cancels itself". (See example below).
Perhaps something non-trivial can be said at least when $\omega$ is parallel (w.r.t the Levi-Civita connection) or harmonic?
Comment: I know that $E$ is a null-Lagrangian, i.e. that every map satisfies its E-L equation. (This follows from the fact that I assumed that $d\omega=0$ and that the variation field is compactly supported in the interior; i.e. the variation do not change the values of the map on boundary $\partial M$). However, it does not mean that the E-L equation is "trivial" i.e. $0=0$. See the example below.
The first task is to find an expression for $\left. \frac{d}{dt} \right|_{t=0} f_t^*\omega$ where $f_t$ is some variation of $f$.
Here are two attempts:
First attempt:
Set $f_t=\phi_t \circ f$, where $\phi_t \in \text{Diff}(N)$. If $\left. \frac{\partial \phi_t}{\partial t} \right|_{t=0}=X$, then $V=X \circ f$.
Then $f_t^* \omega=f^* \phi_t^* \omega$, so $$\left. \frac{d}{dt} \right|_{t=0} f_t^*\omega=f^*L_X\omega,$$ and
$$ \left. \frac{d}{dt} \right|_{t=0} E(f_t) =\int_M f^*L_X\omega=\int_M f^*d(\iota_X \omega).$$
Note that since $V=X \circ f$, then $f(p_1)=f(p_2) \Rightarrow V(p_1)=V(p_2)$, so if $f$ is not injective, not every variation field $V \in \Gamma(f^*TN)$ can be expressed as $V=X \circ f$, for some $X \in \Gamma(TN)$.
I don't know how to pass from here to something of the form $\langle X,\Phi(f) \rangle $.
Second attempt:
We now take $f_t=f \circ \psi_t$, where $\psi_t \in \text{Diff}(M)$. If $\left. \frac{\partial \psi_t}{\partial t} \right|_{t=0} =X \in \Gamma(TM)$, then $V=df \circ X$. (This means that if $df_p$ is not surjective not every $V$ can be realized in this way).
Then $f_t^* \omega= \psi_t^* f^*\omega$, so $$\left. \frac{d}{dt} \right|_{t=0} f_t^*\omega=L_Xf^*\omega,$$ and
$$ \left. \frac{d}{dt} \right|_{t=0} E(f_t) =\int_M L_Xf^*\omega=\int_M d(\iota_X f^*\omega).$$
Again, I don't know how to pass from here to something of the form $\langle X,\Phi(f) \rangle $.
Here is an example for the kind of 'explicit' expression I am looking for:
Let $M=\Omega $ be an open subset of $\mathbb{R}^n$, and $N=\mathbb{R}^n$, both endowed with the Euclidean metrics, and set $\omega= \text{Vol}_{\mathbb{R}^n}$. Then $$ E(f)=\int_{\Omega} f^*\text{Vol}_{\R^n}=\int_{\Omega} \det df \text{Vol}_{\R^n}, $$ and so writing $f_t=f+tV$ for some $V:\Omega \to \R^n$, we obtain
$$\left. \frac{d}{dt} \right|_{t=0} E(f_t)=\int_{\Omega} \left. \frac{d}{dt} \right|_{t=0} \det df_t=\int_{\Omega} \left. \frac{d}{dt} \right|_{t=0} \det (df+tdV)=\int_{\Omega} \langle \Cof df,dV \rangle=\int_{\Omega} \langle \text{div} \Cof df,V \rangle,$$
so the Euler-Lagrange equation is given by $$ \text{div} \Cof df=0,$$ where $\Cof df$ is the cofactor matrix of $df$, and the divergence here is taken row-by-row.
This special case can be generalized to the case where $M,N$ are equidimensional Riemannian manifolds, and $\omega=\text{Vol}_N$ is the Riemannian volume form of $N$. (admittedly, here $\omega$ itself depends upon the target metric $g_N$, which is not exactly the context I presented in this question so far, where the metrics $g_M,g_N$ and the form $\omega$ are chosen independently).
The corresponding Euler-Lagrange equation is $\delta (\Cof df)=0$, where $\delta$ and $\Cof df$ are natural generalizations of the cofactor matrix and the divergence operator to Riemannian settings. For more details about this see this answer of mine.
This in fact gives an infinite family of equations, that any map satisfies-one equation for every two metrics we choose.
The first comment is that you say that you want to get an expression in terms of the metrics, but in fact your functional doesn't depend on the Riemannian structure at all.
Secondly, Lie derivative commutes with pullbacks, so assuming that your vector field is compactly supported in the interior of $M$ you get that $\delta E=0$. This is just a manifestation of the previously mentioned fact that your functional is diffeomorphism invariant. Thus "Euler-Lagrange" equations take the form $0=0$.