$\newcommand{\TM}{\operatorname{TM}}$ $\newcommand{\M}{\mathcal{M}}$ $\newcommand{\N}{\mathcal{N}}$ $\newcommand{\TM}{\operatorname{T\M}}$ $\newcommand{\TN}{\operatorname{T\N}}$ $\newcommand{\TstarM}{T^*\M}$
In a nutshell: Conformal maps form a natural class of maps whose pullback preserve co-closed forms in the right degree. The question is whether or not these are the only examples.
Let $\M,\N$ be $d$-dimensional smooth Riemannian manifolds. Let $f:\M \to \N$ be smooth. Fix $1 \le k < d$.
We say that $f$ is a $k$-order harmonic morphism if pullback by $f$ locally preserves co-closedness of closed and co-closed forms. (For $k=1$ this is the standard notion of "harmonic morphism").
That is, for every open $U \subseteq \N$,
$\omega\in \Omega^k(U)$ is closed and co-closed $\Rightarrow f^*\omega$ is co-closed.
For even $d$, every conformal map is $\frac{d}{2}$-harmonic morphism, as follows immediately from the conformal invariance of the Hodge dual at degree $\frac{d}{2}$.
Is it true that every $k$-order harmonic morphism is conformal or constant?
If it helps, I am ready to assume that $f$ preserves co-closedness of all forms, not just closed ones. (I am also fine with assuming $df$ is everywhere invertible).
For $k=1$, it is known that every harmonic-morphism is "horizontally weakly conformal", which for equidimensional manifolds $\M,\N$ implies weak conformality. (i.e. $f$ is conformal or constant). This is proved by "testing" the condition against many harmonic functions. However, for $k>1$ I don't see how to adapt this argument.
Comment: If $k \neq \frac{d}{2}$, then any conformal $k$-morphism is a homothety: (The conformal factor is constant). Before proving this below, we note this is analogous to the following "pointwise phenomena": A linear map which commutes with the Hodge dual in degree $k<d$ is conformal, and if $k \neq \frac{d}{2}$ it is an isometry.
Proof:
Suppose that $\omega \in \Omega^k(U)$ is closed and co-closed. Then $d(\star \omega)=0 \Rightarrow d(f^*\star \omega)=0$. Since for conformal maps $f^*(\star \omega)=\pm |\det df|^{1-2\frac{k}{d}} \star f^*\omega$, we obtain $d(h \star f^*\omega)=0$ where $h:=|\det df|^{1-2\frac{k}{d}} $. If $f$ is also a $k$-morphism, then $d(\star f^*\omega)=0$. Now, $$ d(h \star f^*\omega)=d(\star f^*\omega)=0 \, \, \text{implies}$$ $$ dh \wedge \star f^*\omega=0 \, \, \text{for every closed and co-closed} \, \, \omega \in \Omega^k(U).$$
Since locally, the closed and co-closed forms make a frame for $\bigwedge^k T^*\N$, and $df$ is an isomorphism, it follows that $dh=0$, so $h$ is constant. Since $k \neq \frac{d}{2}$, it follows that $\det df$ is constant.
