Consider grand canonical partition function $Z = \prod_{k} \left(1-e^{-\beta(\varepsilon_k-\mu)}\right)^{-1}$ for an ideal Bose-gas in a trapping potential with energy levels $\varepsilon_k$ and chemical potential $\mu$. I try to find the fluctuation in particle number $\delta n_k^2$ with this equation:
$$\delta n_k^2 = \frac{1}{\beta^2}\frac{\partial^2}{\partial^2\mu} \log Z_k $$
What I find is the following surprising result:
$$ \delta n_k^2 = \left\langle n_k \right\rangle^2 + \left\langle n_k \right\rangle $$
this means that the number particle fluctuation is greater than mean occupation number, namely: $\delta n_k > \left\langle n_k \right\rangle $ $\forall \varepsilon_k$.
My question: this is clearly unphysical, is there a relative simple way to compute the correct mean particle number fluctuation $\delta n_k^2$ ?