Chapter 1: Strong Coupling Theory and Eliashberg Formalism

Learning Objectives

By the end of this chapter, you will be able to:

• Identify when BCS weak-coupling theory fails and strong-coupling effects dominate
• Understand the Migdal theorem and why vertex corrections can be neglected
• Derive and interpret the Eliashberg equations for gap and renormalization functions
• Extract the electron-phonon spectral function $\alpha^2F(\omega)$ from experimental data
• Apply McMillan and Allen-Dynes formulas to predict $T_c$ from material properties
• Implement numerical solvers for Eliashberg equations in Python
• Understand how strong coupling modifies isotope effect and thermodynamic properties
• Apply these methods to real superconductors (Pb, Nb, Nb$_3$Sn, hydrides)

1. Introduction: Why BCS Is Not Enough

1.1 The Success and Limitations of BCS Theory

BCS theory (1957) provides a remarkably accurate description of conventional superconductors within its regime of validity. Recall that BCS predicts:

• Gap ratio: $2\Delta(0)/k_BT_c = 3.52$ (universal)
• Specific heat jump: $\Delta C/C_n = 1.43$ at $T_c$
• Isotope effect: $T_c \propto M^{-\alpha}$ with $\alpha = 0.5$
• Energy gap that opens smoothly below $T_c$

However, BCS theory rests on several key approximations:

1.2 Strong Coupling Superconductors

For several important materials, these approximations fail spectacularly. Consider the following experimental evidence:

The bold entries show clear deviations from BCS predictions:

• Enhanced gap ratios: $2\Delta(0)/k_BT_c > 3.52$ for Pb and Hg
• Reduced isotope effect: $\alpha \approx 0.1$ for Nb$_3$Sn instead of 0.5
• Large coupling: $\lambda > 1$ indicates strong coupling regime

These discrepancies motivated Eliashberg to develop a more complete theory in 1960 that properly accounts for:

• Retarded (frequency-dependent) electron-phonon interaction
• Energy-dependent gap function $\Delta(\omega)$
• Quasiparticle renormalization $Z(\omega)$
• Coulomb pseudopotential $\mu^*$ screening

2. Migdal-Eliashberg Theory Foundations

2.1 The Migdal Theorem

A natural question arises: if we include retardation effects, shouldn't we also include vertex corrections (higher-order diagrams where electron-phonon vertices are dressed by other phonons)? Migdal (1958) showed these corrections are suppressed by a small parameter:

Physical Interpretation: The electron-phonon interaction involves momentum transfer of order $q \sim k_F$ (the Fermi wavevector). The phonon frequency is $\omega_{\text{ph}} \sim v_s q$ where $v_s$ is the sound velocity. The electron Fermi velocity is $v_F \sim \hbar k_F/m$. The ratio of these velocities gives:

$$\frac{v_s}{v_F} \sim \frac{\omega_{\text{ph}}}{E_F} \sim \left(\frac{m}{M}\right)^{1/2}$$

Vertex corrections involve additional phonon propagators at intermediate energies $\sim \omega_{\text{ph}}$, which are small compared to the electron energy scale $E_F$. Thus, Migdal's theorem allows us to work with renormalized electron propagators without worrying about vertex renormalization.

2.2 Electron Self-Energy from Phonons

In many-body theory, the electron Green's function $G(\mathbf{k}, i\omega_n)$ satisfies the Dyson equation:

$$G(\mathbf{k}, i\omega_n) = \frac{1}{i\omega_n - \xi_{\mathbf{k}} - \Sigma(\mathbf{k}, i\omega_n)}$$

where $\xi_{\mathbf{k}} = \epsilon_{\mathbf{k}} - \mu$ is the energy relative to the chemical potential, and $\Sigma(\mathbf{k}, i\omega_n)$ is the electron self-energy. In the superconducting state, we need to work with the Nambu-Gorkov formalism, where the Green's function becomes a 2×2 matrix:

$$\hat{G}(\mathbf{k}, i\omega_n) = \begin{pmatrix} G(\mathbf{k}, i\omega_n) & F(\mathbf{k}, i\omega_n) \\ F^\dagger(\mathbf{k}, i\omega_n) & -G(-\mathbf{k}, -i\omega_n) \end{pmatrix}$$

The anomalous Green's function $F(\mathbf{k}, i\omega_n)$ describes Cooper pairing. Similarly, the self-energy becomes a matrix:

$$\hat{\Sigma}(\mathbf{k}, i\omega_n) = \begin{pmatrix} \Sigma_N(\mathbf{k}, i\omega_n) & \Sigma_{\text{anom}}(\mathbf{k}, i\omega_n) \\ \Sigma_{\text{anom}}^*(\mathbf{k}, i\omega_n) & -\Sigma_N(-\mathbf{k}, -i\omega_n) \end{pmatrix}$$

where $\Sigma_N$ is the normal self-energy and $\Sigma_{\text{anom}}$ gives the pairing. The full Dyson equation is:

$$\hat{G}^{-1}(\mathbf{k}, i\omega_n) = \hat{G}_0^{-1}(\mathbf{k}, i\omega_n) - \hat{\Sigma}(\mathbf{k}, i\omega_n)$$

2.3 Eliashberg Equations: Derivation

The self-energy from electron-phonon interaction in Migdal approximation is:

$$\Sigma(\mathbf{k}, i\omega_n) = T\sum_{m} \int \frac{d^3k'}{(2\pi)^3} |g_{\mathbf{k}-\mathbf{k}'}|^2 D(\mathbf{k}-\mathbf{k}', i\omega_n - i\omega_m) G(\mathbf{k}', i\omega_m)$$

where $g_{\mathbf{q}}$ is the electron-phonon matrix element and $D(\mathbf{q}, i\nu)$ is the phonon propagator. We now make several simplifications standard in Eliashberg theory:

1. Average over Fermi surface: Replace $\mathbf{k}$-dependent quantities by Fermi surface averages
2. Isotropic approximation: Assume gap depends only on frequency, not momentum direction
3. Define renormalization: Write $\Sigma_N(i\omega_n) = (1 - Z(i\omega_n))i\omega_n + i\omega_n\chi(i\omega_n)$

The key quantities are:

• $Z(i\omega_n)$: Quasiparticle renormalization function (mass renormalization)
• $\phi(i\omega_n) = Z(i\omega_n)\Delta(i\omega_n)$: Renormalized gap function
• $\chi(i\omega_n)$: Band shift (usually neglected)

After considerable algebra, one arrives at the Eliashberg equations:

Here:

• $\omega_n = \pi T(2n+1)$ are fermionic Matsubara frequencies
• $\lambda(i\nu)$ is the electron-phonon interaction kernel (to be defined)
• $\mu^*$ is the Coulomb pseudopotential (repulsive interaction)
• $\omega_c$ is the phonon cutoff energy
• The sums run over $m = 0, \pm 1, \pm 2, \ldots$

The second equation includes Coulomb repulsion with a cutoff: $\mu^*(|\omega_m| < \omega_c)$ is 1 if $|\omega_m| < \omega_c$ and 0 otherwise. This represents the fact that Coulomb repulsion is instantaneous and acts at all energies, but is strongly screened in metals. The effective parameter $\mu^*$ incorporates this screening.

2.4 The Electron-Phonon Kernel

The interaction kernel $\lambda(i\nu)$ is related to the electron-phonon spectral function $\alpha^2F(\omega)$ by:

$$\lambda(i\nu) = 2\int_0^{\omega_{\text{max}}} d\omega \frac{\omega \alpha^2F(\omega)}{\nu^2 + \omega^2}$$

The spectral function $\alpha^2F(\omega)$ contains all the information about electron-phonon coupling:

The total electron-phonon coupling constant is:

$$\lambda = 2\int_0^{\omega_{\text{max}}} \frac{\alpha^2F(\omega)}{\omega} d\omega$$

This dimensionless parameter $\lambda$ quantifies the overall strength of electron-phonon interaction. In BCS theory (weak coupling), $\lambda \ll 1$. In strong-coupling materials like Pb, $\lambda \approx 1.5$.

3. Solving the Eliashberg Equations

3.1 Iterative Solution Method

The Eliashberg equations are coupled nonlinear integral equations. They cannot be solved analytically except in the weak-coupling limit (where they reduce to BCS). The standard approach is iterative:

1. Initialize: Start with BCS-like guess $\phi^{(0)}(i\omega_n) = \Delta_0$, $Z^{(0)}(i\omega_n) = 1$
2. Iterate: Given $\phi^{(k)}$ and $Z^{(k)}$, compute new values from Eliashberg equations: $$Z^{(k+1)}(i\omega_n) = 1 + \frac{\pi T}{\omega_n}\sum_{m} \frac{\text{sign}(\omega_m)}{\sqrt{\omega_m^2 + (\phi^{(k)}/Z^{(k)})^2}}\lambda(i\omega_n - i\omega_m)$$ $$\phi^{(k+1)}(i\omega_n) = \pi T \sum_{m} \frac{\phi^{(k)}(i\omega_m)}{\sqrt{\omega_m^2 + (\phi^{(k)}/Z^{(k)})^2}}[\lambda(i\omega_n - i\omega_m) - \mu^*]$$
3. Converge: Repeat until $|\phi^{(k+1)} - \phi^{(k)}| < \epsilon$ for all $\omega_n$

The critical temperature $T_c$ is found by lowering $T$ until a non-trivial solution $\phi(i\omega_n) \neq 0$ appears.

3.2 Model Spectral Functions

For pedagogical purposes and quick estimates, two model forms for $\alpha^2F(\omega)$ are commonly used:

Real materials have more complex spectra measured by tunneling spectroscopy (see Section 4).

3.3 Python Implementation: Eliashberg Solver

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

class EliashbergSolver:
    """
    Solve Eliashberg equations for strong-coupling superconductors.
    Uses Matsubara formalism with iterative convergence.
    """
    def __init__(self, lambda_ep, omega_D, mu_star, T, n_matsubara=100):
        """
        Parameters:
        -----------
        lambda_ep : float
            Electron-phonon coupling constant
        omega_D : float
            Debye frequency (meV)
        mu_star : float
            Coulomb pseudopotential
        T : float
            Temperature (K)
        n_matsubara : int
            Number of Matsubara frequencies to include
        """
        self.lambda_ep = lambda_ep
        self.omega_D = omega_D
        self.mu_star = mu_star
        self.T = T
        self.n_matsubara = n_matsubara

        # Constants
        self.kB = 0.0861733  # meV/K (Boltzmann constant)

        # Matsubara frequencies (fermionic)
        n = np.arange(-n_matsubara, n_matsubara+1)
        self.omega_n = np.pi * self.kB * self.T * (2*n + 1)  # meV

    def alpha2F_debye(self, omega):
        """Debye model spectral function"""
        if omega < 0 or omega > self.omega_D:
            return 0.0
        return 3 * self.lambda_ep * omega**2 / (2 * self.omega_D**3)

    def compute_kernel(self, nu):
        """Compute electron-phonon kernel λ(iν)"""
        def integrand(omega):
            return omega * self.alpha2F_debye(omega) / (nu**2 + omega**2)

        result, _ = quad(integrand, 0, self.omega_D)
        return 2 * result

    def solve_eliashberg(self, max_iter=100, tol=1e-6):
        """
        Iteratively solve Eliashberg equations.

        Returns:
        --------
        phi : array
            Renormalized gap function φ(iωₙ) (meV)
        Z : array
            Renormalization function Z(iωₙ)
        Delta : array
            Physical gap Δ(iωₙ) = φ(iωₙ)/Z(iωₙ) (meV)
        converged : bool
            Whether iteration converged
        """
        T_meV = self.kB * self.T

        # Initialize with BCS-like values
        phi = np.ones_like(self.omega_n) * T_meV  # Initial gap ~ kB*T
        Z = np.ones_like(self.omega_n)

        # Precompute kernels (symmetric matrix)
        n = len(self.omega_n)
        kernel = np.zeros((n, n))
        for i in range(n):
            for j in range(n):
                nu = self.omega_n[i] - self.omega_n[j]
                kernel[i,j] = self.compute_kernel(nu)

        # Iterative solution
        for iteration in range(max_iter):
            phi_old = phi.copy()
            Z_old = Z.copy()

            # Compute denominator for all frequencies
            Delta = phi / Z
            denom = np.sqrt(self.omega_n**2 + Delta**2)

            # Update Z(iωₙ)
            for i in range(n):
                if abs(self.omega_n[i]) > 1e-10:
                    Z_sum = np.sum(np.sign(self.omega_n) * kernel[i,:] / denom)
                    Z[i] = 1 + (np.pi * T_meV / self.omega_n[i]) * Z_sum
                else:
                    Z[i] = 1.0

            # Update φ(iωₙ)
            for i in range(n):
                # Coulomb term (cutoff at ω_D)
                mu_eff = np.where(np.abs(self.omega_n) < self.omega_D,
                                  self.mu_star, 0.0)

                interaction = kernel[i,:] - mu_eff
                phi_sum = np.sum(phi / denom * interaction)
                phi[i] = np.pi * T_meV * phi_sum

            # Check convergence
            phi_err = np.max(np.abs(phi - phi_old))
            Z_err = np.max(np.abs(Z - Z_old))

            if phi_err < tol and Z_err < tol:
                print(f"Converged in {iteration+1} iterations")
                print(f"φ(0) = {phi[n//2]:.3f} meV")
                print(f"Z(0) = {Z[n//2]:.3f}")
                Delta = phi / Z
                return phi, Z, Delta, True

        print(f"Warning: Not converged after {max_iter} iterations")
        Delta = phi / Z
        return phi, Z, Delta, False

    def find_Tc(self, T_min=1, T_max=50, tol=0.1):
        """
        Find critical temperature by bisection.

        Returns:
        --------
        Tc : float
            Critical temperature (K)
        """
        def has_solution(T):
            """Check if non-trivial solution exists at temperature T"""
            self.T = T
            n = np.arange(-self.n_matsubara, self.n_matsubara+1)
            self.omega_n = np.pi * self.kB * T * (2*n + 1)

            phi, Z, Delta, _ = self.solve_eliashberg(max_iter=50, tol=1e-5)
            # Non-trivial if gap is significant
            return np.max(np.abs(Delta)) > 0.01

        # Bisection search
        T_low, T_high = T_min, T_max

        while T_high - T_low > tol:
            T_mid = (T_low + T_high) / 2

            if has_solution(T_mid):
                T_low = T_mid  # Solution exists, Tc is higher
            else:
                T_high = T_mid  # No solution, Tc is lower

        return (T_low + T_high) / 2


# Example 1: Solve for Pb-like parameters
print("="*60)
print("Example 1: Lead (Pb) - Strong Coupling Superconductor")
print("="*60)

solver_Pb = EliashbergSolver(
    lambda_ep=1.55,   # Strong coupling
    omega_D=8.0,      # Debye frequency ~ 8 meV
    mu_star=0.13,     # Coulomb pseudopotential
    T=7.2,            # Around Tc
    n_matsubara=50
)

phi_Pb, Z_Pb, Delta_Pb, converged = solver_Pb.solve_eliashberg()

# Plot results
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# Gap function
axes[0].plot(solver_Pb.omega_n, Delta_Pb, 'o-', linewidth=2)
axes[0].axhline(y=0, color='k', linestyle='--', alpha=0.3)
axes[0].set_xlabel('Matsubara Frequency ωₙ (meV)', fontsize=11)
axes[0].set_ylabel('Gap Δ(iωₙ) (meV)', fontsize=11)
axes[0].set_title('Gap Function (Pb)', fontsize=12, fontweight='bold')
axes[0].grid(True, alpha=0.3)

# Renormalization
axes[1].plot(solver_Pb.omega_n, Z_Pb, 'o-', color='red', linewidth=2)
axes[1].axhline(y=1, color='k', linestyle='--', alpha=0.3, label='BCS (Z=1)')
axes[1].set_xlabel('Matsubara Frequency ωₙ (meV)', fontsize=11)
axes[1].set_ylabel('Renormalization Z(iωₙ)', fontsize=11)
axes[1].set_title('Mass Renormalization (Pb)', fontsize=12, fontweight='bold')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# φ vs Δ
axes[2].plot(solver_Pb.omega_n, phi_Pb, 'o-', label='φ(iωₙ)', linewidth=2)
axes[2].plot(solver_Pb.omega_n, Delta_Pb, 's-', label='Δ(iωₙ) = φ/Z', linewidth=2)
axes[2].axhline(y=0, color='k', linestyle='--', alpha=0.3)
axes[2].set_xlabel('Matsubara Frequency ωₙ (meV)', fontsize=11)
axes[2].set_ylabel('Gap Functions (meV)', fontsize=11)
axes[2].set_title('Renormalized vs Physical Gap (Pb)', fontsize=12, fontweight='bold')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('eliashberg_Pb.png', dpi=150, bbox_inches='tight')
print("\nPlot saved: eliashberg_Pb.png")

# Calculate gap ratio
Delta_0 = Delta_Pb[len(Delta_Pb)//2]  # Gap at ω=0
gap_ratio = 2 * Delta_0 / (solver_Pb.kB * solver_Pb.T)
print(f"\nGap ratio 2Δ(0)/kBT = {gap_ratio:.2f}")
print(f"BCS prediction: 3.52")
print(f"Enhancement: {gap_ratio/3.52:.2f}x")

4. The Electron-Phonon Spectral Function $\alpha^2F(\omega)$

4.1 Physical Meaning

The spectral function $\alpha^2F(\omega)$ represents the weighted phonon density of states, where the weighting is by the electron-phonon coupling strength:

$$\alpha^2F(\omega) = \sum_{\nu} \alpha^2_{\nu}F_{\nu}(\omega)$$

Here $\alpha_{\nu}$ is the coupling constant for phonon branch $\nu$ and $F_{\nu}(\omega)$ is the phonon density of states for that branch. This function encapsulates:

• Which phonons couple strongly: Peaks indicate phonon modes that contribute most to pairing
• Energy scales: The width shows the range of relevant phonon energies
• Overall strength: The integral $\lambda = 2\int \alpha^2F(\omega)/\omega \, d\omega$ gives total coupling

4.2 Experimental Measurement: Tunneling Spectroscopy

The most direct way to measure $\alpha^2F(\omega)$ is through superconductor-insulator-normal metal (S-I-N) tunneling. The differential conductance $dI/dV$ in such a junction is related to the density of states in the superconductor:

$$\frac{dI}{dV}(V) \propto N_S(eV)$$

where $N_S(E)$ is the superconducting density of states. In BCS theory, this would simply be:

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

However, in strong-coupling materials, phonon structure appears as additional features. The McMillan-Rowell inversion procedure extracts $\alpha^2F(\omega)$ from the measured $d^2I/dV^2$ (second derivative of the I-V curve).

The key relation (simplified) is:

$$\alpha^2F(\omega) \propto \int_{-\infty}^{\infty} \frac{d^2I}{dV^2}(V) K(eV, \omega) dV$$

where $K$ is an inversion kernel. Practically, one sees:

• Sharp features at phonon energies: Peaks in $d^2I/dV^2$ correspond to phonon modes
• Gap enhancement: Stronger coupling leads to larger gap and more structure

4.3 Relation to Phonon DOS

If the electron-phonon coupling were constant across all modes (which it never is), then $\alpha^2F(\omega)$ would be proportional to the phonon density of states $F(\omega)$. In reality:

$$\alpha^2F(\omega) = \langle |g|^2 \rangle_{\omega} F(\omega)$$

where $\langle |g|^2 \rangle_{\omega}$ is the average squared matrix element at frequency $\omega$. Optical phonons often couple more strongly than acoustic phonons due to larger ionic displacements.

4.4 Python Example: Model $\alpha^2F(\omega)$ Functions

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

def alpha2F_einstein(omega, lambda_ep, omega_0):
    """Einstein model: single phonon frequency"""
    return lambda_ep * omega_0 / 2 * (np.abs(omega - omega_0) < 1e-6)

def alpha2F_debye(omega, lambda_ep, omega_D):
    """Debye model: uniform density up to ωD"""
    if omega < 0 or omega > omega_D:
        return 0.0
    return 3 * lambda_ep * omega**2 / (2 * omega_D**3)

def alpha2F_double_peak(omega, lambda_1, omega_1, lambda_2, omega_2, gamma=1.0):
    """
    Two-peak model: acoustic + optical phonons
    Uses Lorentzian broadening for realism
    """
    def lorentzian(w, w0, g):
        return (g/np.pi) / ((w - w0)**2 + g**2)

    return (lambda_1 * omega_1 / 2 * lorentzian(omega, omega_1, gamma) +
            lambda_2 * omega_2 / 2 * lorentzian(omega, omega_2, gamma))

def compute_lambda_total(alpha2F_func, omega_max):
    """Compute total coupling constant λ = 2∫α²F(ω)/ω dω"""
    def integrand(omega):
        if omega < 1e-6:
            return 0.0
        return 2 * alpha2F_func(omega) / omega

    result, _ = quad(integrand, 0, omega_max)
    return result

def compute_omega_log(alpha2F_func, omega_max):
    """Compute logarithmic average frequency"""
    def integrand(omega):
        if omega < 1e-6:
            return 0.0
        a2F = alpha2F_func(omega)
        if a2F < 1e-10:
            return 0.0
        return 2 * a2F / omega * np.log(omega)

    numerator, _ = quad(integrand, 0, omega_max)
    lambda_total = compute_lambda_total(alpha2F_func, omega_max)

    return np.exp(numerator / lambda_total)


# Example 2: Compare different α²F(ω) models
print("\n" + "="*60)
print("Example 2: Electron-Phonon Spectral Functions")
print("="*60)

omega = np.linspace(0, 30, 500)  # meV

# Model 1: Debye (Al-like)
alpha2F_Al = np.array([alpha2F_debye(w, lambda_ep=0.43, omega_D=25) for w in omega])

# Model 2: Two peaks (Pb-like: acoustic + optical)
alpha2F_Pb = np.array([alpha2F_double_peak(w, lambda_1=0.8, omega_1=4,
                                            lambda_2=0.75, omega_2=8, gamma=0.5)
                       for w in omega])

# Model 3: Einstein (simple optical mode)
omega_0 = 10
alpha2F_Einstein = np.zeros_like(omega)
idx = np.argmin(np.abs(omega - omega_0))
alpha2F_Einstein[idx] = 1.5 * omega_0 / 2  # λ = 1.5

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

# Left: Spectral functions
axes[0].plot(omega, alpha2F_Al, label='Debye (Al-like, λ=0.43)', linewidth=2)
axes[0].plot(omega, alpha2F_Pb, label='Two-peak (Pb-like, λ≈1.5)', linewidth=2)
axes[0].stem([omega_0], [alpha2F_Einstein[idx]], linefmt='g-', markerfmt='go',
             basefmt=' ', label='Einstein (ω₀=10 meV, λ=1.5)')
axes[0].set_xlabel('Phonon Frequency ω (meV)', fontsize=11)
axes[0].set_ylabel('α²F(ω)', fontsize=11)
axes[0].set_title('Electron-Phonon Spectral Functions', fontsize=12, fontweight='bold')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(0, 30)

# Right: Cumulative λ(ω) = 2∫₀^ω α²F/ω' dω'
lambda_Al_cumul = []
lambda_Pb_cumul = []

for w_max in omega:
    if w_max < 0.1:
        lambda_Al_cumul.append(0)
        lambda_Pb_cumul.append(0)
        continue

    omega_int = np.linspace(0.1, w_max, 100)
    alpha2F_Al_int = np.array([alpha2F_debye(w, 0.43, 25) for w in omega_int])
    alpha2F_Pb_int = np.array([alpha2F_double_peak(w, 0.8, 4, 0.75, 8, 0.5)
                               for w in omega_int])

    lambda_Al_cumul.append(np.trapz(2 * alpha2F_Al_int / omega_int, omega_int))
    lambda_Pb_cumul.append(np.trapz(2 * alpha2F_Pb_int / omega_int, omega_int))

axes[1].plot(omega, lambda_Al_cumul, label='Al (weak coupling)', linewidth=2)
axes[1].plot(omega, lambda_Pb_cumul, label='Pb (strong coupling)', linewidth=2)
axes[1].axhline(y=1, color='k', linestyle='--', alpha=0.5, label='λ = 1 (BCS limit)')
axes[1].set_xlabel('Integration Upper Limit ω (meV)', fontsize=11)
axes[1].set_ylabel('Cumulative λ(ω)', fontsize=11)
axes[1].set_title('Coupling Constant vs Phonon Cutoff', fontsize=12, fontweight='bold')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim(0, 30)

plt.tight_layout()
plt.savefig('alpha2F_comparison.png', dpi=150, bbox_inches='tight')
print("\nPlot saved: alpha2F_comparison.png")

# Compute characteristic frequencies
print("\nCharacteristic Frequencies:")
print(f"Al (Debye): λ_total = {lambda_Al_cumul[-1]:.3f}")
print(f"Pb (Two-peak): λ_total = {lambda_Pb_cumul[-1]:.3f}")

5. $T_c$ Formulas: From Eliashberg to Practice

5.1 The McMillan Formula

While the full Eliashberg equations must be solved numerically, McMillan (1968) derived a remarkably accurate approximate formula for $T_c$ in terms of $\lambda$ and a characteristic phonon frequency $\omega_{\text{log}}$:

The logarithmic average frequency is defined through:

$$\omega_{\text{log}} = \exp\left[\frac{2}{\lambda}\int_0^{\omega_{\text{max}}} \frac{d\omega}{\omega} \alpha^2F(\omega)\ln\omega\right]$$

This formula is accurate to within 10-20% for most phonon-mediated superconductors with $\lambda < 2$.

5.2 Allen-Dynes Modification

Allen and Dynes (1975) improved McMillan's formula to better account for the shape of $\alpha^2F(\omega)$:

The frequency moments are:

$$\omega_n = \left[\frac{2}{\lambda}\int_0^{\omega_{\text{max}}} d\omega \frac{\alpha^2F(\omega)}{\omega}\omega^n\right]^{1/n}$$

These formulas are widely used in computational materials science to quickly estimate $T_c$ from first-principles calculations of $\alpha^2F(\omega)$.

5.3 The Coulomb Pseudopotential $\mu^*$

The parameter $\mu^*$ represents the screened Coulomb repulsion between electrons. The bare Coulomb interaction $\mu$ is large ($\sim 0.3$-0.5), but screening and retardation effects reduce it to $\mu^* \approx 0.1$-0.15 in most materials.

The Morel-Anderson formula gives:

$$\mu^* = \frac{\mu}{1 + \mu\ln(E_F/\omega_{\text{ph}})}$$

where $E_F$ is the Fermi energy and $\omega_{\text{ph}}$ is a characteristic phonon energy. The logarithm reflects that Coulomb interaction acts at all energies up to $E_F$, while phonon-mediated attraction only acts below $\omega_{\text{ph}}$. The ratio $E_F/\omega_{\text{ph}} \sim 10$-100 provides significant screening.

5.4 Isotope Effect in Strong Coupling

In BCS theory, $T_c \propto \omega_D \propto M^{-1/2}$, giving isotope exponent $\alpha = 0.5$. In strong coupling, the McMillan formula predicts:

$$\alpha = -\frac{d\ln T_c}{d\ln M}$$

Taking the derivative of the McMillan formula (assuming $\lambda$ is mass-independent and $\omega_{\text{log}} \propto M^{-1/2}$):

$$\alpha = \frac{1}{2}\left[1 - \frac{\mu^*(1 + 0.62\lambda)}{1 + \lambda}\right]$$

For $\lambda \gg 1$ (strong coupling), $\alpha \to 0$. This explains why Nb$_3$Sn with $\lambda \approx 2$ has $\alpha \approx 0.1$ instead of 0.5.

5.5 Python Example: McMillan Formula Calculator

import numpy as np
import matplotlib.pyplot as plt

def mcmillan_Tc(lambda_ep, omega_log, mu_star):
    """
    Calculate Tc using McMillan's formula.

    Parameters:
    -----------
    lambda_ep : float
        Electron-phonon coupling constant
    omega_log : float
        Logarithmic average phonon frequency (K)
    mu_star : float
        Coulomb pseudopotential

    Returns:
    --------
    Tc : float
        Critical temperature (K)
    """
    numerator = 1.04 * (1 + lambda_ep)
    denominator = lambda_ep - mu_star * (1 + 0.62 * lambda_ep)

    if denominator <= 0:
        return 0.0  # No superconductivity

    exponent = -numerator / denominator
    Tc = (omega_log / 1.2) * np.exp(exponent)

    return Tc

def isotope_exponent(lambda_ep, mu_star):
    """Calculate isotope exponent from strong-coupling theory"""
    return 0.5 * (1 - mu_star * (1 + 0.62*lambda_ep) / (1 + lambda_ep))


# Example 3: McMillan Tc predictions for various materials
print("\n" + "="*60)
print("Example 3: McMillan Formula - Tc Predictions")
print("="*60)

materials = {
    'Al': {'lambda': 0.43, 'omega_log': 375, 'mu_star': 0.15, 'Tc_exp': 1.18},
    'Sn': {'lambda': 0.72, 'omega_log': 165, 'mu_star': 0.13, 'Tc_exp': 3.72},
    'Pb': {'lambda': 1.55, 'omega_log': 109, 'mu_star': 0.13, 'Tc_exp': 7.20},
    'Hg': {'lambda': 1.60, 'omega_log': 70, 'mu_star': 0.15, 'Tc_exp': 4.15},
    'Nb': {'lambda': 1.05, 'omega_log': 240, 'mu_star': 0.12, 'Tc_exp': 9.25},
    'V₃Si': {'lambda': 1.30, 'omega_log': 420, 'mu_star': 0.14, 'Tc_exp': 17.1},
}

print(f"\n{'Material':<10} {'λ':<6} {'ωlog(K)':<10} {'μ*':<6} {'Tc_pred(K)':<12} {'Tc_exp(K)':<12} {'Error':<10} {'α_iso':<8}")
print("-" * 85)

for name, params in materials.items():
    Tc_pred = mcmillan_Tc(params['lambda'], params['omega_log'], params['mu_star'])
    error = 100 * (Tc_pred - params['Tc_exp']) / params['Tc_exp']
    alpha = isotope_exponent(params['lambda'], params['mu_star'])

    print(f"{name:<10} {params['lambda']:<6.2f} {params['omega_log']:<10.0f} "
          f"{params['mu_star']:<6.2f} {Tc_pred:<12.2f} {params['Tc_exp']:<12.2f} "
          f"{error:>+6.1f}%   {alpha:<8.3f}")

# Example 4: Tc vs λ parameter space
print("\n" + "="*60)
print("Example 4: Tc Dependence on λ and μ*")
print("="*60)

lambda_range = np.linspace(0.1, 3.0, 100)
omega_log = 200  # K

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

# Left: Tc vs λ for different μ*
mu_star_values = [0.0, 0.10, 0.15, 0.20]
for mu_star in mu_star_values:
    Tc_values = [mcmillan_Tc(lam, omega_log, mu_star) for lam in lambda_range]
    axes[0].plot(lambda_range, Tc_values, linewidth=2,
                 label=f'μ* = {mu_star:.2f}')

axes[0].set_xlabel('Coupling Constant λ', fontsize=11)
axes[0].set_ylabel('Critical Temperature Tc (K)', fontsize=11)
axes[0].set_title('McMillan Formula: Tc vs λ (ωlog = 200 K)',
                  fontsize=12, fontweight='bold')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(0, 3)

# Mark material points
materials_plot = {'Al': (0.43, 1.18), 'Pb': (1.55, 7.20), 'Nb': (1.05, 9.25)}
for name, (lam, Tc) in materials_plot.items():
    axes[0].plot(lam, Tc, 'ro', markersize=8)
    axes[0].annotate(name, xy=(lam, Tc), xytext=(5, 5),
                    textcoords='offset points', fontsize=10)

# Right: Isotope effect vs λ
alpha_values = [isotope_exponent(lam, 0.13) for lam in lambda_range]
axes[1].plot(lambda_range, alpha_values, linewidth=2.5, color='purple')
axes[1].axhline(y=0.5, color='k', linestyle='--', alpha=0.5, label='BCS (α=0.5)')
axes[1].set_xlabel('Coupling Constant λ', fontsize=11)
axes[1].set_ylabel('Isotope Exponent α', fontsize=11)
axes[1].set_title('Isotope Effect in Strong Coupling (μ*=0.13)',
                  fontsize=12, fontweight='bold')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim(0, 3)
axes[1].set_ylim(0, 0.6)

plt.tight_layout()
plt.savefig('mcmillan_predictions.png', dpi=150, bbox_inches='tight')
print("\nPlot saved: mcmillan_predictions.png")

print("\nKey Observations:")
print("1. Tc increases strongly with λ, but saturates due to μ* (Coulomb repulsion)")
print("2. Larger μ* significantly suppresses Tc, especially at high λ")
print("3. Isotope exponent α decreases from BCS value 0.5 as λ increases")
print("4. McMillan formula predicts experimental Tc within ~20% for most materials")

6. Applications to Real Superconductors

6.1 Lead (Pb): The Prototypical Strong-Coupling Superconductor

Lead is the textbook example of strong-coupling superconductivity:

• $T_c = 7.20$ K
• $\lambda = 1.55$ (strong coupling)
• $2\Delta(0)/k_BT_c = 4.3$ (enhanced gap)
• $\alpha_{\text{isotope}} = 0.48$ (close to BCS, fortuitously)

The $\alpha^2F(\omega)$ spectrum of Pb shows two main peaks:

• Low-frequency peak (~4 meV): Transverse acoustic phonons
• High-frequency peak (~8 meV): Longitudinal optical phonons

Tunneling spectroscopy on Pb revealed clear phonon structure in $dI/dV$, confirming the strong electron-phonon interaction and validating Eliashberg theory.

6.2 Niobium (Nb): A High-$T_c$ Element

Niobium has the highest $T_c$ among elements ($T_c = 9.25$ K):

• $\lambda \approx 1.05$ (moderate-to-strong coupling)
• High density of states $N(0)$ due to $d$-band
• $\omega_{\text{log}} \approx 240$ K (high phonon frequencies)

Nb is widely used in superconducting RF cavities for particle accelerators due to its high $T_c$ and excellent surface properties.

6.3 A15 Compounds (Nb$_3$Sn, V$_3$Si): Record $T_c$ Before Cuprates

The A15 crystal structure (Cr$_3$Si prototype) hosts several high-$T_c$ superconductors:

Key features:

• Very high $\lambda$: Among the strongest electron-phonon coupling known
• Structural instability: Close to a martensitic transition, softening phonon modes
• Reduced isotope effect: $\alpha \approx 0.1$ for Nb$_3$Sn due to large $\lambda$

6.4 Hydride Superconductors: Pushing the Limits

Recent high-pressure experiments have discovered superconductivity at record-breaking temperatures in compressed hydrides:

These materials achieve extreme $T_c$ values through:

• High phonon frequencies: Light hydrogen atoms give $\omega_D \sim 2000$ K
• Strong coupling: $\lambda \sim 2$-3 from metal-H bond stretching modes
• High density of states: Pressure increases $N(0)$

DFT calculations of $\alpha^2F(\omega)$ for LaH$_{10}$ show strong coupling to H vibrations at 150-200 meV, explaining the high $T_c$. These discoveries renewed interest in room-temperature superconductivity under pressure.

6.5 Python Example: Predicting Hydride $T_c$

import numpy as np
import matplotlib.pyplot as plt

# Example 5: High-Tc hydride superconductors
print("\n" + "="*60)
print("Example 5: Hydride Superconductors - Extreme Strong Coupling")
print("="*60)

# Hydride parameters (from DFT calculations and experiments)
hydrides = {
    'H₃S (155 GPa)': {
        'lambda': 2.0,
        'omega_log': 1550,  # K (high due to light H)
        'mu_star': 0.13,
        'Tc_exp': 203
    },
    'LaH₁₀ (200 GPa)': {
        'lambda': 2.65,
        'omega_log': 1600,
        'mu_star': 0.13,
        'Tc_exp': 260
    },
    'CaH₆ (172 GPa)': {
        'lambda': 2.2,
        'omega_log': 1400,
        'mu_star': 0.13,
        'Tc_exp': 215
    },
}

print(f"\n{'Material':<20} {'λ':<6} {'ωlog(K)':<10} {'Tc_pred(K)':<12} {'Tc_exp(K)':<12} {'Error':<10}")
print("-" * 75)

for name, params in hydrides.items():
    Tc_pred = mcmillan_Tc(params['lambda'], params['omega_log'], params['mu_star'])
    error = 100 * (Tc_pred - params['Tc_exp']) / params['Tc_exp']

    print(f"{name:<20} {params['lambda']:<6.2f} {params['omega_log']:<10.0f} "
          f"{Tc_pred:<12.1f} {params['Tc_exp']:<12.0f} {error:>+6.1f}%")

# Plot: Path to room-temperature SC
fig, ax = plt.subplots(figsize=(12, 7))

# Historical progression
materials_history = [
    ('Hg (1911)', 1911, 4.15),
    ('Nb (1930)', 1930, 9.25),
    ('Nb₃Sn (1954)', 1954, 18.1),
    ('Nb₃Ge (1973)', 1973, 23.2),
    ('La-Ba-Cu-O (1986)', 1986, 35),
    ('Y-Ba-Cu-O (1987)', 1987, 93),
    ('Hg-Ba-Ca-Cu-O (1993)', 1993, 135),
    ('H₃S (2015)', 2015, 203),
    ('LaH₁₀ (2019)', 2019, 260),
]

years = [y for _, y, _ in materials_history]
Tcs = [tc for _, _, tc in materials_history]
names = [n for n, _, _ in materials_history]

# Color code by type
colors = []
for name, _, _ in materials_history:
    if 'Cu' in name:
        colors.append('red')  # Cuprates
    elif 'H' in name and name != 'Hg (1911)':
        colors.append('green')  # Hydrides
    else:
        colors.append('blue')  # Conventional

ax.scatter(years, Tcs, s=200, c=colors, alpha=0.7, edgecolors='black', linewidth=2)

# Annotate points
for i, (name, year, tc) in enumerate(materials_history):
    if i % 2 == 0:
        xytext = (10, 15)
    else:
        xytext = (10, -15)
    ax.annotate(name, xy=(year, tc), xytext=xytext,
                textcoords='offset points', fontsize=9,
                bbox=dict(boxstyle='round,pad=0.3', facecolor=colors[i], alpha=0.3),
                arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=0'))

# Room temperature line
ax.axhline(y=273, color='orange', linestyle='--', linewidth=2,
           label='Room Temperature (273 K)', alpha=0.7)

# Liquid nitrogen line
ax.axhline(y=77, color='cyan', linestyle='--', linewidth=2,
           label='Liquid N₂ (77 K)', alpha=0.7)

ax.set_xlabel('Year of Discovery', fontsize=12)
ax.set_ylabel('Critical Temperature Tc (K)', fontsize=12)
ax.set_title('Historical Progress Toward Room-Temperature Superconductivity',
             fontsize=13, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(1900, 2030)
ax.set_ylim(0, 300)

# Legend for colors
from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor='blue', alpha=0.7, label='Conventional (phonon-mediated)'),
    Patch(facecolor='red', alpha=0.7, label='Cuprates (d-wave, unclear mechanism)'),
    Patch(facecolor='green', alpha=0.7, label='Hydrides (extreme coupling)')
]
ax.legend(handles=legend_elements, loc='upper left', fontsize=10)

plt.tight_layout()
plt.savefig('Tc_history.png', dpi=150, bbox_inches='tight')
print("\nPlot saved: Tc_history.png")

print("\nKey Milestones:")
print("- 1911: Discovery of superconductivity (Hg, 4.15 K)")
print("- 1973: Nb₃Ge reaches 23.2 K (pre-cuprate record)")
print("- 1986: Cuprate era begins, mechanism still debated")
print("- 2015: H₃S breaks 200 K barrier under pressure")
print("- 2019: LaH₁₀ reaches 260 K, ~13°C below room temperature!")

7. Temperature Dependence and Real-Axis Continuation

7.1 Gap Function Temperature Evolution

The Eliashberg equations naturally describe how the gap $\Delta(T)$ evolves with temperature. As $T \to T_c$ from below, the gap vanishes continuously (second-order phase transition):

$$\Delta(T) \approx \Delta(0)\sqrt{1 - T/T_c} \quad \text{near } T_c$$

At $T = 0$, the gap reaches its maximum value $\Delta(0)$. The full temperature dependence must be obtained by solving Eliashberg equations at each $T$.

7.2 From Matsubara to Real Frequencies

Physical quantities like density of states and optical conductivity require the gap on the real frequency axis $\Delta(\omega)$ (real $\omega$), not the imaginary Matsubara axis $\Delta(i\omega_n)$. The continuation process is non-trivial:

$$\Delta(i\omega_n) \xrightarrow{\text{analytic continuation}} \Delta(\omega + i\delta)$$

Methods include:

• Padé approximants: Rational function fitting
• Maximum entropy: Bayesian inference approach
• Direct real-axis solution: Solve Eliashberg equations directly on real axis (more complex due to branch cuts)

7.3 Density of States and Tunneling Conductance

Once $\Delta(\omega)$ is known on the real axis, the density of states is:

$$N_S(\omega) = N(0)\text{Re}\left[\frac{\omega Z(\omega)}{\sqrt{(\omega Z(\omega))^2 - \Delta^2(\omega)}}\right]$$

This generalizes the BCS result and includes phonon structure. Tunneling conductance is proportional to $N_S(\omega)$:

$$\frac{dI}{dV}(V) \propto N_S(eV)$$

7.4 Python Example: Temperature-Dependent Gap

import numpy as np
import matplotlib.pyplot as plt

# Example 6: Temperature dependence of gap function
print("\n" + "="*60)
print("Example 6: Gap Temperature Dependence - Pb vs Al")
print("="*60)

def solve_gap_at_T(lambda_ep, omega_D, mu_star, T, n_matsubara=50, max_iter=100):
    """
    Solve Eliashberg equations at given temperature.
    Returns gap at ω=0.
    """
    if T < 0.1:
        T = 0.1  # Avoid T=0 singularities

    kB = 0.0861733  # meV/K
    T_meV = kB * T

    # Matsubara frequencies
    n = np.arange(-n_matsubara, n_matsubara+1)
    omega_n = np.pi * T_meV * (2*n + 1)

    # Initialize
    phi = np.ones_like(omega_n) * T_meV
    Z = np.ones_like(omega_n)

    # Precompute kernels (Debye model)
    def compute_kernel(nu):
        omega = np.linspace(0, omega_D, 100)
        alpha2F = 3 * lambda_ep * omega**2 / (2 * omega_D**3)
        integrand = omega * alpha2F / (nu**2 + omega**2)
        return 2 * np.trapz(integrand, omega)

    kernel = np.zeros((len(omega_n), len(omega_n)))
    for i in range(len(omega_n)):
        for j in range(len(omega_n)):
            kernel[i,j] = compute_kernel(omega_n[i] - omega_n[j])

    # Iterate
    for iteration in range(max_iter):
        phi_old = phi.copy()
        Delta = phi / Z
        denom = np.sqrt(omega_n**2 + Delta**2)

        # Update Z
        for i in range(len(omega_n)):
            if abs(omega_n[i]) > 1e-10:
                Z[i] = 1 + (np.pi * T_meV / omega_n[i]) * np.sum(
                    np.sign(omega_n) * kernel[i,:] / denom)
            else:
                Z[i] = 1.0

        # Update φ
        for i in range(len(omega_n)):
            mu_eff = np.where(np.abs(omega_n) < omega_D, mu_star, 0.0)
            phi[i] = np.pi * T_meV * np.sum(
                phi / denom * (kernel[i,:] - mu_eff))

        # Check convergence
        if np.max(np.abs(phi - phi_old)) < 1e-6:
            Delta = phi / Z
            return Delta[len(Delta)//2]  # Gap at ω=0

    Delta = phi / Z
    return Delta[len(Delta)//2]


# Compute Δ(T) for Pb and Al
T_range = np.linspace(0.1, 10, 50)

print("\nSolving temperature-dependent gap for Pb...")
Delta_Pb_vs_T = []
for T in T_range:
    if T > 7.5:  # Above Tc for Pb
        Delta_Pb_vs_T.append(0.0)
    else:
        Delta = solve_gap_at_T(lambda_ep=1.55, omega_D=8.0, mu_star=0.13, T=T)
        Delta_Pb_vs_T.append(max(0, Delta))

print("Solving temperature-dependent gap for Al...")
Delta_Al_vs_T = []
for T in T_range:
    if T > 1.5:  # Above Tc for Al
        Delta_Al_vs_T.append(0.0)
    else:
        Delta = solve_gap_at_T(lambda_ep=0.43, omega_D=25.0, mu_star=0.15, T=T)
        Delta_Al_vs_T.append(max(0, Delta))

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

# Left: Absolute gap values
axes[0].plot(T_range, Delta_Pb_vs_T, 'o-', linewidth=2, label='Pb (λ=1.55)')
axes[0].plot(T_range, Delta_Al_vs_T, 's-', linewidth=2, label='Al (λ=0.43)')
axes[0].axvline(x=7.2, color='blue', linestyle='--', alpha=0.5, label='Tc(Pb) = 7.2 K')
axes[0].axvline(x=1.18, color='orange', linestyle='--', alpha=0.5, label='Tc(Al) = 1.18 K')
axes[0].set_xlabel('Temperature T (K)', fontsize=11)
axes[0].set_ylabel('Gap Δ(0,T) (meV)', fontsize=11)
axes[0].set_title('Temperature-Dependent Gap: Pb vs Al', fontsize=12, fontweight='bold')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Right: Normalized gap Δ(T)/Δ(0)
Delta_Pb_0 = Delta_Pb_vs_T[0] if Delta_Pb_vs_T[0] > 0 else 1.0
Delta_Al_0 = Delta_Al_vs_T[0] if Delta_Al_vs_T[0] > 0 else 1.0

Delta_Pb_norm = np.array(Delta_Pb_vs_T) / Delta_Pb_0
Delta_Al_norm = np.array(Delta_Al_vs_T) / Delta_Al_0

# Temperature normalized by Tc
T_Pb_norm = T_range / 7.2
T_Al_norm = T_range / 1.18

# BCS universal curve
T_BCS = np.linspace(0, 1, 100)
Delta_BCS = np.tanh(1.74 * np.sqrt((1/T_BCS - 1)))  # BCS approximation

axes[1].plot(T_Pb_norm, Delta_Pb_norm, 'o-', linewidth=2, label='Pb (strong coupling)')
axes[1].plot(T_Al_norm, Delta_Al_norm, 's-', linewidth=2, label='Al (weak coupling)')
axes[1].plot(T_BCS, Delta_BCS, 'k--', linewidth=2, label='BCS theory', alpha=0.7)
axes[1].set_xlabel('Reduced Temperature T/Tc', fontsize=11)
axes[1].set_ylabel('Normalized Gap Δ(T)/Δ(0)', fontsize=11)
axes[1].set_title('Universal Gap Behavior', fontsize=12, fontweight='bold')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim(0, 1.2)
axes[1].set_ylim(0, 1.2)

plt.tight_layout()
plt.savefig('gap_temperature_dependence.png', dpi=150, bbox_inches='tight')
print("\nPlot saved: gap_temperature_dependence.png")

print("\nResults:")
print(f"Pb: Δ(0) = {Delta_Pb_vs_T[0]:.3f} meV, 2Δ(0)/kBTc = {2*Delta_Pb_vs_T[0]/(0.0861733*7.2):.2f}")
print(f"Al: Δ(0) = {Delta_Al_vs_T[0]:.3f} meV, 2Δ(0)/kBTc = {2*Delta_Al_vs_T[0]/(0.0861733*1.18):.2f}")
print(f"BCS prediction: 2Δ(0)/kBTc = 3.52")

8. Advanced Topics and Extensions

8.1 Anisotropic Eliashberg Theory

The isotropic approximation breaks down in materials with strong anisotropy (e.g., layered cuprates, MgB$_2$). The full anisotropic Eliashberg equations depend on momentum direction $\mathbf{k}$ on the Fermi surface:

$$\Delta(\mathbf{k}, i\omega_n) = \pi T \sum_m \int \frac{d\Omega_{\mathbf{k}'}}{4\pi} \frac{\Delta(\mathbf{k}', i\omega_m)}{...} V(\mathbf{k}, \mathbf{k}', i\omega_n - i\omega_m)$$

where $V(\mathbf{k}, \mathbf{k}', i\nu)$ is the momentum-dependent pairing interaction. This is computationally demanding but essential for:

• MgB$_2$: Two-band superconductor with different gaps on $\sigma$ and $\pi$ bands
• Iron-based superconductors: Multi-band with possible sign changes

8.2 Beyond Migdal's Theorem

In some exotic materials, Migdal's theorem may fail:

• Light atoms + strong coupling: $(m/M)^{1/2}$ not small enough
• Soft phonon modes: Near structural instabilities, $\omega_{\text{ph}} \to 0$
• Non-adiabatic systems: Electron and phonon timescales comparable

Vertex corrections can then become important, requiring more sophisticated many-body techniques (parquet equations, functional renormalization group).

8.3 Multiband Eliashberg Theory

Materials like MgB$_2$ require separate Eliashberg equations for each band:

$$\Delta_i(i\omega_n) = \pi T \sum_{j,m} \lambda_{ij}(...)\Delta_j(i\omega_m)$$

where $\lambda_{ij}$ is the inter-band coupling matrix. This leads to:

• Multiple gaps (e.g., $\Delta_\sigma \approx 7$ meV, $\Delta_\pi \approx 2$ meV in MgB$_2$)
• Enhanced $T_c$ from inter-band pairing

9. Summary and Key Takeaways

10. Exercises and Further Exploration

Exercise 1: Weak vs Strong Coupling

Compare BCS and Eliashberg predictions for a material with $\lambda = 1.2$, $\omega_D = 200$ K, $\mu^* = 0.13$:

1. Calculate $T_c$ using BCS formula: $T_c^{\text{BCS}} = 1.14 \omega_D \exp(-1/\lambda)$
2. Calculate $T_c$ using McMillan formula
3. Compute gap ratio $2\Delta(0)/k_BT_c$ for both theories
4. What is the percentage error if BCS is used?

Exercise 2: Isotope Effect Engineering

Design a hypothetical material to maximize the isotope effect (large $\alpha$):

1. What value of $\lambda$ gives $\alpha = 0.5$ (BCS value) for $\mu^* = 0.13$?
2. Why do A15 compounds have small isotope effects?
3. Could isotope substitution be used to tune $T_c$ in applications?

Exercise 3: Extracting $\alpha^2F(\omega)$ from Data

Given tunneling data (simulated or from literature):

1. Identify phonon peaks in $d^2I/dV^2$ spectrum
2. Construct a model $\alpha^2F(\omega)$ with appropriate peaks
3. Compute $\lambda$ and $\omega_{\text{log}}$
4. Predict $T_c$ and compare with experiment

Exercise 4: Hydride $T_c$ Optimization

Explore parameter space for hydride superconductors:

1. How does $T_c$ scale with pressure (through $\omega_{\text{log}}$ and $\lambda$)?
2. What is the theoretical maximum $T_c$ for $\lambda = 3$, $\omega_{\text{log}} = 2000$ K?
3. Why is high pressure required for hydride superconductivity?

Further Reading

Foundational Papers

• G. M. Eliashberg, Sov. Phys. JETP 11, 696 (1960) - Original Eliashberg theory
• W. L. McMillan, Phys. Rev. 167, 331 (1968) - McMillan $T_c$ formula
• P. B. Allen and R. C. Dynes, Phys. Rev. B 12, 905 (1975) - Improved $T_c$ formula
• A. B. Migdal, Sov. Phys. JETP 7, 996 (1958) - Migdal theorem

Reviews and Textbooks

• J. P. Carbotte, Rev. Mod. Phys. 62, 1027 (1990) - Eliashberg theory review
• P. B. Allen, Handbook of Superconductivity (1999) - Practical guide to electron-phonon coupling
• E. K. H. Salje, Phase Transitions in Ferroelastic and Co-elastic Crystals (1990) - Structural aspects
• M. Tinkham, Introduction to Superconductivity, 2nd ed. (2004), Chapter 3 - Strong-coupling extensions

Computational Methods

• S. Baroni et al., Rev. Mod. Phys. 73, 515 (2001) - DFPT for electron-phonon coupling
• F. Giustino, Rev. Mod. Phys. 89, 015003 (2017) - Modern computational methods
• Quantum ESPRESSO documentation - Practical DFT calculations of $\alpha^2F(\omega)$

Recent Hydride Discoveries

• A. P. Drozdov et al., Nature 525, 73 (2015) - H$_3$S at 203 K
• M. Somayazulu et al., Phys. Rev. Lett. 122, 027001 (2019) - LaH$_{10}$ at 260 K
• I. A. Kruglov et al., Sci. Rep. 9, 10470 (2019) - Theoretical predictions