Bose-Einstein Grand Canonical partition function derivation step 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...
 A: 

"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...

The term
$$
\sum_N \sum_{\{n_\epsilon\}}\left(\ldots\right)
$$
means a sum over all N and for each fixed N a sum over all $n_\epsilon$ consistent with the fixed value of N.
But this is the same as summing over all $n_\epsilon$ regardless of the value of N (i.e., unconstrained summing over $n_\epsilon$. I.e.,
$$
\sum_N \sum_{\{n_\epsilon\}}\left(\ldots\right)
=
\sum_{n_1}\sum_{n_2}\sum_{n_3}{}_{\ldots}\left(\ldots\right)\;,
$$
where the sums over $n_1$, $n_2$, $n_3$, etc are unconstrained.
Consider this Fermionic example. Suppose there are only three single-particle state ($\epsilon$ can be: 1, 2, or 3). Then the different constrained values $n_\epsilon$ for each value of N can be decomposed as:
$$
N=0:\qquad n_1=0,n_2=0,n_3=0
$$
$$
N=1:\qquad n_1=0,n_2=0,n_3=1
$$
$$
N=1:\qquad n_1=0,n_2=1,n_3=0
$$
$$
N=1:\qquad n_1=1,n_2=0,n_3=0
$$
$$
N=2:\qquad n_1=0,n_2=1,n_3=1
$$
$$
N=2:\qquad n_1=1,n_2=0,n_3=1
$$
$$
N=2:\qquad n_1=1,n_2=1,n_3=0
$$
$$
N=3:\qquad n_1=1,n_2=1,n_3=1
$$
And you see there are 8=$2^3$ terms because this is the same as summing over each of the two possible values of $n_1$, $n_2$, and $n_3$ independently (for example, regardless of the values I pick for $n_1$ and $n_2$ I can find one state with $n_3=0$ and one state with $n_3=1$).  
This is true in general not just for the specific example above.
