Chapter 1

Cooper Pairs

The Cooper Instability

In 1956, Leon Cooper demonstrated a remarkable result: the filled Fermi sea is unstable toward the formation of at least one bound pair of electrons, provided there exists any attractive interaction between them, no matter how weak. This is the Cooper instability — a fundamentally many-body effect that does not occur in vacuum, where a minimum threshold attraction is needed for binding in three dimensions.

The key insight is that the Pauli exclusion principle, by filling all states below the Fermi energy $E_F$, restricts the pair to occupy states only above the Fermi surface. This effectively reduces the dimensionality of the problem — the density of states at the Fermi level is finite and constant, mimicking a two-dimensional system where bound states exist for arbitrarily weak attraction.

Why Any Attraction Suffices

In three dimensions, the Schrodinger equation for two particles with an attractive potential $V$ has no bound state unless $V$ exceeds a critical strength. However, when restricted to a thin shell of thickness$\hbar\omega_D$ above the Fermi surface, the effective density of states $N(0)$ is constant, and the problem maps to a 2D-like situation. In 2D, a bound state exists for any $V > 0$.

The Cooper Problem

Consider two electrons added to a filled Fermi sea at $T = 0$. The Fermi sea is inert: all states with $|\mathbf{k}| \leq k_F$ are occupied and blocked by the Pauli principle. The two added electrons interact via an attractive potential $V(\mathbf{r}_1 - \mathbf{r}_2)$ while being forced to occupy states above $k_F$.

We seek a two-particle state with zero total momentum ($\mathbf{K} = \mathbf{k}_1 + \mathbf{k}_2 = 0$), so the two electrons have opposite momenta $\mathbf{k}$ and $-\mathbf{k}$. The pair wavefunction in the center-of-mass frame is:

$$\psi(\mathbf{r}_1, \mathbf{r}_2) = \sum_{\mathbf{k}} g_{\mathbf{k}}\, e^{i\mathbf{k}\cdot(\mathbf{r}_1 - \mathbf{r}_2)}$$

where the sum runs over $|\mathbf{k}| > k_F$ only, and$g_{\mathbf{k}}$ are the pair amplitudes. Substituting into the Schrodinger equation and projecting onto plane waves gives the Cooper equation:

$$(E - 2\epsilon_{\mathbf{k}})\, g_{\mathbf{k}} = \sum_{\mathbf{k}'} V_{\mathbf{k}\mathbf{k}'}\, g_{\mathbf{k}'}$$

where $\epsilon_{\mathbf{k}} = \hbar^2 k^2 / 2m$ is the single-particle energy and $E$ is the total pair energy. For a simplified contact interaction, $V_{\mathbf{k}\mathbf{k}'} = -V$ for$|\epsilon_{\mathbf{k}} - E_F|, |\epsilon_{\mathbf{k}'} - E_F| < \hbar\omega_D$, and zero otherwise.

Solving the Cooper Equation

With the simplified interaction, all $g_{\mathbf{k}}$ within the Debye window are equal, and the self-consistency condition yields:

$$1 = V \int_0^{\hbar\omega_D} \frac{N(0)\, d\epsilon}{2\epsilon + |E_b|}$$

where $E_b = 2E_F - E > 0$ is the binding energy measured from the Fermi level. Evaluating the integral gives:

$$1 = \frac{N(0)V}{2} \ln\!\left(\frac{2\hbar\omega_D + |E_b|}{|E_b|}\right)$$

Cooper Pair Binding Energy

In the weak-coupling limit $N(0)V \ll 1$, the binding energy$|E_b| \ll \hbar\omega_D$, and the logarithm simplifies. The resulting binding energy is:

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

This is a remarkable result with several important features:

●The binding energy is non-analytic in the coupling $V$: it cannot be expanded as a Taylor series around $V = 0$. This means Cooper pairing cannot be found by perturbation theory to any finite order — it is an essentially non-perturbative effect.
●The energy scale is set by the Debye energy $\hbar\omega_D$, confirming the role of phonons in mediating the attraction.
●For typical metals, $N(0)V \sim 0.2\text{--}0.3$ and$\hbar\omega_D \sim 10\text{--}30$ meV, giving$E_b \sim 10^{-4}\text{--}10^{-1}$ meV, far smaller than the Fermi energy.
●A bound state exists for any$V > 0$, confirming the Cooper instability. There is no threshold.

Phonon-Mediated Attraction

The attractive interaction between electrons is mediated by the exchange of virtual phonons. The physical picture is as follows: an electron with momentum$\mathbf{k}$ polarizes the lattice as it moves through the crystal, creating a region of positive charge density. A second electron with momentum$-\mathbf{k}$ is attracted to this positive charge excess.

The key feature is that this is a retarded interaction. The lattice responds on a timescale$\tau_{\text{ph}} \sim 1/\omega_D \sim 10^{-13}$ s, while the electron traverses the distortion in a time$\tau_e \sim 1/\omega_F \sim 10^{-16}$ s. By the time the lattice has relaxed to the polarized configuration, the first electron is far away, and the second electron sees only the attractive potential without the Coulomb repulsion from the first.

Effective Electron-Electron Interaction

In the jellium model with electron-phonon coupling, the effective interaction is:

$$V_{\text{eff}}(\mathbf{q}, \omega) = \frac{4\pi e^2}{q^2 \epsilon(\mathbf{q},\omega)} = V_C(q)\left[1 + \frac{\omega_q^2}{\omega^2 - \omega_q^2}\right]$$

For $|\omega| < \omega_q$ (energy transfer less than the phonon energy), $V_{\text{eff}}$ becomes negative — the interaction is attractive. This overcomes the bare Coulomb repulsion in a narrow energy shell of width $\sim 2\hbar\omega_D$ around the Fermi surface.

The net effective interaction is characterized by the dimensionless coupling constant $\lambda = N(0)V_{\text{ph}} - \mu^*$, where$\mu^* \approx 0.1\text{--}0.15$ is the screened Coulomb pseudopotential. Superconductivity occurs when $\lambda > 0$, i.e., the phonon-mediated attraction exceeds the residual Coulomb repulsion.

Pair Size: The Coherence Length

The spatial extent of the Cooper pair wavefunction defines the BCS coherence length$\xi_0$. It is determined by the uncertainty principle: the pair occupies a momentum shell of width $\delta k \sim \Delta_0 / (\hbar v_F)$around $k_F$, giving a real-space extent:

$$\xi_0 = \frac{\hbar v_F}{\pi \Delta_0}$$

where $v_F$ is the Fermi velocity and $\Delta_0$ is the BCS gap at zero temperature. For conventional superconductors:

Material	$T_c$ (K)	$\Delta_0$ (meV)	$\xi_0$ (nm)	$\xi_0 / a$
Aluminum (Al)	1.18	0.18	~1600	~4000
Niobium (Nb)	9.25	1.55	~38	~115
Lead (Pb)	7.19	1.35	~83	~170

The ratio $\xi_0 / a$ (where $a$ is the lattice constant) is enormous: Cooper pairs in conventional superconductors extend over hundreds to thousands of lattice spacings. This means millions of pairs overlap in any given region — a crucial feature that makes the BCS mean-field theory self-consistent.

Spin-Singlet vs Spin-Triplet Pairing

The total pair wavefunction must be antisymmetric under exchange of the two electrons (Pauli principle). Since the wavefunction factorizes as$\Psi = \psi_{\text{orbital}} \otimes \chi_{\text{spin}}$, we have two possibilities:

Spin-Singlet Pairing

$$\chi_s = \frac{1}{\sqrt{2}}(|\!\uparrow\downarrow\rangle - |\!\downarrow\uparrow\rangle)$$

Spin part is antisymmetric ($S = 0$), so the orbital part must be symmetric: even parity ($\ell = 0, 2, \ldots$). The conventional s-wave ($\ell = 0$) BCS state belongs here. Most conventional superconductors (Al, Nb, Pb, Sn) are s-wave singlet.

Spin-Triplet Pairing

$$\chi_t \in \{|\!\uparrow\uparrow\rangle,\; \frac{|\!\uparrow\downarrow\rangle + |\!\downarrow\uparrow\rangle}{\sqrt{2}},\; |\!\downarrow\downarrow\rangle\}$$

Spin part is symmetric ($S = 1$), so the orbital part must be antisymmetric: odd parity ($\ell = 1, 3, \ldots$). Superfluid$^3$He is a p-wave triplet superfluid. Some unconventional superconductors (e.g., Sr$_2$RuO$_4$) may be triplet.

The symmetry of the gap function $\Delta(\mathbf{k})$ on the Fermi surface reflects the orbital angular momentum: isotropic for s-wave, with nodes for d-wave (as in the cuprate high-$T_c$ superconductors where$\Delta(\mathbf{k}) \propto \cos k_x - \cos k_y$).

Why Pairs Have Zero Total Momentum

Cooper showed that the maximum binding energy is achieved when the pair has zero total momentum, $\mathbf{K} = \mathbf{k}_1 + \mathbf{k}_2 = 0$. The physical argument involves phase space: at $\mathbf{K} = 0$, both electrons in the pair $(\mathbf{k}, -\mathbf{k})$ can scatter to any other pair state $(\mathbf{k}', -\mathbf{k}')$ on the Fermi surface.

For $\mathbf{K} \neq 0$, the constraint$\mathbf{k}_1 + \mathbf{k}_2 = \mathbf{K}$ means that both$\mathbf{k}_1$ and $\mathbf{k}_2$ must lie above the Fermi surface simultaneously. The overlap region of the two displaced Fermi spheres (centered at $\mathbf{K}/2$ and $-\mathbf{K}/2$) shrinks as $|\mathbf{K}|$ increases, reducing the available phase space and hence the binding energy.

Phase Space Argument

The number of available scattering states scales as the overlap volume of two Fermi spheres displaced by $\mathbf{K}$. For$|\mathbf{K}| = 0$, the overlap is maximum (the full Fermi surface shell of width $\hbar\omega_D$). As $|\mathbf{K}|$increases, the overlap decreases, and the effective density of states available for pairing diminishes. Since the binding energy depends exponentially on the density of states, $E_b \propto \exp(-2/N_{\text{eff}}V)$, even a small reduction in $N_{\text{eff}}$ dramatically suppresses binding.

Structure of the Pair Wavefunction

The real-space pair wavefunction for an s-wave Cooper pair takes the form:

$$\psi(r) \sim \frac{\sin(k_F r)}{k_F r}\, e^{-r/\xi_0}$$

This consists of rapid oscillations at the Fermi wavelength$\lambda_F = 2\pi / k_F \sim 0.5$ nm, modulated by an exponential envelope with decay length $\xi_0$. The oscillations reflect the fact that $g_{\mathbf{k}}$ is peaked near $k_F$, while the envelope reflects the narrow spread $\delta k \sim 1/\xi_0$ of the pair amplitudes in momentum space.

The pair amplitude in momentum space, $g_{\mathbf{k}}$, is nonzero only in a thin shell of width $\sim \hbar\omega_D$ above the Fermi surface. The solution of the Cooper equation gives:

$$g_{\mathbf{k}} = \frac{C}{2\epsilon_{\mathbf{k}} + |E_b|}$$

for $0 < \epsilon_{\mathbf{k}} - E_F < \hbar\omega_D$, where$C$ is a normalization constant. The amplitude is largest for states just above $E_F$ and decreases as $1/\epsilon$ toward the Debye cutoff.

Detailed Derivation: The Cooper Equation

We derive the Cooper equation starting from the full two-body Schrodinger equation for two electrons above a filled Fermi sea. The total Hamiltonian is:

$$H = \frac{\mathbf{p}_1^2}{2m} + \frac{\mathbf{p}_2^2}{2m} + V(\mathbf{r}_1 - \mathbf{r}_2)$$

Step 1: Expand in Pair States

For zero total momentum, we write the pair wavefunction as a superposition of plane-wave pair states with opposite momenta:

$$\psi(\mathbf{r}_1, \mathbf{r}_2) = \sum_{|\mathbf{k}| > k_F} g_{\mathbf{k}}\, e^{i\mathbf{k}\cdot\mathbf{r}_1}\, e^{-i\mathbf{k}\cdot\mathbf{r}_2}$$

The restriction $|\mathbf{k}| > k_F$ enforces the Pauli exclusion principle: all states below the Fermi surface are occupied and unavailable.

Step 2: Substitute into Schrodinger Equation

Acting with $H$ on $\psi$, the kinetic energy terms give:

$$\left(\frac{\hbar^2 k^2}{2m} + \frac{\hbar^2 k^2}{2m}\right)g_{\mathbf{k}} = 2\epsilon_{\mathbf{k}}\, g_{\mathbf{k}}$$

The potential term involves the matrix element of $V$ between plane-wave pair states. Multiplying both sides by $e^{-i\mathbf{k}'\cdot\mathbf{r}_1} e^{i\mathbf{k}'\cdot\mathbf{r}_2}$and integrating over $\mathbf{r}_1$ and $\mathbf{r}_2$, we obtain:

$$E\, g_{\mathbf{k}} = 2\epsilon_{\mathbf{k}}\, g_{\mathbf{k}} + \sum_{\mathbf{k}'} V_{\mathbf{k}\mathbf{k}'}\, g_{\mathbf{k}'}$$

where $V_{\mathbf{k}\mathbf{k}'} = \frac{1}{\Omega}\int V(\mathbf{r})\, e^{i(\mathbf{k}-\mathbf{k}')\cdot\mathbf{r}}\, d^3r$is the Fourier transform of the interaction potential and $\Omega$ is the system volume.

Step 3: Simplify with Model Interaction

Cooper's key simplification is to take a separable, constant interaction within a Debye shell:$V_{\mathbf{k}\mathbf{k}'} = -V$ for $E_F < \epsilon_{\mathbf{k}}, \epsilon_{\mathbf{k}'} < E_F + \hbar\omega_D$, and zero otherwise. Then the Cooper equation becomes:

$$(E - 2\epsilon_{\mathbf{k}})\, g_{\mathbf{k}} = -V \sum_{\mathbf{k}'} g_{\mathbf{k}'}$$

Since the right-hand side is independent of $\mathbf{k}$, we can define$\Lambda = \sum_{\mathbf{k}'} g_{\mathbf{k}'}$ and solve for $g_{\mathbf{k}}$:

$$g_{\mathbf{k}} = \frac{-V\Lambda}{E - 2\epsilon_{\mathbf{k}}} = \frac{V\Lambda}{2\epsilon_{\mathbf{k}} - E}$$

Now sum both sides over $\mathbf{k}$ and cancel $\Lambda$:

$$1 = V \sum_{\mathbf{k}} \frac{1}{2\epsilon_{\mathbf{k}} - E}$$

Step 4: Convert Sum to Integral and Solve

Converting the $\mathbf{k}$-sum to an energy integral using the density of states $N(\epsilon) \approx N(0)$ (constant near $E_F$), and measuring energies from $E_F$ so that $\epsilon = \epsilon_{\mathbf{k}} - E_F$ and$E = 2E_F - |E_b|$:

$$1 = N(0)V \int_0^{\hbar\omega_D} \frac{d\epsilon}{2\epsilon + |E_b|}$$

Evaluating the integral:

$$1 = \frac{N(0)V}{2}\left[\ln(2\epsilon + |E_b|)\right]_0^{\hbar\omega_D} = \frac{N(0)V}{2}\ln\!\left(\frac{2\hbar\omega_D + |E_b|}{|E_b|}\right)$$

In the weak-coupling limit $N(0)V \ll 1$, we have $|E_b| \ll \hbar\omega_D$, so $2\hbar\omega_D + |E_b| \approx 2\hbar\omega_D$. Exponentiating:

$$\frac{2\hbar\omega_D}{|E_b|} = e^{2/N(0)V} \quad \Longrightarrow \quad |E_b| = 2\hbar\omega_D\, e^{-2/N(0)V}$$

Derivation: Effective Electron-Electron Interaction

The attractive interaction between electrons arises from virtual phonon exchange. Starting from the Frohlich electron-phonon Hamiltonian, we derive the effective interaction using second-order perturbation theory.

Step 1: The Frohlich Hamiltonian

The electron-phonon interaction takes the form:

$$H_{ep} = \sum_{\mathbf{k},\mathbf{q}} M_{\mathbf{q}} \left(a_{\mathbf{q}} + a_{-\mathbf{q}}^\dagger\right) c_{\mathbf{k}+\mathbf{q}}^\dagger c_{\mathbf{k}}$$

where $M_{\mathbf{q}}$ is the electron-phonon matrix element, $a_{\mathbf{q}}$annihilates a phonon with wavevector $\mathbf{q}$, and $c_{\mathbf{k}}$annihilates an electron with wavevector $\mathbf{k}$.

Step 2: Second-Order Perturbation Theory

Consider an electron in state $\mathbf{k}_1$ scattering to $\mathbf{k}_1 + \mathbf{q}$by emitting a virtual phonon $\mathbf{q}$, which is then absorbed by a second electron scattering from $\mathbf{k}_2$ to $\mathbf{k}_2 - \mathbf{q}$. The second-order matrix element is:

$$V_{\text{eff}}(\mathbf{q}, \omega) = \frac{|M_{\mathbf{q}}|^2}{\omega - \omega_{\mathbf{q}}} + \frac{|M_{\mathbf{q}}|^2}{-\omega - \omega_{\mathbf{q}}} = \frac{2|M_{\mathbf{q}}|^2 \omega_{\mathbf{q}}}{\omega^2 - \omega_{\mathbf{q}}^2}$$

where $\omega = \epsilon_{\mathbf{k}_1} - \epsilon_{\mathbf{k}_1+\mathbf{q}}$ is the energy transfer. The two terms correspond to two time orderings of the phonon exchange process.

Step 3: Condition for Attraction

From the expression $V_{\text{eff}} = 2|M_{\mathbf{q}}|^2\omega_{\mathbf{q}}/(\omega^2 - \omega_{\mathbf{q}}^2)$, the sign depends on the energy transfer relative to the phonon frequency:

●When $|\omega| < \omega_{\mathbf{q}}$: $\omega^2 < \omega_{\mathbf{q}}^2$, so $V_{\text{eff}} < 0$ — the interaction is attractive.
●When $|\omega| > \omega_{\mathbf{q}}$: $\omega^2 > \omega_{\mathbf{q}}^2$, so $V_{\text{eff}} > 0$ — the interaction is repulsive.

For electrons near the Fermi surface, the energy transfer $|\omega|$ is small (of order $k_BT$ or $\Delta$), while $\omega_{\mathbf{q}} \sim \omega_D$. Since $\Delta \ll \hbar\omega_D$, the interaction is attractive in the pairing channel. This is Cooper's key insight: the Fermi sea restricts the available phase space to a thin shell where phonon-mediated attraction dominates.

Derivation: Cooper Pair Size

The real-space size of a Cooper pair follows from a Fourier transform argument applied to the pair amplitude in momentum space.

From k-space to Real Space

The pair wavefunction in real space is:

$$\psi(\mathbf{r}) = \sum_{\mathbf{k}} g_{\mathbf{k}}\, e^{i\mathbf{k}\cdot\mathbf{r}}$$

Since $g_{\mathbf{k}}$ is peaked at $|\mathbf{k}| = k_F$ and has a width $\delta k$ in momentum space, the Fourier transform gives oscillations at wavelength $\lambda_F = 2\pi/k_F$ modulated by an envelope of width $\sim 1/\delta k$.

The momentum spread is determined by the energy range over which $g_{\mathbf{k}}$is appreciable. From $g_{\mathbf{k}} = C/(2\epsilon_{\mathbf{k}} + |E_b|)$, the amplitude falls to half its maximum value when $\epsilon_{\mathbf{k}} \sim |E_b|/2 \sim \Delta_0$. Using $\epsilon_{\mathbf{k}} = \hbar v_F \delta k$ (linearizing near $k_F$):

$$\delta k \sim \frac{\Delta_0}{\hbar v_F} \quad \Longrightarrow \quad \xi_0 \sim \frac{1}{\delta k} = \frac{\hbar v_F}{\Delta_0}$$

The exact BCS calculation gives $\xi_0 = \hbar v_F / (\pi\Delta_0)$, differing from the estimate by a factor of $\pi$. For aluminum: $v_F = 2.03 \times 10^6$ m/s, $\Delta_0 = 0.18$ meV, giving $\xi_0 \approx 1600$ nm — about 4000 lattice constants.

Historical Context

The Discovery of Electron Pairing

Leon Cooper published his landmark paper “Bound Electron Pairs in a Degenerate Fermi Gas” in Physical Review in 1956. At the time, he was a postdoctoral researcher working with John Bardeen at the University of Illinois at Urbana-Champaign. The problem had been considered intractable: how could electrons, which repel each other via the Coulomb interaction, form bound pairs? Cooper's insight was that the filled Fermi sea qualitatively changes the two-body problem. The Pauli exclusion principle restricts the available states to a thin shell above$E_F$, effectively reducing the dimensionality and allowing a bound state for any attractive coupling.

The phonon-mediated attractive interaction had been suggested earlier by Herbert Frohlich (1950) based on the isotope effect discovered independently by Maxwell and by Reynolds, Serin, Wright, and Nesbitt (1950). The isotope effect showed that $T_c \propto M^{-1/2}$, directly implicating lattice vibrations. Cooper's calculation demonstrated that this attraction, combined with the Fermi sea, inevitably leads to pairing — setting the stage for the full BCS theory published one year later in 1957.

Nobel Prize Recognition

John Bardeen, Leon Cooper, and J. Robert Schrieffer received the Nobel Prize in Physics in 1972 “for their jointly developed theory of superconductivity, usually called the BCS theory.” Bardeen remains the only person to have won two Nobel Prizes in Physics (the first in 1956 for the invention of the transistor). Cooper's contribution of electron pairing was the conceptual breakthrough that made the full microscopic theory possible.

Applications of Cooper Pairing

Real-World Applications

1. Superconducting Magnets (MRI and Particle Accelerators)

Cooper pairing in NbTi and Nb$_3$Sn enables zero-resistance current flow in high-field magnets. MRI scanners use superconducting solenoids producing 1.5–7 T fields, while the LHC at CERN uses 1232 NbTi dipole magnets operating at 8.3 T and 1.9 K.

2. Superconducting Quantum Computing

The macroscopic coherence of Cooper pairs allows quantum information to be encoded in the phase difference across Josephson junctions. Transmon qubits in systems from IBM, Google, and others exploit the nonlinear inductance provided by Cooper pair tunneling to create anharmonic quantum oscillators.

3. Superconducting Power Transmission

Cooper pairing eliminates resistive losses in power cables. Projects in cities like Essen (Germany) and Long Island (USA) have demonstrated superconducting cables carrying 5–10 times the current of conventional cables with zero resistive dissipation.

4. Particle Detectors (STJ Detectors)

Superconducting tunnel junctions exploit the small Cooper pair binding energy ($\sim$ meV) for single-photon detection from infrared to X-ray wavelengths. Breaking a Cooper pair creates quasiparticles that are detected as a tunneling current, providing energy resolution orders of magnitude better than semiconductor detectors.

5. Unconventional Pairing in Exotic Materials

Beyond conventional s-wave pairing, d-wave Cooper pairs in cuprate high-$T_c$superconductors and p-wave triplet pairs in Sr$_2$RuO$_4$ and superfluid $ {}^3$He have opened frontiers in topological quantum matter. Topological superconductors with p-wave pairing may host Majorana fermions useful for fault-tolerant quantum computing.

Computational Simulations

Cooper Pair Binding Energy & Wavefunction

Python

script.py129 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
# numpy-only (no scipy needed)

# ──────────────────────────────────────────────
# Cooper Pair: Binding Energy vs Coupling Strength
# and Pair Wavefunction in Real Space
# ──────────────────────────────────────────────

# Physical constants
hbar = 1.0545718e-34   # J·s
kB = 1.380649e-23      # J/K
eV = 1.602176634e-19   # J

# Parameters for a model superconductor
omega_D = 0.03 * eV    # Debye energy ~ 30 meV
E_F = 7.0 * eV         # Fermi energy ~ 7 eV
v_F = 1.5e6            # Fermi velocity (m/s), typical for metals
k_F = E_F / (hbar * v_F)  # Fermi wavevector

# ── Panel 1: Binding energy vs coupling N(0)V ──
N0V_vals = np.linspace(0.05, 0.8, 200)

# Cooper pair binding energy: E_b = 2*hbar*omega_D * exp(-2 / N(0)V)
# (valid in weak-coupling limit)
E_b = 2.0 * omega_D * np.exp(-2.0 / N0V_vals)
E_b_meV = E_b / eV * 1000  # convert to meV

# Also compute the BCS gap Delta_0 for comparison
# Delta_0 = 2*hbar*omega_D * exp(-1 / N(0)V) in simple BCS
Delta_0 = 2.0 * omega_D * np.exp(-1.0 / N0V_vals)
Delta_0_meV = Delta_0 / eV * 1000

fig, axes = plt.subplots(2, 2, figsize=(14, 11))
fig.patch.set_facecolor('#0a0a1a')
for ax in axes.flat:
    ax.set_facecolor('#0a0a1a')
    for spine in ax.spines.values():
        spine.set_edgecolor('#334155')
    ax.tick_params(colors='#94a3b8')

# Panel 1: Binding energy vs N(0)V
axes[0,0].semilogy(N0V_vals, E_b_meV, color='#2dd4bf', linewidth=2.5,
                    label=r'$E_b = 2\hbar\omega_D \exp(-2/N(0)V)$')
axes[0,0].semilogy(N0V_vals, Delta_0_meV, color='#f472b6', linewidth=2,
                    linestyle='--', label=r'$\Delta_0$ (BCS gap)')
axes[0,0].set_xlabel('Coupling strength N(0)V', color='#94a3b8', fontsize=12)
axes[0,0].set_ylabel('Energy (meV)', color='#94a3b8', fontsize=12)
axes[0,0].set_title('Cooper Pair Binding Energy', color='white', fontsize=14)
axes[0,0].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155',
                  labelcolor='white')
axes[0,0].set_ylim(1e-4, 1e3)
axes[0,0].axvline(x=0.3, color='#fbbf24', linestyle=':', alpha=0.6, linewidth=1)
axes[0,0].text(0.32, 1e2, 'Typical\nmetals', color='#fbbf24', fontsize=9)

# Panel 2: Pair wavefunction g(k) near Fermi surface
# In the Cooper problem, g_k ~ 1/(2*epsilon_k - E) for |epsilon_k| < hbar*omega_D
N0V_ex = 0.3
E_bind = 2.0 * omega_D * np.exp(-2.0 / N0V_ex)
epsilon_k = np.linspace(-omega_D * 1.5, omega_D * 1.5, 1000)
# g_k is nonzero only for 0 < epsilon_k < omega_D (states above Fermi sea)
g_k = np.zeros_like(epsilon_k)
mask = (epsilon_k > 0) & (epsilon_k < omega_D)
g_k[mask] = 1.0 / (2.0 * epsilon_k[mask] + E_bind)
g_k /= np.max(g_k) if np.max(g_k) > 0 else 1.0

axes[0,1].fill_between(epsilon_k/eV*1000, 0, g_k, alpha=0.3, color='#2dd4bf')
axes[0,1].plot(epsilon_k/eV*1000, g_k, color='#2dd4bf', linewidth=2)
axes[0,1].axvline(x=0, color='#f472b6', linewidth=1.5, linestyle='--',
                   label=r'$E_F$ (Fermi level)')
axes[0,1].axvline(x=omega_D/eV*1000, color='#fbbf24', linewidth=1, linestyle=':',
                   label=r'$\hbar\omega_D$')
axes[0,1].set_xlabel(r'$\epsilon_k - E_F$ (meV)', color='#94a3b8', fontsize=12)
axes[0,1].set_ylabel(r'$|g_k|$ (normalized)', color='#94a3b8', fontsize=12)
axes[0,1].set_title('Pair Amplitude in k-space', color='white', fontsize=14)
axes[0,1].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155',
                  labelcolor='white')

# Panel 3: Pair wavefunction in real space
# psi(r) = sum_k g_k * exp(ik*r) ~ integral of g(epsilon) * cos(k_F*r + ...)
# For an s-wave pair, psi(r) ~ sin(k_F * r) * exp(-r / xi_0) / (k_F * r)
N0V_vals_r = [0.2, 0.3, 0.5]
r_vals = np.linspace(0.1, 500, 2000)  # in nm

for N0V_r in N0V_vals_r:
    Delta_r = 2.0 * omega_D * np.exp(-1.0 / N0V_r)
    xi_0 = hbar * v_F / (np.pi * Delta_r)  # coherence length in meters
    xi_nm = xi_0 * 1e9  # convert to nm
    # Pair wavefunction envelope
    psi = np.sin(k_F * r_vals * 1e-9) * np.exp(-r_vals / xi_nm) / (k_F * r_vals * 1e-9)
    psi /= np.max(np.abs(psi))
    label = f'N(0)V={N0V_r}, ' + r'$\xi_0$' + f'={xi_nm:.0f} nm'
    axes[1,0].plot(r_vals, np.abs(psi)**2, linewidth=2, label=label)

axes[1,0].set_xlabel('r (nm)', color='#94a3b8', fontsize=12)
axes[1,0].set_ylabel(r'$|\psi(r)|^2$ (normalized)', color='#94a3b8', fontsize=12)
axes[1,0].set_title('Cooper Pair Wavefunction (Real Space)', color='white', fontsize=14)
axes[1,0].legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#334155',
                  labelcolor='white')
axes[1,0].set_xlim(0, 500)

# Panel 4: Coherence length vs coupling
N0V_xi = np.linspace(0.1, 0.8, 200)
Delta_xi = 2.0 * omega_D * np.exp(-1.0 / N0V_xi)
xi_0_vals = hbar * v_F / (np.pi * Delta_xi) * 1e9  # in nm

axes[1,1].semilogy(N0V_xi, xi_0_vals, color='#a78bfa', linewidth=2.5)
axes[1,1].set_xlabel('Coupling strength N(0)V', color='#94a3b8', fontsize=12)
axes[1,1].set_ylabel(r'Coherence length $\xi_0$ (nm)', color='#94a3b8', fontsize=12)
axes[1,1].set_title('Pair Size vs Coupling', color='white', fontsize=14)
axes[1,1].axhline(y=1600, color='#22d3ee', linewidth=1, linestyle=':', alpha=0.7)
axes[1,1].text(0.55, 2000, 'Al (~1600 nm)', color='#22d3ee', fontsize=9)
axes[1,1].axhline(y=38, color='#f472b6', linewidth=1, linestyle=':', alpha=0.7)
axes[1,1].text(0.55, 50, 'Nb (~38 nm)', color='#f472b6', fontsize=9)

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

# Print summary
print("=== Cooper Pair Analysis ===")
for N0V in [0.2, 0.3, 0.4, 0.5]:
    Eb = 2.0 * omega_D * np.exp(-2.0 / N0V) / eV * 1000
    D0 = 2.0 * omega_D * np.exp(-1.0 / N0V) / eV * 1000
    xi = hbar * v_F / (np.pi * D0 * eV / 1000) * 1e9
    print(f"N(0)V = {N0V:.1f}: E_b = {Eb:.4f} meV, Delta_0 = {D0:.4f} meV, xi_0 = {xi:.1f} nm")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Cooper Pair Properties: Al, Nb, Pb

Fortran

program.f90172 lines

program cooper_pair_materials
  implicit none

! ──────────────────────────────────────────────
  ! Cooper Pair Size and Binding Energy
  ! for Typical Superconductors: Al, Nb, Pb
  ! ──────────────────────────────────────────────

integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: n_materials = 3
  integer :: i, j

! Physical constants
  real(dp), parameter :: hbar = 1.0545718d-34    ! J·s
  real(dp), parameter :: kB   = 1.380649d-23     ! J/K
  real(dp), parameter :: eV   = 1.602176634d-19  ! J
  real(dp), parameter :: pi   = 3.14159265358979323846d0

! Material properties arrays
  character(len=4) :: names(n_materials)
  real(dp) :: Tc(n_materials)          ! Critical temperature (K)
  real(dp) :: theta_D(n_materials)     ! Debye temperature (K)
  real(dp) :: v_F(n_materials)         ! Fermi velocity (m/s)
  real(dp) :: N0(n_materials)          ! DOS at Fermi level (states/eV/atom)
  real(dp) :: lambda_ep(n_materials)   ! electron-phonon coupling

! Computed quantities
  real(dp) :: Delta_0, xi_0, omega_D, N0V
  real(dp) :: E_bind, Tc_calc, ratio
  real(dp) :: Delta_meV, xi_nm, Eb_meV

! Temperature-dependent gap
  integer, parameter :: n_temp = 50
  real(dp) :: T, t_red, Delta_T, gap_ratio

! ── Material data ──
  ! Aluminum
  names(1) = 'Al'
  Tc(1) = 1.18d0          ! K
  theta_D(1) = 428.0d0    ! K
  v_F(1) = 2.03d6         ! m/s
  N0(1) = 0.41d0          ! states/eV/spin/atom
  lambda_ep(1) = 0.43d0   ! electron-phonon coupling

! Niobium
  names(2) = 'Nb'
  Tc(2) = 9.25d0          ! K
  theta_D(2) = 275.0d0    ! K
  v_F(2) = 0.61d6         ! m/s
  N0(2) = 1.83d0
  lambda_ep(2) = 1.26d0

! Lead
  names(3) = 'Pb'
  Tc(3) = 7.19d0          ! K
  theta_D(3) = 105.0d0    ! K
  v_F(3) = 1.83d6         ! m/s
  N0(3) = 0.49d0
  lambda_ep(3) = 1.55d0

! ── Header ──
  write(*,'(A)') '============================================================='
  write(*,'(A)') '     Cooper Pair Properties of Conventional Superconductors'
  write(*,'(A)') '============================================================='
  write(*,*)

! ── Main computation loop ──
  do i = 1, n_materials
    write(*,'(A,A4,A)') '─── ', names(i), ' ───────────────────────────────────────'

! Debye energy
    omega_D = kB * theta_D(i)

! BCS gap at T=0: Delta_0 = 1.764 * kB * Tc (weak-coupling BCS)
    Delta_0 = 1.764d0 * kB * Tc(i)
    Delta_meV = Delta_0 / eV * 1.0d3

! Effective coupling N(0)V from BCS relation
    ! Delta_0 = 2*omega_D * exp(-1/N(0)V)
    ! => N(0)V = -1 / ln(Delta_0 / (2*omega_D))
    N0V = -1.0d0 / log(Delta_0 / (2.0d0 * omega_D))

! Cooper pair binding energy (isolated pair above Fermi sea)
    E_bind = 2.0d0 * omega_D * exp(-2.0d0 / N0V)
    Eb_meV = E_bind / eV * 1.0d3

! BCS coherence length
    xi_0 = hbar * v_F(i) / (pi * Delta_0)
    xi_nm = xi_0 * 1.0d9

! 2*Delta_0 / (kB * Tc) ratio (BCS predicts 3.528)
    ratio = 2.0d0 * Delta_0 / (kB * Tc(i))

! Number of pairs in coherence volume (3D)
    ! n_pairs ~ N(0) * Delta_0 * (4/3)*pi*xi_0^3

write(*,'(A,F8.2,A)')    '  Tc            = ', Tc(i), ' K'
    write(*,'(A,F8.1,A)')    '  Theta_D       = ', theta_D(i), ' K'
    write(*,'(A,F8.4,A)')    '  N(0)V         = ', N0V
    write(*,'(A,F8.4,A)')    '  lambda_ep     = ', lambda_ep(i)
    write(*,'(A,ES12.4,A)')  '  Delta_0       = ', Delta_meV, ' meV'
    write(*,'(A,ES12.4,A)')  '  E_bind(Cooper)= ', Eb_meV, ' meV'
    write(*,'(A,F8.1,A)')    '  xi_0          = ', xi_nm, ' nm'
    write(*,'(A,F8.4)')      '  2*Delta/kBTc  = ', ratio
    write(*,'(A,ES12.4,A)')  '  v_F           = ', v_F(i), ' m/s'
    write(*,*)
  end do

! ── Comparison table ──
  write(*,'(A)') '============================================================='
  write(*,'(A)') '  Summary Comparison Table'
  write(*,'(A)') '============================================================='
  write(*,'(A4, A10, A12, A12, A14, A12)') &
    'Mat', 'Tc(K)', 'Delta(meV)', 'xi_0(nm)', 'E_bind(meV)', 'N(0)V'
  write(*,'(A)') '-------------------------------------------------------------'

do i = 1, n_materials
    omega_D = kB * theta_D(i)
    Delta_0 = 1.764d0 * kB * Tc(i)
    Delta_meV = Delta_0 / eV * 1.0d3
    N0V = -1.0d0 / log(Delta_0 / (2.0d0 * omega_D))
    E_bind = 2.0d0 * omega_D * exp(-2.0d0 / N0V)
    Eb_meV = E_bind / eV * 1.0d3
    xi_0 = hbar * v_F(i) / (pi * Delta_0)
    xi_nm = xi_0 * 1.0d9

write(*,'(A4, F10.2, F12.4, F12.1, ES14.4, F12.4)') &
      names(i), Tc(i), Delta_meV, xi_nm, Eb_meV, N0V
  end do

write(*,*)

! ── Temperature dependence of the gap (approximate BCS) ──
  write(*,'(A)') '============================================================='
  write(*,'(A)') '  BCS Gap Temperature Dependence (Aluminum)'
  write(*,'(A)') '============================================================='
  write(*,'(A10, A15, A15)') 'T/Tc', 'Delta/Delta_0', 'Delta (meV)'

Delta_0 = 1.764d0 * kB * Tc(1)
  Delta_meV = Delta_0 / eV * 1.0d3

do j = 0, n_temp
    t_red = dble(j) / dble(n_temp)
    T = t_red * Tc(1)

if (t_red >= 1.0d0) then
      gap_ratio = 0.0d0
    else if (t_red < 0.01d0) then
      gap_ratio = 1.0d0
    else
      ! Approximate BCS gap: Delta(T)/Delta(0) ~ sqrt(1 - (T/Tc)^3) * tanh(...)
      ! Using Muhlschlegel's approximation
      gap_ratio = tanh(1.74d0 * sqrt((1.0d0/t_red) - 1.0d0))
    end if

Delta_T = gap_ratio * Delta_meV

if (mod(j, 5) == 0) then
      write(*,'(F10.4, F15.6, F15.6)') t_red, gap_ratio, Delta_T
    end if
  end do

write(*,*)
  write(*,'(A)') '============================================================='
  write(*,'(A)') '  Key Physics:'
  write(*,'(A)') '  - Cooper pair size >> interparticle spacing (BCS limit)'
  write(*,'(A)') '  - E_bind << Delta_0 (many-body gap exceeds pair binding)'
  write(*,'(A)') '  - xi_0(Al) ~ 1600 nm: pairs overlap enormously'
  write(*,'(A)') '  - xi_0(Pb) ~ 83 nm: stronger coupling, smaller pairs'
  write(*,'(A)') '============================================================='

end program cooper_pair_materials

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Part III Overview Chapter 2: BCS Theory →