← Part I/Free Fermi Gas

1. Free Fermi Gas

Reading time: ~50 minutes | Difficulty: Graduate | Prerequisites: Statistical mechanics, quantum mechanics

Introduction

The free Fermi gas is the foundational model of condensed matter physics. It describes a system of non-interacting fermions — particles obeying the Pauli exclusion principle — confined to a volume $V$. Despite its simplicity, this model correctly predicts the linear electronic specific heat, Pauli paramagnetism, and the qualitative behavior of metals at low temperatures.

The classical Drude model, which treats conduction electrons as a classical ideal gas, fails spectacularly in predicting the electronic contribution to specific heat: it overestimates it by roughly two orders of magnitude at room temperature. The resolution came from Sommerfeld, who applied Fermi-Dirac statistics to the electron gas.

In this chapter, we derive the key properties of the free Fermi gas from first principles, building the foundation for Landau's Fermi liquid theory in subsequent chapters.

Fermi-Dirac Distribution

For a system of non-interacting fermions in thermal equilibrium at temperature $T$ and chemical potential $\mu$, the mean occupation number of a single-particle state with energy $\epsilon$ is given by the Fermi-Dirac distribution function:

$$f(\epsilon) = \frac{1}{e^{(\epsilon - \mu)/k_B T} + 1}$$

This result follows from the grand canonical ensemble. Consider a single-particle state $|\alpha\rangle$ with energy $\epsilon_\alpha$. Since fermions obey the exclusion principle, the occupation number $n_\alpha$ can only take values 0 or 1. The grand partition function for this state is:

$$\mathcal{Z}_\alpha = \sum_{n_\alpha=0}^{1} e^{-\beta(\epsilon_\alpha - \mu)n_\alpha} = 1 + e^{-\beta(\epsilon_\alpha - \mu)}$$

where $\beta = 1/(k_B T)$. The mean occupation number is then:

$$\langle n_\alpha \rangle = -\frac{1}{\beta}\frac{\partial \ln \mathcal{Z}_\alpha}{\partial \epsilon_\alpha} = \frac{e^{-\beta(\epsilon_\alpha - \mu)}}{1 + e^{-\beta(\epsilon_\alpha - \mu)}} = \frac{1}{e^{\beta(\epsilon_\alpha - \mu)} + 1}$$

Key Properties

At $T = 0$: the distribution is a sharp step function. All states with $\epsilon < \mu$ are occupied ($f = 1$), and all states with $\epsilon > \mu$ are empty ($f = 0$).
At $\epsilon = \mu$: the occupation is always $f(\mu) = 1/2$, regardless of temperature.
The distribution is broadened over a width of approximately $4k_BT$ around the chemical potential. Only electrons within this energy window of $\mu$ are thermally active.
In the limit $\epsilon - \mu \gg k_BT$, the Fermi-Dirac distribution reduces to the classical Maxwell-Boltzmann distribution: $f(\epsilon) \approx e^{-(\epsilon - \mu)/k_BT}$.

Detailed Derivation: The Sommerfeld Expansion

Step-by-step derivation of the Sommerfeld expansion

We wish to evaluate integrals of the form $I = \int_0^\infty H(\epsilon)\, f(\epsilon)\, d\epsilon$ where $f(\epsilon)$ is the Fermi-Dirac distribution. Define $\mathcal{H}(\epsilon) = \int_0^\epsilon H(\epsilon')\,d\epsilon'$.

Step 1: Integrate by parts.

$$I = \left[\mathcal{H}(\epsilon)\,f(\epsilon)\right]_0^\infty - \int_0^\infty \mathcal{H}(\epsilon)\,\frac{\partial f}{\partial \epsilon}\,d\epsilon$$

The boundary term vanishes: at $\epsilon = 0$, $\mathcal{H}(0) = 0$; at $\epsilon \to \infty$, $f \to 0$ exponentially. So:

$$I = -\int_0^\infty \mathcal{H}(\epsilon)\,\frac{\partial f}{\partial \epsilon}\,d\epsilon = \int_0^\infty \mathcal{H}(\epsilon)\,\left(-\frac{\partial f}{\partial \epsilon}\right)\,d\epsilon$$

Step 2: Note that $-\partial f/\partial \epsilon$ is sharply peaked around $\epsilon = \mu$ with width $\sim k_BT$. Substitute $x = (\epsilon - \mu)/(k_BT)$:

$$-\frac{\partial f}{\partial \epsilon} = \frac{1}{k_BT}\frac{e^x}{(e^x + 1)^2}$$

Step 3: Taylor expand $\mathcal{H}(\epsilon)$ around $\epsilon = \mu$:

$$\mathcal{H}(\mu + x\,k_BT) = \mathcal{H}(\mu) + x\,k_BT\,\mathcal{H}'(\mu) + \frac{(x\,k_BT)^2}{2!}\mathcal{H}''(\mu) + \cdots$$

Step 4: The integral becomes a sum of moments of the peaked function. Since $e^x/(e^x+1)^2$ is an even function of $x$, odd moments vanish. The even moments are:

$$\int_{-\infty}^{\infty} \frac{x^{2n}\,e^x}{(e^x+1)^2}\,dx = \frac{2(2^{2n-1}-1)\pi^{2n}|B_{2n}|}{(2n)!}$$

For $n = 1$: the integral gives $\pi^2/3$. For $n = 2$: it gives $7\pi^4/15$.

Step 5: Assembling the result with $\mathcal{H}'(\mu) = H(\mu)$ and $\mathcal{H}''(\mu) = H'(\mu)$:

$$\boxed{I = \int_0^\mu H(\epsilon)\,d\epsilon + \frac{\pi^2}{6}(k_BT)^2 H'(\mu) + \frac{7\pi^4}{360}(k_BT)^4 H'''(\mu) + \cdots}$$

Derivation of the chemical potential shift $\mu(T)$

Apply the Sommerfeld expansion to the number equation $N = \int_0^\infty g(\epsilon)\,f(\epsilon)\,d\epsilon$, with $H(\epsilon) = g(\epsilon)$:

$$N = \int_0^\mu g(\epsilon)\,d\epsilon + \frac{\pi^2}{6}(k_BT)^2 g'(\mu) + \cdots$$

At $T = 0$, $N = \int_0^{E_F} g(\epsilon)\,d\epsilon$. For the 3D density of states $g(\epsilon) \propto \sqrt{\epsilon}$, we have $g'(\epsilon) = g(\epsilon)/(2\epsilon)$. Write $\mu = E_F + \delta\mu$ and expand to first order in $\delta\mu$:

$$\int_0^\mu g(\epsilon)\,d\epsilon \approx \int_0^{E_F} g(\epsilon)\,d\epsilon + g(E_F)\,\delta\mu = N + g(E_F)\,\delta\mu$$

Setting the total equal to $N$ and solving for $\delta\mu$:

$$g(E_F)\,\delta\mu + \frac{\pi^2}{6}(k_BT)^2 g'(E_F) = 0$$

$$\delta\mu = -\frac{\pi^2}{6}(k_BT)^2 \frac{g'(E_F)}{g(E_F)} = -\frac{\pi^2}{12}\frac{(k_BT)^2}{E_F}$$

Thus $\mu(T) = E_F\left[1 - \frac{\pi^2}{12}\left(\frac{k_BT}{E_F}\right)^2 + \cdots\right]$, confirming that the chemical potential decreases quadratically with temperature.

Derivation of the electronic specific heat $C_V = \gamma T$

Apply the Sommerfeld expansion to the energy $E = \int_0^\infty \epsilon\,g(\epsilon)\,f(\epsilon)\,d\epsilon$, with $H(\epsilon) = \epsilon\,g(\epsilon)$:

$$E(T) = \int_0^\mu \epsilon\,g(\epsilon)\,d\epsilon + \frac{\pi^2}{6}(k_BT)^2\,\frac{d}{d\epsilon}\left[\epsilon\,g(\epsilon)\right]_{\epsilon=\mu} + \cdots$$

Evaluate the derivative: $\frac{d}{d\epsilon}[\epsilon\,g(\epsilon)] = g(\epsilon) + \epsilon\,g'(\epsilon)$. At $\epsilon = E_F$:

$$\frac{d}{d\epsilon}[\epsilon\,g(\epsilon)]_{E_F} = g(E_F) + E_F\,g'(E_F) = g(E_F) + \frac{g(E_F)}{2} = \frac{3}{2}g(E_F)$$

The first term gives $E_0 + \mu\,g(E_F)\,\delta\mu$. Including the $\delta\mu$ correction (which cancels a corresponding term from the number constraint), the net thermal energy is:

$$E(T) - E_0 = \frac{\pi^2}{6}(k_BT)^2\,g(E_F) + \mathcal{O}(T^4)$$

Differentiating with respect to $T$:

$$\boxed{C_V = \frac{\partial E}{\partial T} = \frac{\pi^2}{3}k_B^2\,g(E_F)\,T = \gamma\,T}$$

Historical Context

The free electron model of metals has a rich history spanning the transition from classical to quantum physics. Paul Drude proposed the first electron theory of metals in 1900, treating conduction electrons as a classical ideal gas. While the Drude model successfully explained Ohm's law and the Wiedemann-Franz law, it predicted an electronic specific heat of $\frac{3}{2}Nk_B$ — roughly 100 times larger than observed.

The resolution came in 1927–1928 when Arnold Sommerfeld applied the newly discovered Fermi-Dirac statistics (independently derived by Enrico Fermi in 1926 and Paul Dirac in 1926) to the electron gas. Sommerfeld showed that at temperatures far below the Fermi temperature ($T_F \sim 10^4$ K for typical metals), only a fraction $\sim T/T_F$ of electrons are thermally active, resolving the specific heat puzzle. The Sommerfeld expansion technique, published in his landmark 1928 paper in the Zeitschrift für Physik, remains a standard tool in condensed matter physics.

Wolfgang Pauli independently applied quantum statistics to explain the weak, temperature-independent paramagnetism of metals (Pauli paramagnetism, 1927), which had been an equally puzzling departure from the Curie law expected classically. Together, the Sommerfeld model and Pauli's result established that quantum statistics fundamentally governs the behavior of electrons in metals, laying the groundwork for the band theory of solids developed by Felix Bloch (1929) and the eventual Fermi liquid theory of Lev Landau (1956).

Applications

Real-World Applications of Free Fermi Gas Theory

1. Low-temperature calorimetry of metals. The linear electronic specific heat $C_V = \gamma T$ is routinely measured in metals at cryogenic temperatures. The Sommerfeld coefficient $\gamma$ is extracted from plots of $C/T$ vs $T^2$ (separating electronic and phonon contributions). Deviations of $\gamma$ from the free-electron value reveal the effective mass enhancement due to electron-electron and electron-phonon interactions.

2. White dwarf stars. The electron degeneracy pressure $P_0 = \frac{2}{5}nE_F$ supports white dwarf stars against gravitational collapse. Chandrasekhar's calculation (1930) of the maximum mass of a white dwarf ($\sim 1.4 M_\odot$) is based directly on the relativistic generalization of the free Fermi gas. This remains one of the most important astrophysical applications of quantum statistics.

3. Thermoelectric devices. The Wiedemann-Franz law, which relates the thermal and electrical conductivities via $\kappa / (\sigma T) = L_0 = \pi^2 k_B^2/(3e^2)$, is a direct consequence of the Sommerfeld model. This universal Lorenz number $L_0$ is used to characterize and design thermoelectric materials for power generation and refrigeration applications.

4. Metallic nanoparticles and quantum dots. In nanoscale metallic structures, the discrete level spacing $\delta \sim E_F / N$ becomes comparable to $k_BT$, and the continuum density of states must be replaced by discrete levels. The crossover between bulk Fermi gas behavior and single-level physics is essential for understanding Coulomb blockade, single-electron transistors, and catalytic properties of metal nanoparticles.

5. Neutron stars and nuclear matter. The free Fermi gas model, applied to neutrons, provides the starting point for understanding the equation of state of dense nuclear matter in neutron star interiors. The neutron degeneracy pressure supports neutron stars up to $\sim 2\text{--}3 M_\odot$, analogously to electron degeneracy in white dwarfs.

Derivation 2: The Sommerfeld Expansion — Full Detail

The Sommerfeld expansion is the central analytical tool for low-temperature properties of Fermi systems. Here we present the complete derivation with all intermediate steps, starting from the generic Fermi integral and arriving at the systematic asymptotic series.

Setup: The Generic Fermi Integral

We consider integrals of the form:

$$I = \int_0^\infty H(\epsilon)\,f(\epsilon)\,d\epsilon$$

where $H(\epsilon)$ is a smooth function (typically the density of states or $\epsilon \cdot g(\epsilon)$) and $f(\epsilon) = 1/(e^{\beta(\epsilon - \mu)} + 1)$ is the Fermi-Dirac distribution with $\beta = 1/(k_BT)$. Define the antiderivative $\mathcal{H}(\epsilon) = \int_0^\epsilon H(\epsilon')\,d\epsilon'$.

Step 1: Integration by Parts

Write $H(\epsilon) = \mathcal{H}'(\epsilon)$ and integrate by parts:

$$I = \bigl[\mathcal{H}(\epsilon)\,f(\epsilon)\bigr]_0^\infty \;-\; \int_0^\infty \mathcal{H}(\epsilon)\,f'(\epsilon)\,d\epsilon$$

The boundary term vanishes because $\mathcal{H}(0) = 0$ and $f(\epsilon) \to 0$ exponentially as $\epsilon \to \infty$. The derivative $f'(\epsilon) = \partial f / \partial \epsilon$ is:

$$f'(\epsilon) = -\frac{\beta\,e^{\beta(\epsilon-\mu)}}{\bigl(e^{\beta(\epsilon-\mu)}+1\bigr)^2}$$

So the integral becomes:

$$I = \int_0^\infty \mathcal{H}(\epsilon)\,\left(-f'(\epsilon)\right)\,d\epsilon$$

The function $-f'(\epsilon)$ is a sharply peaked, bell-shaped function centered at $\epsilon = \mu$ with characteristic width $\sim k_BT$. It is normalized: $\int_0^\infty (-f')\,d\epsilon = f(0) - f(\infty) \approx 1$ for $\mu \gg k_BT$.

Step 2: Taylor Expansion around $\epsilon = \mu$

Since $-f'(\epsilon)$ is peaked at $\mu$, expand $\mathcal{H}(\epsilon)$ in a Taylor series about $\epsilon = \mu$:

$$\mathcal{H}(\epsilon) = \sum_{n=0}^{\infty} \frac{(\epsilon - \mu)^n}{n!}\,\mathcal{H}^{(n)}(\mu)$$

where $\mathcal{H}^{(n)}(\mu) = \left.\frac{d^n \mathcal{H}}{d\epsilon^n}\right|_{\epsilon=\mu}$. Note that $\mathcal{H}^{(0)}(\mu) = \int_0^\mu H(\epsilon)\,d\epsilon$, $\mathcal{H}^{(1)}(\mu) = H(\mu)$, $\mathcal{H}^{(2)}(\mu) = H'(\mu)$, etc.

Step 3: Substitution $x = \beta(\epsilon - \mu)$

Let $x = (\epsilon - \mu)/(k_BT)$, so $d\epsilon = k_BT\,dx$. The kernel transforms as:

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

The lower limit $\epsilon = 0$ maps to $x = -\mu/(k_BT)$. Since $\mu \gg k_BT$ in the degenerate regime, we extend this to $-\infty$ with exponentially small error $\sim e^{-\mu/(k_BT)}$. The integral becomes:

$$I = \sum_{n=0}^{\infty} \frac{(k_BT)^n}{n!}\,\mathcal{H}^{(n)}(\mu) \int_{-\infty}^{\infty} \frac{x^n\,e^x}{(e^x+1)^2}\,dx$$

Step 4: Evaluation of the Moments

Define $\phi(x) = e^x/(e^x+1)^2$. This function is even: $\phi(-x) = \phi(x)$. Therefore:

Odd moments vanish: $\int_{-\infty}^{\infty} x^{2n+1}\,\phi(x)\,dx = 0$ for all $n \geq 0$.
Zeroth moment: $\int_{-\infty}^{\infty} \phi(x)\,dx = 1$.
Second moment: $\int_{-\infty}^{\infty} x^2\,\phi(x)\,dx = \pi^2/3$.
Fourth moment: $\int_{-\infty}^{\infty} x^4\,\phi(x)\,dx = 7\pi^4/15$.

The general even moment is given in terms of Bernoulli numbers $B_{2n}$:

$$\int_{-\infty}^{\infty} x^{2n}\,\phi(x)\,dx = \frac{2(2^{2n-1}-1)\,\pi^{2n}\,|B_{2n}|}{(2n)!}$$

To prove the second moment, note that $\phi(x) = -d/dx\,[1/(e^x+1)]$. Integrating $x^2\,\phi(x)$ by parts twice relates it to $\int_0^\infty x/(e^x+1)\,dx = \pi^2/12$, yielding $\pi^2/3$.

Step 5: Final Result

Collecting the surviving terms ($n = 0$ from the zeroth moment, $n = 2$ from the second moment, $n = 4$ from the fourth moment):

$$I = \mathcal{H}(\mu)\cdot 1 + \frac{(k_BT)^2}{2!}\,\mathcal{H}''(\mu)\cdot\frac{\pi^2}{3} + \frac{(k_BT)^4}{4!}\,\mathcal{H}^{(4)}(\mu)\cdot\frac{7\pi^4}{15} + \cdots$$

Recalling $\mathcal{H}(\mu) = \int_0^\mu H(\epsilon)\,d\epsilon$, $\mathcal{H}''(\mu) = H'(\mu)$, and $\mathcal{H}^{(4)}(\mu) = H'''(\mu)$:

$$\boxed{I = \int_0^\mu H(\epsilon)\,d\epsilon \;+\; \frac{\pi^2}{6}(k_BT)^2\,H'(\mu) \;+\; \frac{7\pi^4}{360}(k_BT)^4\,H'''(\mu) \;+\; \mathcal{O}(T^6)}$$

This is the Sommerfeld expansion. The leading correction is $\mathcal{O}(T^2)$, and each successive term gains two powers of $T$. For metals at room temperature, $k_BT/E_F \sim 10^{-2}$, so the series converges very rapidly and the first correction suffices in almost all practical cases.

Derivation 3: Specific Heat of the Free Electron Gas

The linear-in-$T$ electronic specific heat is one of the most important results of Fermi gas theory. Here we derive it from the Sommerfeld expansion, carefully tracking all cancellations, and compare with the classical prediction.

Step 1: Total Energy via the Sommerfeld Expansion

The total energy of the electron gas is:

$$E(T) = \int_0^\infty \epsilon\,g(\epsilon)\,f(\epsilon)\,d\epsilon$$

This is a Fermi integral with $H(\epsilon) = \epsilon\,g(\epsilon)$. Applying the Sommerfeld expansion:

$$E(T) = \int_0^\mu \epsilon\,g(\epsilon)\,d\epsilon \;+\; \frac{\pi^2}{6}(k_BT)^2\,\frac{d}{d\epsilon}\bigl[\epsilon\,g(\epsilon)\bigr]_{\epsilon=\mu} \;+\; \mathcal{O}(T^4)$$

Step 2: Handling the $\mu$-Dependence

The first term depends on $\mu(T) = E_F + \delta\mu$. Expand about $E_F$:

$$\int_0^\mu \epsilon\,g(\epsilon)\,d\epsilon = \int_0^{E_F} \epsilon\,g(\epsilon)\,d\epsilon \;+\; E_F\,g(E_F)\,\delta\mu \;+\; \mathcal{O}(\delta\mu^2)$$

Similarly, from particle number conservation, $N = \int_0^\infty g(\epsilon)\,f(\epsilon)\,d\epsilon$, the Sommerfeld expansion gives:

$$0 = g(E_F)\,\delta\mu + \frac{\pi^2}{6}(k_BT)^2\,g'(E_F) + \mathcal{O}(T^4)$$

Therefore $\delta\mu = -\frac{\pi^2}{6}(k_BT)^2\,g'(E_F)/g(E_F)$. The correction $E_F\,g(E_F)\,\delta\mu$ in the energy is of order $T^2$ and partially cancels the explicit $T^2$ term.

Step 3: Net Thermal Energy

Combining both $T^2$ contributions (the explicit Sommerfeld correction and the $\delta\mu$ correction), with the derivative $\frac{d}{d\epsilon}[\epsilon\,g(\epsilon)] = g(\epsilon) + \epsilon\,g'(\epsilon)$ evaluated at $E_F$:

$$E(T) - E_0 = \frac{\pi^2}{6}(k_BT)^2\bigl[g(E_F) + E_F\,g'(E_F)\bigr] + E_F\,g(E_F)\,\delta\mu$$

Substituting $\delta\mu$ and simplifying:

$$E(T) - E_0 = \frac{\pi^2}{6}(k_BT)^2\,g(E_F) + \mathcal{O}(T^4)$$

The terms involving $g'(E_F)$ cancel exactly, leaving a remarkably clean result that depends only on the density of states at the Fermi level.

Step 4: Differentiating to Get $C_V$

The electronic specific heat at constant volume is:

$$\boxed{C_V = \frac{\partial E}{\partial T}\bigg|_V = \frac{\pi^2}{3}\,k_B^2\,T\,g(E_F) \;\equiv\; \gamma\,T}$$

where the Sommerfeld coefficient is:

$$\gamma = \frac{\pi^2}{3}\,k_B^2\,g(E_F) = \frac{\pi^2}{2}\,\frac{N\,k_B}{E_F}\,k_B$$

where in the last equality we used the 3D free-electron result $g(E_F) = 3N/(2E_F)$.

Comparison with Classical Dulong-Petit

The classical equipartition theorem predicts that each electron contributes $\frac{3}{2}k_B$ to the specific heat, giving $C_V^{\text{classical}} = \frac{3}{2}Nk_B$. The ratio of the quantum to classical result is:

$$\frac{C_V^{\text{quantum}}}{C_V^{\text{classical}}} = \frac{\pi^2}{3}\,\frac{k_BT}{E_F} \sim \frac{T}{T_F}$$

For copper at room temperature, $T_F \approx 8 \times 10^4$ K, so $C_V^{\text{quantum}}/C_V^{\text{classical}} \approx 0.01$. This factor-of-100 suppression explains why the classical theory overestimated the electronic specific heat so drastically.

Physical picture: At temperatures $T \ll T_F$, only electrons within a shell of thickness $\sim k_BT$ around the Fermi surface can be thermally excited. The fraction of thermally active electrons is $\sim k_BT/E_F$, and each contributes $\sim k_B$ to the specific heat, yielding $C_V \sim Nk_B \cdot (k_BT/E_F) \propto T$.

Derivation 4: Pauli Paramagnetism

The magnetic susceptibility of a free electron gas provides another striking example where quantum statistics resolves a classical paradox. Classically, each magnetic moment $\mu_B$ should contribute $\mu_B^2/(k_BT)$ to the susceptibility (Curie law), predicting a strong, temperature-dependent paramagnetism. Instead, metals exhibit a weak, nearly temperature-independent paramagnetism.

Step 1: Spin-Split Fermi Seas

In an external magnetic field $\mathbf{B}$, the energy of an electron with spin $\sigma = \pm 1$ shifts by:

$$\epsilon_\sigma(\mathbf{k}) = \frac{\hbar^2 k^2}{2m} - \sigma\,\mu_B\,B$$

where $\mu_B = e\hbar/(2m_e)$ is the Bohr magneton (we take $\sigma = +1$ for spin-up, $\sigma = -1$ for spin-down). The spin-up band is shifted down by $\mu_B B$ and the spin-down band is shifted up by $\mu_B B$.

Step 2: Spin-Resolved Populations

The number of spin-up and spin-down electrons at $T = 0$ are:

$$N_\uparrow = \int_0^{E_F + \mu_B B} \frac{1}{2}\,g(\epsilon)\,d\epsilon, \qquad N_\downarrow = \int_0^{E_F - \mu_B B} \frac{1}{2}\,g(\epsilon)\,d\epsilon$$

where $\frac{1}{2}g(\epsilon)$ is the density of states per spin species. The total particle number is constrained: $N_\uparrow + N_\downarrow = N$ (this determines the common chemical potential).

Step 3: Net Magnetization

The magnetization is $M = \mu_B(N_\uparrow - N_\downarrow)$. For weak fields ($\mu_B B \ll E_F$), expand to first order in $B$:

$$N_\uparrow - N_\downarrow = \int_{E_F - \mu_B B}^{E_F + \mu_B B} \frac{1}{2}\,g(\epsilon)\,d\epsilon \approx \frac{1}{2}\,g(E_F)\cdot 2\mu_B B = \mu_B B\,g(E_F)$$

Therefore the magnetization is:

$$M = \mu_B^2\,B\,g(E_F)$$

Step 4: Pauli Susceptibility

The magnetic susceptibility per unit volume is $\chi = \mu_0 M / B$ (in SI units):

$$\boxed{\chi_P = \mu_0\,\mu_B^2\,g(E_F)}$$

Using the 3D free-electron density of states $g(E_F) = 3N/(2E_F) = 3n/(2E_F)$ per unit volume:

$$\chi_P = \frac{3\mu_0\,\mu_B^2\,n}{2E_F} = \frac{3}{2}\,\frac{\mu_0\,\mu_B^2\,n}{k_B T_F}$$

Comparison with Classical Curie Paramagnetism

The classical Curie susceptibility for $N$ spin-1/2 moments is:

$$\chi_{\text{Curie}} = \frac{\mu_0\,\mu_B^2\,n}{k_BT}$$

The Pauli result differs in two essential ways:

Magnitude: $\chi_P / \chi_{\text{Curie}} \sim T/T_F \sim 10^{-2}$ at room temperature. The susceptibility is suppressed because only electrons near the Fermi surface can flip their spins.
Temperature dependence: $\chi_P$ is essentially temperature-independent (up to $\mathcal{O}(T^2/T_F^2)$ corrections from the Sommerfeld expansion), in contrast to the $1/T$ Curie law.

Including finite-temperature corrections via the Sommerfeld expansion: $\chi_P(T) = \mu_0\,\mu_B^2\,g(E_F)\left[1 - \frac{\pi^2}{12}\left(\frac{k_BT}{E_F}\right)^2 + \cdots\right]$, showing the weak quadratic decrease with temperature.

Further Applications of Free Fermi Gas Theory

White Dwarf Stars and Electron Degeneracy Pressure

White dwarf stars are supported against gravitational collapse by the degeneracy pressure of a dense electron gas. The electrons, stripped from fully ionized carbon and oxygen nuclei, form a Fermi gas with density $n \sim 10^{30}$ cm$^{-3}$ and Fermi energy $E_F \sim 0.3$ MeV. The zero-temperature degeneracy pressure is:

$$P_0 = \frac{2}{5}\,n\,E_F = \frac{(3\pi^2)^{2/3}\hbar^2}{5m_e}\,n^{5/3}$$

Balancing this against gravitational pressure yields the mass-radius relation $R \propto M^{-1/3}$: more massive white dwarfs are smaller. When the electron Fermi energy becomes relativistic ($E_F \sim m_e c^2$), the pressure scales as $n^{4/3}$ instead of $n^{5/3}$, and gravitation wins. This sets the Chandrasekhar limit of $M_{\text{Ch}} \approx 1.4\,M_\odot$, above which no stable white dwarf can exist.

The Wiedemann-Franz Law

The Wiedemann-Franz law states that the ratio of thermal conductivity $\kappa$ to electrical conductivity $\sigma$ in a metal is proportional to temperature:

$$\frac{\kappa}{\sigma T} = L_0 = \frac{\pi^2}{3}\left(\frac{k_B}{e}\right)^2 = 2.44 \times 10^{-8}\;\text{W}\,\Omega\,\text{K}^{-2}$$

The universal constant $L_0$ is the Lorenz number. This result emerges naturally from the Sommerfeld model because both heat and charge are carried by electrons near the Fermi surface, and both conductivities share the same scattering time $\tau$. The Drude model also predicted a Lorenz number, but got the numerical prefactor wrong by a factor of $3/\pi^2$ due to incorrect (classical) statistics.

Thermoelectric Effects

The Seebeck effect (generation of a voltage by a temperature gradient) is directly related to the energy dependence of the electronic transport near $E_F$. The Sommerfeld model predicts a thermopower (Seebeck coefficient):

$$S = -\frac{\pi^2 k_B^2 T}{3\,e\,E_F}$$

This is typically a few $\mu$V/K for simple metals. The thermoelectric figure of merit $ZT = S^2 \sigma T / \kappa$ governs the efficiency of thermoelectric devices for power generation and cooling. Modern thermoelectric research aims to engineer the electronic density of states near $E_F$ to maximize $S$ while minimizing thermal conductivity.

Electronic Specific Heat Measurements

At low temperatures, the total specific heat of a metal is $C = \gamma T + \alpha T^3$, where $\gamma T$ is the electronic contribution and $\alpha T^3$ is the Debye phonon contribution. Plotting $C/T$ vs $T^2$ yields a straight line with intercept $\gamma$ and slope $\alpha$:

$$\frac{C}{T} = \gamma + \alpha\,T^2$$

The measured $\gamma$ values for real metals often differ from free-electron predictions. The ratio $\gamma_{\text{exp}}/\gamma_{\text{free}}$ defines the effective mass enhancement $m^*/m_e$, which encodes electron-electron interactions, electron-phonon coupling, and band structure effects. For example, $m^*/m_e \approx 1.3$ for copper and $\approx 13$ for heavy-fermion compounds like CeAl$_3$.

Free Electron Laser

In a free electron laser (FEL), a relativistic electron beam passes through a periodic magnetic field (undulator). The electrons emit coherent radiation whose wavelength is tunable via the beam energy. While the individual electrons are not in a Fermi gas, the physics of electron beams — space-charge effects, beam emittance, and energy spread — is deeply rooted in Fermi gas concepts. The electron beam's phase-space density is limited by the Pauli exclusion principle (the "quantum degeneracy limit"), and understanding the Fermi energy of the beam is essential for next-generation X-ray FELs that approach this limit.

Expanded Historical Context

The Drude Model and Its Limitations (1900)

Paul Drude proposed his electron theory of metals just three years after the discovery of the electron by J.J. Thomson. Drude treated conduction electrons as a classical ideal gas obeying Maxwell-Boltzmann statistics, with a mean free path determined by collisions with the fixed ion cores. The model had notable successes: it explained Ohm's law, derived the Wiedemann-Franz law (albeit with an incorrect numerical prefactor), and predicted the Hall effect.

However, the Drude model had severe failures. It predicted an electronic contribution to the specific heat of $\frac{3}{2}Nk_B$, which was approximately 100 times larger than observed values. It also predicted a strongly temperature-dependent (Curie-law) magnetic susceptibility, in stark contrast to the weak, temperature-independent paramagnetism seen experimentally. These discrepancies remained unresolved for over 25 years.

The Classical Specific Heat Paradox

The specific heat paradox was a central puzzle of early 20th-century physics. By the equipartition theorem, each electron should contribute $\frac{3}{2}k_B$ to the heat capacity. For a metal with one conduction electron per atom, this would add $\frac{3}{2}R$ per mole to the lattice specific heat, for a total of $C \approx \frac{3}{2}R + 3R = \frac{9}{2}R$. Experimentally, the specific heat of metals near room temperature is close to $3R$ (the Dulong-Petit value from lattice vibrations alone), with essentially no measurable electronic contribution.

This paradox was sometimes cited as evidence against the electron theory of metals altogether. Some physicists even questioned whether conduction electrons existed as free particles. The resolution had to wait for quantum statistics.

Sommerfeld's Quantum Theory (1927–1928)

Arnold Sommerfeld, building on the new Fermi-Dirac statistics independently formulated by Enrico Fermi (1926) and Paul Dirac (1926), realized that electrons in metals form a highly degenerate quantum gas. At room temperature, only a fraction $\sim k_BT/E_F \sim 1\%$ of the electrons are thermally excited, because the Pauli exclusion principle prevents most electrons from accessing nearby energy states (they are already occupied).

In his landmark 1928 paper in the Zeitschrift für Physik, Sommerfeld developed the systematic asymptotic expansion (now called the Sommerfeld expansion) and derived the linear electronic specific heat $C_V = \gamma T$. This single result resolved the 27-year-old specific heat paradox. Sommerfeld also computed the correct Lorenz number for the Wiedemann-Franz law, obtaining the quantum-mechanical value $L_0 = \pi^2 k_B^2/(3e^2)$, which agreed with experiment far better than Drude's classical result.

Pauli's Resolution of the Magnetic Paradox (1927)

Independently from Sommerfeld, Wolfgang Pauli applied Fermi-Dirac statistics to the magnetic properties of the electron gas in 1927. He showed that the weak, temperature-independent susceptibility $\chi_P = \mu_0\mu_B^2 g(E_F)$ is a natural consequence of the fact that only electrons near the Fermi surface can re-orient their spins in an applied magnetic field. Electrons deep below the Fermi level cannot flip because the opposite-spin state at the same energy is already occupied.

Pauli's result, together with Sommerfeld's specific heat, provided overwhelming evidence that quantum statistics governs the behavior of electrons in metals. These insights laid the essential groundwork for Felix Bloch's band theory of solids (1929), the Wigner-Seitz cellular method (1933), and ultimately Lev Landau's Fermi liquid theory (1956–1958), which explains why the free electron picture works surprisingly well even in strongly interacting electron systems.

Legacy and Modern Relevance

The free Fermi gas model remains the starting point for nearly all of condensed matter physics. The Sommerfeld expansion technique is used routinely in electronic structure calculations, transport theory, and finite-temperature quantum field theory. The model's predictions — linear specific heat, Pauli paramagnetism, Wiedemann-Franz law, and degeneracy pressure — serve as benchmarks against which all more sophisticated theories are compared. Deviations from free-electron behavior (enhanced effective masses, non-Fermi-liquid scaling, anomalous Lorenz numbers) are among the most actively studied phenomena in modern condensed matter research, from heavy-fermion compounds to high-temperature superconductors to topological semimetals.

Simulation: Fermi-Dirac Distribution

This Python simulation plots the Fermi-Dirac distribution at several temperatures, showing how the sharp step function at $T = 0$ softens as temperature increases. The right panel shows the derivative $-\partial f/\partial \epsilon$, which peaks at $E_F$ and quantifies the thermal broadening.

Python

script.py79 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Physical constants
kB = 8.617333e-5  # eV/K (Boltzmann constant)
EF = 7.0          # Fermi energy in eV (typical for copper)

# Energy range
E = np.linspace(0, 14, 1000)

# Temperatures to plot
temperatures = [0.01, 300, 1000, 3000]  # Kelvin (0.01 ~ T=0)
colors = ['#00ffff', '#4488ff', '#ff8844', '#ff4444']
labels = ['T \u2248 0 K', 'T = 300 K', 'T = 1000 K', 'T = 3000 K']

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5.5))
fig.patch.set_facecolor('#0a0a1a')
ax1.set_facecolor('#0a0a1a')
ax2.set_facecolor('#0a0a1a')

# ---- Left panel: Fermi-Dirac distribution ----
for T, color, label in zip(temperatures, colors, labels):
    beta = 1.0 / (kB * T)
    f = 1.0 / (np.exp(beta * (E - EF)) + 1.0)
    ax1.plot(E, f, color=color, linewidth=2.2, label=label)

ax1.axvline(x=EF, color='#aaaaaa', linestyle='--', alpha=0.6, linewidth=1)
ax1.text(EF + 0.2, 0.5, r'$E_F$', color='#aaaaaa', fontsize=14, va='center')
ax1.set_xlabel('Energy (eV)', fontsize=13, color='white')
ax1.set_ylabel(r'$f(\epsilon)$', fontsize=14, color='white')
ax1.set_title('Fermi-Dirac Distribution', fontsize=15, color='white', fontweight='bold')
ax1.legend(fontsize=11, facecolor='#0a0a1a', edgecolor='#333355',
           labelcolor='white', loc='center right')
ax1.set_xlim(0, 14)
ax1.set_ylim(-0.05, 1.15)
ax1.tick_params(colors='white', labelsize=11)
ax1.grid(True, alpha=0.15, color='#4488ff')
for spine in ax1.spines.values():
    spine.set_color('#333355')

# ---- Right panel: -df/dE (thermal broadening) ----
for T, color, label in zip(temperatures[1:], colors[1:], labels[1:]):
    beta = 1.0 / (kB * T)
    f = 1.0 / (np.exp(beta * (E - EF)) + 1.0)
    dfde = -np.gradient(f, E)
    ax2.plot(E, dfde, color=color, linewidth=2.2, label=label)

ax2.axvline(x=EF, color='#aaaaaa', linestyle='--', alpha=0.6, linewidth=1)
ax2.text(EF + 0.2, ax2.get_ylim()[1] * 0.3 if ax2.get_ylim()[1] > 0 else 0.5,
         r'$E_F$', color='#aaaaaa', fontsize=14, va='center')
ax2.set_xlabel('Energy (eV)', fontsize=13, color='white')
ax2.set_ylabel(r'$-\partial f / \partial \epsilon$', fontsize=14, color='white')
ax2.set_title('Thermal Broadening of Fermi Surface', fontsize=15, color='white', fontweight='bold')
ax2.legend(fontsize=11, facecolor='#0a0a1a', edgecolor='#333355',
           labelcolor='white')
ax2.set_xlim(4, 10)
ax2.tick_params(colors='white', labelsize=11)
ax2.grid(True, alpha=0.15, color='#4488ff')
for spine in ax2.spines.values():
    spine.set_color('#333355')

plt.tight_layout(pad=2.0)
plt.savefig('fermi_dirac.png', dpi=130, bbox_inches='tight', facecolor=fig.get_facecolor())

print("=" * 65)
print("  Fermi-Dirac Distribution at Various Temperatures")
print("=" * 65)
print(f"\n  Fermi energy: E_F = {EF} eV (copper-like)")
print(f"  Thermal energy at 300 K:  k_B T = {kB * 300:.4f} eV")
print(f"  Thermal energy at 1000 K: k_B T = {kB * 1000:.4f} eV")
print(f"  Thermal energy at 3000 K: k_B T = {kB * 3000:.4f} eV")
print(f"\n  Ratio k_B T / E_F at 300 K: {kB * 300 / EF:.4f}")
print(f"  >> Only electrons within ~k_B T of E_F are thermally active")
print(f"\n  At T = 0: sharp step function at E_F")
print(f"  At T > 0: step is broadened over a width ~4 k_B T")
print("\n  Plot saved as fermi_dirac.png")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Density of States in 3D

To compute thermodynamic quantities, we need the density of states $g(\epsilon)$, defined such that $g(\epsilon)\,d\epsilon$ is the number of single-particle states with energies in $[\epsilon, \epsilon + d\epsilon]$.

Derivation

For free electrons in a box of volume $V = L^3$ with periodic boundary conditions, the single-particle states are plane waves with wavevectors:

$$\mathbf{k} = \frac{2\pi}{L}(n_x, n_y, n_z), \quad n_x, n_y, n_z \in \mathbb{Z}$$

Each allowed $\mathbf{k}$-point occupies a volume $(2\pi/L)^3 = (2\pi)^3/V$ in reciprocal space. Including the spin degeneracy factor of 2, the number of states with $|\mathbf{k}| \leq k$ is:

$$N(k) = 2 \cdot \frac{\frac{4}{3}\pi k^3}{(2\pi)^3/V} = \frac{V k^3}{3\pi^2}$$

Using the free-particle dispersion $\epsilon = \hbar^2 k^2/(2m)$, we obtain $k = \sqrt{2m\epsilon}/\hbar$. The density of states is:

$$g(\epsilon) = \frac{dN}{d\epsilon} = \frac{V}{2\pi^2}\left(\frac{2m}{\hbar^2}\right)^{3/2}\sqrt{\epsilon}$$

The characteristic $\sqrt{\epsilon}$ dependence is a hallmark of three-dimensional free-particle systems. In lower dimensions, the density of states has different energy dependences: constant in 2D and $\propto 1/\sqrt{\epsilon}$ in 1D.

Fermi Energy & Fermi Surface

At $T = 0$, all states up to the Fermi energy $E_F$ are filled, and all states above are empty. The Fermi energy is determined by the total number of electrons $N$:

$$N = \int_0^{E_F} g(\epsilon)\,d\epsilon = \frac{V}{3\pi^2}\left(\frac{2mE_F}{\hbar^2}\right)^{3/2}$$

Solving for $E_F$ in terms of the electron density $n = N/V$:

$$E_F = \frac{\hbar^2}{2m}(3\pi^2 n)^{2/3}$$

The corresponding Fermi wavevector and Fermi velocity are:

$$k_F = (3\pi^2 n)^{1/3}, \qquad v_F = \frac{\hbar k_F}{m}$$

In $\mathbf{k}$-space, the occupied states at $T = 0$ fill a sphere of radius $k_F$ called the Fermi sphere, and its surface is the Fermi surface. For typical metals:

$E_F \sim 1$–$10$ eV (e.g., $E_F \approx 7.0$ eV for copper)
$k_F \sim 10^{10}$ m$^{-1}$
$v_F \sim 10^6$ m/s (about 1% of the speed of light)
The Fermi temperature $T_F = E_F/k_B \sim 10^4$–$10^5$ K, far above room temperature

Since $T \ll T_F$ for all practical temperatures, the electron gas in metals is deeply degenerate: the thermal energy $k_BT$ is a tiny perturbation compared to $E_F$. This is the regime where the Sommerfeld expansion is valid.

Ground State Properties

At $T = 0$, the total ground state energy is obtained by summing the energies of all occupied states:

$$E_0 = \int_0^{E_F} \epsilon\, g(\epsilon)\, d\epsilon = \frac{V}{2\pi^2}\left(\frac{2m}{\hbar^2}\right)^{3/2} \cdot \frac{2}{5}E_F^{5/2} = \frac{3}{5}N E_F$$

The average energy per electron at $T = 0$ is therefore $\langle\epsilon\rangle = \frac{3}{5}E_F$. This is a remarkably large energy for a system at zero temperature — it is a purely quantum mechanical effect arising from the Pauli exclusion principle.

The associated degeneracy pressure is:

$$P_0 = -\frac{\partial E_0}{\partial V}\bigg|_N = \frac{2}{3}\frac{E_0}{V} = \frac{2}{5}nE_F$$

This degeneracy pressure is what supports white dwarf stars against gravitational collapse, as shown by Chandrasekhar.

Sommerfeld Expansion

To compute thermodynamic properties at finite temperature, we need to evaluate integrals of the form:

$$I = \int_0^\infty H(\epsilon)\, f(\epsilon)\, d\epsilon$$

where $H(\epsilon)$ is a smooth function and $f(\epsilon)$ is the Fermi-Dirac distribution. The Sommerfeld expansion exploits the fact that $-\partial f/\partial\epsilon$ is sharply peaked near $\mu$ with width $\sim k_BT$.

Integration by parts and Taylor expanding $\mathcal{H}(\epsilon) = \int_0^\epsilon H(\epsilon')\,d\epsilon'$ around $\mu$ yields:

$$I = \int_0^\mu H(\epsilon)\,d\epsilon + \frac{\pi^2}{6}(k_BT)^2 H'(\mu) + \frac{7\pi^4}{360}(k_BT)^4 H'''(\mu) + \cdots$$

The key mathematical identity used is:

$$\int_{-\infty}^{\infty} \frac{x^{2n}}{(e^x + 1)(e^{-x} + 1)}\,dx = \frac{(2^{2n-1} - 1)}{n}\frac{\pi^{2n}|B_{2n}|}{(2n)!} \cdot 2$$

where $B_{2n}$ are Bernoulli numbers. Note that only even powers of $k_BT$ appear in the expansion due to the particle-hole symmetry of $-\partial f/\partial\epsilon$ near the Fermi surface.

Chemical Potential at Finite Temperature

Applying the Sommerfeld expansion to the particle number constraint $N = \int_0^\infty g(\epsilon) f(\epsilon)\,d\epsilon$, we find that the chemical potential shifts with temperature:

$$\mu(T) = E_F\left[1 - \frac{\pi^2}{12}\left(\frac{k_BT}{E_F}\right)^2 + \cdots\right]$$

For metals at room temperature, $k_BT/E_F \sim 10^{-2}$, so the correction is of order $10^{-4}$ — the chemical potential is essentially pinned at $E_F$.

Electronic Specific Heat

The total energy at finite temperature is:

$$E(T) = \int_0^\infty \epsilon\, g(\epsilon)\, f(\epsilon)\, d\epsilon$$

Applying the Sommerfeld expansion with $H(\epsilon) = \epsilon\, g(\epsilon)$:

$$E(T) = E_0 + \frac{\pi^2}{6}(k_BT)^2 g(E_F) + \mathcal{O}(T^4)$$

The electronic specific heat is therefore:

$$C_V = \frac{\partial E}{\partial T} = \frac{\pi^2}{3}k_B^2\, g(E_F)\, T \equiv \gamma T$$

where the Sommerfeld coefficient is:

$$\gamma = \frac{\pi^2}{3}k_B^2\, g(E_F)$$

Physical Interpretation

The linear-in-$T$ specific heat has a beautiful physical interpretation. At temperature $T$, only electrons within an energy shell of width $\sim k_BT$ around the Fermi surface can be thermally excited. The number of such electrons is $\sim g(E_F) \cdot k_BT$, and each gains an energy of order $k_BT$. Therefore:

$$\Delta E \sim g(E_F)(k_BT)^2 \implies C_V \sim k_B^2\, g(E_F)\, T$$

Compared to the classical (Drude) prediction of $C_V^{\text{classical}} = \frac{3}{2}Nk_B$:

$$\frac{C_V}{C_V^{\text{classical}}} = \frac{\pi^2}{3}\frac{T}{T_F} \sim 10^{-2} \quad \text{at room temperature}$$

This factor of $\sim 100$ suppression resolves the classical specific heat puzzle: experiments showed that electrons contribute far less to the specific heat than the equipartition theorem predicted.

Pauli Paramagnetism

When a magnetic field $\mathbf{B}$ is applied to the free electron gas, the single-particle energies are split by the Zeeman interaction:

$$\epsilon_{\mathbf{k},\sigma} = \frac{\hbar^2 k^2}{2m} - \sigma \mu_B B, \quad \sigma = \pm 1$$

where $\mu_B = e\hbar/(2m_e)$ is the Bohr magneton. The spin-up and spin-down Fermi surfaces are displaced by $\pm\mu_B B$, creating a net magnetization.

At $T = 0$, the number of spin-up electrons exceeds the number of spin-down electrons by:

$$\Delta N = \frac{1}{2}g(E_F) \cdot 2\mu_B B = g(E_F)\mu_B B$$

The magnetization is $M = \mu_B \Delta N / V$, giving the Pauli spin susceptibility:

$$\chi_P = \mu_0 \frac{\partial M}{\partial B} = \mu_0 \mu_B^2\, g(E_F)$$

This is temperature-independent to leading order, in contrast to the Curie law $\chi \propto 1/T$ expected for localized magnetic moments. The physical reason is the same as for the specific heat: only electrons near the Fermi surface can flip their spins.

Comparison with Classical Result

Classically (Curie paramagnetism for free spins), one expects:

$$\chi_{\text{Curie}} = \frac{\mu_0 n \mu_B^2}{k_BT}$$

The ratio $\chi_P / \chi_{\text{Curie}} \sim T/T_F \sim 10^{-2}$, showing that Pauli paramagnetism is suppressed by the same factor as the electronic specific heat. Both effects reflect the fundamental constraint that only a fraction $\sim k_BT/E_F$ of electrons participate in thermal or magnetic responses.

Comparison with the Classical Drude Model

The Drude model (1900) treats conduction electrons as a classical ideal gas. While it successfully explains Ohm's law and the Wiedemann-Franz law, it makes several incorrect predictions that are corrected by the quantum Sommerfeld model:

Property	Drude (Classical)	Sommerfeld (Quantum)
Distribution	Maxwell-Boltzmann	Fermi-Dirac
Specific heat	$C_V = \frac{3}{2}Nk_B$	$C_V = \gamma T \propto T$
Spin susceptibility	$\chi \propto 1/T$ (Curie)	$\chi_P = \text{const}$ (Pauli)
Mean velocity	$v_{\text{th}} = \sqrt{3k_BT/m}$	$v_F = \hbar k_F / m$ (T-independent)
Mean free path	$\ell \sim$ 1 Å	$\ell \sim$ 100–1000 Å

The Drude model's success with electrical conductivity and the Wiedemann-Franz law is somewhat accidental: the classical mean velocity is too low by a factor of $\sim\sqrt{T_F/T}$, and the classical mean free path is correspondingly too short, but these errors cancel in the ratio that determines the Lorenz number.

Simulation: Density of States & Specific Heat

This Fortran simulation computes the density of states $g(\epsilon) \propto \sqrt{\epsilon}$, verifies the electron count by numerical integration, and calculates the Sommerfeld coefficient $\gamma$ to show that $C_V = \gamma T$ is two orders of magnitude smaller than the classical Drude prediction.

Fortran

program.f90114 lines

program fermi_gas_dos
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  integer, parameter :: NE = 500, NT = 200

real(dp), parameter :: pi = 3.14159265358979323846_dp
  real(dp), parameter :: kB = 8.617333e-5_dp    ! eV/K
  real(dp), parameter :: EF = 7.0_dp            ! Fermi energy (eV)
  real(dp), parameter :: me = 9.109e-31_dp       ! electron mass (kg)
  real(dp), parameter :: hbar = 1.0546e-34_dp    ! J*s
  real(dp), parameter :: eV = 1.602e-19_dp       ! J per eV

real(dp) :: E, dE, T, dT, f_fd, g_E
  real(dp) :: n_total, E_total, E_total_0, C_V, gamma_th
  real(dp) :: T_arr(NT), CV_arr(NT), gamma_arr(NT)
  real(dp) :: E_arr(NE), g_arr(NE), n_arr(NE)
  real(dp) :: prefactor, beta
  integer :: i, j

! Density of states prefactor: g(E) = A * sqrt(E)
  ! A = V/(2*pi^2) * (2m/hbar^2)^(3/2)  [in 1/eV per unit volume]
  ! For a normalized DOS with n electrons, we use dimensionless units
  prefactor = 1.5_dp / (EF**1.5_dp)  ! normalized so integral from 0 to EF = 1

write(*,'(A)') '=========================================================='
  write(*,'(A)') '  Free Fermi Gas: Density of States & Specific Heat'
  write(*,'(A)') '=========================================================='
  write(*,*)

! --- Compute density of states ---
  dE = 15.0_dp / real(NE, dp)
  write(*,'(A)') '  Density of States g(E) = A * sqrt(E):'
  write(*,'(A)') '  -------------------------------------------'
  write(*,'(A8, A14, A16)') 'E (eV)', 'g(E) (a.u.)', 'n(E,T=0)'

do i = 1, NE
    E = real(i, dp) * dE
    E_arr(i) = E
    g_arr(i) = prefactor * sqrt(E)
    if (E <= EF) then
      n_arr(i) = g_arr(i)
    else
      n_arr(i) = 0.0_dp
    end if
  end do

! Print selected values
  do i = 1, NE, 50
    write(*,'(F8.2, F14.6, F16.6)') E_arr(i), g_arr(i), n_arr(i)
  end do

write(*,*)
  write(*,'(A,F8.4)') '  g(E_F) = ', prefactor * sqrt(EF)

! --- Verify electron count ---
  n_total = 0.0_dp
  do i = 1, NE
    E = E_arr(i)
    if (E <= EF) then
      n_total = n_total + g_arr(i) * dE
    end if
  end do
  write(*,'(A,F10.6)') '  Integral of g(E) from 0 to E_F = ', n_total
  write(*,'(A)') '  (Should be 1.0 for normalized DOS)'

! --- Electronic specific heat vs temperature ---
  write(*,*)
  write(*,'(A)') '  Electronic Specific Heat C_V = gamma * T:'
  write(*,'(A)') '  -------------------------------------------'

! Sommerfeld coefficient: gamma = pi^2/3 * k_B^2 * g(E_F)
  gamma_th = (pi**2 / 3.0_dp) * kB**2 * prefactor * sqrt(EF)
  write(*,'(A,ES12.4,A)') '  Sommerfeld coefficient gamma = ', gamma_th, ' eV/K^2'
  write(*,*)
  write(*,'(A8, A16, A16)') 'T (K)', 'C_V (gamma*T)', 'C_V/C_class'

dT = 2000.0_dp / real(NT, dp)
  do j = 1, NT
    T = real(j, dp) * dT
    T_arr(j) = T
    beta = 1.0_dp / (kB * T)

! Compute total energy by numerical integration
    E_total = 0.0_dp
    E_total_0 = 0.0_dp
    do i = 1, NE
      E = E_arr(i)
      f_fd = 1.0_dp / (exp(beta * (E - EF)) + 1.0_dp)
      E_total = E_total + E * g_arr(i) * f_fd * dE
      if (E <= EF) then
        E_total_0 = E_total_0 + E * g_arr(i) * dE
      end if
    end do

CV_arr(j) = gamma_th * T
    ! Ratio to classical Drude value (3/2 * n * kB)
    gamma_arr(j) = CV_arr(j) / (1.5_dp * kB)
  end do

! Print selected values
  do j = 1, NT, 20
    write(*,'(F8.1, ES16.4, F16.6)') T_arr(j), CV_arr(j), gamma_arr(j)
  end do

write(*,*)
  write(*,'(A)') '  Key results:'
  write(*,'(A,F8.4)') '    C_V/C_classical at 300 K  = ', gamma_arr(30)
  write(*,'(A)') '    >> Electronic C_V is ~1% of classical prediction'
  write(*,'(A)') '    >> This resolves the classical specific heat puzzle!'
  write(*,*)
  write(*,'(A)') '  The Sommerfeld expansion gives C_V = gamma*T at low T,'
  write(*,'(A)') '  showing only electrons within ~k_B*T of E_F contribute.'

end program fermi_gas_dos

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Summary & Key Results

The free Fermi gas model, despite neglecting electron-electron interactions entirely, captures the essential low-temperature physics of metals. The key results are:

The Fermi energy $E_F = \frac{\hbar^2}{2m}(3\pi^2 n)^{2/3}$ sets the characteristic energy scale, with $T_F \sim 10^4$–$10^5$ K for typical metals.

The electronic specific heat is $C_V = \gamma T$, suppressed by a factor of $T/T_F$ relative to the classical prediction.

The Pauli spin susceptibility $\chi_P = \mu_0\mu_B^2 g(E_F)$ is temperature-independent and suppressed by the same factor.

Both effects arise because only electrons within $\sim k_BT$ of the Fermi surface are thermally active — the vast majority are frozen by the Pauli principle.

The Sommerfeld expansion provides a systematic way to compute finite-temperature corrections as a power series in $(k_BT/E_F)^2$.

In the next chapter, we will see how these results are modified when electron-electron interactions are included, leading to the concept of quasiparticles and Landau's Fermi liquid theory.

← Part I Overview Chapter 2: Quasiparticles →