Evaluating the time average over energy 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,

 A: 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.
