Grand Canonical Molecular adsorption onto a surface I am having a hard time applying the grand canonical theory to a simple example. I expose my understanding of the matter, the problem, my attempt of solution, the solution and my question on this solutions; I apologize for the lengthy question and will be very grateful to whoever feels like going through it!
I also add an answer including some ideas and a second solution. 
Preliminars
I follow Kubo, Statistical Mechanics, but having a look around the notation should be standard. An open system is in contact with a reservoir fixing temperature $T$ and chemical potential $\mu$. A microstate of the open system is denoted by $s$; the grand canonical partition function is
$$Z_G(T,V,\mu) = \sum_se^{-\beta E_s + \beta \mu N_s}$$ 
where $s$ denotes each available microstate of the system, $N_s$ the number of particles in that microstate, and $E_s$ the energy of that microstate. This can be related to the canonical partition function $Z$: let $l$ denote a microstate for fixed number of particles, then
$$Z_G = \sum_{N_s=0}^{N_{tot}}\left(\sum_le^{-\beta E_l} \right)e^{\beta\mu N_s} = \sum_{N_s=0}^{N_{tot}}Z(T,V,N_s)\, e^{\beta\mu N_s}$$
This is useful: $Z$ is relatively hard to compute due to the fixed number of particles condition, but $\sum_{N_s=0}^{N_{tot}}$ allows to get rid of this condition. We consider the single particle properties:

  
*
  
*$i$ runs over the single particle possible microstates
  
*$\epsilon_i$ denotes the energy of the state $i$, that is the energy that a single particle has when happens to be in the state $i$
  
*$n_i=$ is the occupation number of the state $i$, that is the number of particles that happen to be in the state $i$. For fermions
  $n_i=0,1$; for bosons $n_i=0,1,2,...$.
  

A microstate of the whole system $s$ is then specified by the sequence of occupation numbers $n_1, n_2, ...$, and
$$N_s=\sum_i n_i, \quad E_s = \sum_i \epsilon_i n_i $$
The canonical partition function is
$$Z=\sum_l e^{-\beta E_l} = \underbrace{\sum_{n_1}\sum_{n_2}...\sum_{n_i}...}_{\text{with the condition } N_s=\sum_i n_i} e^{-\beta E_s}$$
Plugging this in the grand canonical equation 
$$ Z_G = \sum_{N_s=0}^{N_{tot}} \underbrace{\sum_{n_1}\sum_{n_2}...\sum_{n_i}...}_{\text{with the condition } N_s=\sum_i n_i} e^{-\beta E_s} e^{\beta \mu N_s} = \underbrace{\sum_{n_1}\sum_{n_2}...\sum_{n_i}...}_{\text{on all possible values}} e^{-\beta \sum_i \epsilon_i n_i } e^{\beta \mu \sum_i n_i} = \prod _i\sum_{n_i}e^{-\beta(\epsilon_i-\mu)n_i}$$

We define the single state grand canonical partition function
  $$z_{G,i}=\sum_{n_i}e^{-\beta(\epsilon_i-\mu)n_i}$$
  $$Z_G=\prod_i z_{G,i}$$

The problem

We consider a gas in contact with a solid surface (e.g. argon on
  graphene or molecular nitrogen on iron, as in the Haber-Bosch
  synthesis). The gas molecules can be adsorbed at $N$ specific adsorption
  sites while one site can only bind one molecule. The energies of the
  bound and unbound state are $\epsilon$ and 0, respectively. The gas acts as a reservoir fixing $T$ and $\mu$.

How I would procede


*

*The system role is played by the $N$ adsorption sites

*The single particle role is played by one adsorption site

*The site admits two states, empty $i=0$ and full $i=1$

*The corresponding energies are $\epsilon_0=0$ and $\epsilon_1=\epsilon$

*The occupation numbers are $n_0$ = number of empty sites, $n_1$ = number of full sites. They both run from 0 to the total number of available sites, $n_i=0,1,...,N$

*A microstate of the system is determined by $n_0$ and $n_1$ such that $E_s=\sum_i\epsilon_i n_i = n_1 \epsilon$ and $N_s=\sum_i n_i = n_0+n_1=N$.


The grand canonical partition function should then read ($x_i:=e^{-\beta(\epsilon_i-\mu)}$)
$$Z_G=\prod_{i=0}^1\sum_{n_i=0}^N x_i^{n_i} = \prod_{i=0}^{1}\frac{1-x_i^{N+1}}{1-x_i}$$
Which is wrong (evaluating the product). 
What may be wrong


*

*Assigning to the sites the role of "single particle" the total number of particles is fixed, namely $N$, why it should be allowed to change


Question (see answer attempt)

What is wrong with this approach?


Solution A
With no further explanation beyond the fact that the sites are non-interacting the lecturer, this page and this page claim
$$Z_G=z_G^N$$
This $z_G$ is given as
$$z_G=1+e^{-\beta(\epsilon-\mu)}$$
Questions (still open)

  
*
  
*Is the $z_G$ used here the same as the single state grand canonical    partition function $z_{G,i}$ defined above?
  
*Where is $Z_G=z_G^N$ from?

The similar canonical relation $Z=z^N$ for non interacting systems of identical particles goes like this: we start with N distinguishable particles labeled by $j=1,...,N$; $\epsilon_{ij}$ is the $i$-energy level of the $j$-particle. Then
$$Z=\sum_l e^{-\beta E_l} = \sum_{i_1}\sum_{i_2}...\sum_{i_j}...e^{-\beta\sum_{j=1}^N \epsilon_{j  i_j}} = \left( \sum_{i_1} e^{-\beta \epsilon_{1i_1}} \right)...\left( \sum_{i_N} e^{-\beta \epsilon_{Ni_N}} \right)$$
The $j$ subscript can be dropped if the particles are identical, so that
$$z_j=z=\sum_i e^{-\beta \epsilon_i}$$
$$Z=\prod_{j=1}^N z = z^N$$
Question (still open)

The subscript can not be dropped in the relation $Z_G=\prod_i
 z_{G,i}$, as $z_{G,i}$ is an object strictly related to a state $i$,
   so again, how is $Z_G=z_G^N$ obtained?

 A: I'm not sure I will be able to clarify all your doubts, but this is the right approach to tackle this problem.
We consider $N_s$ available adsorption sites, an energy $\epsilon$ for each bound state, chemical potential $\mu$ and temperature $T$.
The grandpartition function $\mathcal Q$ is always expressed as
$$
\mathcal Q= \sum_{N=0}^\infty e^{\beta \mu N} Z_N
$$
in terms of the $N$-particle partition function. The latter is defined by
$$
Z_N=\sum_{\substack{N-\text{particle}\\ \text{states}}} e^{-\beta E(\text{state})}\,.
$$
In our case, the energy of a given $N$-particle state, of the ensemble of bound states, is $N\epsilon$ and there are $\binom{N_s}{N}$ such $N$-particle states, since each of the $N$ bound states can be placed by choosing one site among the $N_s$ sites, without repetition:
$$
Z_N=\binom{N_s}{N}e^{-\beta \epsilon N}\,.
$$
Note that this partition function does not bear a factorized form.
Finally,
$$
\mathcal Q= \sum_{N=0}^\infty \binom{N_s}{N}e^{\beta (\mu-\epsilon) N}
= \sum_{N=0}^{N_s} \binom{N_s}{N}e^{\beta (\mu-\epsilon) N}
=(1+e^{\beta(\mu-\epsilon)})^{N_s}\,,
$$
where in the last step we used the binomial formula
$$
(1+x)^n = \sum_{k=0}^n \binom{n}{k} x^k.
$$
EDIT: One can also reason directly using the grandpartition function as follows. Using the occupation number representation $\{n_{\alpha,k}\}$ for noninteracting particles, with $|\alpha, k\rangle$ labelling single-particle states with energy $\epsilon_\alpha$ and $g_\alpha$-fold degeneracy $k=1,2,\ldots,g_{\alpha}$,
$$
\mathcal Q= \sum_{N=0}^\infty e^{\beta\mu N} 
\sum_{\substack{
\{n_{\alpha, k}\}:\\
\sum_{\alpha, k}n_{\alpha,k}=N}} 
e^{-\beta \sum_{\alpha, k}n_{\alpha,k} \epsilon_\alpha }
=
\sum_{\{n_{\alpha,k}\}} e^{\beta \sum_{\alpha, k}n_{\alpha,k} (\mu-\epsilon_\alpha) } = \prod_{\alpha,k} \sum_{n_{\alpha, k}}e^{\beta(\mu-\epsilon_\alpha)n_{\alpha,k}}\,.
$$
This is a proof that, for noninteracting systems, the grandpartition function always takes the factorized form
$$
\mathcal Q = \prod_{\alpha, k} \left( \sum_{n_{\alpha,k}} e^{\beta(\mu-\epsilon_\alpha)n_{\alpha,k}}\right)\,.
$$
In the case at hand, the single-particle states all have energy $\epsilon$ and have multiplicity $N_s$; in the above notation $\alpha=1$ and $k=1,2,\ldots,N_s$ so
$$
\mathcal Q = \prod_{k=1}^{N_s}(1+e^{\beta(\mu-\epsilon)})=(1+e^{\beta(\mu-\epsilon)})^{N_s}\,.
$$
A: The questions took literally hours to be written and during the writing I may have gained a partial understanding of the problem, which I'll try to expose here.
Question

What is wrong with this approach?

Answer
The choice of the system: a fixed number $N$ of sites does not make a good gran canonical ensamble and $Z_G=\prod_i z_{G,i}$ does not apply.

On the meaning of $z_G$
As given in the solution, $z_G=1+e^{-\beta(\epsilon-\mu)}$. Formally this is precisely
$$\sum_{n=0}^1e^{-\beta(\epsilon-\mu)n}=\sum_{n=0}^1e^{-\beta \epsilon n +\beta \mu n} \sim \sum_n e^{-\beta E_n + \beta \mu n}$$
This expression resembles the one of the full grand canonical partition function and suggests the following interpretation: it describes a system

*

*whose microstate is determined by $n=0,1$

*whose energy when in the microstate $n$ is $\epsilon n$

*such that the number of elements in the system when in the microstate $n$ is precisely $n$
This system is one adsorption site, and the elements in the system are the captured particles. This number is allowed to change (between zero and one), so this is a good gran canonical ensamble. This picture may clarify the situation. On the left the array of $N$ adsorption sites is the system and one site is an element; on the right one site is the system and the captured (or not) particle is the element.

It then makes sense to write, for the system on the right side
$$z_G=\sum_{n=0}^1 e^{-\beta \epsilon n + \beta \mu n }$$
It still has to be clarified how $Z_G=z_G^N$.

Solution B
Kubo, pag. 92. He denotes by $N_s$ the number of full sites, i.e. of captured particles, and calculates the canonical partition function. What should not be done in two wrong albeit attempting ways giving the same result:
Way one
$$Z=z^N=\left(\sum_{i=0}^1 e^{-\beta \epsilon_i}\right)^N = (1+e^{-\beta \epsilon})^N$$
Way two
This time we consider as the system the occupied sites: the number $N_s$ of occupied sites can vary between $0$ and $N$ and the energy of the system in the microstate $s$ is $E_s=\epsilon N_s$, so (using Newton's formula)
$$Z=\sum_s e^{-\beta E_s} = \sum_{N_s=0}^N \underbrace{g_s}_{\text{degeneracy}} e^{-\beta \epsilon N_s} = \sum_{N_s=0}^N \binom{N}{N_s} e^{-\beta \epsilon N_s} = (1+e^{-\beta \epsilon})^N \ $$
What should be done is calculating first the canonical partition function for a given value of $N_s$ between $0$ and $N$. For this fixed value the energy of the system is always $\epsilon N_s$ with degeneracy $\binom{N}{N_s}$:
$$Z(N_s)=\sum_{\text{distributions of } N_s \text{ particles in N boxes}}e^{-\beta \epsilon N_s} = \binom{N}{N_s}e^{-\beta \epsilon N_s}$$
The grand canonical partition function follows (using Newton's formula):
$$Z_G=\sum_{N_s=0}^N Z(N_s) e^{\beta \mu N_s} = \sum_{N_s=0}^N \binom{N}{N_s} \left( e^{-\beta(\epsilon-\mu)} \right)^{N_s} = (1+e^{-\beta(\epsilon-\mu})^N $$
This could be taken as the sought proof that $Z_G=z_G^N$
A: I've found this way to understand (please someone notify me if this is considered a wrong general formula!): 
We know that:
ZG= ZG1 * ZG2 * ... 
Lets consider the case with 1 single-particle energy, but with a degeneracy of g1=2. Then the ZG can be written as:
ZG= ZG1 * ZG1 ( since we have to take into consideration every singe particle state).
So it's actually ZG= (ZG1)^g1 .
You can easily expand to r singe-particle states with degeneracy gr each, to :
ZG= ((ZG1)^g1 )* ((ZG2)^g2) * * * ((ZGr)^gr)
