How to derive the expression for Bose-Einstein distribution variance? 
Can anyone point me to a derivation of this expression? $n_s$ is the number of bosons in a state.
 A: Let $|\mathbf{n}\rangle = |n_0, n_1, n_2, \dots\rangle$ denote a state with $n_i$ particles in the $i^\mathrm{th}$ energy state of the single-particle hamiltonian.  Let $\epsilon_i$ denote the energy of the state with label $i$.  The $n_i$ are "occupation numbers" in the standard terminology.  Note that my notation is such that $i=0,1,2,\dots$ labels distinct energy eigenstates, not energy levels, so, in particular, in this notation one would have $\epsilon_i=\epsilon_j$ when $i\neq j$ if there were a degeneracy in the spectrum of the single-particle Hamiltonian.  
The grand canonical partition function can be written as follows:
$$
  Z = \sum_{\mathbf n} \prod_{i=0}^\infty x_i^{n_i}, \qquad x_i = e^{-\beta(\epsilon_i-\mu)}
$$
where the sum over $\mathbf n$ is a sum over admissible sequences non-negative integers.  For fermions such sequences can only consists of 1's and 0's, but for bosons there is no such restriction.
The population fraction (probability) of the system being in the state $|\mathbf n\rangle$ is therefore
$$
  p(\mathbf n) = \frac{1}{Z}\prod_{i=0}^\infty x_i^{n_i}
$$
So that the mean occupation number is
$$
  \langle n_i\rangle = \sum_{\mathbf n} n_i p(\mathbf n)
$$
A straightforward (but admittedly tedious) computation shows that
$$
  \frac{\partial}{\partial \epsilon_j}\langle n_i\rangle = -\beta\sum_{\mathbf n} (n_in_j-\langle n_j\rangle)p(\mathbf n)
$$
and setting $j=i$ in this expression gives the desired result.
I just went through the computation explicitly myself; it's not obvious (or wasn't to me) that it would work out, but it does.  If you have any trouble, I can post more detail.
