What is wrong with these ways of determining the mean occupation number? Could anyone point out what went wrong in this argument?
Setup:
We have a system with 2 energy levels say with energies $0,e$ respectively.
We consider the grand canonical ensemble for the system occupied by identical spinless fermions. Then there are clearly 4 feasible microstates -- 


*

*no particles, $E=0$;  

*1 particle in $E=0$ level, $E=0$;

*1 particle in $E=e$ level, $E=e$;

*2 particles one in each level, $E=e$


Clearly the corresponding grand canonical function is $$\Xi = 1+f+f\exp(-\beta e)+f^2\exp(-\beta e) $$ where $f\equiv \exp(\mu\beta)$.
Now I wish to find the mean occupation number of a state with energy $e$.
Here is where my arguments go all wrong, please help me: 

Method 1:
$$n(e)=1\times \text{probability of state (3)}+2\times \text{probability of state (4)}\\
= {f\exp(-\beta e)+2f^2\exp(-\beta e)\over \Xi}\\
={1+2f\over 1+f}\cdot{1\over f^{-1}\exp(\beta e)+1}$$
BUT...

Method 2:
If I use the F-D distribution directly, shouldn't I get $$n(e) = {\text{degeneracy of energy level $E=e$}\over f^{-1}\exp(\beta e)+1}\\= 2\cdot{1 \over f^{-1}\exp(\beta e)+1}$$??
What went wrong?
Thanks a lot in advance!!
Christie 
 A: In your Method 1, you forgot that the state (2) also has the occupation number $n=1$ (it has a vanishing energy but you're computing the mean occupation number, not mean energy), so it must be added to the numerator. This changes $1+2f$ in the numerator on the following line to $2+2f$. The simple fraction simplifies to $2$ and you get an agreement with Method 2.
A: There is a difference between energy levels (sometimes called states) and microstates (sometimes also called states).
You want to find how many particles are in energy level $e$, given a system with energy levels $0$ and $e$. Okay. Now, you did the right thing by considering the four distinguishable microstates, assigning likelihoods to them based on temperature, chemical potential, and the energies of the levels within each state.
But, when adding contributions in the numerator of Method 1, you are only supposed to count each microstate with a multiplicity given by how much it populates the energy level of interest. Microstates 1 and 2 have no particle in energy level $e$, so they are excluded. Microstate 3 has a single particle in energy level $e$, so it is counted once. Microstate 4 also has a single particle in energy level $e$, so it should only be counted once. The fact that it has another particle in the ground state doesn't matter - maybe it also describes a boson at energy level $e'$ on the far side of the Moon - that doesn't contribute to the counting of occupations of energy level $e$ by this species of fermion.
The only difference that ground state fermion in microstate 4 makes is in changing the chemical potential of the state, which you already accounted for with the $f^2$. It could also conceivably alter the Boltzmann factor for the state, but we have a normalization where the ground state has a Boltzmann factor of unity, so there is nothing more to consider.
So remove the 2 in method 1.
For Method 2, why the degeneracy factor? That would only be there if there were 2 distinct energy levels $e$. But your fermions are spinless - there is no other quantum degree of freedom for them to differ, so only 1 at most can ever occupy energy level $e$ in any microstate. In fact, you made use of this by not considering all the microstates that involve 2 particles in level $e$. So remove that 2 as well.
