The study of higher order energy functionals was first proposed by Eells and Sampson in 1965 and, later, by Eells and Lemaire in 1983. These functionals provide a natural generalization of the classical energy functional. More precisely, Eells and Sampson suggested the investigation of the so-called $ES-r$-energy functionals $ E_r^{ES}(\varphi) = (1/2)\int_{M}\, |(d^*+d)^r (\varphi)|^2\, dV$, where $r \geq 2 $ and $ \varphi:M \to N$ is a map between two Riemannian manifolds. The initial part of this paper is a short overview on basic definitions, properties, recent developments and open problems concerning the functionals $ E_r^{ES}(\varphi)$ and other, equally interesting, higher order energy functionals $E_r(\varphi)$ which were introduced and studied in various papers by Maeta and other authors. If a critical point $\varphi$ of $E_r^{ES}(\varphi)$ (respectively, $E_r(\varphi)$) is an isometric immersion, then we say that its image is an $ES-r$-harmonic (respectively, $r$-harmonic) submanifold of $N$. We observe that minimal submanifolds are trivially both $ES-r$-harmonic and $r$-harmonic. Therefore, it is natural to say that an $ES-r$-harmonic ($r$-harmonic) submanifold is proper if it is not minimal. In the special case that the ambient space $N$ is the Euclidean space $\mathbb{R}^n$ the notions of $ES-r$-harmonic and $r$-harmonic submanifolds coincide. The Chen-Maeta conjecture is still open: it states that, for all $r \geq2$, any proper, $r$-harmonic submanifold of $\mathbb{R}^n$ is minimal. In the second part of this paper we shall focus on the study of $G = {\rm SO}(p+1) \times {\rm SO}(q+1)$-invariant submanifolds of $\mathbb{R}^n$, $n = p+q+2$. In particular, we shall obtain an explicit description of the relevant Euler-Lagrange equations in the case that $r = 3$ and we shall discuss difficulties and possible developments towards the proof of the Chen-Maeta conjecture for $3$-harmonic $G$-invariant hypersurfaces.