Solution (source code)

= Solution

For noninteracting particles, a many-particle configuration has $E=\sum_jn_j\epsilon_j$ and $N=\sum_jn_j$. Summing independently over allowed occupations therefore factorizes the grand <partition function>:
$$
\Theta=\prod_j\Theta_j,\qquad
\boxed{\Phi=-k_BT\sum_j\ln\Theta_j.}
$$
At fixed $T,V$, differentiating the original state sum gives $\partial_\mu\ln\Theta=\beta\langle N\rangle$ and $\partial_\mu^2\ln\Theta=\beta^2\operatorname{Var}N$, so
$$
\boxed{\langle N\rangle=-\partial_\mu\Phi,\qquad
\operatorname{Var}N=-k_BT\partial_\mu^2\Phi
=k_BT\partial_\mu\langle N\rangle.}
$$
For one level, write $q=e^{-\beta(\epsilon-\mu)}$. <Fermions> allow $n=0,1$, giving $\Theta_F=1+q$ and $\bar n=q/(1+q)$. <Bosons> allow all nonnegative integers, giving $\Theta_B=1/(1-q)$ and $\bar n=q/(1-q)$, with $q<1$. Thus
$$
\boxed{\bar n_F=\frac1{e^{\beta(\epsilon-\mu)}+1},\qquad
\bar n_B=\frac1{e^{\beta(\epsilon-\mu)}-1}.}
$$
The <grand-canonical single-level occupation variance> follows from the same derivative identity applied to one level:
$$
\boxed{\operatorname{Var}n=k_BT\partial_\mu\bar n
=\begin{cases}\bar n(1-\bar n),&\text{fermions},\\
\bar n(1+\bar n),&\text{bosons}.\end{cases}}
$$
In the high-temperature nondegenerate regime, $\bar n\ll1$ and both approach classical Poisson-like fluctuations, $\operatorname{Var}n\simeq\bar n$. At low <temperature>, nearly filled or empty <fermion> levels have strongly suppressed fluctuations; only the Fermi-surface window changes occupation. Highly occupied low-energy <boson> levels can instead have $\operatorname{Var}n\sim\bar n^2$, showing bosonic enhancement. This last statement depends on occupation and reservoir <chemical potential>, not <temperature> alone; grand-canonical condensate fluctuations must not be confused with a fixed-total-number ensemble.