Energy density, pressure and temperature for massive neutrinos in cosmology

Question

I want to be able to numerically compute the mean energy density and pressure for a massive neutrino species in cosmology, at any given scale factor $a$. These are given in terms of the distribution function $f(a, p)$ as

$$ \rho(a) = 2\int \frac{4\pi p^2 dp}{(2\pi)^3} f(a,p)\sqrt{m^2 + p^2}\,, $$ $$ P(a) = 2\int \frac{4\pi p^2 dp}{(2\pi)^3} f(a,p)\frac{p^2}{3\sqrt{m^2 + p^2}}\,, $$

where $m$ is some mass. In the early Universe, typical momenta were much larger than the mass, $\sqrt{\langle p^2\rangle}\gg m$ (the neutrinos are relativistic), in which case the Fermi-Dirac distribution of the neutrinos was

$$ f(a, p) = \biggl[\exp\biggl(\frac{p}{T(a)}\biggr) + 1\biggr]^{-1} $$

with some temperature $T(a)$. At later times, the neutrinos keep this distribution but with a decreasing $T(a)$.

I know the temperature at early times, and so I can compute $\rho(a)$ and $P(a)$ at these times. As long as the neutrinos remain relativistic, the temperature goes like $T(a)\propto a^{-1}$. When they are completely non-relativistic, I think this changes to $T(a)\propto a^{-2}$, and so there exist some transition period in which $T(a)$ goes from one behavior to the other. How can I calculate this behavior, so that I can define $f$ at any time $a$ and ultimately compute $\rho(a)$ and $P(a)$ for all $a$, without treating the neutrinos as just ultra-relativistic or completely non-relativistic?

Answer 1

It turns out that $T(a)\propto a$ holds for all $a$, even at non-relativistic times. This can e.g. be seen from computing the number density $n(a)$ (simply the integral over $f$ without any additional weight), from which one find $n\propto T^3$. After decoupling, conservation of neutrino numbers imply $n\propto a^{-3}$ regardless of the neutrinos be relativistic or not, and so $T\propto a^{-1}$.

The confusion about $T\propto a^{-2}$ stems from the fact that this relation holds for a non-relativistic species in thermal equilibrium, which the neutrinos are not after the time of decoupling