Compute fluctuation in particle number in ideal Bose Gas

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$ ?

This is a completely correct result. If there is no condensate, then $\langle n_k\rangle\ll1$ for almost all the states. Then the total fluctuation in particle number would be, $$\langle \Delta N^2\rangle\sim N$$ So, $$\frac{\sqrt{\langle \Delta N^2\rangle}}{N}\sim\frac{1}{\sqrt{N}}\rightarrow0\quad\text{as}\quad N\rightarrow\infty$$ So, it does not give any absurd or weird results. If we have a condensate, then the fluctuations in the ground state, where all the particles crowd into, would be of $\mathcal{O}(1)$.