Fermi-Dirac statistics

Have you ever wondered how Fermi-Dirac statistics arises? I may have wondered during my bachelor degree studies, but now I only remember derivation done by the means of combinatorics. In this post let us allow Fermi-Dirac distribution to emerge from simulation of a simple particle system.


  1. Let the system be composed of \( 100 \) energy levels.
  2. Each energy level has energy equal to its index. In other words \( E_i = i \) with \( i = 1, 2, \ldots, 100 \).
  3. Let there be \( 50 \) particles jumping from one energy level to another.
  4. Particle may jump to an energy level if it is empty. Namely, single particle can occupy single energy level.
  5. Probability of the jump from energy level \( i \) to energy level \( j \) is given by the Boltzmann factor:

\begin{equation} P_{i \rightarrow j} = \min\left[1, \exp\left( \frac{E_i - E_j}{k_B T} \right) \right]. \end{equation}

This is a rather simple model, which effectively mirrors the implementation of the Ising model. Boltzmann factor being equivalent to the acceptance ratio in the Metropolis algorithm for the Ising model.

Interactive app

All is fine and interesting, but we would like to have average occupation of the energy levels to follow Fermi-Dirac distribution function:

\begin{equation} F\left(E_i\right) = \frac{1}{1+\exp\left(\frac{E_i - E_F}{k_B T}\right)} . \end{equation}

In the above \( E_F \) is Fermi energy, which marks the inflection point of the distribution. In our simulation it will be at \( E_i = 50 \). Will it work? Check it yourself using the app below (black curve represents \( F\left(E\right) \) from the above).