Solution (source code)

= Solution

Substitute $x=pc/(k_BT)$ in the <Fermi-Dirac distribution> energy <integral>. Then
$$
\epsilon=\frac{4\pi k_B^4}{h^3c^3}T^4\int_0^\infty\frac{x^3}{e^x+1}\,dx.
$$
The <integral> converges, is positive, and is independent of temperature. Thus \b[$\epsilon=\alpha T^4$], with $\alpha$ equal to the displayed constant <integral>. The given normalization is for the specified degeneracy; further internal degrees of freedom multiply it.