Chapter 2

BCS Theory

The BCS Ground State

The Bardeen-Cooper-Schrieffer (BCS) theory, published in 1957, provides a microscopic explanation of superconductivity. The central insight is that an arbitrarily weak attractive interaction between electrons near the Fermi surface leads to the formation of Cooper pairs and a macroscopic quantum condensate. The BCS ground state is a coherent superposition of pair-occupied and pair-unoccupied states:

$$|\text{BCS}\rangle = \prod_{\mathbf{k}} \left(u_k + v_k \, c_{\mathbf{k}\uparrow}^\dagger c_{-\mathbf{k}\downarrow}^\dagger\right)|0\rangle$$

Here $c_{\mathbf{k}\uparrow}^\dagger$ creates an electron with momentum $\mathbf{k}$ and spin up, while $c_{-\mathbf{k}\downarrow}^\dagger$ creates its time-reversed partner. The variational parameters $u_k$ and $v_k$ are the probability amplitudes for the pair state $(\mathbf{k}\uparrow, -\mathbf{k}\downarrow)$ to be empty or occupied, respectively.

The state $|\text{BCS}\rangle$ is not an eigenstate of the particle number operator — it is a superposition of states with different particle numbers. This spontaneous breaking of the $U(1)$ gauge symmetry is intimately connected to the emergence of long-range order and the Meissner effect.

Coherence Factors

Minimizing the expectation value $\langle\text{BCS}|H|\text{BCS}\rangle$ with respect to $u_k$ and $v_k$ subject to the normalization constraint$|u_k|^2 + |v_k|^2 = 1$ yields the BCS coherence factors:

$$|v_k|^2 = \frac{1}{2}\left(1 - \frac{\xi_k}{E_k}\right), \qquad |u_k|^2 = \frac{1}{2}\left(1 + \frac{\xi_k}{E_k}\right)$$

where $\xi_k = \varepsilon_k - \mu$ is the single-particle energy measured from the Fermi level and $E_k$ is the quasiparticle energy. At $T = 0$, $|v_k|^2$ is a smeared step function that transitions from 1 to 0 over an energy range of order $\Delta$ around the Fermi surface, replacing the sharp Fermi-Dirac distribution of the normal state.

Physical Interpretation

The coherence factors appear in transition matrix elements. Type-I processes (e.g., electromagnetic absorption) involve $(u_k v_{k'} - v_k u_{k'})^2$, while type-II processes (e.g., NMR relaxation) involve $(u_k v_{k'} + v_k u_{k'})^2$. The two types predict qualitatively different behaviors below $T_c$, providing a critical experimental test of BCS theory.

The BCS Gap Equation

The self-consistency condition on the order parameter $\Delta$ leads to the fundamental BCS gap equation. With a simplified pairing interaction $V_{kk'} = -V$ for states within $\hbar\omega_D$ of the Fermi surface:

$$\Delta = V \sum_{\mathbf{k}} \frac{\Delta}{2E_k} \tanh\!\left(\frac{E_k}{2k_BT}\right)$$

Converting the sum to an integral over the density of states $N(0)$ at the Fermi level:

$$1 = N(0)V \int_0^{\hbar\omega_D} \frac{d\xi}{\sqrt{\xi^2 + \Delta^2}} \tanh\!\left(\frac{\sqrt{\xi^2 + \Delta^2}}{2k_BT}\right)$$

This integral equation determines $\Delta(T)$ at every temperature. It has a non-trivial solution for $\Delta > 0$ only below a critical temperature $T_c$.

Zero-Temperature Gap

At $T = 0$, $\tanh(E_k / 2k_BT) \to 1$ and the gap equation reduces to:

$$1 = N(0)V \int_0^{\hbar\omega_D} \frac{d\xi}{\sqrt{\xi^2 + \Delta_0^2}} = N(0)V \sinh^{-1}\!\left(\frac{\hbar\omega_D}{\Delta_0}\right)$$

In the weak-coupling limit $N(0)V \ll 1$, this gives the celebrated BCS result:

$$\Delta_0 = 2\hbar\omega_D \exp\!\left(-\frac{1}{N(0)V}\right)$$

Critical Temperature

Setting $\Delta \to 0$ in the gap equation determines $T_c$:

$$1 = N(0)V \int_0^{\hbar\omega_D} \frac{d\xi}{\xi} \tanh\!\left(\frac{\xi}{2k_BT_c}\right)$$

This yields $k_BT_c = 1.134\,\hbar\omega_D\,e^{-1/N(0)V}$, and the universal BCS ratio:

$$\frac{\Delta_0}{k_BT_c} = \frac{\pi}{e^\gamma} \approx 1.764 \qquad \text{or equivalently} \qquad \frac{2\Delta_0}{k_BT_c} \approx 3.528$$

Quasiparticle Spectrum

The BCS Hamiltonian can be diagonalized by a Bogoliubov transformation, introducing quasiparticle operators $\gamma_{\mathbf{k}\sigma}^\dagger$. The quasiparticle excitation spectrum is:

$$E_k = \sqrt{\xi_k^2 + \Delta^2}$$

The minimum excitation energy is $\Delta$, occurring at the Fermi surface where $\xi_k = 0$. This energy gap is directly observable in tunneling experiments (Giaever tunneling), infrared absorption, and specific heat measurements.

Quasiparticle Density of States

The superconducting density of states follows from $N_s(E)\,dE = N(0)\,d\xi$:

$$N_s(E) = N(0) \frac{|E|}{\sqrt{E^2 - \Delta^2}} \quad \text{for } |E| > \Delta, \qquad N_s(E) = 0 \quad \text{for } |E| < \Delta$$

The DOS diverges as $1/\sqrt{E^2 - \Delta^2}$ at the gap edges, producing the characteristic coherence peaks observed in tunneling spectroscopy. This singular DOS was first confirmed by Giaever’s tunneling experiments in 1960.

Bogoliubov Transformation

The canonical transformation $\gamma_{\mathbf{k}\uparrow} = u_k c_{\mathbf{k}\uparrow} - v_k c_{-\mathbf{k}\downarrow}^\dagger$ diagonalizes the BCS Hamiltonian. The quasiparticles are superpositions of electrons and holes, and the BCS ground state is the quasiparticle vacuum:$\gamma_{\mathbf{k}\sigma}|\text{BCS}\rangle = 0$ for all $\mathbf{k}, \sigma$.

BCS Thermodynamics

Specific Heat Jump

At $T_c$, the electronic specific heat exhibits a discontinuous jump. BCS theory predicts a universal ratio for the specific heat discontinuity:

$$\frac{\Delta C}{\gamma T_c} = \frac{C_s(T_c) - C_n(T_c)}{\gamma T_c} = 1.43$$

where $\gamma = \frac{2}{3}\pi^2 N(0)k_B^2$ is the Sommerfeld coefficient. This jump is a hallmark of the second-order phase transition at $T_c$ and has been confirmed experimentally in many conventional superconductors.

Condensation Energy

The condensation energy is the free energy difference between the normal and superconducting states at $T = 0$:

$$F_n(0) - F_s(0) = \frac{1}{2}N(0)\Delta_0^2 = \frac{B_c^2(0)}{2\mu_0}$$

This connects the microscopic gap $\Delta_0$ to the macroscopic thermodynamic critical field $B_c(0)$. For typical conventional superconductors, the condensation energy is only about $10^{-8}$ eV per atom — a remarkably small energy that nevertheless produces dramatic macroscopic effects.

Entropy and Free Energy

The electronic entropy $S_s = -2k_B \sum_{\mathbf{k}} [f_k \ln f_k + (1-f_k)\ln(1-f_k)]$, where $f_k = [e^{E_k/k_BT}+1]^{-1}$, is exponentially suppressed at low $T$as $S_s \propto e^{-\Delta/k_BT}$. The thermodynamic critical field follows from$B_c^2(T)/2\mu_0 = F_n(T) - F_s(T)$, with $F_s = F_n$ at $T_c$.

Low-Temperature Specific Heat

Well below $T_c$, the electronic specific heat is dominated by thermally excited quasiparticles across the gap:

$$C_{es} \propto \left(\frac{\Delta_0}{k_BT}\right)^{5/2} e^{-\Delta_0/k_BT} \qquad (T \ll T_c)$$

This exponential suppression (compared to the linear $C_n = \gamma T$ of the normal state) provides direct evidence for the energy gap and can be used to extract $\Delta_0$from calorimetric measurements.

Universal BCS Predictions

A remarkable feature of BCS theory in the weak-coupling limit is that several thermodynamic ratios are universal — independent of material parameters. These dimensionless numbers provide stringent tests of the theory:

Gap-to-$T_c$ ratio$2\Delta_0 / k_BT_c = 3.528$

Specific heat jump$\Delta C / \gamma T_c = 1.43$

Critical field slope$[T_c/B_c(0)][dB_c/dT]_{T_c} = -1.74$

Deviation parameter$D(t) = B_c(T)/B_c(0) - (1-t^2) \approx 0$

Deviations indicate strong-coupling effects: lead has $2\Delta_0/k_BT_c \approx 4.38$, well above the weak-coupling value, as captured by Eliashberg theory.

Detailed Derivation: BCS Variational Calculation

We derive the BCS coherence factors and gap equation from the variational minimization of the ground state energy. The BCS trial wavefunction is:

$$|\text{BCS}\rangle = \prod_{\mathbf{k}} \left(u_k + v_k\, c_{\mathbf{k}\uparrow}^\dagger c_{-\mathbf{k}\downarrow}^\dagger\right)|0\rangle$$

Step 1: Compute the Expectation Value of H

The reduced BCS Hamiltonian (keeping only pairing terms) is:

$$H = \sum_{\mathbf{k}\sigma} \xi_k\, c_{\mathbf{k}\sigma}^\dagger c_{\mathbf{k}\sigma} + \sum_{\mathbf{k},\mathbf{k}'} V_{\mathbf{k}\mathbf{k}'}\, c_{\mathbf{k}\uparrow}^\dagger c_{-\mathbf{k}\downarrow}^\dagger c_{-\mathbf{k}'\downarrow} c_{\mathbf{k}'\uparrow}$$

Using the expectation values $\langle c_{\mathbf{k}\uparrow}^\dagger c_{\mathbf{k}\uparrow}\rangle = |v_k|^2$and $\langle c_{\mathbf{k}\uparrow}^\dagger c_{-\mathbf{k}\downarrow}^\dagger\rangle = u_k^* v_k$(from the BCS state), the ground state energy is:

$$E_{\text{BCS}} = 2\sum_{\mathbf{k}} \xi_k |v_k|^2 + \sum_{\mathbf{k},\mathbf{k}'} V_{\mathbf{k}\mathbf{k}'}\, u_k^* v_k\, u_{k'} v_{k'}^*$$

Step 2: Parametrize with an Angle

The normalization constraint $|u_k|^2 + |v_k|^2 = 1$ is automatically satisfied by writing $u_k = \cos(\theta_k/2)$ and $v_k = \sin(\theta_k/2)$. Then:

$$|v_k|^2 = \frac{1}{2}(1 - \cos\theta_k), \qquad u_k v_k = \frac{1}{2}\sin\theta_k$$

Substituting:

$$E = \sum_{\mathbf{k}} \xi_k(1 - \cos\theta_k) + \frac{1}{4}\sum_{\mathbf{k},\mathbf{k}'} V_{\mathbf{k}\mathbf{k}'}\sin\theta_k\sin\theta_{k'}$$

Step 3: Minimize with Respect to $\theta_k$

Setting $\partial E / \partial \theta_k = 0$:

$$\xi_k \sin\theta_k + \frac{1}{2}\cos\theta_k \sum_{\mathbf{k}'} V_{\mathbf{k}\mathbf{k}'}\sin\theta_{k'} = 0$$

Define the gap function $\Delta_k = -\frac{1}{2}\sum_{\mathbf{k}'} V_{\mathbf{k}\mathbf{k}'}\sin\theta_{k'}$. Then $\xi_k\sin\theta_k = \Delta_k\cos\theta_k$, giving:

$$\tan\theta_k = \frac{\Delta_k}{\xi_k} \quad \Longrightarrow \quad \sin\theta_k = \frac{\Delta_k}{E_k}, \quad \cos\theta_k = \frac{\xi_k}{E_k}$$

where $E_k = \sqrt{\xi_k^2 + \Delta_k^2}$ is the quasiparticle energy. From these we recover the BCS coherence factors:

$$|v_k|^2 = \frac{1}{2}\left(1 - \frac{\xi_k}{E_k}\right), \qquad |u_k|^2 = \frac{1}{2}\left(1 + \frac{\xi_k}{E_k}\right)$$

Step 4: Self-Consistent Gap Equation

Substituting $\sin\theta_{k'} = \Delta_{k'}/E_{k'}$ back into the definition of$\Delta_k$:

$$\Delta_k = -\sum_{\mathbf{k}'} V_{\mathbf{k}\mathbf{k}'} \frac{\Delta_{k'}}{2E_{k'}}$$

With the BCS model interaction $V_{\mathbf{k}\mathbf{k}'} = -V$ (constant and attractive within $\hbar\omega_D$ of $E_F$), the gap becomes $\mathbf{k}$-independent ($\Delta_k = \Delta$), and canceling $\Delta$:

$$1 = V\sum_{\mathbf{k}} \frac{1}{2\sqrt{\xi_k^2 + \Delta^2}} = N(0)V\int_0^{\hbar\omega_D} \frac{d\xi}{\sqrt{\xi^2 + \Delta^2}}$$

This integral evaluates to $N(0)V\,\text{arcsinh}(\hbar\omega_D/\Delta)$. In the weak-coupling limit $\Delta \ll \hbar\omega_D$, we use $\text{arcsinh}(x) \approx \ln(2x)$for large $x$, giving $\Delta_0 = 2\hbar\omega_D\,e^{-1/N(0)V}$.

Derivation: Bogoliubov-Valatin Transformation

The BCS Hamiltonian can be diagonalized exactly by introducing quasiparticle operators via the Bogoliubov-Valatin canonical transformation.

Step 1: Define Quasiparticle Operators

Introduce new fermionic operators as linear combinations of electron creation and annihilation operators:

$$\gamma_{\mathbf{k}\uparrow} = u_k c_{\mathbf{k}\uparrow} - v_k c_{-\mathbf{k}\downarrow}^\dagger$$

$$\gamma_{-\mathbf{k}\downarrow}^\dagger = u_k c_{-\mathbf{k}\downarrow}^\dagger + v_k c_{\mathbf{k}\uparrow}$$

These satisfy canonical anticommutation relations$\{\gamma_{\mathbf{k}\sigma}, \gamma_{\mathbf{k}'\sigma'}^\dagger\} = \delta_{\mathbf{k}\mathbf{k}'}\delta_{\sigma\sigma'}$provided $|u_k|^2 + |v_k|^2 = 1$.

Step 2: Invert the Transformation

Solving for the original electron operators:

$$c_{\mathbf{k}\uparrow} = u_k \gamma_{\mathbf{k}\uparrow} + v_k \gamma_{-\mathbf{k}\downarrow}^\dagger$$

$$c_{-\mathbf{k}\downarrow}^\dagger = -v_k \gamma_{\mathbf{k}\uparrow} + u_k \gamma_{-\mathbf{k}\downarrow}^\dagger$$

Substituting into the mean-field BCS Hamiltonian and choosing $u_k, v_k$ to eliminate the off-diagonal (pair-creating) terms yields the condition$2\xi_k u_k v_k = \Delta(u_k^2 - v_k^2)$, which gives the standard BCS coherence factors. The diagonalized Hamiltonian is:

$$H_{\text{BCS}} = E_{\text{gs}} + \sum_{\mathbf{k}\sigma} E_k\, \gamma_{\mathbf{k}\sigma}^\dagger \gamma_{\mathbf{k}\sigma}$$

where $E_k = \sqrt{\xi_k^2 + \Delta^2}$ and $E_{\text{gs}}$ is the BCS ground state energy. The quasiparticle excitations are gapped: the minimum energy to create an excitation is $\Delta$, occurring at $\xi_k = 0$ (the Fermi surface).

Derivation: The BCS Critical Temperature

The critical temperature $T_c$ is determined by the condition that the gap vanishes: $\Delta(T_c) = 0$. We derive the universal BCS ratio$\Delta_0/k_BT_c$ from the finite-temperature gap equation.

Step 1: Finite-Temperature Gap Equation

At finite temperature, the gap equation acquires a Fermi function factor:

$$1 = N(0)V\int_0^{\hbar\omega_D} \frac{d\xi}{\sqrt{\xi^2 + \Delta^2}}\tanh\!\left(\frac{\sqrt{\xi^2 + \Delta^2}}{2k_BT}\right)$$

At $T = T_c$, $\Delta \to 0$, so the equation becomes:

$$1 = N(0)V\int_0^{\hbar\omega_D} \frac{d\xi}{\xi}\tanh\!\left(\frac{\xi}{2k_BT_c}\right)$$

Step 2: Evaluate the Integral

Using the substitution $x = \xi/(2k_BT_c)$ and the result$\int_0^X \frac{\tanh x}{x}\,dx \approx \ln(1.134\,X)$ for $X \gg 1$(derived from the integral identity involving the Euler-Mascheroni constant $\gamma_E$and the relation $\ln(4e^{\gamma_E}/\pi) = \ln(1.134\ldots)$):

$$\frac{1}{N(0)V} = \ln\!\left(\frac{1.134\,\hbar\omega_D}{k_BT_c}\right)$$

Therefore:

$$k_BT_c = 1.134\,\hbar\omega_D\,e^{-1/N(0)V}$$

Step 3: The Universal Ratio

Dividing the zero-temperature gap $\Delta_0 = 2\hbar\omega_D\,e^{-1/N(0)V}$by $k_BT_c = 1.134\,\hbar\omega_D\,e^{-1/N(0)V}$:

$$\frac{\Delta_0}{k_BT_c} = \frac{2}{1.134} = \frac{\pi}{e^{\gamma_E}} \approx 1.764$$

where $\gamma_E \approx 0.5772$ is the Euler-Mascheroni constant. This gives the celebrated BCS ratio $2\Delta_0/(k_BT_c) \approx 3.528$, which is universal in the weak-coupling limit: it depends on no material parameters whatsoever. Experimental values for conventional superconductors range from 3.5 (Al, Sn, In) to 4.4 (Pb, Hg), with deviations indicating strong-coupling effects.

Historical Context

The Road to BCS

Superconductivity was discovered by Heike Kamerlingh Onnes in 1911 in mercury at 4.2 K, yet a microscopic explanation eluded physicists for 46 years. Many of the greatest theorists of the 20th century — including Einstein, Feynman, Bohr, Heisenberg, Bloch, and Landau — attempted and failed to explain the phenomenon. John Bardeen, already a Nobel laureate for the transistor, recognized that a successful theory required three ingredients: (1) an attractive interaction between electrons (provided by phonons, as shown by Frohlich in 1950), (2) a mechanism for pair formation in the presence of the Fermi sea (Cooper's 1956 result), and (3) a many-body ground state incorporating all pairs coherently.

The third ingredient was provided by J. Robert Schrieffer, then a graduate student, who wrote down the variational wavefunction on a New York City subway ride in January 1957. The complete BCS paper, “Theory of Superconductivity,” was published in Physical Review in December 1957 and immediately recognized as a landmark achievement. It explained the energy gap, the Meissner effect, the isotope effect, the specific heat jump at $T_c$, coherence factors, and the London penetration depth — all from a single microscopic framework.

Experimental Confirmations

The BCS prediction of the energy gap was dramatically confirmed by Ivar Giaever's tunneling experiments in 1960, for which he shared the 1973 Nobel Prize. The simultaneous explanation of the contrasting behavior in ultrasonic attenuation (type-I coherence factor, dropping rapidly below$T_c$) and NMR relaxation (type-II coherence factor, showing the Hebel-Slichter peak just below $T_c$) provided a uniquely stringent test that no competing theory could match. The Eliashberg extension (1960) generalized BCS to strong-coupling superconductors, resolving the quantitative discrepancies in Pb and Hg.

Applications of BCS Theory

Real-World Applications

1. Materials Design for Superconductors

BCS theory guides the search for new superconducting materials through the McMillan-Allen-Dynes formula, which predicts $T_c$ from the electron-phonon spectral function $\alpha^2F(\omega)$. This has led to the optimization of A15 compounds (Nb$_3$Sn, Nb$_3$Ge), the discovery of MgB$_2$($T_c = 39$ K), and the prediction of high-$T_c$ superconductivity in hydrogen-rich compounds under pressure (LaH$_{10}$ at 250 K).

2. Tunneling Spectroscopy

The BCS density of states $N_s(E) = N(0)|E|/\sqrt{E^2 - \Delta^2}$ is directly measured by scanning tunneling microscopy (STM) and planar tunnel junctions. This technique is used to characterize superconducting gaps, identify gap symmetry (s-wave vs d-wave), and detect subgap states due to impurities or vortex cores.

3. Superconducting Electronics (SFQ Logic)

Single flux quantum (SFQ) digital logic exploits the BCS energy gap for ultra-fast, low-power computation. Clock speeds exceed 100 GHz with power dissipation per gate orders of magnitude below CMOS. Applications include digital signal processing, precision timing, and readout electronics for quantum computers.

4. Superconducting Detectors

The BCS energy gap ($\sim$1 meV) sets the energy threshold for breaking Cooper pairs, enabling photon detection with energy resolution $\delta E \sim$ eV. Transition-edge sensors (TES), kinetic inductance detectors (KIDs), and superconducting nanowire single-photon detectors (SNSPDs) are used in astronomy (CMB measurements), quantum communication, and dark matter searches.

5. Nuclear and Astrophysics

The BCS pairing mechanism applies beyond electrons in metals. Nucleon pairing in atomic nuclei (proton-proton and neutron-neutron Cooper pairs) explains nuclear superfluidity and even-odd mass differences. In neutron stars, neutron pairing in the$ {}^1S_0$ and $ {}^3P_2$ channels produces superfluid phases that affect cooling rates, glitch dynamics, and gravitational wave emission.

Python Simulation: BCS Gap Equation

Self-consistent numerical solution of the BCS gap equation. Left panel shows the temperature dependence of the superconducting gap $\Delta(T)/\Delta_0$. Right panel shows the quasiparticle density of states with the characteristic coherence peaks at $E = \pm\Delta_0$.

Self-Consistent BCS Gap and Quasiparticle DOS

Python

Solves the BCS gap equation numerically and plots the quasiparticle density of states

script.py128 lines

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

# ──────────────────────────────────────────────
# Self-consistent BCS Gap Equation: Delta(T)
# and Quasiparticle Density of States
# ──────────────────────────────────────────────

# Parameters
omega_D = 1.0        # Debye frequency (energy units)
N0V = 0.3            # Dimensionless coupling N(0)V
kB = 1.0             # Boltzmann constant

# Zero-temperature gap from BCS result
Delta_0 = 2.0 * omega_D * np.exp(-1.0 / N0V)
Tc = Delta_0 / 1.764  # BCS relation

print(f"BCS Parameters:")
print(f"  Debye energy: omega_D = {omega_D:.3f}")
print(f"  Coupling: N(0)V = {N0V:.3f}")
print(f"  Zero-T gap: Delta_0 = {Delta_0:.6f}")
print(f"  Critical temperature: T_c = {Tc:.6f}")
print(f"  2*Delta_0 / (kB*Tc) = {2*Delta_0/(kB*Tc):.3f} (BCS: 3.528)")

def solve_gap(T, Delta_init, tol=1e-8, max_iter=300):
    """Solve self-consistent BCS gap equation at temperature T using numerical integration."""
    if T >= Tc * 0.999:
        return 0.0
    Delta = max(Delta_init, 1e-12)
    xi = np.linspace(0, omega_D, 2000)
    for _ in range(max_iter):
        Ek = np.sqrt(xi**2 + Delta**2)
        if T < 1e-15:
            integrand = 1.0 / Ek
        else:
            beta_Ek = np.minimum(Ek / (kB * T), 500.0)
            integrand = np.tanh(0.5 * beta_Ek) / Ek
        rhs = N0V * np.trapezoid(integrand, xi)
        if rhs <= 1.0:
            return 0.0
        # Fixed-point iteration: Delta_new so that integral matches
        Delta_new = Delta * np.sqrt(rhs)
        if abs(Delta_new - Delta) < tol * Delta_0:
            return Delta_new
        Delta = 0.7 * Delta_new + 0.3 * Delta  # damped iteration
    return Delta

# Compute Delta(T) over temperature range
N_temps = 200
T_ratio = np.linspace(0.0, 1.05, N_temps)
T_vals = T_ratio * Tc
Delta_vals = np.zeros(N_temps)

Delta_current = Delta_0
for i, T in enumerate(T_vals):
    if T / Tc > 1.0:
        Delta_vals[i] = 0.0
    else:
        Delta_current = solve_gap(T, Delta_current)
        Delta_vals[i] = Delta_current

# BCS approximate formula for comparison
T_bcs = np.linspace(0, 0.999, 200)
Delta_bcs_approx = Delta_0 * np.sqrt(np.maximum(np.cos(0.5 * np.pi * T_bcs**2), 0))

# Quasiparticle DOS at T = 0
E_range = np.linspace(-3 * Delta_0, 3 * Delta_0, 1000)
dos_normal = np.ones_like(E_range)  # N(E)/N(0) = 1 in normal state
E_abs = np.abs(E_range)
dos_sc = np.where(E_abs > Delta_0 * 1.001,
                  E_abs / np.sqrt(E_abs**2 - Delta_0**2),
                  0.0)

# ── Plotting ──
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0a0a1a')

for ax in axes:
    ax.set_facecolor('#0a0a1a')
    for spine in ax.spines.values():
        spine.set_edgecolor('#334155')

# Left: Delta(T) vs T/Tc
axes[0].plot(T_ratio, Delta_vals / Delta_0, color='#2dd4bf', linewidth=2.5,
             label=r'Self-consistent $\Delta(T)$')
axes[0].plot(T_bcs, Delta_bcs_approx / Delta_0, '--', color='#818cf8',
             linewidth=1.5, alpha=0.7, label='BCS approx.')
axes[0].axvline(x=1.0, color='#ef4444', linewidth=1, linestyle=':', alpha=0.6,
                label=r'$T_c$')
axes[0].axhline(y=1.0, color='#64748b', linewidth=0.5, linestyle='--', alpha=0.4)
axes[0].set_xlabel(r'$T / T_c$', color='#94a3b8', fontsize=13)
axes[0].set_ylabel(r'$\Delta(T) / \Delta_0$', color='#94a3b8', fontsize=13)
axes[0].set_title('BCS Gap vs Temperature', color='white', fontsize=14, fontweight='bold')
axes[0].set_xlim(0, 1.1)
axes[0].set_ylim(-0.05, 1.15)
axes[0].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[0].tick_params(colors='#94a3b8')
axes[0].fill_between(T_ratio, 0, Delta_vals / Delta_0, alpha=0.15, color='#2dd4bf')

# Right: Quasiparticle DOS
axes[1].plot(E_range / Delta_0, dos_normal, '--', color='#64748b', linewidth=1.5,
             label='Normal state', alpha=0.7)
axes[1].plot(E_range / Delta_0, dos_sc, color='#f59e0b', linewidth=2,
             label='Superconducting')
axes[1].fill_between(E_range / Delta_0, 0, dos_sc, alpha=0.15, color='#f59e0b')
axes[1].axvline(x=-1.0, color='#ef4444', linewidth=1, linestyle=':', alpha=0.5)
axes[1].axvline(x=1.0, color='#ef4444', linewidth=1, linestyle=':', alpha=0.5,
                label=r'$\pm\Delta_0$')
axes[1].set_xlabel(r'$E / \Delta_0$', color='#94a3b8', fontsize=13)
axes[1].set_ylabel(r'$N_s(E) / N(0)$', color='#94a3b8', fontsize=13)
axes[1].set_title('Quasiparticle Density of States', color='white', fontsize=14,
                  fontweight='bold')
axes[1].set_xlim(-3, 3)
axes[1].set_ylim(0, 5)
axes[1].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[1].tick_params(colors='#94a3b8')

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

print(f"\nGap values at selected temperatures:")
for t_frac in [0.0, 0.25, 0.5, 0.75, 0.9, 0.95, 1.0]:
    idx = np.argmin(np.abs(T_ratio - t_frac))
    print(f"  T/Tc = {t_frac:.2f}: Delta/Delta_0 = {Delta_vals[idx]/Delta_0:.4f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Implementation: BCS Thermodynamics

High-performance Fortran computation of BCS thermodynamic quantities: specific heat, entropy, and free energy as functions of temperature, with comparison between superconducting and normal states.

BCS Thermodynamic Properties

Fortran

Computes specific heat, entropy, and free energy vs temperature in the BCS framework

program.f90163 lines

program bcs_thermodynamics
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  integer, parameter :: N_temp = 300
  integer, parameter :: N_xi = 2000

real(dp), parameter :: pi = 3.14159265358979323846_dp
  real(dp), parameter :: omega_D = 1.0_dp
  real(dp), parameter :: N0V = 0.30_dp
  real(dp) :: Delta_0, Tc, kB
  real(dp) :: T, dT, xi, dxi, Ek, beta_Ek
  real(dp) :: Delta, Delta_old, rhs
  real(dp) :: S_sc, S_n, E_sc, E_n, fermi_f
  real(dp) :: Delta_arr(N_temp), T_arr(N_temp)
  real(dp) :: entropy_sc(N_temp), entropy_n(N_temp)
  real(dp) :: sp_heat_sc(N_temp), sp_heat_n(N_temp)
  real(dp) :: free_en_sc(N_temp), free_en_n(N_temp)
  integer :: i, j, iter

kB = 1.0_dp
  Delta_0 = 2.0_dp * omega_D * exp(-1.0_dp / N0V)
  Tc = Delta_0 / 1.764_dp

write(*,'(A)') '=== BCS Thermodynamics (Fortran) ==='
  write(*,'(A,F10.6)') '  Delta_0 = ', Delta_0
  write(*,'(A,F10.6)') '  T_c     = ', Tc
  write(*,'(A)') ''

dxi = omega_D / real(N_xi, dp)

! ── Solve gap equation at each temperature ──
  Delta = Delta_0
  do i = 1, N_temp
    T_arr(i) = real(i - 1, dp) / real(N_temp - 1, dp) * 1.05_dp * Tc
    T = T_arr(i)

if (T / Tc > 0.9999_dp) then
      Delta_arr(i) = 0.0_dp
      Delta = 0.0_dp
      cycle
    end if

if (T < 1.0e-12_dp) then
      Delta_arr(i) = Delta_0
      Delta = Delta_0
      cycle
    end if

! Iterate to self-consistency
    Delta = max(Delta, 1.0e-8_dp)
    do iter = 1, 500
      Delta_old = Delta
      rhs = 0.0_dp
      do j = 1, N_xi
        xi = (real(j, dp) - 0.5_dp) * dxi
        Ek = sqrt(xi**2 + Delta**2)
        beta_Ek = Ek / (kB * T)
        if (beta_Ek > 500.0_dp) then
          rhs = rhs + dxi / Ek
        else
          rhs = rhs + dxi * tanh(0.5_dp * beta_Ek) / Ek
        end if
      end do
      rhs = N0V * rhs

if (rhs <= 1.0_dp) then
        Delta = 0.0_dp
        exit
      end if
      Delta = Delta_old * sqrt(rhs)
      if (abs(Delta - Delta_old) < 1.0e-12_dp * Delta_0) exit
    end do
    Delta_arr(i) = Delta
  end do

! ── Compute thermodynamic quantities ──
  do i = 1, N_temp
    T = T_arr(i)
    Delta = Delta_arr(i)

! Superconducting entropy: S_sc = -2*N(0) * sum_k [f*ln(f) + (1-f)*ln(1-f)]
    S_sc = 0.0_dp
    E_sc = 0.0_dp
    S_n = 0.0_dp
    E_n = 0.0_dp

if (T < 1.0e-12_dp) then
      entropy_sc(i) = 0.0_dp
      entropy_n(i) = 0.0_dp
      free_en_sc(i) = -0.5_dp * Delta_0**2
      free_en_n(i) = 0.0_dp
      sp_heat_sc(i) = 0.0_dp
      sp_heat_n(i) = 0.0_dp
      cycle
    end if

do j = 1, N_xi
      xi = (real(j, dp) - 0.5_dp) * dxi
      Ek = sqrt(xi**2 + Delta**2)
      beta_Ek = Ek / (kB * T)

! Superconducting state
      if (beta_Ek < 500.0_dp) then
        fermi_f = 1.0_dp / (exp(beta_Ek) + 1.0_dp)
        if (fermi_f > 1.0e-15_dp .and. fermi_f < 1.0_dp - 1.0e-15_dp) then
          S_sc = S_sc - 2.0_dp * dxi * (fermi_f * log(fermi_f) + &
                 (1.0_dp - fermi_f) * log(1.0_dp - fermi_f))
        end if
        E_sc = E_sc + 2.0_dp * dxi * Ek * fermi_f
      end if

! Normal state (Delta = 0, so Ek = |xi|)
      beta_Ek = abs(xi) / (kB * T)
      if (beta_Ek < 500.0_dp) then
        fermi_f = 1.0_dp / (exp(beta_Ek) + 1.0_dp)
        if (fermi_f > 1.0e-15_dp .and. fermi_f < 1.0_dp - 1.0e-15_dp) then
          S_n = S_n - 2.0_dp * dxi * (fermi_f * log(fermi_f) + &
                (1.0_dp - fermi_f) * log(1.0_dp - fermi_f))
        end if
        E_n = E_n + 2.0_dp * dxi * abs(xi) * fermi_f
      end if
    end do

entropy_sc(i) = S_sc * kB
    entropy_n(i) = S_n * kB

! Free energy: F = E - TS - condensation energy offset
    free_en_sc(i) = E_sc - T * S_sc * kB - 0.5_dp * Delta**2
    free_en_n(i) = E_n - T * S_n * kB
  end do

! Specific heat by numerical differentiation: C = T * dS/dT
  sp_heat_sc(1) = 0.0_dp
  sp_heat_n(1) = 0.0_dp
  do i = 2, N_temp
    dT = T_arr(i) - T_arr(i-1)
    if (dT > 1.0e-15_dp) then
      sp_heat_sc(i) = T_arr(i) * (entropy_sc(i) - entropy_sc(i-1)) / dT
      sp_heat_n(i) = T_arr(i) * (entropy_n(i) - entropy_n(i-1)) / dT
    else
      sp_heat_sc(i) = 0.0_dp
      sp_heat_n(i) = 0.0_dp
    end if
  end do

! ── Output results ──
  write(*,'(A)') '  T/Tc     Delta/D0   S_sc/S_n   C_sc/C_n    F_sc-F_n'
  write(*,'(A)') ' ------   --------   --------   --------   ----------'
  do i = 1, N_temp, 20
    T = T_arr(i)
    if (T / Tc > 1.02_dp) cycle
    if (entropy_n(i) > 1.0e-10_dp) then
      write(*,'(F7.3,3X,F8.5,3X,F8.4,3X,F8.4,3X,F10.6)') &
        T/Tc, Delta_arr(i)/Delta_0, entropy_sc(i)/max(entropy_n(i),1e-15_dp), &
        sp_heat_sc(i)/max(sp_heat_n(i),1e-15_dp), free_en_sc(i)-free_en_n(i)
    end if
  end do

write(*,'(A)') ''
  write(*,'(A,F10.6)') '  Condensation energy: ', 0.5_dp * Delta_0**2
  write(*,'(A,F8.4)')  '  2*Delta_0/(kB*Tc):   ', 2.0_dp*Delta_0/(kB*Tc)

end program bcs_thermodynamics

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Experimental Verification

Tunneling Spectroscopy

Giaever’s tunneling experiments (1960, Nobel Prize 1973) provided the most direct confirmation of the BCS density of states. In a normal-metal/insulator/superconductor (NIS) junction, the differential conductance is:

$$\frac{dI}{dV} \propto N_s(eV) = \text{Re}\left[\frac{|eV|}{\sqrt{(eV)^2 - \Delta^2}}\right]$$

The sharp coherence peaks at $eV = \pm\Delta$ and the hard gap for $|eV| < \Delta$match the BCS prediction with remarkable precision.

Isotope Effect

BCS theory predicts $T_c \propto M^{-\alpha}$ with $\alpha_{\text{BCS}} = 1/2$, since $\omega_D \propto M^{-1/2}$. This was confirmed for Hg, Sn, and Pb, providing key evidence for phonon-mediated pairing.

Coherence Factor Tests

BCS coherence factors predict contrasting behavior for type-I processes (ultrasonic attenuation drops rapidly below $T_c$) and type-II processes (NMR relaxation shows the Hebel-Slichter peak just below $T_c$). The simultaneous description of both by BCS theory firmly established the coherent nature of the condensate.

Beyond Weak-Coupling BCS

The Eliashberg extension replaces the constant BCS pairing interaction with the full retarded electron-phonon spectral function $\alpha^2 F(\omega)$. The McMillan formula gives an approximate $T_c$:

$$T_c = \frac{\omega_{\log}}{1.2}\exp\!\left[-\frac{1.04(1 + \lambda)}{\lambda - \mu^*(1 + 0.62\lambda)}\right]$$

where $\lambda$ is the electron-phonon coupling constant and$\mu^*$ is the Coulomb pseudopotential. Strong-coupling corrections modify the universal BCS ratios: for Pb ($\lambda \approx 1.55$),$2\Delta_0/k_BT_c \approx 4.38$ and$\Delta C/\gamma T_c \approx 2.71$, both exceeding the weak-coupling values.

← Chapter 1: Cooper Pairs Chapter 3: Ginzburg-Landau →