Solution (source code)

= Solution

For one species with two internal states, the <relic-neutrino energy density> obtained from the frozen <Fermi-Dirac distribution> is
$$
\rho_\nu=2\int\frac{d^3p}{(2\pi)^3}
\frac{\sqrt{p^2+m_\nu^2}}{e^{ap/T_\nu}+1}
=\frac1{\pi^2}\int_0^\infty
\frac{p^2\sqrt{p^2+m_\nu^2}}{e^{ap/T_\nu}+1}\,dp.
$$
In the relativistic limit, set $x=ap/T_\nu$ and use the standard <Riemann zeta function> integral
$$
\int_0^\infty\frac{x^3}{e^x+1}\,dx=\frac{7\pi^4}{120}.
$$
This gives
$$
\boxed{\rho_\nu=\frac{7\pi^2}{120}\left(\frac{T_\nu}{a}\right)^4}.
$$