Relation between Canonical and Grand Canonical Ensemble for a two levels system For a two-level system of classical non interacting distinguishable paricles, labelling energy levels as $+\epsilon, -\epsilon$, in the canonical ensemble I found for the partition function Q$=2\cosh(\beta\epsilon)$ and so the Gibbs free energy $A=-kTN\ln(2\cosh(\beta\epsilon))$ , $S=-(\frac{\partial A}{\partial T})_{N,V}$ and $\mu=(\frac{\partial A}{\partial N})_{V,T}$ so the fugacity $z=$Q$^{-1}$
But in the grand canonical frame I have $Q=\sum_{N=0}^\infty(z$Q$)^N$ which i can sum iff $z$Q$<1$. After introducing the $q$ potential I used the expression for the entropy in the GC formalism obtaining: $$S_{GC}=\frac{\epsilon}{T}\frac{2z\sinh(\beta\epsilon)}{1-2z\cosh(\beta\epsilon)}-k\ln(1-2z\cosh(\beta\epsilon))-Nk\ln(z)$$to compare with: $$S_C=-(\frac{\partial A}{\partial T})_{N,V}=Nk(\ln(2\cosh(\beta\epsilon))-\beta\tanh(\beta\epsilon))$$ In which way I can employ the thermodynamic limit to establish the equivalence between these two expressions?
 A: The way it usually works is by showing that the relative fluctuations of the variable that becomes a random variable when changing ensemble ($N$ in this case) vanish for large system sizes. This then shows that to a single value of the chemical potential, one can assign in very good approximation a single value to the random variable $N$ that corresponds to the grand canonical average value.
To to this the following relations are helpful:
\begin{equation}
\langle N\rangle_{\mu, \beta} = z\left(\frac{\partial \ln Q}{\partial z}\right)_{\beta}; \: Var(N) = \langle N^2\rangle_{\mu, \beta}-\langle N\rangle^2_{\mu, \beta} = z\left(\frac{\partial \langle N\rangle_{\mu, \beta}}{\partial z}\right)_{\beta}
\end{equation}
Now, if indeed $zQ < 1$ (${\rm Q}$ being the single particle partition function) then we have $Q = (1-{\rm Q}z)^{-1}$.
It follows that $\ln Q = -\ln(1-{\rm Q}z)$:
\begin{equation}
\langle N\rangle_{\mu, \beta} = \frac{z{\rm Q}}{1-{\rm Q}z}; \: Var(N) = z\frac{{\rm Q}(1-z{\rm Q})+z{\rm Q}^2}{(1-z{\rm Q})^2} = \frac{\langle N\rangle_{\mu, \beta}}{(1-z{\rm Q})}
\end{equation}
It follows that:
\begin{equation}
\frac{\sqrt{Var(N)}}{\langle N\rangle} = \frac{1}{\sqrt{(1-z{\rm Q})\langle N \rangle}}
\end{equation}
The problem is that for $zQ < 1$, it is always the case that $1/(1-z{\rm Q}) > \langle N \rangle$ so that $(1-z{\rm Q})\langle N \rangle < 1$.
Overall, the relative fluctuations will never vanish no matter the system size for the system you are considering. There is no equivalence between the ensembles.
This is understandable because a necessary condition for ensemble equivalence to appear is that the overall probability of the random variable (here $N$) has a bell-shape and one then observes how the waist of this curve varies as the system size changes.
For the system you are looking into, the probability function is always decreasing (exponentially so) such that no particular value is preferred in the overall sum making the grand partition function.
To get a system very close to yours but with equivalence of the grand canonical and canonical ensembles you need to add a little, but very important, ingredient: indistinguishability of the particles. This will put a $1/N!$ in each term of the sum and yield a grand partition function that is $Q = e^{z{\rm Q}}$ (no constraint anymore for the convergence of the sum by the way). Its logarithm is of course $\ln Q = z{\rm Q}$. We then get that $\langle N \rangle = z{\rm Q}$ and $Var(N) = z{\rm Q}$. Which gives directly 
\begin{equation}
\frac{\sqrt{Var(N)}}{\langle N\rangle} = \frac{1}{\sqrt{\langle N \rangle}}
\end{equation}
In that case, there exists a most probable value of $N$ that ends up dominating the sum in the grand partition function so that $Q \approx (z{\rm Q})^{\langle N \rangle}$ as the system tends to the thermodynamic limit.
