Confusion about using single particle or $N$ particle partition function in Boltzmann probability in canonical ensemble Suppose we have a canonical ensemble, where $N$ particles have been divided among $\epsilon_i$ energy levels, each with degeneracy $g_i$. The partition function for a single particle is given by : $$Z_{sp}=\sum_{i}^{r} g_i e^{-\beta\epsilon_i}$$
There are $r$ total energy levels here.
Anyway, the partition function for all the $N$ particles can be found using the same formula, by checking all the possible combinations and values of total energy due to each of these particles, which would be a long and tedious process. However, we can write the $N$ particle partition function as follows :
$$Z_{N}=\prod_{i}^{N}(Z_{sp})_i$$
Now we know that the probability of the system being in a particular 'energy level' is :
$$P(\epsilon_i)=\frac{g_i e^{-\beta \epsilon_i}}{Z}$$
If we cared about a particular state, and not the energy level, then we would just drop the $g_i$ term in the numerator I suppose. My question is, what exactly is $Z$ in this example? Is it $Z_{sp}$ or $Z_N$ ? According to an example in my book, it should be the single-particle partition function. However my system consists of $N$ particles, so shouldn't we consider the $N$ particle partition function instead ?
Suppose my particles are bosons or classical particles, in the sense that there is no limitation on the number of particles in a state, we can say the following :
$$P(\epsilon_i)= \frac{n_i}{N}$$
Hence, we are finding the number of particles with energy $\epsilon_i$ divided by the total number of particles. Hence we can write :
$$n_i=N\frac{g_i e^{-\beta \epsilon_i}}{Z}$$
But now, shouldn't we consider the partition function of all these $N$ particles?
According to my book, the number of particles in a particular energy level is the product of the total number of particles and the probability of a single particle in that level. Because of this, they use the single-particle partition function. This seems a little wrong to me. Since we are talking about $N$ particles, shouldn't we just use the $N$ particle partition function instead
? The probability of the system in a particular energy level is given using the $N$ particle partition function, since there are $N$ particles in the system. So, by that argument, shouldn't number of particles in a particular energy levels also be given using the $N$ particle.
So, if I want to find the number of particles in a particular energy level, in a system of $N$ particles, what should I use :
$$n_i=N\frac{g_i e^{-\beta \epsilon_i}}{Z_{sp}} \space \space or\space\space \frac{g_i e^{-\beta \epsilon_i}}{Z_{N}}$$
Any help on understanding this concept would be highly appreciated. Thanks !
 A: When we go to $N$ particles, degeneracy factors are produced automatically anyway. So I am going to dispense with the headache of having degeneracy factors already with one particle. So just use $\epsilon_1, \dots, \epsilon_R$. If $R > r$, that means some of the $\epsilon_i$ are the same.
The single particle system obeys
\begin{equation}
Z_{\mathrm{sp}} = \sum_{i = 1}^R e^{-\beta \epsilon_i}, \quad P(\epsilon_i) = \frac{e^{-\beta \epsilon_i}}{Z_{\mathrm{sp}}}
\end{equation}
so now we have to talk about how the multi-particle system is defined. If it is defined as $N$ non-interacting copies of the single particle system where there is no restriction on levels being filled multiple times (i.e. the particles are already known to be distinguishable), then and only then do we get to say
\begin{equation}
Z_N = \prod_{i = 1}^N Z_{\mathrm{sp}} = Z_{\mathrm{sp}}^N.
\end{equation}
So now what is the probability of a given energy? Well we have to remember that almost all of the possible values for the energy of the multi-particle system are not elements of $(\epsilon_1, \dots \epsilon_R)$. Instead, we can consider the microstates labelled by $(n_1, \dots, n_R)$ such that $n_1 + \dots + n_R = N$. Then, the probability of one of them being realized is
\begin{equation}
P(n_1, \dots, n_R) = \frac{e^{-\beta(\epsilon_1 n_1 + \dots + \epsilon_R n_R)}}{Z_N}.
\end{equation}
Is this also the probability that an energy of $E = n_1 \epsilon_1 + \dots + n_R \epsilon_R$ will be measured? Not necessarily if the $\epsilon_i$ are regularly spaced. To get $P(E)$, we need to sum the above over all possible sets of $n_i$ such that $n_1 \epsilon_1 + \dots + n_R \epsilon_R = E$.
We finally come to your equation $P(\epsilon_i) = n_i / N$. I think that is either wrong or reliant on some confusing definitions. As with the total energy, the number of particles that have energy $\epsilon_i$ is random because it depends on the microstate. Its expectation will be given by
\begin{equation}
N_i = Z_N^{-1} \sum_{n_1 + \dots + n_R = N} n_i \binom{N}{n_1, \dots, n_R} e^{-\beta(n_1 \epsilon_1 + \dots + n_R \epsilon_R)}.
\end{equation}
A: Assume you have a 2-level atom, with the energy of the lower level 0 (electron in ground state) and the energy of the upper level E (electron in excited state). In thermodynamic equilibrium at temperature T, the probability of the atom to be in the ground state is according to the Boltzmann distribution
$$p_1=\frac{1}{1+e^{-E/(kT)}}$$
and for the upper state
$$p_2=\frac{e^{-E/(kT)}}{1+e^{-E/(kT)}}$$
Note that $p_1+p_2 =1$.
This means if you have a large number N of atoms, you will find about $N_1=N\cdot p_1$ in the lower state and $N_2=N\cdot p_2$ in the upper state. If the upper state is degenerate, $N_2$ will be shared between all the degenerate states, correspondingly for $N_1$ if the lower state is degenerate (see also my first answer to this SE question in this context).
Having said this, the condition for this is the assumption of thermodynamic equilibrium to hold, which implies that the levels are populated and de-populated by collisions only, i.e. the collisional time scale is assumed to be much shorter than any other time scales. In most practical cases this is actually not the case. The quantum mechanical decay times for dipole-allowed transitions are practically always muvh shorter than the collisional time scale, so calculations based on the Boltzmann distribution will be badly wrong in this case.  Only for sufficiently highly excited states or dipole-forbidden transitions do collisions become dominant.
So you always have to be careful with applying results based on the theory of Thermodynamics. It is dangerous taking some equations from textbooks before having analyzed the physical problem at hand in detail and made sure that they are actually applicable to the problem.
