= Solution
Set $s=a(t)/a(t_{\mathrm{dec}})$. Free propagation in the homogeneous expanding universe redshifts frequencies by $\nu=\nu_{\mathrm{dec}}/s$ and dilutes physical number density by $s^{-3}$. Number conservation in corresponding frequency intervals gives
$$
n_\nu(t)d\nu=s^{-3}n_{\nu_{\mathrm{dec}}}(t_{\mathrm{dec}})d\nu_{\mathrm{dec}},\qquad d\nu_{\mathrm{dec}}=s\,d\nu.
$$
Inserting the initial Planck spectrum, the factor $s^{-2}(s\nu)^2$ becomes $\nu^2$, while its exponent is $hs\nu/(kT_{\mathrm{dec}})$. Hence
$$
\boxed{n_\nu(t)=\frac{8\pi\nu^2}{c^3}\frac1{e^{h\nu/(kT(t))}-1},\qquad
T(t)=\frac{a(t_{\mathrm{dec}})}{a(t)}T_{\mathrm{dec}}.}
$$
The shape is preserved even without collisions maintaining <thermal equilibrium>; it follows from redshift of an initially Planckian distribution.
For the total <photon number density> put $x=h\nu/(kT)$. For the <energy density> include the energy $h\nu$ of each <photon>. Then
$$
\boxed{n=\alpha T^3,\qquad \epsilon=\xi T^4,}
$$
where $\alpha=(8\pi/c^3)(k/h)^3\int_0^\infty x^2/(e^x-1)\,dx$ and $\xi=(8\pi h/c^3)(k/h)^4\int_0^\infty x^3/(e^x-1)\,dx$. Both integrals are finite positive constants. Thus $n\propto a^{-3}$ and $\epsilon\propto a^{-4}$, with the extra power supplied by <photon>-energy redshift.
Back to article page
