$\newcommand{\id}{\operatorname{Id}}$
$\newcommand{\TM}{\operatorname{TM}}$
$\newcommand{\Hom}{\operatorname{Hom}}$
$\newcommand{\M}{\mathcal{M}}$
$\newcommand{\N}{\mathcal{N}}$
$\newcommand{\tr}{\operatorname{tr}}$
$\newcommand{\TM}{\operatorname{T\M}}$
$\newcommand{\TN}{\operatorname{T\N}}$
$\newcommand{\TstarM}{T^*\M}$
Edit: I narrowed the focus of the question.
Summary: I suspect $\delta(df \wedge df)=0$ is not an E-L equation because it pus too many constraints on $f$. Can this heuristic be formalized?
Let $\M,\N$ be $d$-dimensional oriented Riemannian manifolds, and let $f:\M \to \N$ be smooth; $df \in \Omega^1(\M,f^*{\TN})$, and for $1 < k \le d$ let $\bigwedge^k df\in \Omega^k\big(\M,\Lambda_k(f^*{\TN})\big)$ be the induced map.
Notation:
$\nabla^{\Lambda_k(f^*{\TN})}$ is the induced connection on $\Lambda_k(f^*{\TN})$, $\delta_{\nabla^{\Lambda_k(f^*{\TN})}} $ is the adjoint of $d_{\nabla^{\Lambda_k(f^*{\TN})}}$.
Question: Is $\delta_{\nabla^{\Lambda_k(f^*{\TN})}} \big( \bigwedge^k df \big) =0$ an Euler-Lagrange equation of some functional $E(f)$?
I suspect the answer is negative.
Edit: Can this be proved using known results on inverse problems in calculus of variations? e.g if we can show the equation in the Euclidean setting ($M=N=\mathbb{R}^d$) is not an E-L eq, then we are done.
(I am not aware of much work on such inverse problems in general Riemannian settings, but if the Euclidean can be decided, it's enough).
Heuristic:
$\delta_{\nabla^{\Lambda_k(f^*{\TN})}} \big( \bigwedge^k df \big) =0$ puts too many constraints on $f$:
$\delta_{\nabla^{\Lambda_k(f^*{\TN})}} \big( \bigwedge^k df \big) \in \Omega^{k-1}(\M,\Lambda_k(f^*{\TN}))$, so locally the "equation" is a system of ${d \choose k-1}{d \choose k}$ scalar equations.
However, it seems to me an Euler-Lagrange eq of a functional cannot consist of more than $\dim(f^*\TN)=d$ equations, since (roughly) this is the number of degrees of freedom we have in choosing the variation field $V\in\Gamma(f^*\TN)$.
More explicitly, since $dE(V)$ is linear in $V$, it should always be in the form of $dE(V)=\langle A(f),V \rangle$, where $A(f) \in \Gamma(f^*\TN)$ so the E-L eq should be in the form of $A(f)=0$.
Of course, this argument gives only an upper bound to the number of "independent" constraints – degeneracies can occur (e.g Null-Lagrangians).
So, suppose ${d \choose k-1}{d \choose k}>d$. Is it possible $\delta_{\nabla^{\Lambda_k(f^*{\TN})}} \big( \bigwedge^k df \big) =0$ is degenerate and reduces to no more than $d$ independent eqs (locally)?
Can we prove this is not happening? (If we can, than this shows our equation is not an E-L eq of any functional).
Analysis of the borderline cases:
We now turn to see what happens when ${d \choose k-1}{d \choose k}=d$. This happens (for $1 \le k \le d$) iff $k=1$ or $k=d$. The former corresponds to harmonicity. Suppose that $k=d$. We will try to understand if our equation is an E-L equation in this case.
For concreteness, let's work with $k=d=2$.
I proved in my answer that the E-L equation of the functional
$$ E_2(f)=\frac{1}{2}\int_{M} \| \bigwedge^2 df\|^2 \text{Vol}_{M}, \, \text{is}$$
$$h_{f^*\TN} \big( \tr_{\TM}\big( df \otimes \delta_{\nabla^{\Lambda_2(f^*{\TN})}}(\bigwedge^2 df)\big) \bigg)=0.$$
Since $h_{f^*\TN}$ is injective in this case, the E-L eq is equivalent to
$$\tr_{\TM}\big( df \otimes \delta_{\nabla^{\Lambda_2(f^*{\TN})}}(\bigwedge^2 df)\big)=0. \tag{1}$$
When restricting the discussion to immersions, this reduces to $ \delta_{\nabla^{\Lambda_2(f^*{\TN})}}(\bigwedge^2 df)=0$.
It is not clear what happens in the general case; If we don't assume $f$ is an immersion, are the equations equivalent?
(If not, perhaps there is another way to realize $\delta_{\nabla^{\Lambda_2(f^*{\TN})}}(\bigwedge^2 df)=0$ as an E-L eq).
The problem is that we cannot conclude (at least not immediately) that at
every point where $df$ is not invertible, $\delta_{\nabla^{\Lambda_2(f^*{\TN})}}(\bigwedge^2 df)=0$. (We know $\bigwedge^2 df=0$ at the point, but $\delta$ is a differential operator, it sees beyond the pointwise behaviour).
Best Answer
Consider the following functional:
$$ E_k(f)=\frac{1}{2}\int_{M} \| \bigwedge^k df\|^2 \text{Vol}_{M}.$$
Proof:
Let $\phi$ be a map $\M \to \N$, and $\phi_t:\M \to \N$ as smooth family, where $\phi_0=\phi$ and $\frac{\partial \phi_t}{\partial t}|_{t=0}:=V \in \Gamma(\phi^*(\TN))$. Then
$$ \frac{d}{dt}|_{t=0}E(\phi_t)=\frac{1}{2}\int_{\M}\frac{\partial{}}{dt}|_{t=0} \| d\phi_t \wedge d\phi_t \|^2 \text{Vol}_{\M}= \int_{\M} \langle d\phi \wedge d\phi, \nabla_{\frac{\partial{}}{dt}} (d\phi_t \wedge d\phi_t)|_{t=0}\rangle \text{Vol}_{\M}. \tag{1}$$
It is well-known that $\nabla_{\frac{\partial{}}{dt}} d\phi_t|_{t=0}=\nabla^{\phi^*(TN)}V \in \Gamma(T^*\M \otimes \phi^*(\TN))$.
Now,
$$\bigg(\nabla_{\frac{\partial{}}{dt}} (d\phi_t \wedge d\phi_t)|_{t=0}\bigg)(X,Y)=$$ $$(\nabla_{\frac{\partial{}}{dt}} d\phi_t|_{t=0})(X) \wedge d\phi(Y)+d\phi(X) \wedge (\nabla_{\frac{\partial{}}{dt}} d\phi_t|_{t=0})(Y)=$$
$$ \nabla V(X)\wedge d\phi(Y)+d\phi(X) \wedge \nabla V(Y)= $$
$$\big(\nabla V \wedge d\phi+d\phi \wedge \nabla V\big) (X,Y), \tag{2}$$
where $\nabla V \wedge d\phi+d\phi \wedge \nabla V \in \Omega^2\Big(\M,\Lambda_2 \big(\phi^*T\N\big)\Big)$ is defined by the last equality.
Thus, we have obtained
$$ \nabla_{\frac{\partial{}}{dt}} (d\phi_t \wedge d\phi_t)|_{t=0}= \nabla V \wedge d\phi+d\phi \wedge \nabla V. \tag{3}$$
Define also $\xi=V \wedge d\phi \in \Omega^1\Big(\M,\Lambda_2 \big(\phi^*T\N\big)\Big)$.
Lemma: $d_{{\nabla}^{\Lambda_2(\phi^*T\N)}}(\xi)=\nabla V \wedge d\phi+d\phi \wedge \nabla V$.
Assuming the lemma, we combine equations $(1),(3)$ and get
$$ \frac{d}{dt}|_{t=0}E(\phi_t)=\int_{\M} \langle d\phi \wedge d\phi, d_{{\nabla}^{\Lambda_2(\phi^*T\N)}}(\xi)\rangle \text{Vol}_{\M}=\int_{\M} \langle \delta_{\nabla^{\Lambda_2(\phi^*{\TN})}}(d\phi \wedge d\phi),\xi\rangle \text{Vol}_{\M}. \tag{5}$$
To find the exact $E-L$ equations, one further step needs to be taken: $$V \to \langle \delta_{\nabla^{\Lambda_2(\phi^*{\TN})}}(d\phi \wedge d\phi),\xi\rangle=\langle \delta_{\nabla^{\Lambda_2(\phi^*{\TN})}}(d\phi \wedge d\phi),V \wedge d\phi\rangle$$ is a linear functional in $V$, so it can be expressed as $V \to \langle V, A(\phi) \rangle_{\phi^*\TN}$,
where $A(\phi) \in \Gamma(\phi^*\TN)$. The E-L equation is $A(\phi)=0$.
We now turn to finding an explicit expression for this "representation":
The corresponding pointwise linear algebra situation is this:
We have two oriented $d$-dimensional inner product spaces $V,W$, together with maps $A \in \Hom(V,W),B \in \Hom(V,\Lambda_2(W))$, and we look for a bilinear map $$\psi: \Hom(V,W) \times \Hom(V,\Lambda_2(W)) \to W,$$ satisfying
$$ \langle w \wedge A,B \rangle_{\Hom(V,\Lambda_2(W))}=\langle w, \psi(A,B) \rangle_W \, \text{ for every $w \in W$}$$
Proposition: With the notation as above, $\psi(A,B)=h_W\big(\tr_{V} (A \otimes B)\big)$ where $h_W:W \otimes \Lambda_2(W) \to W$ is defined by the linear extension of
$$ \tilde w \otimes (w_1 \wedge w_2) \to \langle \tilde w,w_2 \rangle w_1-\langle \tilde w,w_1 \rangle w_2. \tag{6}$$
Note $A \otimes B \in V^* \otimes V^* \otimes W \otimes \Lambda_2(W)$, so $\tr_{V} (A \otimes B) \in W \otimes \Lambda_2(W)$.
Proof:
It suffices to prove this for $A,B$ "pure" tensors, i.e $A=\alpha \otimes \tilde w,B=\beta \otimes (w_1 \wedge w_2)$, where $\alpha,\beta \in V^*,\tilde w,w_1,w_2 \in W$.
Now, on the one hand
$$ \langle w \wedge A,B \rangle_{\Hom(V,\Lambda_2(W))}= \langle \alpha \otimes (w \wedge \tilde w) ,\beta \otimes (w_1 \wedge w_2) \rangle_{\Hom(V,\Lambda_2(W))}=$$
$$ \langle \alpha , \beta \rangle_{V^*} \langle w \wedge \tilde w ,w_1 \wedge w_2 \rangle_{\Lambda_2(W)}. $$
On the other hand
$$ \tr_{V} (A \otimes B)= \langle \alpha , \beta \rangle_{V^*} \tilde w \otimes (w_1 \wedge w_2).$$
Thus, it's enough to show
$$ \langle w \wedge \tilde w ,w_1 \wedge w_2 \rangle_{\Lambda_2(W)}=\langle w , h_W\big(\tilde w \otimes (w_1 \wedge w_2)\big) \rangle_W,$$
but this nows follows directly form the definition of the induced inner product on $\Lambda_2(W)$, and the definition of $h_W$ (see $(6)$).
Using the above proposition, we deduce that $$ A(\phi)=h_{\phi^*TN}\bigg(\tr_{\TM}\big(d\phi \otimes \delta_{\nabla^{\Lambda_2(\phi^*{\TN})}}(d\phi \wedge d\phi)\big)\bigg).$$