Learning Objectives
- Implement BCS gap equation simulations
- Visualize magnetic field penetration using London equations
- Create interactive superconductor phase diagrams
- Build a superconductor materials database and analyzer
- Simulate vortex behavior in Type II superconductors
5.1 Setup and Prerequisites
Required Libraries
# Install required packages
# pip install numpy matplotlib scipy pandas
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize
import pandas as pd
# Set plotting style
plt.style.use('seaborn-v0_8-whitegrid')
plt.rcParams['figure.figsize'] = [10, 6]
plt.rcParams['font.size'] = 12
5.2 BCS Gap Equation Simulation
The BCS theory predicts how the superconducting energy gap varies with temperature. Let's implement this numerically.
Theory
The BCS gap equation at finite temperature is:
1 = N(0)V ∫₀^ωD dε / √(ε² + Δ²) × tanh(√(ε² + Δ²) / 2kBT)
We can solve this self-consistently to find Δ(T).
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize
def bcs_gap_equation(Delta, T, Tc, Delta0):
"""
BCS gap equation for numerical solution.
Returns the difference that should be zero at the solution.
"""
if Delta <= 0 or T >= Tc:
return Delta # Return Delta itself to indicate no solution
# Reduced temperature
t = T / Tc
# Use the approximate BCS formula for the integral
# This is more stable than direct numerical integration
kB = 1.0 # Using natural units where kB*Tc = 1
# Number of integration points
N = 1000
omega_D = 10 * Delta0 # Debye cutoff (>> Delta)
epsilon = np.linspace(0.001, omega_D, N)
E = np.sqrt(epsilon**2 + Delta**2)
# BCS kernel
if T > 0:
kernel = np.tanh(E / (2 * T)) / E
else:
kernel = 1.0 / E
# Numerical integration
integral = np.trapz(kernel, epsilon)
# Gap equation: should equal the coupling strength
# We normalize so that at T=0, Delta = Delta0
coupling = Delta0 / np.trapz(1.0 / np.sqrt(epsilon**2 + Delta0**2), epsilon)
return coupling * integral - 1.0
def solve_bcs_gap(Tc, Delta0, num_points=50):
"""
Solve BCS gap equation for a range of temperatures.
Returns T and Delta arrays.
"""
T_values = np.linspace(0.01 * Tc, 0.999 * Tc, num_points)
Delta_values = []
for T in T_values:
# Initial guess: use approximate formula
t = T / Tc
Delta_approx = Delta0 * np.sqrt(1 - t) if t < 1 else 0
try:
# Find the root of the gap equation
result = optimize.brentq(
lambda D: bcs_gap_equation(D, T, Tc, Delta0),
1e-10, Delta0 * 1.5
)
Delta_values.append(result)
except:
# Use approximate formula if numerical solution fails
Delta_values.append(Delta_approx)
return np.array(T_values), np.array(Delta_values)
# Simulation parameters
Tc = 9.3 # K (Niobium)
kB = 8.617e-5 # eV/K
Delta0 = 1.76 * kB * Tc # BCS prediction for Delta(0)
# Use simpler approximate BCS formula for stability
def bcs_gap_approx(T, Tc, Delta0):
"""Approximate BCS gap temperature dependence."""
t = T / Tc
if t >= 1:
return 0
# Approximation valid over full temperature range
return Delta0 * np.tanh(1.74 * np.sqrt((1/t - 1)))
# Generate data
T_range = np.linspace(0.1, Tc, 200)
Delta_approx = [bcs_gap_approx(T, Tc, Delta0) for T in T_range]
# Plotting
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# Left: Gap vs Temperature
ax1 = axes[0]
ax1.plot(T_range, np.array(Delta_approx) * 1000, 'b-', linewidth=2,
label='BCS prediction')
ax1.axvline(x=Tc, color='r', linestyle='--', label=f'Tc = {Tc} K')
ax1.axhline(y=Delta0 * 1000, color='gray', linestyle=':', alpha=0.5,
label=f'Δ(0) = {Delta0*1000:.3f} meV')
ax1.set_xlabel('Temperature (K)', fontsize=12)
ax1.set_ylabel('Energy Gap Δ (meV)', fontsize=12)
ax1.set_title('BCS Energy Gap vs Temperature (Niobium)', fontsize=14)
ax1.legend(fontsize=10)
ax1.set_xlim(0, Tc * 1.1)
ax1.set_ylim(0, Delta0 * 1000 * 1.2)
ax1.grid(True, alpha=0.3)
# Right: Normalized gap
ax2 = axes[1]
t_norm = T_range / Tc
delta_norm = np.array(Delta_approx) / Delta0
ax2.plot(t_norm, delta_norm, 'b-', linewidth=2, label='BCS approximation')
# Add weak-coupling BCS limit
ax2.plot([0, 1], [1, 0], 'g--', alpha=0.5, linewidth=1.5,
label='Linear approximation')
ax2.set_xlabel('T / Tc', fontsize=12)
ax2.set_ylabel('Δ(T) / Δ(0)', fontsize=12)
ax2.set_title('Normalized BCS Gap', fontsize=14)
ax2.legend(fontsize=10)
ax2.set_xlim(0, 1.1)
ax2.set_ylim(0, 1.1)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Print key values
print("BCS Theory Predictions for Niobium:")
print(f" Critical temperature Tc = {Tc} K")
print(f" Zero-temperature gap Δ(0) = {Delta0*1000:.3f} meV")
print(f" BCS ratio 2Δ(0)/kB*Tc = {2*Delta0/(kB*Tc):.3f} (BCS: 3.52)")
5.3 London Penetration Depth Simulation
The London equations describe how magnetic fields penetrate into superconductors. Let's simulate the field decay at the surface.
import numpy as np
import matplotlib.pyplot as plt
class LondonPenetration:
"""
Simulate magnetic field penetration using London theory.
"""
def __init__(self, lambda_L, thickness=None):
"""
Initialize with London penetration depth.
Parameters:
-----------
lambda_L : float
London penetration depth in nm
thickness : float, optional
Film thickness in nm (None for bulk)
"""
self.lambda_L = lambda_L
self.thickness = thickness
def field_profile_bulk(self, x, B0=1.0):
"""
Magnetic field profile in bulk superconductor.
B(x) = B0 * exp(-x/lambda_L)
"""
return B0 * np.exp(-x / self.lambda_L)
def field_profile_film(self, x, B0=1.0):
"""
Magnetic field profile in thin film.
For a film of thickness d with field applied on both sides.
"""
if self.thickness is None:
return self.field_profile_bulk(x, B0)
d = self.thickness
lam = self.lambda_L
# Solution for thin film: cosh profile
return B0 * np.cosh((x - d/2) / lam) / np.cosh(d / (2*lam))
def current_density(self, x, B0=1.0, mu0=4*np.pi*1e-7):
"""
Screening current density from Maxwell equation.
J = -dB/dx / mu0
"""
return B0 / (self.lambda_L * mu0) * np.exp(-x / self.lambda_L)
# Create instances for different materials
materials = {
'Aluminum': LondonPenetration(50),
'Niobium': LondonPenetration(40),
'YBCO': LondonPenetration(150),
'NbTi': LondonPenetration(300),
}
# Spatial range
x = np.linspace(0, 1000, 1000) # nm
# Create comprehensive visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Plot 1: Field profiles comparison
ax1 = axes[0, 0]
colors = plt.cm.viridis(np.linspace(0.2, 0.8, len(materials)))
for (name, sc), color in zip(materials.items(), colors):
B = sc.field_profile_bulk(x)
ax1.plot(x, B, linewidth=2, color=color,
label=f'{name} (λ={sc.lambda_L} nm)')
ax1.axhline(y=1/np.e, color='gray', linestyle='--', alpha=0.5)
ax1.text(800, 1/np.e + 0.05, f'B/B₀ = 1/e', fontsize=10, color='gray')
ax1.set_xlabel('Distance from surface (nm)', fontsize=12)
ax1.set_ylabel('B/B₀', fontsize=12)
ax1.set_title('Magnetic Field Penetration in Bulk Superconductors', fontsize=14)
ax1.legend(fontsize=10)
ax1.set_xlim(0, 1000)
ax1.set_ylim(0, 1.1)
ax1.grid(True, alpha=0.3)
# Plot 2: Current density profiles
ax2 = axes[0, 1]
for (name, sc), color in zip(materials.items(), colors):
# Normalized current density
J = np.exp(-x / sc.lambda_L) / sc.lambda_L
ax2.plot(x, J * 1000, linewidth=2, color=color, label=name)
ax2.set_xlabel('Distance from surface (nm)', fontsize=12)
ax2.set_ylabel('Screening current (arb. units)', fontsize=12)
ax2.set_title('Screening Current Density', fontsize=14)
ax2.legend(fontsize=10)
ax2.set_xlim(0, 500)
ax2.grid(True, alpha=0.3)
# Plot 3: Thin film effects
ax3 = axes[1, 0]
thicknesses = [50, 100, 200, 500] # nm
colors_film = plt.cm.plasma(np.linspace(0.2, 0.8, len(thicknesses)))
lambda_L = 100 # nm
for d, color in zip(thicknesses, colors_film):
film = LondonPenetration(lambda_L, thickness=d)
x_film = np.linspace(0, d, 200)
B = film.field_profile_film(x_film)
ax3.plot(x_film, B, linewidth=2, color=color, label=f'd = {d} nm')
ax3.set_xlabel('Position in film (nm)', fontsize=12)
ax3.set_ylabel('B/B₀', fontsize=12)
ax3.set_title(f'Field Profile in Thin Films (λ = {lambda_L} nm)', fontsize=14)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)
# Plot 4: Temperature dependence of λ
ax4 = axes[1, 1]
def lambda_vs_T(T, Tc, lambda_0):
"""Two-fluid model: λ(T) = λ(0) / sqrt(1 - (T/Tc)^4)"""
t = T / Tc
if t >= 1:
return np.inf
return lambda_0 / np.sqrt(1 - t**4)
# Niobium example
Tc = 9.3
lambda_0 = 40
T_range = np.linspace(0.1, 0.99 * Tc, 100)
lambda_T = [lambda_vs_T(T, Tc, lambda_0) for T in T_range]
ax4.plot(T_range / Tc, np.array(lambda_T) / lambda_0, 'b-', linewidth=2)
ax4.axhline(y=1, color='gray', linestyle=':', alpha=0.5)
ax4.axvline(x=1, color='r', linestyle='--', alpha=0.5, label='Tc')
ax4.set_xlabel('T / Tc', fontsize=12)
ax4.set_ylabel('λ(T) / λ(0)', fontsize=12)
ax4.set_title('Temperature Dependence of Penetration Depth', fontsize=14)
ax4.set_xlim(0, 1.05)
ax4.set_ylim(0.9, 5)
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
5.4 Superconductor Materials Database
Let's build a comprehensive database of superconducting materials and analyze their properties.
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
# Create superconductor database
sc_data = {
'Material': ['Al', 'Hg', 'Pb', 'Nb', 'NbTi', 'Nb3Sn', 'MgB2',
'YBCO', 'BSCCO-2223', 'Hg-1223', 'FeSe', 'H3S'],
'Tc_K': [1.2, 4.2, 7.2, 9.3, 10, 18.3, 39,
93, 110, 133, 8, 203],
'Hc2_T': [0.01, 0.04, 0.08, 0.4, 15, 24, 16,
100, 60, 100, 17, 70],
'Type': ['I', 'I', 'I', 'II', 'II', 'II', 'II',
'II', 'II', 'II', 'II', 'II'],
'Class': ['Element', 'Element', 'Element', 'Element', 'Alloy', 'A15', 'Compound',
'Cuprate', 'Cuprate', 'Cuprate', 'Iron-based', 'Hydride'],
'Year': [1933, 1911, 1913, 1930, 1962, 1954, 2001,
1987, 1988, 1993, 2008, 2015],
'Practical': [False, False, False, True, True, True, True,
True, True, False, False, False],
'Lambda_nm': [50, 40, 39, 40, 300, 80, 140,
150, 200, 200, 400, 50],
'Xi_nm': [1600, 170, 83, 38, 4, 3, 5,
1.5, 2, 1.5, 2, 1]
}
df = pd.DataFrame(sc_data)
# Calculate Ginzburg-Landau parameter
df['Kappa'] = df['Lambda_nm'] / df['Xi_nm']
# Display database
print("Superconductor Database:")
print(df.to_string(index=False))
print()
# Analysis and visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Plot 1: Tc by material class
ax1 = axes[0, 0]
class_colors = {'Element': 'blue', 'Alloy': 'green', 'A15': 'orange',
'Compound': 'purple', 'Cuprate': 'red', 'Iron-based': 'brown',
'Hydride': 'magenta'}
for mat_class in df['Class'].unique():
subset = df[df['Class'] == mat_class]
ax1.scatter(subset['Year'], subset['Tc_K'], s=100,
c=class_colors[mat_class], label=mat_class, alpha=0.7)
ax1.axhline(y=77, color='cyan', linestyle='--', alpha=0.7, linewidth=2)
ax1.text(2018, 80, 'LN₂ (77K)', fontsize=10, color='cyan')
ax1.set_xlabel('Discovery Year', fontsize=12)
ax1.set_ylabel('Critical Temperature (K)', fontsize=12)
ax1.set_title('Evolution of Superconductor Tc', fontsize=14)
ax1.legend(fontsize=9, loc='upper left')
ax1.grid(True, alpha=0.3)
# Plot 2: Tc vs Hc2
ax2 = axes[0, 1]
for mat_class in df['Class'].unique():
subset = df[df['Class'] == mat_class]
ax2.scatter(subset['Tc_K'], subset['Hc2_T'], s=100,
c=class_colors[mat_class], label=mat_class, alpha=0.7)
ax2.set_xlabel('Critical Temperature Tc (K)', fontsize=12)
ax2.set_ylabel('Upper Critical Field Hc2 (T)', fontsize=12)
ax2.set_title('Tc vs Hc2', fontsize=14)
ax2.set_xscale('log')
ax2.set_yscale('log')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3, which='both')
# Plot 3: Kappa values (Type I vs II boundary)
ax3 = axes[1, 0]
colors_kappa = ['blue' if k < 0.71 else 'red' for k in df['Kappa']]
bars = ax3.bar(df['Material'], df['Kappa'], color=colors_kappa, alpha=0.7)
ax3.axhline(y=1/np.sqrt(2), color='green', linestyle='--', linewidth=2,
label=f'κ = 1/√2 ≈ 0.71')
ax3.set_ylabel('κ = λ/ξ', fontsize=12)
ax3.set_title('Ginzburg-Landau Parameter κ', fontsize=14)
ax3.set_yscale('log')
ax3.legend(fontsize=10)
plt.setp(ax3.xaxis.get_majorticklabels(), rotation=45, ha='right')
# Plot 4: Practical vs Research materials
ax4 = axes[1, 1]
practical = df[df['Practical'] == True]
research = df[df['Practical'] == False]
ax4.scatter(practical['Tc_K'], practical['Hc2_T'], s=150, c='green',
marker='s', label=f'Practical ({len(practical)})', alpha=0.7)
ax4.scatter(research['Tc_K'], research['Hc2_T'], s=150, c='red',
marker='o', label=f'Research ({len(research)})', alpha=0.7)
# Add material names
for _, row in df.iterrows():
ax4.annotate(row['Material'], (row['Tc_K'], row['Hc2_T']),
xytext=(5, 5), textcoords='offset points', fontsize=8)
ax4.set_xlabel('Critical Temperature Tc (K)', fontsize=12)
ax4.set_ylabel('Upper Critical Field Hc2 (T)', fontsize=12)
ax4.set_title('Practical vs Research Materials', fontsize=14)
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Statistical summary
print("\nStatistical Summary by Class:")
print(df.groupby('Class')[['Tc_K', 'Hc2_T', 'Kappa']].mean().round(2))
5.5 Vortex State Simulation
In Type II superconductors, magnetic flux penetrates in quantized vortices. Let's simulate a simple vortex lattice.
import numpy as np
import matplotlib.pyplot as plt
class VortexSimulation:
"""
Simulate vortex states in Type II superconductors.
"""
def __init__(self, lambda_L=100, xi=5, Phi_0=2.07e-15):
"""
Initialize vortex simulation.
Parameters:
-----------
lambda_L : float
Penetration depth (nm)
xi : float
Coherence length (nm)
Phi_0 : float
Flux quantum (Wb)
"""
self.lambda_L = lambda_L
self.xi = xi
self.Phi_0 = Phi_0
self.kappa = lambda_L / xi
def single_vortex_field(self, x, y, x0=0, y0=0):
"""
Magnetic field of a single vortex.
Approximate solution using modified Bessel function.
"""
r = np.sqrt((x - x0)**2 + (y - y0)**2)
# Avoid singularity at center
r = np.maximum(r, self.xi)
# Field profile (approximate)
B = np.exp(-r / self.lambda_L) / (2 * np.pi * self.lambda_L**2)
return B
def vortex_lattice(self, x, y, B_avg, lattice_type='triangular'):
"""
Create a vortex lattice for given average field.
"""
# Vortex spacing from average field
# B_avg = Phi_0 / (area per vortex)
# For triangular: a = sqrt(2*Phi_0 / (sqrt(3) * B_avg))
a = np.sqrt(2 * self.Phi_0 / (np.sqrt(3) * B_avg)) * 1e9 # Convert to nm
# Generate vortex positions
vortex_positions = []
n_range = 5 # Number of vortices in each direction
for i in range(-n_range, n_range + 1):
for j in range(-n_range, n_range + 1):
if lattice_type == 'triangular':
x0 = a * (i + 0.5 * (j % 2))
y0 = a * j * np.sqrt(3) / 2
else: # square
x0 = a * i
y0 = a * j
vortex_positions.append((x0, y0))
# Calculate total field
B_total = np.zeros_like(x)
for x0, y0 in vortex_positions:
B_total += self.single_vortex_field(x, y, x0, y0)
return B_total, vortex_positions, a
# Create simulation
sim = VortexSimulation(lambda_L=150, xi=2) # YBCO-like parameters
# Create spatial grid
L = 500 # nm
N = 200
x = np.linspace(-L, L, N)
y = np.linspace(-L, L, N)
X, Y = np.meshgrid(x, y)
# Different applied fields (in Tesla)
fields = [0.5, 1.0, 2.0]
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for ax, B_applied in zip(axes, fields):
B, vortex_pos, spacing = sim.vortex_lattice(X, Y, B_applied)
# Plot field intensity
im = ax.contourf(X, Y, B, levels=50, cmap='hot')
plt.colorbar(im, ax=ax, label='B (arb. units)')
# Mark vortex centers
vx = [p[0] for p in vortex_pos if abs(p[0]) < L and abs(p[1]) < L]
vy = [p[1] for p in vortex_pos if abs(p[0]) < L and abs(p[1]) < L]
ax.scatter(vx, vy, c='cyan', s=20, marker='o', alpha=0.5)
ax.set_xlabel('x (nm)', fontsize=11)
ax.set_ylabel('y (nm)', fontsize=11)
ax.set_title(f'B = {B_applied} T\nVortex spacing ≈ {spacing:.0f} nm', fontsize=12)
ax.set_aspect('equal')
plt.suptitle('Vortex Lattice in Type II Superconductor', fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
# Additional visualization: Order parameter near vortex core
fig, ax = plt.subplots(figsize=(10, 6))
r = np.linspace(0, 100, 200) # nm
# Order parameter recovery
psi = np.tanh(r / (np.sqrt(2) * sim.xi))
# Magnetic field decay
B_normalized = np.exp(-r / sim.lambda_L)
ax.plot(r, psi**2, 'b-', linewidth=2, label='|ψ|² (order parameter)')
ax.plot(r, B_normalized, 'r-', linewidth=2, label='B/B_core (magnetic field)')
ax.axvline(x=sim.xi, color='blue', linestyle='--', alpha=0.5,
label=f'ξ = {sim.xi} nm')
ax.axvline(x=sim.lambda_L, color='red', linestyle='--', alpha=0.5,
label=f'λ = {sim.lambda_L} nm')
ax.set_xlabel('Distance from vortex center (nm)', fontsize=12)
ax.set_ylabel('Normalized value', fontsize=12)
ax.set_title('Structure of a Single Vortex', fontsize=14)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 100)
plt.tight_layout()
plt.show()
print(f"Vortex parameters:")
print(f" Penetration depth λ = {sim.lambda_L} nm")
print(f" Coherence length ξ = {sim.xi} nm")
print(f" GL parameter κ = {sim.kappa:.1f}")
print(f" Flux quantum Φ₀ = {sim.Phi_0:.2e} Wb")
5.6 Complete Simulation Project
Let's combine everything into a comprehensive superconductor analysis tool.
import numpy as np
import matplotlib.pyplot as plt
class SuperconductorAnalyzer:
"""
Complete analyzer for superconductor properties.
"""
def __init__(self, name, Tc, lambda_0, xi_0, Hc2_0=None):
"""
Initialize superconductor with basic parameters.
"""
self.name = name
self.Tc = Tc # K
self.lambda_0 = lambda_0 # nm
self.xi_0 = xi_0 # nm
self.kappa = lambda_0 / xi_0
self.type = 'I' if self.kappa < 1/np.sqrt(2) else 'II'
# BCS constants
self.kB = 8.617e-5 # eV/K
self.Delta_0 = 1.76 * self.kB * Tc # BCS gap
# Critical field
if Hc2_0 is None:
# Estimate from GL theory
Phi_0 = 2.07e-15 # Wb
self.Hc2_0 = Phi_0 / (2 * np.pi * (xi_0 * 1e-9)**2)
else:
self.Hc2_0 = Hc2_0
def gap(self, T):
"""BCS gap at temperature T."""
if T >= self.Tc:
return 0
t = T / self.Tc
return self.Delta_0 * np.tanh(1.74 * np.sqrt(1/t - 1))
def lambda_L(self, T):
"""Penetration depth at temperature T."""
if T >= self.Tc:
return np.inf
t = T / self.Tc
return self.lambda_0 / np.sqrt(1 - t**4)
def xi(self, T):
"""Coherence length at temperature T."""
if T >= self.Tc:
return np.inf
t = T / self.Tc
return self.xi_0 / np.sqrt(1 - t**4)
def Hc2(self, T):
"""Upper critical field at temperature T."""
if T >= self.Tc:
return 0
return self.Hc2_0 * (1 - (T/self.Tc)**2)
def full_analysis(self):
"""Generate comprehensive analysis plots."""
T_range = np.linspace(0.1, self.Tc * 1.1, 200)
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
fig.suptitle(f'Superconductor Analysis: {self.name}', fontsize=16)
# 1. Energy gap
ax1 = axes[0, 0]
gap_vals = [self.gap(T) * 1000 for T in T_range]
ax1.plot(T_range, gap_vals, 'b-', linewidth=2)
ax1.axvline(x=self.Tc, color='r', linestyle='--', alpha=0.5)
ax1.set_xlabel('T (K)')
ax1.set_ylabel('Δ (meV)')
ax1.set_title('Energy Gap')
ax1.grid(True, alpha=0.3)
# 2. Penetration depth
ax2 = axes[0, 1]
lambda_vals = [self.lambda_L(T) if T < self.Tc else np.nan
for T in T_range]
ax2.plot(T_range, lambda_vals, 'g-', linewidth=2)
ax2.axvline(x=self.Tc, color='r', linestyle='--', alpha=0.5)
ax2.set_xlabel('T (K)')
ax2.set_ylabel('λ (nm)')
ax2.set_title('Penetration Depth')
ax2.set_ylim(0, self.lambda_0 * 5)
ax2.grid(True, alpha=0.3)
# 3. Coherence length
ax3 = axes[0, 2]
xi_vals = [self.xi(T) if T < self.Tc else np.nan for T in T_range]
ax3.plot(T_range, xi_vals, 'm-', linewidth=2)
ax3.axvline(x=self.Tc, color='r', linestyle='--', alpha=0.5)
ax3.set_xlabel('T (K)')
ax3.set_ylabel('ξ (nm)')
ax3.set_title('Coherence Length')
ax3.set_ylim(0, self.xi_0 * 5)
ax3.grid(True, alpha=0.3)
# 4. Critical field
ax4 = axes[1, 0]
Hc2_vals = [self.Hc2(T) for T in T_range]
ax4.plot(T_range, Hc2_vals, 'r-', linewidth=2)
ax4.axvline(x=self.Tc, color='r', linestyle='--', alpha=0.5)
ax4.set_xlabel('T (K)')
ax4.set_ylabel('Hc2 (T)')
ax4.set_title('Upper Critical Field')
ax4.grid(True, alpha=0.3)
# 5. Phase diagram
ax5 = axes[1, 1]
T_phase = np.linspace(0, self.Tc, 100)
H_phase = [self.Hc2(T) for T in T_phase]
ax5.plot(T_phase, H_phase, 'k-', linewidth=2)
ax5.fill_between(T_phase, 0, H_phase, alpha=0.3, color='cyan',
label='Superconducting')
ax5.fill_between(T_phase, H_phase, max(H_phase)*1.2, alpha=0.3,
color='lightcoral', label='Normal')
ax5.set_xlabel('T (K)')
ax5.set_ylabel('H (T)')
ax5.set_title('Phase Diagram')
ax5.legend()
ax5.grid(True, alpha=0.3)
# 6. Info panel
ax6 = axes[1, 2]
ax6.axis('off')
info_text = f"""
Material: {self.name}
Type: {self.type}
Critical Parameters:
Tc = {self.Tc} K
Hc2(0) = {self.Hc2_0:.2f} T
Δ(0) = {self.Delta_0*1000:.3f} meV
Length Scales:
λ(0) = {self.lambda_0} nm
ξ(0) = {self.xi_0} nm
κ = {self.kappa:.2f}
BCS Predictions:
2Δ/kBTc = {2*self.Delta_0/(self.kB*self.Tc):.2f}
"""
ax6.text(0.1, 0.9, info_text, transform=ax6.transAxes,
fontfamily='monospace', fontsize=11, verticalalignment='top')
ax6.set_title('Material Properties')
plt.tight_layout()
plt.show()
# Analyze different superconductors
materials_to_analyze = [
('Niobium', 9.3, 40, 38, 0.4),
('YBCO', 93, 150, 1.5, 100),
('MgB2', 39, 140, 5, 16),
]
for params in materials_to_analyze:
sc = SuperconductorAnalyzer(*params)
print(f"\n{'='*50}")
print(f"Analyzing: {sc.name}")
print(f"{'='*50}")
sc.full_analysis()
Summary
Key Takeaways
- Implemented BCS gap equation showing temperature dependence of energy gap
- Simulated London penetration with material-specific parameters
- Built a superconductor database for comparative analysis
- Visualized vortex states in Type II superconductors
- Created a comprehensive analyzer combining all aspects
Exercises
Exercise 1: Extend the Database
Add 5 more superconductors to the database with their actual parameters. Include both low-Tc and high-Tc materials. Regenerate the comparison plots.
Exercise 2: Temperature Sweep
Modify the vortex simulation to show how the vortex lattice spacing changes from T = 0.1Tc to T = 0.9Tc at constant applied field.
Exercise 3: Critical Current
Implement a function that calculates the critical current density Jc(T, H) using the scaling law: Jc ∝ [1 - (T/Tc)]² × [1 - (H/Hc2)]². Plot Jc as a 2D colormap.
Further Reading
- M. Tinkham, "Introduction to Superconductivity" - Comprehensive textbook
- J. Bardeen, L. Cooper, J. Schrieffer, Phys. Rev. 108, 1175 (1957) - Original BCS paper
- SuperConducting QUantum Computing resources at IBM/Google