For more info see the article equations 37
Edit: The $\varepsilon ^3 $ has vanished due to time average. But how to get the 4th order?
Let us define some function for scalar field $$\phi= \sum_{k=1}^\infty\varepsilon^k \phi_k \tag{1}=\varepsilon^1 \phi_1+\varepsilon^2 \phi_2 +\varepsilon^3 \phi_3 ...$$ position scaling \begin{equation}\tag{2} \zeta^i=\varepsilon x^i \end{equation} time scaling \begin{equation}\tag{3} \tau=\omega(\varepsilon) t \end{equation} \begin{equation} \Delta S -S +S^3=0\,,\quad S=p_1\sqrt{\lambda} \,. \end{equation} \begin{equation} p_3=\frac{1}{\lambda^2\sqrt{\lambda}}\left[ \left(\frac{1}{24}\lambda^2-\frac{1}{6}\lambda g_2^2+\frac{5}{8}g_5 -\frac{7}{4}g_2g_4+\frac{35}{27}g_2^4\right)Z -\frac{1}{54}\lambda g_2^2S(32+19S^2) \right] \end{equation} the solutions for equation (1) can be written as, \begin{eqnarray} \phi_1&=&p_1\cos \tau \tag{4}\\ \phi_2&=&\frac16 g_2p_1^2\left(\cos(2\tau)-3\right) \tag{5}\\ \phi_3&=&p_3\cos \tau+\frac{1}{72}(4g_2^2-3\lambda)p_1^3\cos(3\tau)\tag{6} \end{eqnarray} My intention is to determine the time average energy density
The energy corresponding can be written as \begin{equation} E = \int d^{D}x\,{\cal E}\,, \quad {\cal{E}} = \frac{1}{2} \left(\partial_t \phi\right)^2 + \frac{1}{2} \left(\partial_i \phi\right)^2 +U(\phi)\,, \end{equation} where ${\cal{E}}$ denotes the energy density. \begin{equation} E = \frac{1}{\varepsilon^D}\int d^{D}\zeta\,{\cal E}\,, \quad{\rm where}\quad {\cal{E}} = \frac{1}{2} (1-\varepsilon^2)\left(\partial_\tau \phi\right)^2 + \varepsilon^2\frac{1}{2} \left(\partial_i \phi\right)^2 +U(\phi)\,. \end{equation} Because of the periodic time dependence we shall compute the energy density averaged over a period, \begin{equation} \bar{E}=\frac{1}{2\pi}\int_0^{2\pi}d\tau E\,, \end{equation} The bar over a quantity will denote its time average. Using the results of the $\varepsilon$ expansion, Eqs. (4-6) the time averaged energy density, $\bar{{\cal E}}$, up to fourth order in $\varepsilon$ can be written as (My problem is here that how will this term arise) \begin{eqnarray} \bar{{\cal E}}&=&\frac{\varepsilon^2}{2\lambda}S^2 -\frac{\varepsilon^4}{216\lambda^3}\biggl[ \lambda S^2(64g_2^2+27\lambda)(S^2+2) -54\lambda^2(\nabla S)^2 \nonumber\\ &&-SZ(135g_5-378g_2g_4+280g_2^4-36\lambda g_2^2+9\lambda^2) \biggl] \,. \end{eqnarray} For more info see the article equations 37
If you have problem to get the question then ask me please, I will edit for you. Thanks advance,
Look,
$$ {\cal{E}} = \frac{1}{2} \left(\partial_t \phi\right)^2 + \frac{1}{2} \left(\partial_i \phi\right)^2 +U(\phi)$$
Now fill in the expansion
$$\phi= \varepsilon^1 \phi_1+\varepsilon^2 \phi_2 +\varepsilon^3 \phi_3 + \ldots$$
this should give
$${\cal{E}} = \frac{1}{2} \left(\varepsilon^1 \partial_t\phi_1+\varepsilon^2 \partial_t\phi_2 +\varepsilon^3 \partial_t\phi_3 + \ldots \right)^2 + \frac{1}{2} \left(\varepsilon^1 \partial_i\phi_1+\varepsilon^2 \partial_i\phi_2 +\varepsilon^3 \partial_i\phi_3 + \ldots\right)^2 +\frac{1}{8}\left(\varepsilon^1 \phi_1+\varepsilon^2 \phi_2 +\varepsilon^3 \phi_3 + \ldots\right)^2\left( \varepsilon^1 \phi_1+\varepsilon^2 \phi_2 +\varepsilon^3 \phi_3 + \ldots -2\right)^2$$
This is to second order in $\varepsilon$
$${\cal{E}}=\varepsilon^2 \left[\frac{1}{2}(\partial_t \phi_1)^2+\frac{1}{2}(\partial_i \phi_1)^2+\frac{1}{2}(\phi_1)^2\right] + O(\varepsilon^3)$$
Filling in $\phi_1 = p_1 \cos \tau$, you get
$${\cal{E}}=\varepsilon^2 \left[\frac{1}{2}(p_1 \omega(\varepsilon)\sin\tau)^2+\frac{1}{2}(\partial_ip_1 \cos\tau)^2+\frac{1}{2}(p_1 \cos\tau)^2\right] + O(\varepsilon^3)$$
or
$${\cal{E}}=\varepsilon^2 \left[\frac{1}{2}p_1^2+\frac{1}{2}(\partial_ip_1 \cos\tau)^2\right] + O(\varepsilon^3)$$
With the rescaling of the spatial coordinates introduced in the paper, the second term is proportional to $\varepsilon^2$ and thus as a whole to $\varepsilon^4$. Then using $S=p_1\sqrt{\lambda}$ we get
$${\cal{E}}=\varepsilon^2 \frac{1}{2\lambda}S^2 + O(\varepsilon^3)$$
So, if you continue working diligently in the same fashion, you should be able to find out what the fourth order correction is.