The total grand canonical partition function is $$\mathcal{Z} = \sum_{all\ states}{e^{-\beta(E-N\mu)}} = \sum_{N=0}^\infty\sum_{\{E\}}{e^{-\beta(E-N\mu)}}$$
For Bose-Einstein or Fermi-Dirac, the energy eigenstates are countable since $E=\sum_{\epsilon}{n_\epsilon \epsilon}$, then
$$\mathcal{Z} = \sum_{N=0}^{\infty}\sum_{\{n_\epsilon\}}e^{-\beta\sum_{\epsilon}{n_\epsilon(\epsilon-\mu)}} = \sum_{N=0}^{\infty}\sum_{\{n_\epsilon\}}\prod_\epsilon e^{-\beta n_\epsilon(\epsilon-\mu)}= (\sum_{n_0}e^{-\beta{n_0(\epsilon_0-\mu)}})(\sum_{n_1}e^{-\beta{n_1(\epsilon_1-\mu)}})...$$
This method is used in Pathria's Statistical Mechanics book (page 133). I am having trouble following the argument behind how the last step occurs. The book states the following:
"The double summation (first over the numbers $n_\epsilon$ constrained by a fixed value of the total number N, and then over all possible values of N) is equivalent to a summation over all possible values of the numbers $n_\epsilon$ independent of one another"
Why is this true? I have been stuck on this argument for hours...