Fermi–Dirac Distribution

Purpose

For an independent fermionic mode of one-particle energy $\epsilon_i$, the grand-canonical equilibrium occupation is

$$f_i \equiv \langle \hat n_i\rangle = \frac{1}{ e^{\beta(\epsilon_i-\mu)}+1 }, \qquad \beta = \frac{1}{k_{\mathrm B}T}.$$

The value $f_i$ is a mean occupation, not an allowed outcome of one number measurement. A complete fermionic mode has

$$n_i\in\{0,1\}, \qquad 0\leq f_i\leq1.$$

The canonical derivation and physical interpretation are at Fermi–Dirac Statistics. This card collects the forms most useful in calculations and consistency checks.

At a glance

Quantity	Formula
Dimensionless energy	$x_i=\beta(\epsilon_i-\mu)$
Mean occupation	$f_i=(e^{x_i}+1)^{-1}$
Empty-mode probability	$P_i(0)=1-f_i$
Occupied-mode probability	$P_i(1)=f_i$
One-mode partition factor	$\mathcal Z_i=1+e^{-x_i}$
Number variance	$\operatorname{Var}(n_i)=f_i(1-f_i)$
Energy derivative	$\partial f/\partial\epsilon=-\beta f(1-f)$
Chemical-potential derivative	$\partial f/\partial\mu=\beta f(1-f)$
Particle–hole identity	$f(\mu+\delta)=1-f(\mu-\delta)$
Zero-temperature limit	$f(\epsilon)\to\Theta(\mu-\epsilon)$
Dilute classical limit	$f(\epsilon)\simeq e^{-\beta(\epsilon-\mu)}$

Unless a degeneracy factor is stated explicitly, $i$ labels one complete spin-orbital or other complete one-particle mode.

One-mode probability law

For the quadratic grand-energy contribution

$$\hat K_i = (\epsilon_i-\mu)\hat n_i,$$

the two allowed number states have weights

$$P_i(n) = \frac{ e^{-\beta(\epsilon_i-\mu)n} }{ 1+e^{-\beta(\epsilon_i-\mu)} }, \qquad n=0,1.$$

Consequently,

$$P_i(1) = \frac{ e^{-\beta(\epsilon_i-\mu)} }{ 1+e^{-\beta(\epsilon_i-\mu)} } = f_i,$$

and

$$P_i(0) = \frac{1}{ 1+e^{-\beta(\epsilon_i-\mu)} } = 1-f_i.$$

Only for a binary fermionic mode does the mean occupation equal the probability of the occupied outcome. For a degenerate energy level containing $g_i$ independent modes,

$$\langle N_i\rangle = g_i f_i,$$

so the level can contain more than one fermion even though each complete mode cannot.

Equivalent forms

The following forms are algebraically identical:

$$f(\epsilon) = \frac{1}{e^{\beta(\epsilon-\mu)}+1} = \frac{ e^{-\beta(\epsilon-\mu)} }{ 1+e^{-\beta(\epsilon-\mu)} } = \frac12 \left[ 1- \tanh \left( \frac{\beta(\epsilon-\mu)}{2} \right) \right].$$

With fugacity

$$z=e^{\beta\mu},$$

one may write

$$f(\epsilon) = \frac{ ze^{-\beta\epsilon} }{ 1+ze^{-\beta\epsilon} }.$$

The first form is best for analytic work. The second avoids overflow when $\beta(\epsilon-\mu)$ is large and positive. In numerical code, use a stable logistic implementation or branch on the sign of $x=\beta(\epsilon-\mu)$ rather than evaluating $e^x$ blindly.

Symmetry and derivatives

Relative to the chemical potential,

$$f(-x) = 1-f(x), \qquad f(x) = \frac{1}{e^x+1}.$$

Equivalently,

$$f(\mu+\delta) = 1-f(\mu-\delta).$$

This identity is a property of the Fermi function. It does not imply that a material has particle–hole-symmetric bands, density of states, or interactions.

Differentiation gives

$$\frac{\partial f}{\partial\epsilon} = - \beta f(1-f), \qquad \frac{\partial f}{\partial\mu} = \beta f(1-f).$$

At fixed $\mu$,

$$\frac{\partial f}{\partial T} = \frac{\epsilon-\mu}{ k_{\mathrm B}T^2 } f(1-f).$$

The sign is physically useful: heating depletes modes below $\mu$ and fills modes above $\mu$ when the chemical potential is held fixed. At fixed particle number, however, $\mu$ generally depends on $T$, and the total derivative must include $d\mu/dT$.

Thermal window

The derivative kernel

$$w_T(\epsilon) \equiv - \frac{\partial f}{\partial\epsilon} = \frac{1}{ 4k_{\mathrm B}T } \operatorname{sech}^2 \left( \frac{\epsilon-\mu}{ 2k_{\mathrm B}T } \right)$$

is nonnegative and centered at $\epsilon=\mu$. When the spectrum may be extended far beyond the thermal window,

$$\int_{-\infty}^{\infty} w_T(\epsilon)\,d\epsilon = 1.$$

Its maximum is

$$w_T(\mu) = \frac{1}{4k_{\mathrm B}T},$$

and its full width at half maximum is

$$\Delta\epsilon_{\mathrm{FWHM}} = 4 \operatorname{arcosh} \sqrt2\, k_{\mathrm B}T \simeq 3.53\,k_{\mathrm B}T.$$

Thus low-temperature response integrals weighted by $-\partial f/\partial\epsilon$ sample an energy shell only a few $k_{\mathrm B}T$ wide around $\mu$.

Fluctuations and response

Because $\hat n_i^2=\hat n_i$,

$$\operatorname{Var}(n_i) = \langle \hat n_i^2\rangle - \langle \hat n_i\rangle^2 = f_i(1-f_i).$$

The variance is largest at $\epsilon_i=\mu$:

$$f_i=\frac12, \qquad \operatorname{Var}(n_i)=\frac14.$$

Combining the variance with the derivative identity gives the one-mode fluctuation–response relation

$$\frac{\partial f_i}{\partial\mu} = \beta \operatorname{Var}(n_i).$$

For statistically independent modes in a grand-canonical ideal gas,

$$\langle N\rangle = \sum_i f_i, \qquad \operatorname{Var}(N) = \sum_i f_i(1-f_i) = k_{\mathrm B}T \frac{\partial\langle N\rangle}{\partial\mu}.$$

The final equality uses fixed $T$ and volume. It does not apply unchanged in a finite canonical ensemble with exactly fixed $N$, where occupation numbers are correlated by the number constraint.

Zero-temperature limit

For fixed $\mu$ and $\epsilon\neq\mu$,

$$\lim_{T\to0^+} f(\epsilon) = \begin{cases} 1, & \epsilon<\mu,\\ 0, & \epsilon>\mu. \end{cases}$$

At every nonzero temperature,

$$f(\mu)=\frac12.$$

The value assigned to the limiting step function exactly at its jump does not affect continuum integrals. For an ideal Fermi gas at fixed density,

$$\lim_{T\to0^+}\mu(T) = \epsilon_{\mathrm F}.$$

The occupied zero-temperature region is the Fermi sea. A Fermi surface is the boundary in momentum space satisfying $\epsilon_{\mathbf k}=\epsilon_{\mathrm F}$, not a boundary in ordinary position space.

Dilute classical limit

When

$$\epsilon-\mu \gg k_{\mathrm B}T,$$

the mode is weakly occupied and

$$f(\epsilon) \simeq e^{-\beta(\epsilon-\mu)}.$$

More precisely, if

$$q \equiv ze^{-\beta\epsilon} <1,$$

then

$$f(\epsilon) = \frac{q}{1+q} = \sum_{\ell=1}^{\infty} (-1)^{\ell-1}q^\ell.$$

The first term is the Maxwell–Boltzmann occupation. The alternating higher terms encode fermionic corrections. Low mean occupation is the relevant criterion; high temperature by itself is not enough if the density is raised at the same time.

Sums and density-of-states integrals

For a diagonal noninteracting Hamiltonian

$$\hat H = \sum_i \epsilon_i \hat n_i, \qquad \hat N = \sum_i \hat n_i,$$

the basic grand-canonical sums are

$$\langle N\rangle = \sum_i f_i, \qquad U = \sum_i \epsilon_i f_i,$$

and

$$\Omega = - k_{\mathrm B}T \sum_i \ln \left( 1+ e^{-\beta(\epsilon_i-\mu)} \right).$$

The entropy is

$$S = - k_{\mathrm B} \sum_i \left[ f_i\ln f_i + (1-f_i)\ln(1-f_i) \right],$$

with $0\ln0\equiv0$.

If $D(\epsilon)$ denotes the total one-particle density of states, including the physical volume and every degeneracy not already included in the mode label, then

$$N = \int D(\epsilon) f(\epsilon) \,d\epsilon,$$ $$U = \int \epsilon D(\epsilon) f(\epsilon) \,d\epsilon,$$

and

$$\Omega = - k_{\mathrm B}T \int D(\epsilon) \ln \left[ 1+ e^{-\beta(\epsilon-\mu)} \right] \,d\epsilon.$$

At fixed $N$, the number equation determines $\mu(T)$. Do not insert $\mu=\epsilon_{\mathrm F}$ at nonzero temperature unless the approximation being used makes that replacement consistently.

Low-temperature expansion

For a smooth function $\phi(\epsilon)$ and a chemical potential sufficiently far from spectral edges,

$$\begin{aligned} \int_{\epsilon_0}^{\infty} \phi(\epsilon)f(\epsilon)\,d\epsilon ={}& \int_{\epsilon_0}^{\mu} \phi(\epsilon)\,d\epsilon \\ &+ \frac{\pi^2}{6} (k_{\mathrm B}T)^2 \phi'(\mu) + O(T^4). \end{aligned}$$

This is the leading Sommerfeld expansion. For fixed particle number, set $\phi=D$ in the number equation. If $\epsilon_{\mathrm F}$ is the zero-temperature chemical potential, then

$$\mu(T) = \epsilon_{\mathrm F} - \frac{\pi^2}{6} (k_{\mathrm B}T)^2 \frac{ D'(\epsilon_{\mathrm F}) }{ D(\epsilon_{\mathrm F}) } + O(T^4).$$

For a three-dimensional free gas, $D(\epsilon)\propto\sqrt{\epsilon}$, so

$$\mu(T) = \epsilon_{\mathrm F} \left[ 1 - \frac{\pi^2}{12} \left( \frac{T}{T_{\mathrm F}} \right)^2 + O \left( \frac{T^4}{T_{\mathrm F}^4} \right) \right].$$

These formulas require a smooth density of states near the chemical potential. They can fail near band edges, van Hove singularities, narrow levels, or a gap.

Worked numerical check

Suppose a mode lies two thermal energies above the chemical potential:

$$\epsilon-\mu = 2k_{\mathrm B}T.$$

Then

$$f(\epsilon) = \frac{1}{e^2+1} \simeq 0.1192.$$

By the particle–hole identity, a mode two thermal energies below $\mu$ has

$$f(\mu-2k_{\mathrm B}T) = 1-0.1192 \simeq 0.8808.$$

This pair is a quick check on signs in numerical implementations.

Assumptions and scope

The formula is exact when:

• the state is a grand-canonical equilibrium state;
• the Hamiltonian is diagonal in independent fermionic modes;
• $\mu$ couples to a conserved particle number or charge;
• each mode label includes all quantum numbers needed to make its occupation binary.

It can also be used for controlled quasiparticles with renormalized energies. In that setting, the distribution describes quasiparticle occupations and does not automatically equal the momentum distribution of bare particles.

The bare Fermi–Dirac form is generally insufficient for:

• strongly correlated states without an independent-mode description;
• superconducting or superfluid pairing, where Bogoliubov coherence factors enter;
• driven steady states and other nonequilibrium distributions;
• finite canonical systems whose exact particle number correlates modes;
• open systems not equilibrated to a bath with the stated $T$ and $\mu$.

Common mistakes

• Treating $f_i$ as a fractional eigenvalue of $\hat n_i$ rather than an ensemble mean.
• Applying the one-particle cap to an energy level while omitting its spin or other degeneracy.
• Multiplying by a degeneracy factor already included in $D(\epsilon)$.
• Assuming $\mu=\epsilon_{\mathrm F}$ at every temperature.
• Treating $\mu$ as constrained below the ground-state energy; that convergence condition belongs to an ideal Bose gas, not a finite fermionic mode.
• Replacing a discrete sum by a continuum integral when the level spacing is comparable to $k_{\mathrm B}T$.
• Using the $T\to0$ step before evaluating a quantity controlled by the thermal shell.
• Inferring particle–hole symmetry of a system from the identity $f(-x)=1-f(x)$.
• Applying grand-canonical fluctuation formulas to an exactly fixed-$N$ ensemble.
• Evaluating exponentials in a numerically unstable form.

Exercises

1. Fluctuation–response identity

Differentiate the Fermi function with respect to $\mu$ and show that the result equals $\beta\operatorname{Var}(n)$ for one ideal mode.

Solution

Let $x=\beta(\epsilon-\mu)$. Then

$$\frac{df}{dx} = - \frac{e^x}{(e^x+1)^2} = - f(1-f),$$

and $\partial x/\partial\mu=-\beta$. Therefore

$$\frac{\partial f}{\partial\mu} = \beta f(1-f).$$

Since $n^2=n$ for $n=0,1$,

$$\operatorname{Var}(n) = \langle n\rangle-\langle n\rangle^2 = f(1-f),$$

which proves the identity.

2. Locate a prescribed occupation

Find $\epsilon-\mu$ when $f=0.9$.

Solution

Invert the distribution:

$$e^{\beta(\epsilon-\mu)} = \frac{1-f}{f}.$$

For $f=0.9$,

$$\epsilon-\mu = k_{\mathrm B}T \ln \left( \frac{0.1}{0.9} \right) = - k_{\mathrm B}T\ln9.$$

The negative sign is sensible because a mode with occupation greater than one half lies below the chemical potential.

3. Width of the thermal shell

Derive the full width at half maximum of $-\partial f/\partial\epsilon$.

Solution

Writing

$$w_T(\epsilon) = \frac{1}{4k_{\mathrm B}T} \operatorname{sech}^2 y, \qquad y = \frac{\epsilon-\mu}{2k_{\mathrm B}T},$$

the half-maximum condition is

$$\operatorname{sech}^2 y = \frac12.$$

Hence

$$\cosh y = \sqrt2, \qquad \lvert y\rvert = \operatorname{arcosh}\sqrt2.$$

The separation between the two solutions is therefore

$$\Delta\epsilon_{\mathrm{FWHM}} = 4 \operatorname{arcosh}\sqrt2\, k_{\mathrm B}T \simeq 3.53\,k_{\mathrm B}T.$$

4. Fixed-number chemical-potential shift

Use the leading Sommerfeld expansion to derive the low-temperature shift of $\mu$ for a smooth density of states at fixed $N$.

Solution

At zero temperature,

$$N = \int_{\epsilon_0}^{\epsilon_{\mathrm F}} D(\epsilon)\,d\epsilon.$$

At low temperature,

$$N = \int_{\epsilon_0}^{\mu} D(\epsilon)\,d\epsilon + \frac{\pi^2}{6} (k_{\mathrm B}T)^2 D'(\mu) + O(T^4).$$

Write $\mu=\epsilon_{\mathrm F}+\delta\mu$ with $\delta\mu=O(T^2)$. Expanding to this order gives

$$0 = D(\epsilon_{\mathrm F}) \delta\mu + \frac{\pi^2}{6} (k_{\mathrm B}T)^2 D'(\epsilon_{\mathrm F}).$$

Thus

$$\delta\mu = - \frac{\pi^2}{6} (k_{\mathrm B}T)^2 \frac{ D'(\epsilon_{\mathrm F}) }{ D(\epsilon_{\mathrm F}) } + O(T^4).$$

The sign depends on the local slope of the density of states; it is negative for the three-dimensional free-particle density of states.

Canonical explanations

• Fermi–Dirac Statistics
• Occupation Numbers
• Ideal Fermi Gas
• Degenerate Fermi Gas
• Fermi Momentum and Fermi Energy
• Fermi Surface
• Sommerfeld Expansion
• Grand-Canonical Ensemble
• Fermionic Fock Space

Related lookup pages

• Fermi Gas Formula Sheet
• Quantum Gas Formula Sheet
• Bose–Einstein Distribution
• Fermi Gas Model Card
• Statistical Mechanics Checklist

References

• R. K. Pathria and P. D. Beale, Statistical Mechanics, 4th ed., Academic Press, 2021.
• M. Kardar, Statistical Physics of Particles, Cambridge University Press, 2007.
• A. L. Fetter and J. D. Walecka, Quantum Theory of Many-Particle Systems, Dover, 2003.
• N. W. Ashcroft and N. D. Mermin, Solid State Physics, Holt, Rinehart and Winston, 1976.