Chapter 4: BCS Theory In Depth

Learning Objectives

• Understand the Fermi liquid framework and quasiparticle description of normal metals
• Derive the Cooper instability and show bound state formation for weak attractive interactions
• Construct the Cooper pair wave function and analyze pairing symmetries (singlet vs triplet, s-wave)
• Master the BCS variational ground state and coherence factors $u_k$, $v_k$
• Solve the BCS gap equation numerically and verify universal ratios
• Calculate temperature-dependent gap $\Delta(T)$ and understand phase transition behavior
• Derive the superconducting density of states with coherence peaks
• Analyze tunneling spectroscopy as experimental verification of BCS theory
• Implement comprehensive Python simulations of BCS phenomena

4.1 Normal Metal Fermi Liquid Theory

The Fermi Sea at Zero Temperature

Before understanding superconductivity, we must review the normal metallic state. In a non-interacting electron gas (Sommerfeld model), electrons fill single-particle states up to the Fermi energy $E_F$:

where the single-particle energy is:

The Fermi Surface

The Fermi surface is the boundary in momentum space separating occupied and unoccupied states:

The density of states at the Fermi surface per spin is:

Landau Fermi Liquid Theory

When electron-electron interactions are included (Coulomb repulsion, screening), Landau's Fermi liquid theory tells us that:

• Low-energy excitations are quasiparticles (dressed electrons with renormalized mass $m^*$)
• Quasiparticles have one-to-one correspondence with non-interacting electrons
• The Fermi surface survives (though shape may be modified)
• Lifetime of quasiparticles $\tau \propto (\epsilon - E_F)^{-2}$ (long-lived near Fermi surface)

Fermi-Dirac Distribution at Finite Temperature

At temperature $T > 0$, the occupation probability is:

where $\mu \approx E_F$ is the chemical potential. The distribution has a thermal width $\sim k_B T$ around $E_F$.

import numpy as np
import matplotlib.pyplot as plt

# Fermi-Dirac distribution
def fermi_dirac(epsilon, mu, T, k_B=1.0):
    """
    Calculate Fermi-Dirac distribution.

    Parameters:
    -----------
    epsilon : array
        Energy values
    mu : float
        Chemical potential (Fermi energy)
    T : float
        Temperature (in energy units if k_B=1)
    k_B : float
        Boltzmann constant (default 1.0)
    """
    if T == 0:
        return (epsilon < mu).astype(float)
    return 1.0 / (np.exp((epsilon - mu) / (k_B * T)) + 1.0)

# Energy range around Fermi surface
epsilon = np.linspace(-2, 2, 1000)  # In units of E_F
E_F = 0.0  # Set Fermi energy to zero
k_B = 1.0

# Different temperatures
temperatures = [0.01, 0.05, 0.1, 0.2, 0.5]

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Left: Fermi-Dirac distributions
ax1 = axes[0]
for T in temperatures:
    f = fermi_dirac(epsilon, E_F, T, k_B)
    label = f'T = {T:.2f} E_F/k_B' if T > 0 else 'T = 0'
    ax1.plot(epsilon, f, linewidth=2, label=label)

ax1.axvline(x=E_F, color='k', linestyle='--', linewidth=1.5, alpha=0.5, label='E_F')
ax1.axhline(y=0.5, color='gray', linestyle=':', linewidth=1, alpha=0.5)
ax1.set_xlabel(r'Energy $(\epsilon - E_F) / E_F$', fontsize=13)
ax1.set_ylabel(r'Occupation $f(\epsilon)$', fontsize=13)
ax1.set_title('Fermi-Dirac Distribution', fontsize=14, pad=15)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-2, 2)
ax1.set_ylim(-0.05, 1.05)

# Right: Thermal derivative -df/dε (density of thermally excited quasiparticles)
ax2 = axes[1]
for T in temperatures[1:]:  # Skip T=0
    f = fermi_dirac(epsilon, E_F, T, k_B)
    df_deps = -np.gradient(f, epsilon)  # -df/dε
    label = f'T = {T:.2f} E_F/k_B'
    ax2.plot(epsilon, df_deps, linewidth=2, label=label)

ax2.axvline(x=E_F, color='k', linestyle='--', linewidth=1.5, alpha=0.5)
ax2.set_xlabel(r'Energy $(\epsilon - E_F) / E_F$', fontsize=13)
ax2.set_ylabel(r'$-\partial f/\partial \epsilon$', fontsize=13)
ax2.set_title('Thermal Excitation Density', fontsize=14, pad=15)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(-1, 1)

# Add annotation
ax2.annotate('Thermal width\n~ k_B T', xy=(0.2, 1.5), xytext=(0.7, 2.5),
            fontsize=11, arrowprops=dict(arrowstyle='->', color='red', lw=1.5),
            bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.5))

plt.tight_layout()
plt.show()

# Print thermal width information
print("Thermal Excitation Characteristics:")
print("-" * 60)
for T in [0.01, 0.1, 0.5]:
    thermal_width = 4.4 * k_B * T  # FWHM of -df/dε
    print(f"T = {T:.2f} E_F/k_B  →  Thermal width ≈ {thermal_width:.3f} E_F")

4.2 Cooper Instability: Bound States Above the Fermi Sea

The Cooper Problem Setup

In 1956, Leon Cooper posed a crucial question: Can two electrons with opposite momenta and opposite spins form a bound state when added above a filled Fermi sea, even with weak attractive interaction?

Key constraints:

• Two electrons with momenta $\mathbf{k}\uparrow$ and $-\mathbf{k}\downarrow$ (zero total momentum, spin singlet)
• Both electrons have energies above $E_F$ (Pauli principle blocks states below)
• Weak attractive interaction $V > 0$ (mediated by phonons)
• Interaction limited to energy shell $|\epsilon_{\mathbf{k}}| < \hbar\omega_D$ (Debye cutoff)

Two-Particle Schrödinger Equation

The pair wave function in momentum space $\phi(\mathbf{k})$ satisfies:

where $E_{\text{pair}}$ is the pair energy measured from $2E_F$ (binding energy if negative).

Contact Interaction Approximation

For simplicity, Cooper used a contact interaction:

This captures the essential physics: attractive pairing near the Fermi surface within a Debye energy window.

Cooper's Derivation

Summing over all momenta in the interaction shell:

Multiplying both sides by $(2\epsilon_{\mathbf{k}} - E_{\text{pair}})$ and summing over $\mathbf{k}$:

Defining $A = \sum_{\mathbf{k}} \phi(\mathbf{k})$ (non-zero for bound state):

Integral Form and Solution

Converting the sum to an integral with density of states $N(0)$:

where we assume $E_{\text{pair}} < 0$ (bound state). Evaluating:

For weak coupling ($VN(0) \ll 1$) and $|E_{\text{pair}}| \ll \hbar\omega_D$:

import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

# Cooper pair binding energy calculation
def cooper_equation(E_pair, N0_V, omega_D):
    """
    Cooper's equation: 1 = (V*N(0)/2) * ln[(2*omega_D - E_pair) / (-E_pair)]

    Parameters:
    -----------
    E_pair : float
        Pair binding energy (negative for bound state)
    N0_V : float
        Dimensionless coupling N(0)*V
    omega_D : float
        Debye energy cutoff
    """
    if E_pair >= 0:
        return 1e10  # No bound state

    lhs = 1.0
    rhs = (N0_V / 2.0) * np.log((2 * omega_D - E_pair) / (-E_pair))
    return lhs - rhs

def cooper_binding_analytical(N0_V, omega_D):
    """Analytical weak-coupling result."""
    return -2 * omega_D * np.exp(-2.0 / N0_V)

# Parameters
omega_D = 1.0  # Normalized Debye energy
N0_V_values = np.linspace(0.1, 1.0, 100)

# Calculate binding energies
E_pair_numerical = []
E_pair_analytical = []

for N0_V in N0_V_values:
    # Numerical solution
    E_guess = -0.01 * omega_D
    try:
        E_sol = fsolve(cooper_equation, E_guess, args=(N0_V, omega_D))[0]
        E_pair_numerical.append(abs(E_sol))
    except:
        E_pair_numerical.append(np.nan)

    # Analytical approximation
    E_anal = abs(cooper_binding_analytical(N0_V, omega_D))
    E_pair_analytical.append(E_anal)

# Plot
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Left: Log plot of binding energy vs coupling
ax1 = axes[0]
ax1.semilogy(N0_V_values, E_pair_analytical, 'b-', linewidth=2.5,
             label='Analytical: $2\hbar\omega_D e^{-2/N(0)V}$')
ax1.semilogy(N0_V_values, E_pair_numerical, 'ro', markersize=4, alpha=0.6,
             label='Numerical solution')
ax1.set_xlabel(r'Coupling strength $N(0)V$', fontsize=13)
ax1.set_ylabel(r'Binding energy $|E_{\rm pair}| / \hbar\omega_D$', fontsize=13)
ax1.set_title('Cooper Pair Binding Energy', fontsize=14, pad=15)
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3, which='both')
ax1.set_xlim(0.1, 1.0)

# Add annotation
ax1.text(0.5, 1e-2, r'Non-perturbative: $\sim e^{-2/N(0)V}$',
         fontsize=12, bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

# Right: Pair size (coherence length)
ax2 = axes[1]
# ξ_0 ~ ħv_F / |E_pair| (in units of ħv_F / ħω_D)
v_F = 1.0  # Normalized Fermi velocity
xi_0 = v_F / np.array(E_pair_analytical)

ax2.semilogy(N0_V_values, xi_0, 'g-', linewidth=2.5)
ax2.set_xlabel(r'Coupling strength $N(0)V$', fontsize=13)
ax2.set_ylabel(r'Coherence length $\xi_0 / (\hbar v_F / \hbar\omega_D)$', fontsize=13)
ax2.set_title('Cooper Pair Size', fontsize=14, pad=15)
ax2.grid(True, alpha=0.3, which='both')
ax2.set_xlim(0.1, 1.0)

# Add annotation
ax2.text(0.3, 1e2, 'Weak coupling:\nLarge pairs',
         fontsize=11, bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.5))

plt.tight_layout()
plt.show()

# Print example values
print("Cooper Pair Properties:")
print("-" * 70)
print(f"{'N(0)V':<10} {'|E_pair|/ℏω_D':<20} {'ξ_0 (units of ℏv_F/ℏω_D)':<30}")
print("-" * 70)
for N0_V in [0.2, 0.3, 0.4, 0.5]:
    E_b = abs(cooper_binding_analytical(N0_V, omega_D))
    xi = v_F / E_b
    print(f"{N0_V:<10.2f} {E_b:<20.6e} {xi:<30.1f}")

4.3 Cooper Pair Wave Function and Pairing Symmetry

Momentum Space Pairing

The Cooper pair wave function in momentum space couples states $(\mathbf{k}\uparrow, -\mathbf{k}\downarrow)$:

where $g(\mathbf{k})$ is the pair amplitude. Key property: zero total momentum ($\mathbf{k} + (-\mathbf{k}) = 0$).

Spin Singlet vs Triplet Pairing

The two-particle spin wave function can be either:

Singlet (Total spin $S=0$, antisymmetric):

This is the pairing in conventional BCS superconductors.

Triplet (Total spin $S=1$, symmetric):

Triplet pairing occurs in unconventional superconductors (e.g., Sr₂RuO₄, superfluid ³He).

Pairing Symmetry: s-wave, p-wave, d-wave

The spatial part of the pair wave function has angular momentum $L$. Combined with spin symmetry (Pauli principle requires overall antisymmetry):

Real-Space Cooper Pair Size

In real space, the pair wave function is:

For s-wave pairing with $g(\mathbf{k}) \approx \text{const}$ near $k_F$, the pair decays over the BCS coherence length:

This is much larger than the interatomic spacing ($\sim 0.3$ nm), meaning Cooper pairs heavily overlap.

4.4 BCS Ground State Wave Function

The BCS Variational Ansatz

Bardeen, Cooper, and Schrieffer (1957) generalized Cooper's result to a many-body ground state. They proposed a trial wave function that creates pairs coherently across all momentum states:

where:

• $|0\rangle$ is the vacuum state (no electrons)
• $c_{\mathbf{k}\sigma}^\dagger$ creates an electron with momentum $\mathbf{k}$ and spin $\sigma$
• $u_{\mathbf{k}}$ and $v_{\mathbf{k}}$ are complex amplitudes (coherence factors)
• Normalization: $|u_{\mathbf{k}}|^2 + |v_{\mathbf{k}}|^2 = 1$

Physical Interpretation

For each momentum pair $(\mathbf{k}\uparrow, -\mathbf{k}\downarrow)$:

• $|u_{\mathbf{k}}|^2$ = probability that the pair is empty
• $|v_{\mathbf{k}}|^2$ = probability that the pair is occupied
• The state is a coherent quantum superposition of different pair configurations

Unlike a classical system, $|\Psi_{\text{BCS}}\rangle$ does not have a definite particle number. It describes a macroscopic quantum state with definite phase.

Occupation Number

The expectation value of the occupation operator:

This smoothly interpolates from 1 (filled) below $E_F$ to 0 (empty) above $E_F$, with transition width $\sim \Delta$ around the Fermi surface.

Coherence Factors

The coherence factors are determined by minimizing the ground-state energy. The result (derived in next section) is:

where $\epsilon_{\mathbf{k}} = \hbar^2k^2/2m - E_F$ and $E_{\mathbf{k}} = \sqrt{\epsilon_{\mathbf{k}}^2 + \Delta^2}$ is the BCS quasiparticle energy.

Condensation Energy

The ground-state energy difference between superconducting and normal states (at $T=0$):

This positive condensation energy stabilizes the superconducting state.

import numpy as np
import matplotlib.pyplot as plt

# BCS coherence factors and occupation
def bcs_coherence_factors(epsilon, Delta):
    """
    Calculate BCS coherence factors u_k^2 and v_k^2.

    Parameters:
    -----------
    epsilon : array
        Single-particle energy measured from Fermi surface
    Delta : float
        BCS gap parameter
    """
    E_k = np.sqrt(epsilon**2 + Delta**2)
    u_k_sq = 0.5 * (1.0 + epsilon / E_k)
    v_k_sq = 0.5 * (1.0 - epsilon / E_k)
    u_k = np.sqrt(u_k_sq)
    v_k = np.sqrt(v_k_sq)
    return u_k_sq, v_k_sq, u_k, v_k, E_k

# Energy range
epsilon = np.linspace(-5, 5, 1000)  # In units of Delta
Delta = 1.0

u_sq, v_sq, u, v, E_k = bcs_coherence_factors(epsilon, Delta)

# Create figure
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Top-left: Coherence factors u_k^2 and v_k^2
ax1 = axes[0, 0]
ax1.plot(epsilon, u_sq, 'r-', linewidth=2.5, label=r'$u_k^2$ (empty pair probability)')
ax1.plot(epsilon, v_sq, 'b-', linewidth=2.5, label=r'$v_k^2$ (occupied pair probability)')
ax1.axvline(x=0, color='k', linestyle='--', linewidth=1.5, alpha=0.5)
ax1.axhline(y=0.5, color='gray', linestyle=':', linewidth=1, alpha=0.5)
ax1.fill_between(epsilon, 0, 1, where=(np.abs(epsilon) < Delta),
                  alpha=0.15, color='yellow', label='Pairing region (width ~ Δ)')
ax1.set_xlabel(r'$\epsilon_k / \Delta$', fontsize=12)
ax1.set_ylabel('Probability', fontsize=12)
ax1.set_title(r'BCS Coherence Factors $u_k^2$ and $v_k^2$', fontsize=13, pad=10)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-5, 5)
ax1.set_ylim(-0.05, 1.05)

# Top-right: Quasiparticle energy E_k
ax2 = axes[0, 1]
ax2.plot(epsilon, epsilon, 'k--', linewidth=1.5, alpha=0.5, label='Normal state')
ax2.plot(epsilon, E_k, 'b-', linewidth=2.5, label='BCS quasiparticle')
ax2.plot(epsilon, -E_k, 'b-', linewidth=2.5)
ax2.axhline(y=Delta, color='r', linestyle=':', linewidth=2, label=r'Gap $\pm\Delta$')
ax2.axhline(y=-Delta, color='r', linestyle=':', linewidth=2)
ax2.fill_between(epsilon, -Delta, Delta, alpha=0.2, color='yellow', label='Energy gap')
ax2.set_xlabel(r'$\epsilon_k / \Delta$', fontsize=12)
ax2.set_ylabel(r'$E_k / \Delta$', fontsize=12)
ax2.set_title('BCS Quasiparticle Dispersion', fontsize=13, pad=10)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(-5, 5)
ax2.set_ylim(-6, 6)

# Bottom-left: Product u_k * v_k (anomalous average)
ax3 = axes[1, 0]
uv_product = u * v
ax3.plot(epsilon, uv_product, 'g-', linewidth=2.5)
ax3.axvline(x=0, color='k', linestyle='--', linewidth=1.5, alpha=0.5)
ax3.axhline(y=0, color='gray', linestyle=':', linewidth=1, alpha=0.5)
max_uv = 0.5  # Maximum at epsilon=0
ax3.axhline(y=max_uv, color='r', linestyle=':', linewidth=1.5, alpha=0.7,
            label=f'Max = {max_uv:.2f}')
ax3.set_xlabel(r'$\epsilon_k / \Delta$', fontsize=12)
ax3.set_ylabel(r'$u_k v_k$', fontsize=12)
ax3.set_title(r'Anomalous Average $\langle c_{-k\downarrow} c_{k\uparrow}\rangle \propto u_k v_k$',
              fontsize=13, pad=10)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)
ax3.set_xlim(-5, 5)
ax3.set_ylim(-0.05, 0.55)

# Bottom-right: Comparison with Fermi-Dirac at T=0
ax4 = axes[1, 1]
f_fermi = (epsilon < 0).astype(float)  # Step function at T=0
ax4.plot(epsilon, f_fermi, 'k--', linewidth=2.5, alpha=0.7, label='Normal state (T=0)')
ax4.plot(epsilon, v_sq, 'b-', linewidth=2.5, label='BCS state')
ax4.fill_between(epsilon, f_fermi, v_sq, where=(f_fermi > v_sq),
                  alpha=0.3, color='blue', label='Hole excitations')
ax4.fill_between(epsilon, f_fermi, v_sq, where=(f_fermi < v_sq),
                  alpha=0.3, color='red', label='Particle excitations')
ax4.axvline(x=0, color='k', linestyle='--', linewidth=1.5, alpha=0.5)
ax4.set_xlabel(r'$\epsilon_k / \Delta$', fontsize=12)
ax4.set_ylabel(r'Occupation $\langle n_k \rangle$', fontsize=12)
ax4.set_title('BCS vs Normal State Occupation', fontsize=13, pad=10)
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)
ax4.set_xlim(-5, 5)
ax4.set_ylim(-0.05, 1.05)

plt.tight_layout()
plt.show()

# Print coherence factor values
print("BCS Coherence Factors at Key Energies:")
print("-" * 80)
print(f"{'ε_k/Δ':<10} {'u_k²':<12} {'v_k²':<12} {'u_k v_k':<12} {'E_k/Δ':<12}")
print("-" * 80)
for eps_val in [-3, -2, -1, -0.5, 0, 0.5, 1, 2, 3]:
    eps_arr = np.array([eps_val])
    u2, v2, u_val, v_val, E = bcs_coherence_factors(eps_arr * Delta, Delta)
    print(f"{eps_val:<10.1f} {u2[0]:<12.4f} {v2[0]:<12.4f} {u_val[0]*v_val[0]:<12.4f} {E[0]/Delta:<12.4f}")

4.5 BCS Gap Equation: Self-Consistent Determination of Δ

Mean-Field Hamiltonian

The BCS model Hamiltonian with attractive pairing interaction:

The mean-field approximation decouples the four-fermion term by defining the gap parameter:

Pair Amplitude in BCS State

In the BCS ground state, the anomalous average is:

where we used $u_{\mathbf{k}} v_{\mathbf{k}} = \frac{1}{2}\sqrt{1 - (\epsilon_{\mathbf{k}}/E_{\mathbf{k}})^2} = \frac{\Delta}{2E_{\mathbf{k}}}$.

Gap Equation at Zero Temperature

Substituting into the definition of $\Delta$:

Assuming constant density of states $N(0)$ and converting to integral:

Evaluating the integral:

For weak coupling ($\Delta_0 \ll \hbar\omega_D$), $\sinh^{-1}(x) \approx \ln(2x)$:

Finite Temperature Gap Equation

At temperature $T > 0$, thermal population of quasiparticles modifies the gap equation:

This must be solved numerically for $\Delta(T)$.

Critical Temperature

At $T = T_c$, the gap vanishes: $\Delta(T_c) = 0$. The gap equation becomes:

This yields:

Universal BCS Ratio

This universal ratio is a hallmark prediction of BCS theory, verified in many conventional superconductors (e.g., Al: 3.4, Sn: 3.5, Pb: 4.3).

import numpy as np
from scipy.integrate import quad
from scipy.optimize import brentq, fsolve
import matplotlib.pyplot as plt

# BCS gap equation solver
class BCSGapSolver:
    def __init__(self, N0_V, omega_D):
        """
        Initialize BCS gap solver.

        Parameters:
        -----------
        N0_V : float
            Dimensionless coupling constant N(0)*V
        omega_D : float
            Debye energy cutoff
        """
        self.N0_V = N0_V
        self.omega_D = omega_D
        self.k_B = 1.0  # Set Boltzmann constant to 1 (energy units)

        # Calculate T_c and Delta_0
        self.Delta_0 = 2.0 * omega_D * np.exp(-1.0 / N0_V)
        self.T_c = self.calculate_Tc()

    def calculate_Tc(self):
        """Calculate critical temperature from linearized gap equation."""
        # 1 = N(0)V * int_0^omega_D (dε/ε) * tanh(ε/2kT_c)
        def tc_equation(T):
            if T <= 0:
                return 1e10
            integrand = lambda eps: np.tanh(eps / (2 * self.k_B * T)) / eps
            integral, _ = quad(integrand, 1e-10, self.omega_D)
            return 1.0 - self.N0_V * integral

        # Find T_c
        T_guess = 0.5 * self.Delta_0 / self.k_B
        try:
            T_c = brentq(tc_equation, 0.01 * T_guess, 2.0 * T_guess)
            return T_c
        except:
            return 1.134 * self.omega_D * np.exp(-1.0 / self.N0_V)  # Analytical

    def gap_equation(self, Delta, T):
        """
        BCS gap equation at temperature T.
        Returns: 1 - RHS of gap equation (zero when satisfied)
        """
        if Delta < 0:
            return 1e10
        if Delta < 1e-12:
            return -1.0  # No solution

        def integrand(epsilon):
            E_k = np.sqrt(epsilon**2 + Delta**2)
            if T == 0:
                return 1.0 / E_k
            else:
                return np.tanh(E_k / (2 * self.k_B * T)) / E_k

        integral, _ = quad(integrand, 0, self.omega_D)
        return 1.0 - self.N0_V * integral

    def solve_gap(self, T):
        """Solve for gap Delta at temperature T."""
        if T >= self.T_c:
            return 0.0

        # Initial guess: weak-coupling approximation near T_c
        if T > 0.9 * self.T_c:
            Delta_guess = self.Delta_0 * np.sqrt(1 - (T / self.T_c)**2)
        else:
            Delta_guess = self.Delta_0 * (1 - 0.5 * (T / self.T_c)**2)

        try:
            Delta_sol = brentq(self.gap_equation, 0, 2*self.Delta_0, args=(T,))
            return Delta_sol
        except:
            try:
                result = fsolve(self.gap_equation, Delta_guess, args=(T,))
                if result[0] > 0:
                    return result[0]
                else:
                    return 0.0
            except:
                return 0.0

# Parameters for typical BCS superconductor
N0_V = 0.25
omega_D = 1.0  # Normalized

solver = BCSGapSolver(N0_V, omega_D)

print("BCS Theory Parameters:")
print("=" * 70)
print(f"Coupling constant N(0)V = {N0_V}")
print(f"Debye energy ℏω_D = {omega_D:.4f} (normalized)")
print(f"Zero-temperature gap Δ_0 = {solver.Delta_0:.6f}")
print(f"Critical temperature k_B T_c = {solver.T_c:.6f}")
print(f"Ratio 2Δ_0 / k_B T_c = {2*solver.Delta_0 / solver.T_c:.3f} (BCS: 3.53)")
print("=" * 70)

# Solve gap equation over temperature range
T_range = np.linspace(0, 1.2 * solver.T_c, 150)
Delta_T = np.array([solver.solve_gap(T) for T in T_range])

# BCS weak-coupling approximation near T_c
T_norm = T_range / solver.T_c
Delta_approx = solver.Delta_0 * np.sqrt(np.maximum(1 - T_norm**2, 0))

# Plot gap vs temperature
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Left: Absolute gap
ax1 = axes[0]
ax1.plot(T_norm, Delta_T / solver.Delta_0, 'b-', linewidth=2.5,
         label='Numerical BCS solution')
ax1.plot(T_norm, Delta_approx / solver.Delta_0, 'r--', linewidth=2,
         label=r'Approximation: $\Delta(T) \propto \sqrt{1-(T/T_c)^2}$')
ax1.axvline(x=1.0, color='gray', linestyle=':', linewidth=1.5, alpha=0.7,
            label=r'$T_c$')
ax1.fill_between(T_norm, 0, Delta_T / solver.Delta_0, alpha=0.15, color='blue')
ax1.set_xlabel(r'Temperature $T / T_c$', fontsize=13)
ax1.set_ylabel(r'Gap $\Delta(T) / \Delta_0$', fontsize=13)
ax1.set_title('Temperature Dependence of BCS Gap', fontsize=14, pad=15)
ax1.legend(fontsize=11, loc='upper right')
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 1.2)
ax1.set_ylim(0, 1.1)

# Right: Low-temperature exponential behavior
ax2 = axes[1]
T_low = T_range[T_range < 0.5 * solver.T_c]
Delta_low = Delta_T[T_range < 0.5 * solver.T_c]
ax2.plot(T_low / solver.T_c, Delta_low / solver.Delta_0, 'b-', linewidth=2.5,
         label='BCS gap')
ax2.axhline(y=1.0, color='k', linestyle='--', linewidth=1.5, alpha=0.5,
            label=r'$\Delta_0$ (T=0)')
ax2.set_xlabel(r'Temperature $T / T_c$', fontsize=13)
ax2.set_ylabel(r'Gap $\Delta(T) / \Delta_0$', fontsize=13)
ax2.set_title('Low-Temperature Gap (Nearly Constant)', fontsize=14, pad=15)
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, 0.5)
ax2.set_ylim(0.95, 1.01)

plt.tight_layout()
plt.show()

# Print gap values at selected temperatures
print("\nGap Values at Key Temperatures:")
print("-" * 70)
print(f"{'T/T_c':<10} {'Δ(T)/Δ_0':<15} {'Δ(T) (abs)':<15}")
print("-" * 70)
for T_ratio in [0.0, 0.2, 0.4, 0.6, 0.8, 0.9, 0.95, 1.0]:
    T = T_ratio * solver.T_c
    Delta = solver.solve_gap(T)
    print(f"{T_ratio:<10.2f} {Delta/solver.Delta_0:<15.4f} {Delta:<15.6f}")

4.6 Temperature Dependence of the Gap

Behavior Near Critical Temperature

Close to $T_c$, the gap vanishes as:

This square-root singularity is characteristic of a second-order phase transition.

Low Temperature Behavior

For $T \ll T_c$, the gap remains nearly constant:

The exponential factor reflects thermally excited quasiparticles above the gap.

Full Numerical Solution

The full temperature dependence requires numerical solution of:

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

# Full numerical solution for Delta(T)
def solve_delta_T_full(N0_V, omega_D, num_points=200):
    """
    Solve BCS gap equation for full temperature range.

    Returns:
    --------
    T_array : array
        Temperature values
    Delta_array : array
        Gap values at each temperature
    T_c : float
        Critical temperature
    Delta_0 : float
        Zero-temperature gap
    """
    k_B = 1.0
    Delta_0 = 2.0 * omega_D * np.exp(-1.0 / N0_V)

    # Find T_c
    def tc_equation(T):
        integrand = lambda eps: np.tanh(eps / (2 * k_B * T)) / eps
        integral, _ = quad(integrand, 1e-10, omega_D)
        return 1.0 - N0_V * integral

    T_c = brentq(tc_equation, 0.01 * Delta_0, 2.0 * Delta_0)

    # Solve for Delta(T)
    T_array = np.linspace(0, T_c * 0.999, num_points)
    Delta_array = np.zeros_like(T_array)

    for i, T in enumerate(T_array):
        if T == 0:
            Delta_array[i] = Delta_0
        else:
            def gap_eq(Delta):
                if Delta < 1e-12:
                    return -1.0
                integrand = lambda eps: np.tanh(np.sqrt(eps**2 + Delta**2) / (2*k_B*T)) / np.sqrt(eps**2 + Delta**2)
                integral, _ = quad(integrand, 0, omega_D)
                return 1.0 - N0_V * integral

            try:
                Delta_guess = Delta_0 * np.sqrt(1 - (T/T_c)**2)
                Delta_sol = brentq(gap_eq, 0, 2*Delta_0)
                Delta_array[i] = Delta_sol
            except:
                Delta_array[i] = 0.0

    return T_array, Delta_array, T_c, Delta_0

# Solve for different coupling strengths
N0_V_values = [0.2, 0.25, 0.3, 0.4]
omega_D = 1.0

plt.figure(figsize=(12, 7))

for N0_V in N0_V_values:
    T_arr, Delta_arr, Tc, D0 = solve_delta_T_full(N0_V, omega_D, num_points=150)

    # Plot normalized gap
    plt.plot(T_arr / Tc, Delta_arr / D0, linewidth=2.5,
             label=f'N(0)V = {N0_V:.2f}, 2Δ₀/k_BT_c = {2*D0/Tc:.2f}')

# Add universal BCS curve
T_universal = np.linspace(0, 1, 100)
Delta_universal = np.sqrt(np.maximum(1 - T_universal**2, 0))
plt.plot(T_universal, Delta_universal, 'k--', linewidth=2, alpha=0.5,
         label='Universal BCS (approx)')

plt.axvline(x=1.0, color='gray', linestyle=':', linewidth=1.5, alpha=0.7)
plt.axhline(y=1.0, color='gray', linestyle=':', linewidth=1, alpha=0.5)

plt.xlabel(r'Reduced Temperature $T / T_c$', fontsize=14)
plt.ylabel(r'Reduced Gap $\Delta(T) / \Delta_0$', fontsize=14)
plt.title('Universal BCS Gap Function', fontsize=15, pad=15)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.xlim(0, 1.05)
plt.ylim(0, 1.1)

# Add annotation for critical behavior
plt.annotate(r'$\Delta(T) \sim \sqrt{1-T/T_c}$ near $T_c$',
            xy=(0.95, 0.25), xytext=(0.6, 0.4),
            fontsize=12, arrowprops=dict(arrowstyle='->', color='red', lw=1.5),
            bbox=dict(boxstyle='round', facecolor='pink', alpha=0.3))

plt.tight_layout()
plt.show()

# Compare analytical approximations
print("Temperature Dependence Approximations:")
print("=" * 70)
print("Near T_c: Δ(T) ≈ Δ_0 √(1 - T/T_c)")
print("Low T:   Δ(T) ≈ Δ_0 [1 - √(2πk_BT/Δ_0) exp(-Δ_0/k_BT)]")
print("=" * 70)

4.7 Density of States and Coherence Peaks

Normal State Density of States

In the normal metallic state, the DOS near the Fermi surface is approximately constant:

Superconducting Density of States

In the BCS state, the quasiparticle DOS is:

This can be derived from the dispersion $E_{\mathbf{k}} = \sqrt{\epsilon_{\mathbf{k}}^2 + \Delta^2}$ using:

Key Features

• Energy gap: $N_{\text{SC}}(E) = 0$ for $|E| < \Delta$ (no single-particle excitations below gap)
• Coherence peaks: $N_{\text{SC}}(E) \to \infty$ as $E \to \Delta^+$ (square-root singularity)
• High-energy limit: $N_{\text{SC}}(E) \to N(0)$ for $E \gg \Delta$ (recovers normal state)

Density of States Formula

For $|E| > \Delta$:

Near the gap edge ($E = \Delta + \delta$, $\delta \ll \Delta$):

This square-root divergence is a hallmark of BCS superconductivity.

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter1d

# Superconducting density of states
def dos_bcs(E, Delta, N0=1.0):
    """
    BCS density of states.

    Parameters:
    -----------
    E : array
        Energy values
    Delta : float
        Gap parameter
    N0 : float
        Normal state DOS at Fermi level
    """
    dos = np.zeros_like(E)
    mask = np.abs(E) >= Delta
    dos[mask] = N0 * np.abs(E[mask]) / np.sqrt(E[mask]**2 - Delta**2)
    return dos

# Energy range
E = np.linspace(-4, 4, 2000)  # In units of Delta
Delta = 1.0
N0 = 1.0

# Calculate DOS
dos_sc = dos_bcs(E, Delta, N0)
dos_normal = N0 * np.ones_like(E)

# Create figure
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Top-left: BCS DOS with coherence peaks
ax1 = axes[0, 0]
ax1.plot(E, dos_normal, 'k--', linewidth=2, alpha=0.6, label='Normal state')
ax1.plot(E, dos_sc, 'b-', linewidth=2.5, label='Superconducting state')
ax1.axvline(x=Delta, color='r', linestyle=':', linewidth=2, alpha=0.7)
ax1.axvline(x=-Delta, color='r', linestyle=':', linewidth=2, alpha=0.7)
ax1.fill_between(E, 0, dos_sc.max() * 1.2, where=(np.abs(E) < Delta),
                  alpha=0.2, color='yellow', label='Energy gap')
ax1.set_xlabel(r'Energy $E / \Delta$', fontsize=12)
ax1.set_ylabel(r'DOS $N(E) / N(0)$', fontsize=12)
ax1.set_title('BCS Density of States', fontsize=13, pad=10)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-4, 4)
ax1.set_ylim(0, 5)

# Add annotation
ax1.annotate('Coherence\npeak', xy=(Delta, 4.0), xytext=(2, 3.5),
            fontsize=11, arrowprops=dict(arrowstyle='->', color='red', lw=1.5),
            bbox=dict(boxstyle='round', facecolor='pink', alpha=0.3))

# Top-right: Zoom near gap edge
ax2 = axes[0, 1]
E_zoom = E[(E > 0.9*Delta) & (E < 2*Delta)]
dos_zoom = dos_bcs(E_zoom, Delta, N0)
ax2.plot(E_zoom, dos_zoom, 'b-', linewidth=2.5)
ax2.axvline(x=Delta, color='r', linestyle='--', linewidth=2, label=r'Gap edge $\Delta$')
ax2.set_xlabel(r'Energy $E / \Delta$', fontsize=12)
ax2.set_ylabel(r'DOS $N(E) / N(0)$', fontsize=12)
ax2.set_title('Coherence Peak Detail', fontsize=13, pad=10)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0.9*Delta, 2*Delta)

# Bottom-left: Temperature broadening effect
ax3 = axes[1, 0]
temperatures = [0.01, 0.05, 0.1, 0.2]
for kT in temperatures:
    # Thermal broadening via convolution with derivative of Fermi function
    sigma_pts = int(kT / (E[1] - E[0]) * 2)  # Convert to array points
    if sigma_pts > 0:
        dos_broadened = gaussian_filter1d(dos_sc, sigma=sigma_pts)
    else:
        dos_broadened = dos_sc

    label = f'k_B T = {kT:.2f} Δ' if kT > 0 else 'T = 0'
    ax3.plot(E, dos_broadened, linewidth=2, label=label)

ax3.axvline(x=Delta, color='k', linestyle=':', linewidth=1.5, alpha=0.5)
ax3.axvline(x=-Delta, color='k', linestyle=':', linewidth=1.5, alpha=0.5)
ax3.set_xlabel(r'Energy $E / \Delta$', fontsize=12)
ax3.set_ylabel(r'DOS $N(E) / N(0)$', fontsize=12)
ax3.set_title('Temperature Broadening Effect', fontsize=13, pad=10)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)
ax3.set_xlim(-3, 3)
ax3.set_ylim(0, 4)

# Bottom-right: Integrated DOS (check particle conservation)
ax4 = axes[1, 1]
# Integrate DOS from -infinity to E
E_int = E[E >= 0]
integrated_dos = np.zeros_like(E_int)
for i, e in enumerate(E_int):
    mask = (E >= -e) & (E <= e)
    integrated_dos[i] = np.trapz(dos_sc[mask], E[mask])

# Normalization: should equal 2*N(0)*E for normal state
integrated_normal = 2 * N0 * E_int

ax4.plot(E_int, integrated_normal, 'k--', linewidth=2, alpha=0.6, label='Normal state')
ax4.plot(E_int, integrated_dos, 'b-', linewidth=2.5, label='Superconducting')
ax4.set_xlabel(r'Integration limit $E / \Delta$', fontsize=12)
ax4.set_ylabel(r'Integrated DOS', fontsize=12)
ax4.set_title('Particle Number Conservation', fontsize=13, pad=10)
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)
ax4.set_xlim(0, 4)

plt.tight_layout()
plt.show()

# Print DOS values near gap
print("BCS Density of States Near Gap Edge:")
print("-" * 60)
print(f"{'E/Δ':<10} {'N_SC(E)/N(0)':<20}")
print("-" * 60)
for E_val in [0.5, 0.9, 0.99, 1.0, 1.01, 1.1, 1.5, 2.0, 3.0]:
    dos_val = dos_bcs(np.array([E_val]), Delta, N0)[0]
    if dos_val > 0:
        print(f"{E_val:<10.2f} {dos_val:<20.4f}")
    else:
        print(f"{E_val:<10.2f} {'0 (in gap)':<20}")

4.8 Tunneling Spectroscopy: Experimental Verification

Superconductor-Insulator-Normal Metal (SIN) Junction

In a tunnel junction between a superconductor (S) and a normal metal (N) separated by a thin insulator (I), the differential conductance is proportional to the superconductor's DOS:

where $V$ is the applied bias voltage.

SIN Tunneling Characteristics

• At $|eV| < \Delta$: No current (gap region)
• At $|eV| = \Delta$: Sharp conductance peak (coherence peak)
• At $|eV| > \Delta$: Current flows, conductance decreases with energy

Superconductor-Insulator-Superconductor (SIS) Junction

For two identical superconductors with gap $\Delta$:

At $T = 0$:

• No current for $|eV| < 2\Delta$ (sum of two gaps)
• Sharp onset at $|eV| = 2\Delta$
• Strong peaks at $|eV| = 2\Delta$ (convolution of coherence peaks)

Experimental Significance

Tunneling spectroscopy provides direct measurement of:

• Gap magnitude $\Delta(T)$ as function of temperature
• Coherence peak structure (BCS vs non-BCS behavior)
• Gap anisotropy (for unconventional superconductors)
• Gap symmetry (s-wave, d-wave, etc.)

import numpy as np
import matplotlib.pyplot as plt

# SIN and SIS tunneling conductance
def sin_conductance(V, Delta, T, N0=1.0):
    """
    SIN junction differential conductance.
    Proportional to superconductor DOS at energy eV.

    Parameters:
    -----------
    V : array
        Bias voltage (in units of Delta/e)
    Delta : float
        Gap parameter
    T : float
        Temperature
    N0 : float
        Normal state DOS
    """
    E = V  # Energy = eV (in units where e=1)
    dos = dos_bcs(E, Delta, N0)

    # Add thermal broadening
    if T > 0:
        from scipy.ndimage import gaussian_filter1d
        sigma_pts = int(T / (V[1] - V[0]) * 2)
        if sigma_pts > 0:
            dos = gaussian_filter1d(dos, sigma=sigma_pts)

    return dos

def sis_conductance(V, Delta, T):
    """
    SIS junction differential conductance at T=0.
    Convolution of two superconductor DOS functions.
    """
    if T > 0:
        # Full thermal calculation (complex, simplified here)
        return sin_conductance(V, Delta, T)  # Approximation

    # At T=0: step at |V| = 2*Delta/e
    conductance = np.zeros_like(V)
    mask = np.abs(V) >= 2*Delta

    # Approximate convolution peak structure
    for i, v in enumerate(V):
        if abs(v) >= 2*Delta:
            # Simplified: peak at threshold
            E1_range = np.linspace(-10*Delta, 10*Delta, 500)
            E2 = E1_range - v
            dos1 = dos_bcs(E1_range, Delta, 1.0)
            dos2 = dos_bcs(E2, Delta, 1.0)
            conductance[i] = np.trapz(dos1 * dos2, E1_range) / (20*Delta)

    return conductance

# Voltage range
V = np.linspace(-3, 3, 1000)  # In units of Delta/e
Delta = 1.0
temperatures = [0.0, 0.05, 0.1, 0.2]

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Left: SIN tunneling
ax1 = axes[0]
for T in temperatures:
    dIdV = sin_conductance(V, Delta, T)
    label = f'T = {T:.2f} Δ/k_B' if T > 0 else 'T = 0'
    ax1.plot(V, dIdV, linewidth=2, label=label)

ax1.axvline(x=Delta, color='r', linestyle=':', linewidth=1.5, alpha=0.7, label='Gap ±Δ')
ax1.axvline(x=-Delta, color='r', linestyle=':', linewidth=1.5, alpha=0.7)
ax1.fill_between(V, 0, ax1.get_ylim()[1], where=(np.abs(V) < Delta),
                  alpha=0.15, color='yellow')
ax1.set_xlabel(r'Bias Voltage $eV / \Delta$', fontsize=13)
ax1.set_ylabel(r'$dI/dV$ (arb. units)', fontsize=13)
ax1.set_title('SIN Junction Tunneling Spectroscopy', fontsize=14, pad=15)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-3, 3)
ax1.set_ylim(0, 5)

# Right: SIS tunneling
ax2 = axes[1]
dIdV_sis = sis_conductance(V, Delta, T=0)
ax2.plot(V, dIdV_sis, 'b-', linewidth=2.5, label='T = 0')

# Add thermal broadened version
dIdV_sis_T = sis_conductance(V, Delta, T=0.1)
ax2.plot(V, dIdV_sis_T, 'r-', linewidth=2, label='T = 0.1 Δ/k_B')

ax2.axvline(x=2*Delta, color='g', linestyle='--', linewidth=2, alpha=0.7,
            label='Threshold ±2Δ')
ax2.axvline(x=-2*Delta, color='g', linestyle='--', linewidth=2, alpha=0.7)
ax2.fill_between(V, 0, ax2.get_ylim()[1], where=(np.abs(V) < 2*Delta),
                  alpha=0.15, color='lightblue')
ax2.set_xlabel(r'Bias Voltage $eV / \Delta$', fontsize=13)
ax2.set_ylabel(r'$dI/dV$ (arb. units)', fontsize=13)
ax2.set_title('SIS Junction Tunneling Spectroscopy', fontsize=14, pad=15)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(-3, 3)

plt.tight_layout()
plt.show()

print("Tunneling Spectroscopy Summary:")
print("=" * 70)
print("SIN Junction:")
print("  - Gap appears at |eV| = Δ")
print("  - Coherence peaks at gap edges")
print("  - Direct measurement of N_SC(E)")
print("\nSIS Junction:")
print("  - Current threshold at |eV| = 2Δ (sum of two gaps)")
print("  - Strong peaks at threshold from coherence peak convolution")
print("  - More pronounced structure than SIN")
print("=" * 70)

4.9 Specific Heat Jump at Critical Temperature

Thermodynamic Signature of Phase Transition

The superconducting transition at $T_c$ is a second-order phase transition, characterized by a discontinuous jump in the specific heat.

Normal State Electronic Specific Heat

In the normal metallic state (Fermi liquid):

where the Sommerfeld coefficient is:

Superconducting State Specific Heat

Below $T_c$, the electronic specific heat has two regimes:

Low temperature ($T \ll T_c$):

This exponential suppression reflects the energy cost of thermally exciting quasiparticles across the gap.

Near $T_c$:

Universal Jump Ratio

The BCS prediction for the specific heat discontinuity at $T_c$:

This universal ratio is confirmed in many conventional superconductors (Al: 1.45, Sn: 1.60, Nb: 1.87).

import numpy as np
import matplotlib.pyplot as plt

# Specific heat calculations
def specific_heat_normal(T, gamma):
    """Normal state: C_n = gamma * T"""
    return gamma * T

def specific_heat_bcs_full(T, T_c, Delta_0, gamma):
    """
    BCS specific heat (approximate form).

    For T << T_c: exponential suppression
    Near T_c: discontinuous jump
    """
    if T >= T_c:
        return gamma * T

    if T == 0:
        return 0.0

    # Temperature-dependent gap (simplified)
    Delta_T = Delta_0 * np.sqrt(np.maximum(1 - (T / T_c)**2, 0))

    # BCS specific heat (approximation)
    # Full formula involves derivatives of gap equation
    # Here we use phenomenological form
    if T < 0.5 * T_c:
        # Exponential regime
        k_B = 1.0
        x = Delta_T / (k_B * T)
        C_exp = gamma * T_c * 8.5 * x**2 * np.exp(-x)
        return C_exp
    else:
        # Near T_c: use jump formula
        # Linear interpolation to jump value 1.43 * gamma * T_c
        t_reduced = T / T_c
        C_jump = gamma * T * (1.0 + 0.43 * (1 - t_reduced)**0.5)
        return C_jump

# Parameters
k_B = 1.0
T_c = 1.0  # Normalized
Delta_0 = 1.76 * k_B * T_c  # BCS ratio
gamma = 1.0  # Normalized

# Temperature range
T = np.linspace(0.01, 2*T_c, 500)

# Calculate specific heats
C_normal = np.array([specific_heat_normal(t, gamma) for t in T])
C_sc = np.array([specific_heat_bcs_full(t, T_c, Delta_0, gamma) for t in T])

# Normalized by gamma*T_c
C_n_norm = C_normal / (gamma * T_c)
C_sc_norm = C_sc / (gamma * T_c)

# Create plots
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Left: Full temperature range
ax1 = axes[0]
ax1.plot(T/T_c, C_n_norm, 'k--', linewidth=2, alpha=0.7, label='Normal state')
ax1.plot(T/T_c, C_sc_norm, 'b-', linewidth=2.5, label='Superconducting state')

# Mark the jump at T_c
jump_height = 1.43
ax1.plot([1.0, 1.0], [1.0, jump_height], 'ro-', linewidth=3, markersize=8,
         label=f'Jump: ΔC/C_n = {jump_height - 1.0:.2f}')

ax1.axvline(x=1.0, color='gray', linestyle=':', linewidth=1.5, alpha=0.7)
ax1.fill_between(T/T_c, 0, C_sc_norm, where=(T < T_c), alpha=0.15, color='blue')

ax1.set_xlabel(r'Temperature $T / T_c$', fontsize=13)
ax1.set_ylabel(r'Specific Heat $C / \gamma T_c$', fontsize=13)
ax1.set_title('BCS Electronic Specific Heat', fontsize=14, pad=15)
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 2)
ax1.set_ylim(0, 2.5)

# Add annotations
ax1.text(0.3, 0.2, r'$C_{SC} \sim e^{-\Delta/k_BT}$', fontsize=12, color='blue')
ax1.text(1.5, 1.6, r'$C_n = \gamma T$', fontsize=12, color='black')

# Right: Low-temperature exponential behavior (log scale)
ax2 = axes[1]
T_low = T[T < 0.6 * T_c]
C_sc_low = np.array([specific_heat_bcs_full(t, T_c, Delta_0, gamma) for t in T_low])
C_sc_low_norm = C_sc_low / (gamma * T_c)

# Avoid log(0)
C_sc_low_norm[C_sc_low_norm == 0] = 1e-10

ax2.semilogy(T_low/T_c, C_sc_low_norm, 'b-', linewidth=2.5, label='BCS (numerical)')

# Analytical exponential fit
Delta_eff = Delta_0
analytical_exp = 8.5 * (Delta_eff / (k_B * T_low))**2 * np.exp(-Delta_eff / (k_B * T_low))
ax2.semilogy(T_low/T_c, analytical_exp, 'r--', linewidth=2,
             label=r'$\propto (\Delta/k_BT)^2 e^{-\Delta/k_BT}$')

ax2.set_xlabel(r'Temperature $T / T_c$', fontsize=13)
ax2.set_ylabel(r'Specific Heat $C / \gamma T_c$ (log)', fontsize=13)
ax2.set_title('Low-Temperature Exponential Suppression', fontsize=14, pad=15)
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3, which='both')
ax2.set_xlim(0, 0.6)
ax2.set_ylim(1e-4, 1e1)

plt.tight_layout()
plt.show()

# Print specific heat values
print("Specific Heat at Key Temperatures:")
print("=" * 70)
print(f"{'T/T_c':<10} {'C_SC/γT_c':<15} {'C_n/γT_c':<15} {'Ratio C_SC/C_n':<15}")
print("=" * 70)
for T_ratio in [0.1, 0.3, 0.5, 0.7, 0.9, 0.95, 1.0]:
    t_val = T_ratio * T_c
    C_sc_val = specific_heat_bcs_full(t_val, T_c, Delta_0, gamma) / (gamma * T_c)
    C_n_val = specific_heat_normal(t_val, gamma) / (gamma * T_c)
    if C_n_val > 0:
        ratio = C_sc_val / C_n_val
    else:
        ratio = 0
    print(f"{T_ratio:<10.2f} {C_sc_val:<15.4f} {C_n_val:<15.4f} {ratio:<15.4f}")

print("\nBCS Prediction: ΔC/C_n(T_c) = 1.43")
print("Experimental values: Al (1.45), Sn (1.60), Nb (1.87)")

Summary

Practice Problems

Further Reading

• J. Bardeen, L. N. Cooper, and J. R. Schrieffer, "Theory of Superconductivity," Phys. Rev. 108, 1175 (1957) — The original BCS paper (Nobel Prize 1972)
• L. N. Cooper, "Bound Electron Pairs in a Degenerate Fermi Gas," Phys. Rev. 104, 1189 (1956) — Cooper problem
• M. Tinkham, "Introduction to Superconductivity," 2nd ed. (Dover, 2004) — Comprehensive graduate textbook
• J. R. Schrieffer, "Theory of Superconductivity," (Perseus Books, 1999) — Detailed theoretical treatment by BCS co-author
• A. A. Abrikosov, "Fundamentals of the Theory of Metals," (North-Holland, 1988) — Advanced many-body theory
• P. G. de Gennes, "Superconductivity of Metals and Alloys," (Westview Press, 1999) — Classic reference with BdG formalism
• I. Giaever, "Energy Gap in Superconductors Measured by Electron Tunneling," Phys. Rev. Lett. 5, 147 (1960) — Experimental tunneling verification