Solution (source code)

= Solution

For bosons the occupation number is $n_{\mathbf p}=0,1,2,\ldots$. Summing the <geometric series> for the single-mode <grand canonical ensemble> gives the <Bose-Einstein distribution>
$$
\boxed{\langle N(\mathbf p)\rangle
=\frac1{e^{\beta(E_{\mathbf p}-\mu)}-1}}.
$$