fermi-dirac distribution

# 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: - $E$ is the electron energy level (eV). - $E_F$ is the Fermi energy or Fermi level (eV). - $k_B T$ is the thermal energy ($0.025852\text{ eV}$ at $300\text{ K}$). --- ## 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}$)**: - For $E < E_F$: $e^{(E - E_F)/k_B T} = e^{-\infty} = 0 \implies f_{\text{FD}}(E) = 1$. - For $E > E_F$: $e^{(E - E_F)/k_B T} = e^{+\infty} = \infty \implies f_{\text{FD}}(E) = 0$. - At $T = 0\text{ K}$, the distribution collapses into a sharp step function. All energy states up to $E_F$ are 100% filled, while all states above $E_F$ are 100% empty. 2. **Finite Temperature ($T > 0\text{ K}$)**: - Exactly at $E = E_F$: $f_{\text{FD}}(E_F) = \frac{1}{1 + e^0} = 0.5$ (50% occupation probability regardless of $T$). - Thermal excitation creates an "energy transition window" of width $\approx 4 k_B T$ centered at $E_F$. Electrons below $E_F$ are thermally excited into empty states above $E_F$. 3. **High Energy Tail ($E - E_F \ge 3 k_B T$)**: - The exponential term dominates: $e^{(E - E_F)/k_B T} \gg 1$. - The distribution reduces to the non-degenerate **Maxwell–Boltzmann approximation**: $$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: - $N_c = 2 \left( \frac{2\pi m_n^* k_B T}{h^2} \right)^{3/2}$ is the effective density of states in the conduction band ($2.86 \times 10^{19}\text{ cm}^{-3}$ for Si at $300\text{ K}$). - $F_{1/2}(\eta_c)$ is the **Complete Fermi–Dirac Integral of order 1/2**: $$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: - $A_1 = \frac{1}{\sqrt{8}} \approx 0.353553$ - $A_2 = -\left( \frac{3}{16} - \frac{\sqrt{3}}{9} \right) \approx -0.004950$ - $A_3 = 0.000148$ - $A_4 = -0.00000489$ ### 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. ```python 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.