The axion fields keeps oscillating around \(\theta _a=0\)  until nowadays and its mass is almost constant, \(\dot{m_a} \simeq 0\) . It is adequate to apply the WKB approximation provided \(m_a\gg \big \lbrace \dot{m}_{a}/{m}_{a},\dot{H}/H\big \rbrace \) , i.e. provided the oscillations are much faster than both the expansion rate of the Universe and the rate at which the mass changes, which is the case here. The adiabatic approximation predicts the existence of a conserved quantity [1]}, [2]}, [3]}, an adiabatic invariant \(N_a\)  that can be interpreted as the comoving number of non-relativistic axion quanta, defined for a generic axion mass \(m_a\)  as
\(N_a\equiv \frac{\rho _a a^3}{m_a}= {\rm const.}\, , \qquad \text{with}\qquad {\rho _{a}} \equiv \frac{1}{2} {\dot{\theta }}_a^2 {f_a^2}+\frac{1}{2} m_{a}^2{f_a^2} \theta _a^2\,,\) 
