The problem
Let's have axion lagrangian (without aF∧F terms) in FRLW Universe (gμν=diag(1,−R2(t),−R2(t),−R2(t))):
L=12gμν∂μa∂νa−m2a(t)2a2
Corresponding EOM for spatially homogeneous field reads
¨a+3H˙a+m2a(t)a=0,H=˙RR
Corresponding solution for large times and for adiabatically changing mass |˙ma|<m2a reads
a(t)≈a0(R(tin)R(t))32cos(∫ma(t)dt)
Corresponding energy density reads (here I've inserted the solution (0) and used relation H<<ma(t), as typical H∼1t, while ma∼t2)
ρa=T00=12((∂0a)2+m2aa2)≈m2a(t)a20(R(tin)R(t))3≈Cam2a(T)T3,
where I've used relations R(t)∼√t and t∼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):
ρa≈˜Cama(T)T3,
It could be obtained from the following way: from EOM ¨a+3H˙a+m2aa=0 we can obtain EOM for ρ=12˙a2+12m2aa2,
˙ρa=˙a¨a+˙mamaa2+˙aam2a=−3˙a2H+˙mamaa2
Let's average over an oscillations, ⟨m2aa2⟩=⟨˙a2⟩=ρa (this action implies conditions H<<ma(t),˙maa<<m2a(t)), and thus
˙ρa=−3Hρa+˙mamaρa,
and the solution reads
ρa≈(R(tin)R(t))3ma(t)ρa(tin)=˜CaT3ma(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,˙mama<<ma, so it seems that they must coincide.
Can you explain me where is the mistake in derivation of expression (1)?