I am studying the distribution functions in statistical mechanics (Maxwell-Boltzmann, Bose-Einstein, and Fermi-Dirac). Are these distribution functions give the number of particles in an energy level or the probability of seeing a particle at an energy level? If they give the probability they should be less than one...
I am confused Could anyone please help me to clarify the concept?
The way in which they are derived makes it clear that they give the expected number of particles in a given energy level. If we look at the partition functions first, \begin{align} Z_F &= \sum_{n = 0, 1} e^{-n \beta (\epsilon - \mu)} = 1 + e^{-\beta (\epsilon - \mu)} \\ Z_B &= \sum_{n = 0}^\infty e^{-n \beta (\epsilon - \mu)} = (1 - e^{-\beta (\epsilon - \mu)})^{-1} \end{align} because the occupancy is $0$ or $1$ for a fermion while it is unrestricted for a boson. Partition functions are defined so that they normalize the Gibbs factor. So \begin{equation} e^{-n \beta (\epsilon - \mu)} / Z_B, \quad e^{-n \beta (\epsilon - \mu)} / Z_F \end{equation} are probabilities for any $n$ which you can check. You then use these probabilities in the formula for expected value: \begin{align} f_F(\epsilon) &= \sum_{n = 0, 1} n e^{-n\beta(\epsilon - \mu)} / Z_F = (e^{\beta(\epsilon - \mu)} + 1)^{-1} \\ f_B(\epsilon) &= \sum_{n = 0}^\infty n e^{-n\beta(\epsilon - \mu)} / Z_B = (e^{\beta(\epsilon - \mu)} - 1)^{-1}. \end{align}
Yes you are right, they cannot be probabilities as probabilities need to be less than 1 and Bose-Einstein Distribution clearly violates this. These formulae give the expected number of particles at a given energy state. In the Grand Canonical Ensemble you can have exchange of particles between different energy states and therefore, you can calculate the average number of particles per energy state.
Note that you can always talk about probabilities. Take the generic probability distribution in a Grand Canonical Ensemble:
\begin{align} p_{i} = \frac{e^{-\beta(E_{i} - \mu N_{i})}}{\mathcal{Z}} \end{align}
where the subscript $i$ refers to the microstate. We can then take our microstate to be a specific energy level. We then note that the total energy of this microstate will be precisely the number of particles occupying that state! We then arrive at the following formula:
\begin{align} p_{i} = \frac{e^{-\beta N_i(\epsilon_{i} - \mu)}}{\mathcal{Z}} \end{align}
where $\epsilon_i$ refers to the energy level of that state. This gives the probability distribution of finding $N_i$ particles at a state which has energy $E_{i}$. You can then use this formula to determine the average number of particles at a given energy state and derive the different types of distribution as Connor explained. Hope this helps.
To add to the answer by @ConnorBehan: since for non-interacting Fermions the state is either occupied or not, one could treat occupation number as distributed according to Bernoulli distribution, whose mean is equal to the occupation probability. Then Fermi function is both the mean/average and the probability.
This does not hold for the bosons. It also does not exactly hold when we try to include interactions, but scattering processes would be still often treated as if the initial means occupations were independent probabilities of states being occupied or empty. Thus, the probability of a Coulomb scattering process $\mathbf{k},\mathbf{k}'\longrightarrow\mathbf{k}+\mathbf{q}, \mathbf{k}'-\mathbf{q}$ would be given by a Fermi Golden rule expression like $$ w_{\mathbf{k},\mathbf{k}'\longrightarrow\mathbf{k}+\mathbf{q}, \mathbf{k}'-\mathbf{q}}= \frac{2\pi}{\hbar}\left|V(\mathbf{k},\mathbf{k}';\mathbf{k}+\mathbf{q}, \mathbf{k}'-\mathbf{q})\right|^2 f(\epsilon_\mathbf{k})f(\epsilon_{\mathbf{k}'}) \left[1-f(\epsilon_{\mathbf{k}+\mathbf{q}})\right]\left[1-f(\epsilon_{\mathbf{k}'-\mathbf{q}})\right]\times\\\delta\left(\epsilon_\mathbf{k}+\epsilon_{\mathbf{k}'}-\epsilon_{\mathbf{k}+\mathbf{q}}-\epsilon_{\mathbf{k}'-\mathbf{q}})\right) $$
