Is this expression for the Laplacian of conformal maps between Riemannian manifolds known? $\newcommand{\M}{\mathcal{M}}$
$\newcommand{\N}{\mathcal{N}}$
$\newcommand{\Hom}{\operatorname{Hom}}$
$\newcommand{\tr}{\operatorname{tr}}$
$\newcommand{\TM}{\operatorname{T\M}}$
$\newcommand{\TN}{\operatorname{T\N}}$
$\newcommand{\sAverage}[1]{\langle#1\rangle} $
$\newcommand{\IP}[2]{\sAverage{#1,#2}}$
$\newcommand{\Cof}{\text{Cof}}$
Lately I derived the following equation (details are at the end):

Theorem:
    Let $(\M,g),(\N,h)$ be smooth $d$-dimensional oriented Riemannian manifolds, 
Let $\,f:\M \to \N\,$ be a smooth conformal map. Then
$$ \det df \cdot  \delta(df) + (\frac{2}{d}-1)     \tr_{g}\big(df \otimes d \det(df)  \big) =0. \tag{1} $$

($\delta$ is the adjoint of the pullback connection on $f^*{\TN}$. The equation is an equality of sections of $f^*{\TN}$.)
Equation $(1)$ has a few immediate (well-known) corollaries:


*

*In dimension $d=2$, conformal mappings are harmonic.

*Homotheties are harmonic (in every dimension).

*In dimension $d > 2$, conformal harmonic maps are homotheties (assuming $\M$ is connected): 
Let $e_i$ be an orthonormal frame for $\TM$. Then  $$ 0=\tr_{g} \big( df \otimes d \det(df) \big)= \sum_{i=1}^d  d \det(df)(e_i) \cdot df(e_i) .$$
Since $f$ is an immersion, $df(e_i)$ is a basis of $f^*\TN$.


(Summary: A conformal map is harmonic if and only if $d=2$ or it is a homothety).
Questions:
1. Is equation $(1)$ known? (Perhaps in a Euclidean setting? Or in coordinates?)
2. Is equation $(1)$ an Euler-Lagrange's equation of some functional?

Sketch of derivation:
My colleagues and I showed the following: 

Theorem(2):
    Let $(\M,g),(\N,h)$ be smooth $d$-dimensional oriented Riemannian manifolds, 
Let $\,f:\M \to \N\,$ be an arbitrary smooth map. Then
$$ \delta (\Cof df)=0. \tag{2}$$ 
  where $\Cof df:= (-1)^{d-1} \star_{f^*TN}^{d-1} (\wedge^{d-1} df) \star_{TM}^1.$ (It's a section  of $T^*\M \otimes f^*{\TN}$).

(For a proof see Proposition 3.4 here). 
Now, if $f$ is conformal, then
$$ \Cof df=(\det df)^{1-\frac{2}{d}} df. \tag{3}$$
Plug $(3)$ into $(2)$  (the rest is Leibniz rule + simplification).
 A: Well, for conformal maps equation $(1)$ is merely $d$-harmonicity in disguise:)
The equation is
$$
\delta\big((\det df)^{1-\frac{2}{d}} df\big)=0. \tag{1}
$$
Since for conformal maps, $\det df=\|df\|^d$ up to a constant, we equivalently get
$$
\delta\big(\|df\|^{d-2} df\big)=0. \tag{2}
$$
(Which is trivial of course, since conformal maps are source symmetries of the $d$-energy).
The funny part is that while  equation $(2)$ is an Euler-Lagrange's equation, equation $(1)$ is not (for $d \neq 2$). 
The moral here is that if we start from an E-L equation, write it differently for specific class of maps (the conformal maps in this example) then the resulted equation may not be an E-L equation.
Let's prove equation $(1)$ is not E-L:
Indeed, if it were an E-L equation (of an isometrically invariant smooth functional), this would imply the existence of a smooth map $ h:M_d \to \mathbb{R}$, such that $dh_X(V)=\langle (\det X)^r X,V\rangle $.
This is not possible since $d^2h\neq 0$:
$d^2h=\sum_{ij}(\det X)^rx_{ij}dx^{ij} \Rightarrow d^h=\sum_{ijks} \frac{\partial((\det X)^rx_{ij}) }{\partial x_{sk}} dx^{sk} \wedge dx^{ij}=r \sum_{ijks} x_{ij} \frac{\partial \det X }{\partial x_{sk}} dx^{sk} \wedge dx^{ij}$ which is zero iff 
$$x_{ij} \frac{\partial \det X }{\partial x_{sk}} =x_{sk} \frac{\partial \det X }{\partial x_{ij}} , $$ which is false. (The RHS does not depend on $x_{ij}$ while the LHS is dependent of it). 
