Solution (source code)

= Solution

Let $V$ be proper volume, let $g_i$ count the populated internal states, and take $\mu_i=0$ and $\epsilon=pc$. Spherical integration in momentum space gives
$$
\frac{N_i}{V}=\frac{g_i}{2\pi^2}\left(\frac{k_{\rm B}T}{\hbar c}\right)^3
\int_0^\infty\frac{x^2\,dx}{e^x\mp1},\qquad x=\frac{pc}{k_{\rm B}T}.
$$
Using the <fermion-to-boson thermal integral ratio>, the minus-denominator <boson> integral is $2\zeta(3)$ and the plus-denominator <fermion> integral is $3\zeta(3)/2$. Hence
$$
\boxed{N_B=\frac{g_B\zeta(3)}{\pi^2}V\left(\frac{k_{\rm B}T}{\hbar c}\right)^3,\qquad
N_F=\frac34\frac{g_F\zeta(3)}{\pi^2}V\left(\frac{k_{\rm B}T}{\hbar c}\right)^3.}
$$
The energy integral has one extra factor $pc=k_{\rm B}Tx$:
$$
\frac{U_i}{V}=\frac{g_i}{2\pi^2}\frac{(k_{\rm B}T)^4}{(\hbar c)^3}
\int_0^\infty\frac{x^3\,dx}{e^x\mp1}.
$$
The supplied distribution identity gives the <fermion> integral as $(1-2^{-3})\pi^4/15=7\pi^4/120$. Therefore
$$
\boxed{U_B=\frac{g_B\pi^2}{30}V\frac{(k_{\rm B}T)^4}{(\hbar c)^3},\qquad
U_F=\frac78\frac{g_F\pi^2}{30}V\frac{(k_{\rm B}T)^4}{(\hbar c)^3}.}
$$
For an isotropic <ultrarelativistic> gas, momentum flux gives $P_i=\tfrac13\int p v f_i\,d^3p$. Since $v=c$ and $pc=\epsilon$,
$$
\boxed{P_i=\frac{U_i}{3V}.}
$$
<Photons> have $g_\gamma=2$, <electrons> and <Positrons> each have $g=2$, and each <neutrino> or antineutrino species has one populated <helicity> in the standard thermal bath.