Query about derivation for the grand partition function for an ideal fermi gas I am having trouble understanding the following derivation:

I am having a bit of trouble understanding what it means to sum over all states. 
So this is how I am interpreting the above.
We have a gas of $N$ particles in different states $n_1, n_2, n_3,...$ with energies $Ɛ_1, Ɛ_2, Ɛ_3,...$ etc. 
In the above derivation $E$ is just the total energy, given by $$\sum_{i=1}^{i=\infty} {n_iƐ_i}, $$
and $N$ is just the sum of the occupancies $n_1 + n_2 + n_3 + ...$.
My issue might be more with the mathematical notation than the physics but when we have a sum like the one above, the different terms in the sum are given by the indice $i$, i.e $$\sum_{i=1}^{i=\infty} {n_iƐ_i} ,$$ = $n_1Ɛ_1 + n_2Ɛ_2 + n_3Ɛ_3 +...$.
What are the different terms in the sum given by line 2 in the above derivation? 
I can't even write the second term in the above sum as I don't see what the summing indice is. 
Secondly, I don't understand what the logic is going from line 3 to line 4 is.
If someone understands the above derivation can you please explain it.
 A: The index $i$ in $n_i$ and $\epsilon_i$ refers to the label of the one-particle states. I.e., the Hamiltonian $\hat H_N$ of the non-interacting gas of $N$ particles is written as a sum of $N$ one-particle Hamiltonian $\hat h_{\alpha}$:
$$
\hat H_N= \sum_{\alpha=1}^N \hat h_{\alpha}
$$
and the eigenstates of $\hat h_{\alpha}$ are the states $\left| i\right>$ such that
$$
\hat h_{\alpha}\left| i\right>=\epsilon_i \left| i\right>.
$$
In general, the number of the one-particle states is infinite ($i=1, \dots, \infty$).
$N$-particle states are obtained as symmetrized or anti-symmetrized tensor products of $N$ of such one-particle states. The corresponding $N$-particle energy, $E_N$, is just the sum of the one particle energy for each state, each multiplied by the number of particles which are in that state  (the occupation numbers $n_i$):
$$
E_N= \sum_{i=1}^{\infty} n_i \epsilon_i,
$$
with the constraint $\sum_{i=1}^{\infty} n_i=N$. Therefore, one could say that a single $N$_particle state is uniquely identified by the infinite set of occupation numbers $\left\{ n_i \right\}$, or equivalently by the ordered sequence $(n_1,n_2,n_3,\dots)$, always with the constraint $\sum_{i=1} n_i=N$.
Now, let's go to grand-canonical statistical mechanics. Your starting expression is a possible way to slightly simplify the derivation, by starting immediately with a sum ove all the states, without any constraint on the number of particles. It is a way to get almost immediately the correct final result for the gran-canonical partition function. 
It could be interesting and probably pedagogically more useful to start with the expression for the gran canonical partition function, written as:
$$
\mathcal{Z} = \sum_{N=0}^{\infty} e^{\beta \mu N} Q_N ~~~~~~~~~~~[1]
$$
where $Q_N$ is the canonical partition function for a system of $N$ particles.
$$
Q_N=\sum_{\left\{ n_i \right\}; \sum_{i=1} n_i=N   } e^{-\beta E_N(\left\{ n_i \right\})}.
$$
The constraint $ \sum_{i=1} n_i=N $ can be absorbed into a Konecker's delta, in order to rewrite $Q_N$ as an unconstrained the sum over all possible values of the occupation numbers:
$$
Q_N=\sum_{\left\{ n_i \right\}  } e^{-\beta E_N(\left\{ n_i \right\})}\delta_{N,\sum_{i=1} n_i}~~~~~~~~~~[2]
$$.
Using this expression [2], the gran canonical partition function [1] can be rewritten as:
$$
\mathcal{Z} = \sum_{N=0}^{\infty} e^{\beta \mu N}\sum_{\left\{ n_i \right\}  } e^{-\beta E_N(\left\{ n_i \right\})}\delta_{N,\sum_{i=1} n_i}=
\sum_{N=0}^{\infty} \sum_{\left\{ n_i \right\}  } e^{\beta \mu N}e^{-\beta E_N(\left\{ n_i \right\})}\delta_{N,\sum_{i=1} n_i}=
\sum_{N=0}^{\infty} \sum_{\left\{ n_i \right\}  } e^{\beta \mu \sum_{i=1}^N n_i}e^{-\beta E_N(\left\{ n_i \right\})}\delta_{N,\sum_{i=1} n_i}=
\sum_{N=0}^{\infty} \sum_{\left\{ n_i \right\}  } e^{-\beta \sum_{i=1}^{\infty} n_i (\epsilon_i-\mu)}\delta_{N,\sum_{i=1} n_i}
$$ 
In the last formula, it is possible to invert the order of summation, leaving an external, unconstrained sum over all the possible choices of occupation numbers, while the sum over $N$ of the Kronecker's delta, for each choice of of the occupation numbers $\left\{ n_i \right\} $, is equal to $1$ (only one value of $N$ is equal to $\sum_i n_i $ ).
Therefore, we have arrived to the unconstrained sum over the occupation numbers.
Its explicit meaning is 
$$
\sum_{\left\{ n_i \right\}  }=\sum_{n_1}\sum_{n_2}\sum_{n_3}\dots.
$$
So that, if each term in the sum is factorized [1] gets factorized as 
$$
\left(\sum_{n_1} e^{-\beta n_1 (\epsilon_1-\mu)}\right)\left(\sum_{n_2} e^{-\beta n_2 (\epsilon_2-\mu)}\right)\left(\sum_{n_3} e^{-\beta n_3 (\epsilon_3-\mu)}\right)\dots.
$$
At this point the exact calculation depends on the statistics (two values of each occupation number for fermions, infinite values for bosons), but it is straightforward.
