Axion field energy density

Question

The problem

Let's have axion lagrangian (without $a F \wedge F $ terms) in FRLW Universe ($g_{\mu \nu} = diag(1, -R^{2}(t), -R^{2}(t), -R^{2}(t))$):
$$
L = \frac{1}{2}g_{\mu \nu}\partial^{\mu}a\partial^{\nu}a - \frac{m^{2}_{a}(t)}{2}a^{2}
$$

Corresponding EOM for spatially homogeneous field reads
$$
\ddot{a} + 3H\dot{a} + m^{2}_{a}(t)a = 0, \quad H = \frac{\dot{R}}{R}
$$

Corresponding solution for large times and for adiabatically changing mass $|\dot{m}_{a}| < m_{a}^{2}$ reads
$$
\tag 0 a(t) \approx a_{0}\left(\frac{R(t_{\text{in}})}{R(t)}\right)^{\frac{3}{2}}cos(\int m_{a}(t)dt)
$$

Corresponding energy density reads (here I've inserted the solution $(0)$ and used relation $H << m_{a}(t)$, as typical $H \sim \frac{1}{t}$, while $m_{a} \sim t^{2}$)
$$
\tag 1 \rho_{a} = T_{00} = \frac{1}{2}\left((\partial_{0}a)^{2} + m_{a}^{2}a^{2} \right) \approx m_{a}^{2}(t)a_{0}^{2}\left( \frac{R(t_{\text{in}} )}{R(t)}\right)^{3} \approx C_{a}m_{a}^{2}(T)T^{3},
$$

where I've used relations $R(t) \sim \sqrt{t}$ and $t \sim T^{-2}$.

But recently I've found the true expression for axion energy density (for example, in Rubakov's "Introduction to the Theory of the Early Universe: Hot Big Bang Theory", chapter 9):
$$
\tag 2 \rho_{a} \approx \tilde{C}_{a}m_{a}(T)T^{3},
$$

It could be obtained from the following way: from EOM $\ddot{a} + 3H \dot{a} + m_{a}^{2}a = 0$ we can obtain EOM for $\rho = \frac{1}{2}\dot{a}^{2} + \frac{1}{2}m_{a}^{2}a^{2}$,
$$
\dot{\rho}_{a} = \dot{a}\ddot{a} + \dot{m}_{a}m_{a}a^{2} + \dot{a}am_{a}^{2} = -3\dot{a}^{2}H + \dot{m}_{a}m_{a}a^{2}
$$

Let's average over an oscillations, $\langle m_{a}^{2}a^{2}\rangle = \langle \dot{a}^{2}\rangle = \rho_{a}$ (this action implies conditions $H << m_{a}(t), \dot{m_{a}}{a} << m_{a}^{2}(t)$), and thus
$$
\dot{\rho}_{a} = -3H\rho_{a} + \frac{\dot{m}_{a}}{m_{a}}\rho_{a},
$$

and the solution reads
$$
\rho_{a} \approx \left(\frac{R(t_{\text{in}})}{R(t)}\right)^{3}m_{a}(t)\rho_{a}(t_{\text{in}}) = \tilde{C}_{a}T^{3}m_{a}(T)
$$

The question

Expressions $(1)$ and $(2)$ for energy density for large times don't coincide; in particular, they state different behaviours of energy density with time (due to time dependence of axion mass). But both expressions are derived by taking into account conditions $H, \frac{\dot{m_{a}}}{m_{a}} << m_{a}$, so it seems that they must coincide.

Can you explain me where is the mistake in derivation of expression $(1)$?