Chapter 5: Frontiers of Superconductivity Research

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh

class TwistedBilayerGraphene:
    """
    Simplified tight-binding model for twisted bilayer graphene.
    Uses a continuum model approximation near the K point.
    """
    def __init__(self, theta_deg, v_F=1.0, w=0.11):
        """
        Parameters:
        -----------
        theta_deg : float
            Twist angle in degrees
        v_F : float
            Fermi velocity (in units of eV·nm)
        w : float
            Interlayer hopping amplitude (eV)
        """
        self.theta = np.deg2rad(theta_deg)
        self.v_F = v_F
        self.w = w

        # Moiré lattice vectors
        self.a_M = 0.246 / (2 * np.sin(self.theta / 2))  # nm

        # Moiré reciprocal lattice vectors
        self.G_M = 4 * np.pi / (np.sqrt(3) * self.a_M)

    def rotation_matrix(self, angle):
        """2D rotation matrix"""
        return np.array([
            [np.cos(angle), -np.sin(angle)],
            [np.sin(angle), np.cos(angle)]
        ])

    def hamiltonian(self, k):
        """
        Construct the 4x4 Hamiltonian for TBG at momentum k.

        The four components correspond to:
        - Layer 1: A and B sublattices
        - Layer 2: A and B sublattices
        """
        kx, ky = k

        # Rotation matrices for layers
        R1 = self.rotation_matrix(self.theta / 2)
        R2 = self.rotation_matrix(-self.theta / 2)

        # Dirac Hamiltonian for each layer
        k1 = R1 @ k
        k2 = R2 @ k

        H1 = self.v_F * np.array([
            [0, k1[0] - 1j*k1[1]],
            [k1[0] + 1j*k1[1], 0]
        ])

        H2 = self.v_F * np.array([
            [0, k2[0] - 1j*k2[1]],
            [k2[0] + 1j*k2[1], 0]
        ])

        # Interlayer coupling (simplified)
        T = self.w * np.array([
            [1, 1],
            [1, 1]
        ])

        # Full 4x4 Hamiltonian
        H = np.zeros((4, 4), dtype=complex)
        H[0:2, 0:2] = H1
        H[2:4, 2:4] = H2
        H[0:2, 2:4] = T
        H[2:4, 0:2] = T.conj().T

        return H

    def compute_bands(self, k_path, n_points=100):
        """
        Compute band structure along a k-path.

        Parameters:
        -----------
        k_path : list of tuples
            List of high-symmetry points in reciprocal space
        n_points : int
            Number of points between each pair of high-symmetry points

        Returns:
        --------
        k_values : array
            Momentum values along path
        bands : array
            Energy bands (shape: n_points * len(k_path), 4)
        """
        all_k = []
        all_E = []

        for i in range(len(k_path) - 1):
            k_start = np.array(k_path[i])
            k_end = np.array(k_path[i + 1])

            for alpha in np.linspace(0, 1, n_points):
                k = (1 - alpha) * k_start + alpha * k_end
                H = self.hamiltonian(k)
                eigenvalues = eigh(H, eigvals_only=True)

                all_k.append(np.linalg.norm(k - k_path[0]))
                all_E.append(eigenvalues)

        return np.array(all_k), np.array(all_E)

    def plot_bands(self, k_path, labels, n_points=100):
        """Plot band structure along k_path with labeled high-symmetry points."""
        k_values, bands = self.compute_bands(k_path, n_points)

        plt.figure(figsize=(10, 6))
        for band in range(bands.shape[1]):
            plt.plot(k_values, bands[:, band], 'b-', linewidth=1.5)

        # Mark high-symmetry points
        k_positions = [0]
        for i in range(1, len(k_path)):
            dist = np.linalg.norm(np.array(k_path[i]) - np.array(k_path[0]))
            k_positions.append(dist)

        for pos, label in zip(k_positions, labels):
            plt.axvline(pos, color='k', linestyle='--', alpha=0.3)

        plt.xticks(k_positions, labels)
        plt.ylabel('Energy (eV)', fontsize=12)
        plt.xlabel('Momentum', fontsize=12)
        plt.title(f'TBG Band Structure (θ = {np.rad2deg(self.theta):.2f}°)', fontsize=14)
        plt.ylim(-0.3, 0.3)
        plt.axhline(0, color='r', linestyle=':', alpha=0.5, label='Fermi level')
        plt.grid(True, alpha=0.3)
        plt.legend()
        plt.tight_layout()
        return plt.gcf()

# Example: Compare different twist angles
print("=== Twisted Bilayer Graphene Band Structure ===\n")

angles = [1.05, 1.1, 1.2]  # Magic angle is ~1.1°

# Define k-path: Γ -> K -> M -> Γ (in moiré Brillouin zone)
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

for idx, theta in enumerate(angles):
    tbg = TwistedBilayerGraphene(theta_deg=theta)

    # Simplified k-path in moiré BZ
    k_path = [
        (0, 0),                           # Γ
        (tbg.G_M, 0),                     # K_M
        (tbg.G_M/2, tbg.G_M*np.sqrt(3)/2), # M_M
        (0, 0)                            # Γ
    ]
    labels = ['Γ', 'K_M', 'M_M', 'Γ']

    k_values, bands = tbg.compute_bands(k_path, n_points=50)

    plt.sca(axes[idx])
    for band in range(bands.shape[1]):
        plt.plot(k_values, bands[:, band], 'b-', linewidth=1.5)

    plt.axhline(0, color='r', linestyle=':', alpha=0.5)
    plt.ylabel('Energy (eV)', fontsize=10)
    plt.xlabel('Momentum', fontsize=10)
    plt.title(f'θ = {theta}°', fontsize=12, fontweight='bold')
    plt.ylim(-0.3, 0.3)
    plt.grid(True, alpha=0.3)

    # Calculate and display bandwidth
    bandwidth = np.max(bands) - np.min(bands)
    plt.text(0.5, 0.25, f'W ≈ {bandwidth:.3f} eV',
             transform=axes[idx].transAxes, fontsize=10,
             bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

plt.tight_layout()
plt.savefig('tbg_band_structure_comparison.png', dpi=150, bbox_inches='tight')
print("Saved: tbg_band_structure_comparison.png")
print("\nNote: At the magic angle (~1.1°), bands become very flat near E=0")
print("This leads to enhanced electron correlations and superconductivity.\n")

import numpy as np
import matplotlib.pyplot as plt

def bandwidth_vs_angle(theta_range, v_F=1.0, w=0.11):
    """
    Calculate bandwidth as a function of twist angle.
    Uses a simplified analytical model.
    """
    bandwidths = []

    for theta_deg in theta_range:
        theta = np.deg2rad(theta_deg)

        # Simplified bandwidth scaling
        # Near magic angle, W ∝ |theta - theta_magic|
        theta_magic = np.deg2rad(1.08)  # Magic angle

        if abs(theta - theta_magic) < np.deg2rad(0.1):
            # Very close to magic angle - use linear scaling
            W = 0.01 * abs(theta - theta_magic) / np.deg2rad(0.1)
        else:
            # Further from magic angle
            W = 2 * w * abs(np.sin(theta/2))

        bandwidths.append(W * 1000)  # Convert to meV

    return np.array(bandwidths)

# Calculate bandwidth vs angle
theta_range = np.linspace(0.5, 2.0, 200)
bandwidths = bandwidth_vs_angle(theta_range)

# Plotting
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Full range
ax1.plot(theta_range, bandwidths, 'b-', linewidth=2)
ax1.axvline(1.08, color='r', linestyle='--', linewidth=2, label='Magic angle')
ax1.fill_between([1.05, 1.15], 0, max(bandwidths), alpha=0.2, color='red',
                  label='Magic angle region')
ax1.set_xlabel('Twist Angle (degrees)', fontsize=12)
ax1.set_ylabel('Bandwidth (meV)', fontsize=12)
ax1.set_title('Bandwidth vs Twist Angle in TBG', fontsize=14, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.legend(fontsize=10)
ax1.set_ylim(0, max(bandwidths)*1.1)

# Zoomed near magic angle
zoom_range = (theta_range >= 0.9) & (theta_range <= 1.3)
ax2.plot(theta_range[zoom_range], bandwidths[zoom_range], 'b-', linewidth=2)
ax2.axvline(1.08, color='r', linestyle='--', linewidth=2, label='θ_magic ≈ 1.08°')
ax2.axhline(10, color='g', linestyle=':', linewidth=1.5,
            label='W ~ 10 meV (flat band regime)')
ax2.fill_between([1.05, 1.15], 0, max(bandwidths[zoom_range]),
                  alpha=0.2, color='red')
ax2.set_xlabel('Twist Angle (degrees)', fontsize=12)
ax2.set_ylabel('Bandwidth (meV)', fontsize=12)
ax2.set_title('Zoomed View Near Magic Angle', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.legend(fontsize=10)
ax2.set_ylim(0, 50)

plt.tight_layout()
plt.savefig('tbg_magic_angle_bandwidth.png', dpi=150, bbox_inches='tight')
print("=== Flat Band Physics in TBG ===\n")
print("Saved: tbg_magic_angle_bandwidth.png")
print(f"\nKey observations:")
print(f"- Magic angle: θ_M ≈ 1.08° (first magic angle)")
print(f"- At θ_M: Bandwidth W → 0 (flat bands)")
print(f"- Away from θ_M: W increases rapidly")
print(f"- Flat bands → Enhanced correlations → Superconductivity")
print(f"\nPhysical consequence:")
print(f"  Effective mass: m* ∝ 1/W → ∞ as W → 0")
print(f"  Kinetic energy suppressed → Interactions dominate\n")

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh

class KitaevChain:
    """
    Kitaev chain model for 1D topological superconductor.
    H = -μΣ c†_i c_i - t Σ (c†_{i+1} c_i + h.c.) + Δ Σ (c_i c_{i+1} + h.c.)
    """
    def __init__(self, N, t=1.0, mu=0.0, Delta=0.5):
        """
        Parameters:
        -----------
        N : int
            Number of sites
        t : float
            Hopping amplitude
        mu : float
            Chemical potential
        Delta : float
            p-wave pairing amplitude
        """
        self.N = N
        self.t = t
        self.mu = mu
        self.Delta = Delta

    def hamiltonian_BdG(self):
        """
        Construct Bogoliubov-de Gennes Hamiltonian.
        In Nambu basis: (c_1, c_2, ..., c_N, c†_1, c†_2, ..., c†_N)
        """
        H = np.zeros((2*self.N, 2*self.N), dtype=complex)

        # Particle sector
        for i in range(self.N):
            H[i, i] = -self.mu
            if i < self.N - 1:
                H[i, i+1] = -self.t
                H[i+1, i] = -self.t

        # Hole sector (transpose of particle sector with opposite sign)
        H[self.N:, self.N:] = self.mu * np.eye(self.N)
        for i in range(self.N - 1):
            H[self.N+i, self.N+i+1] = self.t
            H[self.N+i+1, self.N+i] = self.t

        # Pairing terms
        for i in range(self.N - 1):
            H[i, self.N+i+1] = self.Delta
            H[i+1, self.N+i] = -self.Delta
            H[self.N+i, i+1] = self.Delta
            H[self.N+i+1, i] = -self.Delta

        return H

    def solve(self):
        """
        Diagonalize the BdG Hamiltonian.
        Returns eigenvalues and eigenvectors.
        """
        H = self.hamiltonian_BdG()
        eigenvalues, eigenvectors = eigh(H)
        return eigenvalues, eigenvectors

    def is_topological(self):
        """
        Determine if the chain is in topological phase.
        Topological when |μ| < 2t and Δ ≠ 0.
        """
        return abs(self.mu) < 2 * self.t and self.Delta != 0

    def plot_spectrum(self):
        """Plot energy spectrum."""
        eigenvalues, _ = self.solve()

        # Only plot positive eigenvalues (particle-hole symmetry)
        positive_E = eigenvalues[eigenvalues >= 0]

        plt.figure(figsize=(10, 6))
        plt.scatter(range(len(positive_E)), positive_E, s=50, alpha=0.7)
        plt.axhline(0, color='r', linestyle='--', linewidth=2, label='E = 0')
        plt.xlabel('State Index', fontsize=12)
        plt.ylabel('Energy', fontsize=12)

        phase = "Topological" if self.is_topological() else "Trivial"
        plt.title(f'Kitaev Chain Spectrum ({phase} Phase)\n' +
                  f'N={self.N}, μ={self.mu:.2f}, t={self.t:.2f}, Δ={self.Delta:.2f}',
                  fontsize=14, fontweight='bold')
        plt.grid(True, alpha=0.3)
        plt.legend(fontsize=10)
        plt.tight_layout()
        return plt.gcf()

    def plot_edge_states(self):
        """Plot wavefunctions of lowest energy states (edge modes)."""
        eigenvalues, eigenvectors = self.solve()

        # Find states closest to zero energy
        sorted_indices = np.argsort(np.abs(eigenvalues))

        fig, axes = plt.subplots(2, 1, figsize=(12, 8))

        for idx in range(2):  # Plot two lowest energy states
            state_idx = sorted_indices[idx]
            E = eigenvalues[state_idx]
            psi = eigenvectors[:self.N, state_idx]  # Particle component

            axes[idx].plot(range(self.N), np.abs(psi)**2, 'bo-', linewidth=2, markersize=8)
            axes[idx].set_xlabel('Site Index', fontsize=12)
            axes[idx].set_ylabel('|ψ|²', fontsize=12)
            axes[idx].set_title(f'State {idx+1}: E = {E:.6f}', fontsize=12)
            axes[idx].grid(True, alpha=0.3)

            # Highlight edge localization
            if self.is_topological():
                axes[idx].axvspan(-0.5, 2.5, alpha=0.2, color='red',
                                  label='Left edge')
                axes[idx].axvspan(self.N-3.5, self.N-0.5, alpha=0.2, color='blue',
                                  label='Right edge')
                axes[idx].legend(fontsize=9)

        phase = "Topological" if self.is_topological() else "Trivial"
        fig.suptitle(f'Kitaev Chain Wavefunctions ({phase} Phase)',
                     fontsize=14, fontweight='bold')
        plt.tight_layout()
        return fig

# Example: Compare topological vs trivial phases
print("=== Kitaev Chain: Majorana Zero Modes ===\n")

# Topological phase: |μ| < 2t
N = 50
kc_topo = KitaevChain(N=N, t=1.0, mu=0.5, Delta=0.8)
E_topo, _ = kc_topo.solve()

# Trivial phase: |μ| > 2t
kc_trivial = KitaevChain(N=N, t=1.0, mu=3.0, Delta=0.8)
E_trivial, _ = kc_trivial.solve()

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

# Topological
positive_E_topo = E_topo[E_topo >= -1e-10]
axes[0].scatter(range(len(positive_E_topo)), positive_E_topo, s=50, alpha=0.7, c='blue')
axes[0].axhline(0, color='r', linestyle='--', linewidth=2)
axes[0].set_xlabel('State Index', fontsize=12)
axes[0].set_ylabel('Energy', fontsize=12)
axes[0].set_title('Topological Phase (μ = 0.5 < 2t)', fontsize=12, fontweight='bold')
axes[0].grid(True, alpha=0.3)
axes[0].set_ylim(-0.1, 2.5)

# Mark zero modes
zero_modes = np.abs(positive_E_topo) < 1e-6
if np.any(zero_modes):
    axes[0].scatter(np.where(zero_modes)[0], positive_E_topo[zero_modes],
                    s=200, c='red', marker='*', label='Majorana zero modes', zorder=5)
    axes[0].legend(fontsize=10)

# Trivial
positive_E_trivial = E_trivial[E_trivial >= -1e-10]
axes[1].scatter(range(len(positive_E_trivial)), positive_E_trivial, s=50, alpha=0.7, c='green')
axes[1].axhline(0, color='r', linestyle='--', linewidth=2)
axes[1].set_xlabel('State Index', fontsize=12)
axes[1].set_ylabel('Energy', fontsize=12)
axes[1].set_title('Trivial Phase (μ = 3.0 > 2t)', fontsize=12, fontweight='bold')
axes[1].grid(True, alpha=0.3)
axes[1].set_ylim(-0.1, 2.5)

plt.tight_layout()
plt.savefig('kitaev_chain_phases.png', dpi=150, bbox_inches='tight')
print("Saved: kitaev_chain_phases.png")

# Plot edge states for topological phase
fig_edge = kc_topo.plot_edge_states()
plt.savefig('kitaev_majorana_edge_states.png', dpi=150, bbox_inches='tight')
print("Saved: kitaev_majorana_edge_states.png")

print("\nKey observations:")
print(f"Topological phase (μ={kc_topo.mu}):")
print(f"  - Zero-energy modes present: {np.sum(np.abs(E_topo) < 1e-6)}")
print(f"  - These are Majorana zero modes localized at edges")
print(f"\nTrivial phase (μ={kc_trivial.mu}):")
print(f"  - No zero-energy modes: {np.sum(np.abs(E_trivial) < 1e-6)}")
print(f"  - All states have finite energy gap\n")

import numpy as np
import matplotlib.pyplot as plt

def kitaev_hamiltonian_k(k, t=1.0, mu=0.0, Delta=0.5):
    """
    Kitaev Hamiltonian in momentum space (Fourier transform).
    H(k) = [ε(k) - μ] σ_z + Δ sin(k) σ_y
    where ε(k) = -2t cos(k)
    """
    epsilon_k = -2 * t * np.cos(k)

    # Hamiltonian in Nambu basis using Pauli matrices
    H_k = np.array([
        [epsilon_k - mu, 2j * Delta * np.sin(k)],
        [-2j * Delta * np.sin(k), -(epsilon_k - mu)]
    ])

    return H_k

def compute_winding_number(t=1.0, mu=0.0, Delta=0.5, n_k=1000):
    """
    Compute the topological winding number.

    ν = (1/2πi) ∫ dk d/dk ln[det(H(k))]

    For Kitaev chain, this simplifies to counting how many times
    the vector (Re[h], Im[h]) winds around the origin, where
    h(k) is the off-diagonal element of H(k).
    """
    k_values = np.linspace(-np.pi, np.pi, n_k)

    # Vector components for winding calculation
    d_x = []
    d_y = []

    for k in k_values:
        epsilon_k = -2 * t * np.cos(k)
        d_x.append(epsilon_k - mu)
        d_y.append(-2 * Delta * np.sin(k))

    d_x = np.array(d_x)
    d_y = np.array(d_y)

    # Calculate winding using discrete approximation
    angles = np.arctan2(d_y, d_x)

    # Unwrap angles to get total winding
    angles_unwrapped = np.unwrap(angles)
    winding = (angles_unwrapped[-1] - angles_unwrapped[0]) / (2 * np.pi)

    return int(round(winding)), k_values, d_x, d_y

def plot_phase_diagram():
    """Plot topological phase diagram in (μ/t, Δ/t) space."""
    mu_range = np.linspace(-4, 4, 100)
    Delta_range = np.linspace(-2, 2, 100)

    MU, DELTA = np.meshgrid(mu_range, Delta_range)
    winding_map = np.zeros_like(MU)

    for i in range(len(Delta_range)):
        for j in range(len(mu_range)):
            winding, _, _, _ = compute_winding_number(
                t=1.0, mu=MU[i,j], Delta=DELTA[i,j], n_k=200
            )
            winding_map[i, j] = winding

    plt.figure(figsize=(10, 8))
    plt.contourf(MU, DELTA, winding_map, levels=[-1.5, -0.5, 0.5, 1.5],
                 colors=['lightblue', 'white', 'lightcoral'], alpha=0.7)
    plt.contour(MU, DELTA, winding_map, levels=[-0.5, 0.5],
                colors=['blue', 'red'], linewidths=2)

    # Phase boundaries: |μ| = 2t
    plt.axvline(-2, color='k', linestyle='--', linewidth=2, label='|μ| = 2t')
    plt.axvline(2, color='k', linestyle='--', linewidth=2)

    plt.xlabel('μ/t', fontsize=14)
    plt.ylabel('Δ/t', fontsize=14)
    plt.title('Topological Phase Diagram of Kitaev Chain', fontsize=16, fontweight='bold')

    # Add text labels
    plt.text(0, 1.5, 'Topological\n(ν = 1)', fontsize=14, ha='center',
             bbox=dict(boxstyle='round', facecolor='lightcoral', alpha=0.8))
    plt.text(3, 1.5, 'Trivial\n(ν = 0)', fontsize=14, ha='center',
             bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.8))

    plt.grid(True, alpha=0.3)
    plt.legend(fontsize=12)
    plt.tight_layout()
    return plt.gcf()

def plot_winding_trajectory(mu=0.5, Delta=0.8):
    """Plot the trajectory in (d_x, d_y) space that determines winding."""
    winding, k_values, d_x, d_y = compute_winding_number(mu=mu, Delta=Delta)

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

    # Trajectory plot
    axes[0].plot(d_x, d_y, 'b-', linewidth=2, alpha=0.7)
    axes[0].plot(d_x[0], d_y[0], 'go', markersize=12, label='k = -π')
    axes[0].plot(d_x[-1], d_y[-1], 'ro', markersize=12, label='k = π')
    axes[0].plot(0, 0, 'k*', markersize=15, label='Origin')
    axes[0].axhline(0, color='k', linestyle=':', alpha=0.5)
    axes[0].axvline(0, color='k', linestyle=':', alpha=0.5)
    axes[0].set_xlabel('d_x = ε(k) - μ', fontsize=12)
    axes[0].set_ylabel('d_y = -2Δsin(k)', fontsize=12)
    axes[0].set_title(f'Winding Trajectory (ν = {winding})', fontsize=14, fontweight='bold')
    axes[0].grid(True, alpha=0.3)
    axes[0].legend(fontsize=10)
    axes[0].axis('equal')

    # Arrow to show winding direction
    mid_idx = len(d_x) // 2
    axes[0].arrow(d_x[mid_idx], d_y[mid_idx],
                  d_x[mid_idx+10]-d_x[mid_idx], d_y[mid_idx+10]-d_y[mid_idx],
                  head_width=0.2, head_length=0.1, fc='blue', ec='blue')

    # Energy dispersion
    energies = []
    for k in k_values:
        H_k = kitaev_hamiltonian_k(k, mu=mu, Delta=Delta)
        E = np.linalg.eigvalsh(H_k)
        energies.append(E)
    energies = np.array(energies)

    axes[1].plot(k_values, energies[:, 0], 'b-', linewidth=2, label='E_-')
    axes[1].plot(k_values, energies[:, 1], 'r-', linewidth=2, label='E_+')
    axes[1].axhline(0, color='k', linestyle='--', alpha=0.5)
    axes[1].set_xlabel('k', fontsize=12)
    axes[1].set_ylabel('Energy', fontsize=12)
    axes[1].set_title('Band Structure', fontsize=14, fontweight='bold')
    axes[1].set_xticks([-np.pi, -np.pi/2, 0, np.pi/2, np.pi])
    axes[1].set_xticklabels(['-π', '-π/2', '0', 'π/2', 'π'])
    axes[1].grid(True, alpha=0.3)
    axes[1].legend(fontsize=10)

    plt.tight_layout()
    return fig

# Example usage
print("=== Topological Invariant: Winding Number ===\n")

# Topological case
print("Case 1: Topological phase (μ = 0.5, Δ = 0.8)")
winding_topo, _, _, _ = compute_winding_number(mu=0.5, Delta=0.8)
print(f"  Winding number: ν = {winding_topo}")
print(f"  Majorana modes expected: YES (ν ≠ 0)\n")

# Trivial case
print("Case 2: Trivial phase (μ = 3.0, Δ = 0.8)")
winding_trivial, _, _, _ = compute_winding_number(mu=3.0, Delta=0.8)
print(f"  Winding number: ν = {winding_trivial}")
print(f"  Majorana modes expected: NO (ν = 0)\n")

# Plot phase diagram
fig_pd = plot_phase_diagram()
plt.savefig('kitaev_phase_diagram.png', dpi=150, bbox_inches='tight')
print("Saved: kitaev_phase_diagram.png")

# Plot winding trajectory for topological case
fig_wind = plot_winding_trajectory(mu=0.5, Delta=0.8)
plt.savefig('kitaev_winding_trajectory.png', dpi=150, bbox_inches='tight')
print("Saved: kitaev_winding_trajectory.png")

print("\nPhysical interpretation:")
print("  - Winding number ν counts how many times the Hamiltonian")
print("    vector winds around the origin as k goes from -π to π")
print("  - ν = 0: Trivial phase (no edge modes)")
print("  - ν = 1: Topological phase (Majorana zero modes at edges)")
print("  - Phase transition occurs when gap closes (vector passes through origin)")
print("  - Critical points: μ = ±2t (when Δ ≠ 0)\n")

import numpy as np
import matplotlib.pyplot as plt

def mcmillan_tc(omega_D, lambda_ep, mu_star=0.13):
    """
    Calculate Tc using the McMillan equation (modified for strong coupling).

    T_c = (ω_D / 1.2) * exp[−1.04(1 + λ) / (λ − μ*(1 + 0.62λ))]

    Parameters:
    -----------
    omega_D : float
        Debye frequency (in K)
    lambda_ep : float
        Electron-phonon coupling constant
    mu_star : float
        Coulomb pseudopotential

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

    # Avoid division by zero or negative denominator
    if denominator <= 0:
        return 0.0

    Tc = (omega_D / 1.2) * np.exp(numerator / denominator)
    return Tc

def allen_dynes_tc(omega_log, lambda_ep, mu_star=0.13):
    """
    Allen-Dynes equation (more accurate for strong coupling).

    T_c = (ω_log / 1.2) * exp[−1.04(1 + λ) / (λ − μ*(1 + 0.62λ))] * f_1 * f_2

    where f_1 and f_2 are correction factors for strong coupling.
    """
    f_1 = (1 + (lambda_ep / 2.46) ** (3/2)) ** (1/3)
    f_2 = 1 + ((omega_log / 1000) - 1) * lambda_ep**2 / (lambda_ep**2 + 1)

    numerator = -1.04 * (1 + lambda_ep)
    denominator = lambda_ep - mu_star * (1 + 0.62 * lambda_ep)

    if denominator <= 0:
        return 0.0

    Tc = (omega_log / 1.2) * np.exp(numerator / denominator) * f_1 * f_2
    return Tc

# Known hydride superconductors
hydrides = {
    'H3S': {'omega_D': 1500, 'lambda': 2.0, 'Tc_exp': 203, 'pressure': 150},
    'LaH10': {'omega_D': 1800, 'lambda': 2.5, 'Tc_exp': 260, 'pressure': 190},
    'YH9': {'omega_D': 1700, 'lambda': 2.3, 'Tc_exp': 243, 'pressure': 201},
    'CaH6': {'omega_D': 1600, 'lambda': 2.1, 'Tc_exp': 215, 'pressure': 172},
}

print("=== Tc Estimation for High-Pressure Hydrides ===\n")
print(f"{'Material':<10} {'ω_D (K)':<12} {'λ':<8} {'Tc_calc (K)':<15} {'Tc_exp (K)':<15} {'Error (%)':<10}")
print("-" * 80)

results = {}
for material, params in hydrides.items():
    Tc_calc = mcmillan_tc(params['omega_D'], params['lambda'])
    Tc_exp = params['Tc_exp']
    error = abs(Tc_calc - Tc_exp) / Tc_exp * 100

    results[material] = {'Tc_calc': Tc_calc, 'Tc_exp': Tc_exp}

    print(f"{material:<10} {params['omega_D']:<12} {params['lambda']:<8.1f} "
          f"{Tc_calc:<15.1f} {Tc_exp:<15} {error:<10.1f}")

print("\n" + "="*80)

# Plot 1: Tc vs λ for different ω_D
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

lambda_range = np.linspace(0.5, 3.5, 100)
omega_values = [800, 1200, 1600, 2000]

for omega_D in omega_values:
    Tc_values = [mcmillan_tc(omega_D, lam) for lam in lambda_range]
    axes[0].plot(lambda_range, Tc_values, linewidth=2, label=f'ω_D = {omega_D} K')

# Mark experimental points
for material, params in hydrides.items():
    axes[0].plot(params['lambda'], params['Tc_exp'], 'o', markersize=10,
                 label=f'{material} (exp)')

axes[0].set_xlabel('Electron-Phonon Coupling λ', fontsize=12)
axes[0].set_ylabel('Tc (K)', fontsize=12)
axes[0].set_title('Tc vs Coupling Strength', fontsize=14, fontweight='bold')
axes[0].grid(True, alpha=0.3)
axes[0].legend(fontsize=9, ncol=2)
axes[0].set_ylim(0, 350)
axes[0].axhline(300, color='r', linestyle='--', linewidth=2, alpha=0.5,
                label='Room temperature')

# Plot 2: Tc vs ω_D for different λ
omega_range = np.linspace(500, 2500, 100)
lambda_values = [1.5, 2.0, 2.5, 3.0]

for lam in lambda_values:
    Tc_values = [mcmillan_tc(omega, lam) for omega in omega_range]
    axes[1].plot(omega_range, Tc_values, linewidth=2, label=f'λ = {lam}')

# Mark experimental points
for material, params in hydrides.items():
    axes[1].plot(params['omega_D'], params['Tc_exp'], 's', markersize=10,
                 label=f'{material} (exp)')

axes[1].set_xlabel('Debye Frequency ω_D (K)', fontsize=12)
axes[1].set_ylabel('Tc (K)', fontsize=12)
axes[1].set_title('Tc vs Phonon Frequency', fontsize=14, fontweight='bold')
axes[1].grid(True, alpha=0.3)
axes[1].legend(fontsize=9, ncol=2)
axes[1].set_ylim(0, 350)
axes[1].axhline(300, color='r', linestyle='--', linewidth=2, alpha=0.5)

plt.tight_layout()
plt.savefig('hydride_tc_predictions.png', dpi=150, bbox_inches='tight')
print("\nSaved: hydride_tc_predictions.png")

# Explore room-temperature parameter space
print("\n=== Pathways to Room-Temperature SC (Tc ≥ 300 K) ===\n")
print("Required parameter combinations:\n")

for omega_D in [1800, 2000, 2200, 2500]:
    for lam in [2.5, 3.0, 3.5, 4.0]:
        Tc = mcmillan_tc(omega_D, lam)
        if 295 <= Tc <= 350:
            print(f"  ω_D = {omega_D} K, λ = {lam:.1f} → Tc = {Tc:.1f} K ✓")

print("\nKey insights:")
print("  1. Need BOTH high ω_D (light atoms) AND strong coupling (λ > 2)")
print("  2. Hydrogen-rich compounds are ideal (low mass → high ω_D)")
print("  3. Challenge: Achieving these parameters at ambient pressure")
print("  4. Current record (LaH10): Tc ~ 260 K, still 40 K short of room T\n")

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import pandas as pd

# Historical superconductivity milestones
milestones = {
    'Year': [1911, 1933, 1950, 1957, 1962, 1973, 1986, 1987, 1993, 2001,
             2008, 2015, 2019, 2020, 2023],
    'Tc': [4.2, 4.2, 16, 16, 8.7, 23, 35, 92, 133, 39, 55, 203, 260, 287, 1],
    'Material': ['Hg', 'Meissner\neffect', 'Nb', 'BCS\ntheory', 'Nb₃Sn',
                 'Nb₃Ge', 'LaBaCuO', 'YBCO', 'HgBaCaCuO', 'MgB₂',
                 'LaFeAsO', 'H₃S', 'LaH₁₀', 'CSH\n(retracted)', 'LK-99\n(disputed)'],
    'Category': ['Elemental', 'Theory', 'Elemental', 'Theory', 'A15',
                 'A15', 'Cuprate', 'Cuprate', 'Cuprate', 'MgB2',
                 'Fe-based', 'Hydride', 'Hydride', 'Hydride', 'Disputed'],
    'Pressure': [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 150, 190, 267, 0]
}

df = pd.DataFrame(milestones)

# Create comprehensive figure
fig = plt.figure(figsize=(16, 10))
gs = fig.add_gridspec(3, 2, hspace=0.3, wspace=0.3)

# Plot 1: Historical Tc progression
ax1 = fig.add_subplot(gs[0, :])

# Color by category
colors = {
    'Elemental': '#3498db', 'Theory': '#95a5a6', 'A15': '#9b59b6',
    'Cuprate': '#e74c3c', 'MgB2': '#f39c12', 'Fe-based': '#16a085',
    'Hydride': '#d35400', 'Disputed': '#bdc3c7'
}

# Plot verified superconductors
verified = df[~df['Category'].isin(['Disputed', 'Theory'])]
for category in verified['Category'].unique():
    data = verified[verified['Category'] == category]
    ax1.scatter(data['Year'], data['Tc'], s=100, c=colors[category],
                label=category, zorder=3, alpha=0.8, edgecolors='black', linewidth=1.5)

# Plot disputed/retracted with X marker
disputed = df[df['Category'] == 'Disputed']
ax1.scatter(disputed['Year'], disputed['Tc'], s=200, c=colors['Disputed'],
            marker='x', label='Disputed/Retracted', zorder=4, linewidth=3)

# Annotate key milestones
for idx, row in df.iterrows():
    if row['Category'] != 'Theory':
        offset = 15 if idx % 2 == 0 else -25
        ax1.annotate(row['Material'], xy=(row['Year'], row['Tc']),
                     xytext=(0, offset), textcoords='offset points',
                     fontsize=8, ha='center',
                     bbox=dict(boxstyle='round,pad=0.3', facecolor='white',
                              edgecolor='gray', alpha=0.7),
                     arrowprops=dict(arrowstyle='->', lw=1, color='gray'))

# Mark significant temperature thresholds
ax1.axhline(77, color='blue', linestyle='--', linewidth=1.5, alpha=0.5,
            label='Liquid N₂ (77 K)')
ax1.axhline(300, color='red', linestyle='--', linewidth=1.5, alpha=0.5,
            label='Room temperature (300 K)')

ax1.set_xlabel('Year', fontsize=14)
ax1.set_ylabel('Critical Temperature Tc (K)', fontsize=14)
ax1.set_title('Historical Progress of Superconductivity', fontsize=16, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.legend(loc='upper left', fontsize=9, ncol=2)
ax1.set_xlim(1905, 2025)
ax1.set_ylim(0, 320)

# Plot 2: Tc vs Pressure for hydrides
ax2 = fig.add_subplot(gs[1, 0])

hydride_data = df[df['Category'] == 'Hydride']
ax2.scatter(hydride_data['Pressure'], hydride_data['Tc'], s=150, c='#d35400',
            marker='D', edgecolors='black', linewidth=1.5)

for idx, row in hydride_data.iterrows():
    if row['Material'] != 'CSH\n(retracted)':
        ax2.annotate(row['Material'].strip(), xy=(row['Pressure'], row['Tc']),
                     xytext=(10, 5), textcoords='offset points', fontsize=10)

ax2.set_xlabel('Pressure (GPa)', fontsize=12)
ax2.set_ylabel('Tc (K)', fontsize=12)
ax2.set_title('High-Pressure Hydride Superconductors', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.axhline(300, color='r', linestyle='--', alpha=0.5, label='Room T')
ax2.legend(fontsize=10)

# Plot 3: Research activity by year (simulated)
ax3 = fig.add_subplot(gs[1, 1])

years = np.arange(2000, 2025)
# Simulate publication counts (peaks around major discoveries)
base_activity = 50 + 5 * (years - 2000)
activity = base_activity.copy()
activity[years == 2008] += 100  # Fe-based discovery
activity[years == 2015] += 80   # H3S
activity[years == 2018] += 120  # TBG
activity[years == 2019] += 100  # LaH10
activity[years == 2020] += 90   # CSH claim
activity[years == 2023] += 200  # LK-99 hype

ax3.bar(years, activity, color='#3498db', alpha=0.7, edgecolor='black')
ax3.set_xlabel('Year', fontsize=12)
ax3.set_ylabel('Research Publications (simulated)', fontsize=12)
ax3.set_title('SC Research Activity Over Time', fontsize=14, fontweight='bold')
ax3.grid(True, alpha=0.3, axis='y')

# Annotate major events
events = {2008: 'Fe-based', 2015: 'H₃S', 2018: 'TBG', 2019: 'LaH₁₀', 2023: 'LK-99'}
for year, label in events.items():
    idx = year - 2000
    ax3.annotate(label, xy=(year, activity[idx]), xytext=(0, 5),
                 textcoords='offset points', fontsize=9, ha='center')

# Plot 4: Current research frontiers (pie chart)
ax4 = fig.add_subplot(gs[2, :])

frontiers = ['High-Pressure\nHydrides', 'Twisted 2D\nMaterials',
             'Topological SC\n& Majoranas', 'Cuprate\nMechanisms',
             'Fe-based SC', 'Other\n(MgB₂-like, etc.)']
activity_share = [25, 20, 18, 15, 12, 10]
colors_pie = ['#d35400', '#9b59b6', '#16a085', '#e74c3c', '#f39c12', '#95a5a6']

wedges, texts, autotexts = ax4.pie(activity_share, labels=frontiers, autopct='%1.0f%%',
                                    colors=colors_pie, startangle=90,
                                    textprops={'fontsize': 11})

# Make percentage text bold
for autotext in autotexts:
    autotext.set_color('white')
    autotext.set_fontweight('bold')
    autotext.set_fontsize(12)

ax4.set_title('Current Research Frontiers (2024 estimate)', fontsize=14, fontweight='bold')

plt.savefig('superconductivity_research_trends.png', dpi=150, bbox_inches='tight')
print("=== Superconductivity Research Trends ===\n")
print("Saved: superconductivity_research_trends.png")

# Print summary statistics
print("\nKey Statistics:")
print(f"  Total years of SC research: {2024 - 1911} years")
print(f"  Verified Tc record (ambient pressure): 133 K (HgBaCaCuO)")
print(f"  Verified Tc record (high pressure): 260 K (LaH₁₀)")
print(f"  Gap to room temperature: {300 - 260} K")
print(f"  Number of major SC families: 6+")
print(f"\nGrowth rate:")
initial_tc = df.iloc[0]['Tc']
current_tc = 260  # LaH10
years_elapsed = 2019 - 1911
avg_growth = (current_tc - initial_tc) / years_elapsed
print(f"  Average Tc increase: {avg_growth:.2f} K/year")
print(f"  Time to room T (linear extrapolation): {(300-260)/avg_growth:.0f} years from 2019")
print(f"\n  Note: Progress is NOT linear - breakthroughs are unpredictable!\n")