Home Knowledge Base fermions

Fermi–Dirac Statistics: Quantum Electron Distributions, Fermi Level Kinetics, and Degenerate Semiconductor Physics

Executive Overview

Fermi–Dirac (FD) statistics governs the thermodynamic equilibrium and energy state occupation of identical half-integer spin quantum particles ($s = 1/2, 3/2, \dots$), known as fermions. In solid-state physics, electrons and holes are fermions subject to the Pauli Exclusion Principle, which dictates that no two identical fermions can occupy the exact same quantum state simultaneously. Fermi–Dirac statistics is the foundational quantum framework for semiconductor physics, governing electron and hole concentrations, Fermi level positioning ($E_F$), Quasi-Fermi level splitting ($E_{Fn}, E_{Fp}$) under optical or electrical bias, threshold voltage engineering ($V_t$) in metal-gate FinFETs/GAAFETs, contact barrier kinetics, and degenerate transport regimes in modern sub-5 nm integrated circuits. This article provides a rigorous mathematical derivation, complete analytical formulations (including Fermi–Dirac integral approximations), Python simulation models, and semiconductor device engineering applications.


Quantum Derivation & Grand Canonical Ensemble

Pauli Exclusion & Anti-Symmetric Wave Functions

For a system of $N$ identical fermions, the total quantum mechanical wave function $\Psi(\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N)$ must be strictly anti-symmetric under particle exchange:

$$\Psi(\dots, \mathbf{r}_i, \dots, \mathbf{r}_j, \dots) = -\Psi(\dots, \mathbf{r}_j, \dots, \mathbf{r}_i, \dots)$$

If two fermions occupy the same spatial and spin quantum state ($\mathbf{r}_i = \mathbf{r}_j$), then $\Psi = -\Psi \implies \Psi = 0$. Thus, the occupation number $n_i$ of any single-particle state $i$ is restricted to binary values:

$$n_i \in \{0, 1\}$$

Grand Canonical Partition Function

Consider a single-particle state $i$ with energy $\epsilon_i$ in thermal and particle equilibrium with a reservoir at temperature $T$ ($\beta = 1 / (k_B T)$) and Fermi energy $E_F$ (chemical potential $\mu = E_F$). The grand partition function $\Xi_i$ sums over allowable occupation numbers $n_i = 0$ and $n_i = 1$:

$$\Xi_i = \sum_{n_i \in \{0, 1\}} e^{-\beta n_i (\epsilon_i - E_F)} = 1 + e^{-\beta (\epsilon_i - E_F)}$$

The grand potential contribution is $\Phi_i = -k_B T \ln \Xi_i = -k_B T \ln \left( 1 + e^{-\beta (\epsilon_i - E_F)} \right)$. The mean occupation probability $f_{\text{FD}}(\epsilon_i) = \langle n_i \rangle$ is derived via partial differentiation:

$$f_{\text{FD}}(\epsilon_i) = -\frac{\partial \Phi_i}{\partial E_F} = \frac{e^{-\beta (\epsilon_i - E_F)}}{1 + e^{-\beta (\epsilon_i - E_F)}}$$

Dividing the numerator and denominator by $e^{-\beta (\epsilon_i - E_F)}$ yields the Fermi–Dirac distribution function:

$$f_{\text{FD}}(E) = \frac{1}{1 + e^{(E - E_F) / k_B T}}$$

Where:


Temperature Dependent Kinetics of $f_{\text{FD}}(E)$

The Fermi–Dirac distribution exhibits distinct behavior across temperature regimes:

1. Absolute Zero Limit ($T \to 0\text{ K}$):

2. Finite Temperature ($T > 0\text{ K}$):

3. High Energy Tail ($E - E_F \ge 3 k_B T$):

$$f_{\text{FD}}(E) \approx e^{-(E - E_F)/k_B T}$$
   f_FD(E) Occupation Probability
   1.0 |=========\  (T = 0 K Step Function)
       |          \
   0.5 |-----------\----------  at E = E_F
       |            \ (T = 300 K Thermal Smearing ~4 k_B T)
   0.0 +-------------+------------------------> Energy E
                    E_F

Carrier Densities & Fermi–Dirac Integrals

Conduction Band Electron Density $n$

In a 3D bulk semiconductor with parabolic band edge $E_c$ and effective mass $m_n^*$, the density of states is:

$$N_c(E) = \frac{1}{2\pi^2} \left( \frac{2m_n^*}{\hbar^2} \right)^{3/2} \sqrt{E - E_c} \quad (E \ge E_c)$$

The total conduction band electron concentration $n$ is calculated by integrating $N_c(E) f_{\text{FD}}(E)$:

$$n = \int_{E_c}^{\infty} N_c(E) f_{\text{FD}}(E) dE = \frac{1}{2\pi^2} \left( \frac{2m_n^*}{\hbar^2} \right)^{3/2} \int_{E_c}^{\infty} \frac{\sqrt{E - E_c}}{1 + e^{(E - E_F)/k_B T}} dE$$

Defining dimensionless variables $\eta_c = \frac{E_F - E_c}{k_B T}$ and $x = \frac{E - E_c}{k_B T}$:

$$n = N_c \cdot F_{1/2}(\eta_c)$$

Where:

$$F_{1/2}(\eta) = \frac{2}{\sqrt{\pi}} \int_{0}^{\infty} \frac{x^{1/2}}{1 + e^{x - \eta}} dx$$

Valence Band Hole Density $p$

Similarly, the hole density $p$ in the valence band (edge $E_v$, effective mass $m_p^*$) with hole occupation $1 - f_{\text{FD}}(E)$ is:

$$p = N_v \cdot F_{1/2}(\eta_v)$$

Where $\eta_v = \frac{E_v - E_F}{k_B T}$ and $N_v = 2 \left( \frac{2\pi m_p^* k_B T}{h^2} \right)^{3/2}$ ($3.10 \times 10^{19}\text{ cm}^{-3}$ for Si at $300\text{ K}$).


Analytical Approximations for $F_{1/2}(\eta)$

Because $F_{1/2}(\eta)$ cannot be solved analytically in closed form, explicit analytical approximations are required for TCAD solvers and device modeling:

1. Non-Degenerate Limit ($\eta \ll -2$, $E_c - E_F \gg 2 k_B T$)

When $E_F$ lies deep inside the bandgap ($> 2 k_B T$ below $E_c$), $e^{x - \eta} \gg 1$, yielding:

$$F_{1/2}(\eta) \approx \frac{2}{\sqrt{\pi}} \int_{0}^{\infty} x^{1/2} e^{-(x-\eta)} dx = e^{\eta}$$
$$n \approx N_c e^{\eta_c} = N_c e^{-(E_c - E_F)/k_B T}$$

2. Joyce–Dixon Approximation

To extract the Fermi level position $\eta_c$ continuously across non-degenerate and moderately degenerate regimes ($n / N_c \le 5$):

$$\eta_c = \ln\left( \frac{n}{N_c} \right) + \sum_{m=1}^{4} A_m \left( \frac{n}{N_c} \right)^m$$

Where the Joyce–Dixon coefficients are:

3. Bednarczyk–Bednarczyk / Nilsson Approximations

For ultra-high accuracy across all regimes ($\eta \in [-\infty, +\infty]$) with relative error $< 0.4\%$:

$$F_{1/2}(\eta) \approx \left[ e^{-\eta} + \frac{3\sqrt{\pi}}{4} (\eta + 2.13 + (\eta - 2.13)^2 + 9.6)^{-3/8} \right]^{-1}$$

Non-Equilibrium Quasi-Fermi Levels

Under external excitation—such as optical illumination, forward bias in a p-n junction, or high electric field transport—the electron and hole populations deviate from thermal equilibrium ($n \cdot p \ne n_i^2$).

While a single Fermi level $E_F$ is no longer defined, electrons and holes within their respective bands thermalize rapidly ($\sim 100\text{ fs}$) via intraband carrier-carrier scattering to separate quasi-equilibrium distributions characterized by Quasi-Fermi levels:

$$f_n(E) = \frac{1}{1 + e^{(E - E_{Fn})/k_B T}} \implies n = N_c F_{1/2}\left(\frac{E_{Fn} - E_c}{k_B T}\right)$$
$$f_p(E) = \frac{1}{1 + e^{(E_{Fp} - E)/k_B T}} \implies p = N_v F_{1/2}\left(\frac{E_v - E_{Fp}}{k_B T}\right)$$

The product of non-equilibrium carrier concentrations scales exponentially with the Quasi-Fermi level separation:

$$n \cdot p = n_i^2 \exp\left( \frac{E_{Fn} - E_{Fp}}{k_B T} \right)$$

This splitting $\Delta E_F = E_{Fn} - E_{Fp} = q V_a$ defines the internal electrochemical potential difference across a forward-biased junction ($V_a$).


Quantitative Python Model: Fermi Level & Occupation Solver

The following Python program computes $F_{1/2}(\eta)$, compares Maxwell–Boltzmann vs Fermi–Dirac occupation, and extracts $E_F$ across donor doping concentrations ($10^{14}$ to $10^{21}\text{ cm}^{-3}$) in silicon.

import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt

# Physical Constants
k_B = 8.617333262145e-5  # eV/K
q = 1.602176634e-19      # C
T = 300.0                # K
kBT = k_B * T            # eV (~0.02585 eV)

N_c_Si = 2.86e19         # cm^-3 (Silicon Conduction Band DOS at 300K)
N_v_Si = 3.10e19         # cm^-3 (Silicon Valence Band DOS at 300K)
E_g_Si = 1.12            # eV

def F_half_exact(eta):
    """Calculates exact Complete Fermi-Dirac Integral F_{1/2}(eta)."""
    integrand = lambda x: np.sqrt(x) / (1.0 + np.exp(x - eta))
    val, _ = quad(integrand, 0, 100)
    return (2.0 / np.sqrt(np.pi)) * val

def joyce_dixon_eta(r):
    """Joyce-Dixon approximation for eta = (E_F - E_c) / kBT from r = n / N_c."""
    A1 = 1.0 / np.sqrt(8.0)
    A2 = -(3.0/16.0 - np.sqrt(3.0)/9.0)
    A3 = 0.000148
    A4 = -0.00000489
    return np.log(r) + A1*r + A2*(r**2) + A3*(r**3) + A4*(r**4)

# Doping Sweep (N_D from 1e14 to 1e21 cm^-3)
N_D_array = np.logspace(14, 21, 100)
E_F_mb = []
E_F_jd = []

for N_D in N_D_array:
    # Maxwell-Boltzmann
    eta_mb = np.log(N_D / N_c_Si)
    E_F_mb.append(eta_mb * kBT)
    
    # Joyce-Dixon Fermi-Dirac
    r = N_D / N_c_Si
    eta_jd = joyce_dixon_eta(r)
    E_F_jd.append(eta_jd * kBT)

E_F_mb = np.array(E_F_mb)
E_F_jd = np.array(E_F_jd)

print("==================================================================")
print("FERMI-DIRAC VS MAXWELL-BOLTZMANN FERMI LEVEL POSITION (E_F - E_c)")
print("==================================================================")
test_dopings = [1e15, 1e18, 1e19, 5e19, 1e20, 5e20]
for nd in test_dopings:
    mb_val = np.log(nd / N_c_Si) * kBT
    jd_val = joyce_dixon_eta(nd / N_c_Si) * kBT
    diff = jd_val - mb_val
    print(f"N_D = {nd:8.1e} cm^-3 | MB: {mb_val:+.4f} eV | FD (JD): {jd_val:+.4f} eV | Error: {diff*1000:6.1f} meV")
print("==================================================================")

Semiconductor Engineering Applications

1. Threshold Voltage ($V_t$) Engineering in High-k Metal Gate (HKMG) FinFETs: In sub-5 nm FinFETs, the threshold voltage $V_t$ is set by adjusting the metal gate work function $\Phi_m$. Because the metal electrode's Fermi level $E_{F,m}$ determines the surface potential $\psi_s$ via $q\psi_s = \Phi_m - \chi_{\text{Si}} - (E_c - E_F)_{\text{bulk}}$, precise alignment of $E_{F,m}$ relative to the silicon conduction/valence band edges enables symmetric $V_t$ tuning for nFET and pFET devices.

2. Heavy Doping & Degenerate Source/Drain Contacts: In advanced source/drain contacts ($N_D > 10^{20}\text{ cm}^{-3}$), the Fermi level enters the conduction band ($E_F > E_c$, $\eta_c > 0$). MB statistics underestimates contact resistance $R_c$ by failing to account for Pauli blocking of incoming tunneling electrons. Fermi–Dirac statistics is mandatory for modeling field emission (tunneling) through Schottky barriers.

3. Solar Cell Open-Circuit Voltage ($V_{oc}$): The maximum achievable open-circuit voltage in silicon heterojunction solar cells is constrained by Quasi-Fermi level splitting: $$q V_{oc} = E_{Fn} - E_{Fp} = E_g - k_B T \ln\left( \frac{N_c N_v}{n \cdot p} \right)$$ Maximizing passivation reduces surface recombination, maintaining wide Quasi-Fermi level separation under solar illumination.


References

1. Joyce, W. B., & Dixon, R. W. (1977). "Analytic approximations for the Fermi energy of an ideal Fermi-Dirac gas." Applied Physics Letters, 31(5), 354–356. 2. Sze, S. M., & Ng, K. K. (2006). Physics of Semiconductor Devices (3rd ed.). John Wiley & Sons. 3. Blakemore, J. S. (1987). Semiconductor Statistics. Dover Publications. 4. Pierret, R. F. (1996). Semiconductor Device Fundamentals. Addison-Wesley.

fermi-dirac distributionfermi level kineticspauli exclusion principle fermionsfermi-dirac integral f12joyce-dixon approximationquasi-fermi levels non-equilibriumdegenerate semiconductor transport

Explore 500+ Semiconductor & AI Topics

From EUV lithography to CUDA optimization — search the full knowledge base or chat with our AI assistant.