Solution (source code)

= Solution

The <Bose-Einstein distribution> gives one occupation for the nondegenerate ground level and $M$ identical contributions from the excited levels:
$$
\boxed{\langle N_0\rangle=\frac1{e^{-\beta\mu}-1},
\qquad
\langle N_\epsilon\rangle=\frac{M}{e^{\beta(\epsilon-\mu)}-1}.}
$$