Chapter 1: Ginzburg-Landau Theory

Learning Objectives

• Understand the fundamental concepts of GL theory as a phenomenological theory
• Explain the physical meaning of the order parameter ψ
• Derive the GL free energy functional and understand the meaning of each term
• Derive the first and second GL equations
• Understand the definitions and physical interpretations of coherence length ξ and penetration depth λ
• Explain the GL parameter κ = λ/ξ and the classification of Type I/II superconductors
• Understand how the sign of surface energy determines the superconductor type
• Visualize GL theory predictions using Python

1.1 Introduction to Phenomenological Theory

Historical Background of GL Theory

In 1950, Soviet physicists Vitaly Ginzburg and Lev Landau proposed a phenomenological theory of superconductivity. Since this theory was developed before the establishment of BCS theory (1957), it aimed to describe the macroscopic properties of superconductivity rather than explaining the microscopic mechanism.

Landau's Theory of Phase Transitions

GL theory is based on Landau's general theory of second-order phase transitions. In this theory, phase transitions are characterized by an order parameter:

• Normal state (T > Tc): Order parameter = 0
• Superconducting state (T < Tc): Order parameter ≠ 0

Near the critical temperature Tc, the order parameter is small, and the free energy can be expanded as a power series in terms of it.

1.2 The Concept of Order Parameter

Complex Scalar Field ψ(r)

In GL theory, the superconducting state is described by a complex order parameter ψ(r):

ψ(r) = |ψ(r)| exp(iθ(r))

where:

• |ψ(r)|² is proportional to the density ns of superconducting electron pairs (Cooper pairs)
• θ(r) is the phase, related to supercurrent flow
• r is the spatial coordinate

Normalization

The order parameter is usually normalized by its bulk equilibrium value:

|ψ∞|² = ns0

where ns0 is the density of superconducting electron pairs in a uniform superconductor.

1.3 Ginzburg-Landau Free Energy Functional

Construction of Free Energy

The GL free energy functional is written as the excess of superconducting free energy density over the normal state:

Fs - Fn = ∫ dV [ α|ψ|² + (β/2)|ψ|⁴ + (1/2m*)|(−iℏ∇ − q*A)ψ|² + B²/(2μ0) ]

Physical Meaning of Each Term

Temperature Dependence

The coefficient α is temperature-dependent and varies linearly near the critical temperature:

α(T) = α0(T - Tc)

T > Tc: α > 0 → ψ = 0 is minimum (normal state)
T < Tc: α < 0 → ψ ≠ 0 is minimum (superconducting state)

Python Visualization of Free Energy

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Parameters
alpha_0 = 1.0
beta = 1.0
Tc = 10.0  # Critical temperature (arbitrary units)

# Temperature range
temperatures = [12, 10, 8, 5]  # T > Tc, T = Tc, T < Tc

# Order parameter range
psi = np.linspace(-2, 2, 400)

# Plot free energy landscape
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
axes = axes.flatten()

for i, T in enumerate(temperatures):
    alpha = alpha_0 * (T - Tc)

    # Free energy density (excluding field terms)
    f = alpha * psi**2 + (beta/2) * psi**4

    axes[i].plot(psi, f, 'b-', linewidth=2)
    axes[i].axhline(y=0, color='k', linestyle='--', alpha=0.3)
    axes[i].axvline(x=0, color='k', linestyle='--', alpha=0.3)
    axes[i].set_xlabel('Order parameter ψ', fontsize=11)
    axes[i].set_ylabel('Free energy density f', fontsize=11)
    axes[i].set_title(f'T = {T} K (α = {alpha:.1f})', fontsize=12)
    axes[i].grid(True, alpha=0.3)

    # Mark minima
    if T > Tc:
        axes[i].plot(0, 0, 'ro', markersize=10, label='Minimum (normal state)')
    elif T == Tc:
        axes[i].plot(0, 0, 'yo', markersize=10, label='Minimum (critical point)')
    else:
        psi_eq = np.sqrt(-alpha / beta)
        f_min = alpha * psi_eq**2 + (beta/2) * psi_eq**4
        axes[i].plot([psi_eq, -psi_eq], [f_min, f_min], 'ro', markersize=10,
                    label='Minimum (SC state)')
        axes[i].plot(0, 0, 'go', markersize=8, alpha=0.5, label='Local maximum')

    axes[i].legend(fontsize=9)
    axes[i].set_ylim(-1, 3)

plt.tight_layout()
plt.show()

# Print equilibrium values
print("\nEquilibrium order parameter:")
for T in temperatures:
    alpha = alpha_0 * (T - Tc)
    if T <= Tc:
        psi_eq = np.sqrt(-alpha / beta)
        print(f"T = {T} K: |ψ| = {psi_eq:.3f}")
    else:
        print(f"T = {T} K: |ψ| = 0 (normal state)")

1.4 First and Second GL Equations

Variational Principle

At equilibrium, the free energy is minimized. Applying the variational principle yields two GL equations:

First GL Equation (Order Parameter Equation)

αψ + β|ψ|²ψ + (1/2m*)(−iℏ∇ − q*A)²ψ = 0

This is obtained from the variation δF/δψ* = 0.

Second GL Equation (Current Density Equation)

j = (q*/2m*i) [ψ*(−iℏ∇ − q*A)ψ − ψ(−iℏ∇ − q*A)*ψ*]
  = (q*ℏ/m*) |ψ|² [∇θ − (q*/ℏ)A]

This gives the supercurrent density, obtained from variation with respect to vector potential A.

Relationship to London Equations

Assuming a uniform order parameter (|ψ| = const), the second GL equation reduces to the London equations:

j = −(ns q*²/m*) A
∇ × j = −(ns q*²/m*) B

1.5 Coherence Length (ξ)

Definition and Physical Meaning

The GL coherence length ξ(T) represents the characteristic length scale over which the order parameter can vary spatially:

ξ²(T) = ℏ² / (2m*|α(T)|)
      = ℏ² / (2m*α0(Tc - T))

Physical Interpretation

• ξ characterizes the distance over which the order parameter changes from the normal state value (ψ = 0) to the superconducting equilibrium value
• ξ increases as temperature rises and diverges at Tc
• Typical values: ξ(0) ≈ 1-100 nm (material dependent)

Order Parameter Profile Near Boundaries

Consider a boundary between a superconductor and vacuum. In one dimension without magnetic field, the first GL equation becomes:

−(ℏ²/2m*) d²ψ/dx² + αψ + β|ψ|²ψ = 0

Boundary conditions: ψ → ψ∞ as x → ∞, dψ/dx = 0 at x = 0 (boundary condition)

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

# Parameters
hbar = 1.0545718e-34  # Reduced Planck constant (J·s)
m_star = 2 * 9.10938e-31  # Cooper pair mass (2 * electron mass)
alpha_0 = 1e-10  # Arbitrary units
Tc = 10.0
T = 5.0
alpha = alpha_0 * (T - Tc)
beta = 1e-11

# Coherence length
xi = np.sqrt(hbar**2 / (2 * m_star * abs(alpha)))

# Equilibrium order parameter
psi_infty = np.sqrt(-alpha / beta)

# Dimensionless GL equation
# Let ψ = ψ∞ · f(x/ξ), where x/ξ = s
# Then: d²f/ds² = f - f³

def gl_1d(y, s):
    """
    y[0] = f (dimensionless order parameter)
    y[1] = df/ds
    """
    f, df = y
    d2f = f - f**3  # Dimensionless GL equation
    return [df, d2f]

# Solve from bulk (s = large) back to boundary (s = 0)
s_span = np.linspace(10, 0, 1000)
y0 = [1.0, 0.0]  # Initial conditions: f = 1, df/ds = 0 at s = large

solution = odeint(gl_1d, y0, s_span)
f_profile = solution[:, 0]

# Reverse to go from boundary to bulk
s_plot = s_span[::-1]
f_plot = f_profile[::-1]

# Convert to real coordinates
x_plot = s_plot * xi * 1e9  # Convert to nm
psi_plot = f_plot * psi_infty

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

# Plot 1: Dimensionless profile
axes[0].plot(s_plot, f_plot, 'b-', linewidth=2)
axes[0].axhline(y=1, color='r', linestyle='--', label='f = 1 (bulk value)')
axes[0].axhline(y=0, color='k', linestyle='--', alpha=0.3)
axes[0].set_xlabel('x/ξ (dimensionless distance)', fontsize=12)
axes[0].set_ylabel('f = ψ/ψ∞ (dimensionless order parameter)', fontsize=12)
axes[0].set_title('Order Parameter Profile Near Boundary', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(0, 5)

# Shading to show healing length
axes[0].fill_between(s_plot, 0, f_plot, where=(s_plot <= 2),
                     alpha=0.3, color='cyan', label='ξ region')

# Plot 2: Temperature dependence of ξ
T_range = np.linspace(0.1, Tc*0.99, 100)
alpha_T = alpha_0 * (T_range - Tc)
xi_T = np.sqrt(hbar**2 / (2 * m_star * abs(alpha_T))) * 1e9  # nm

axes[1].plot(T_range, xi_T, 'r-', linewidth=2)
axes[1].axvline(x=Tc, color='b', linestyle='--', label=f'Tc = {Tc} K')
axes[1].set_xlabel('Temperature T (K)', fontsize=12)
axes[1].set_ylabel('Coherence length ξ(T) (nm)', fontsize=12)
axes[1].set_title('Temperature-Dependent Coherence Length', fontsize=14)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
axes[1].set_ylim(0, 100)

plt.tight_layout()
plt.show()

print(f"\nCoherence length ξ(T={T} K) = {xi*1e9:.2f} nm")
print(f"Equilibrium order parameter ψ∞ = {psi_infty:.3e}")

1.6 Penetration Depth (λ)

Definition and Physical Meaning

The GL penetration depth λ(T) represents the characteristic depth to which a magnetic field can penetrate into a superconductor:

λ²(T) = m* β / (μ0 q*² |ψ∞|²)
      = m* / (μ0 q*² ns)

where ns = |ψ∞|² / β is the superconducting electron pair density.

Physical Interpretation

• λ is the exponential decay length of magnetic field from surface to interior
• λ is equivalent to the London penetration depth
• Typical values: λ(0) ≈ 20-500 nm (material dependent)

Mathematical Description of Meissner Effect

For a semi-infinite superconductor (x > 0) with parallel applied field B(0) = B0, the internal field is:

B(x) = B0 exp(−x/λ)

import numpy as np
import matplotlib.pyplot as plt

# Parameters
lambda_L = 50e-9  # London penetration depth (m)
B0 = 0.1  # Applied field (T)

# Distance from surface
x = np.linspace(0, 300e-9, 500)  # 0 to 300 nm

# Magnetic field profile
B = B0 * np.exp(-x / lambda_L)

# Screening current density (from Ampere's law)
# ∇ × B = μ0 j  →  dB/dx = −μ0 j_y
mu0 = 4 * np.pi * 1e-7  # Vacuum permeability
j = (B0 / (mu0 * lambda_L)) * np.exp(-x / lambda_L)

# Convert to more convenient units
x_nm = x * 1e9  # nm
lambda_nm = lambda_L * 1e9  # nm
j_MA = j * 1e-6  # MA/m²

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

# Plot 1: Magnetic field penetration
axes[0].plot(x_nm, B / B0, 'b-', linewidth=2)
axes[0].axhline(y=np.exp(-1), color='r', linestyle='--',
               label=f'B/B₀ = 1/e (x = λ = {lambda_nm:.0f} nm)')
axes[0].axvline(x=lambda_nm, color='r', linestyle='--', alpha=0.5)
axes[0].fill_between(x_nm, 0, B / B0, where=(x_nm <= lambda_nm),
                    alpha=0.3, color='cyan', label='Penetration region')
axes[0].set_xlabel('Distance from surface x (nm)', fontsize=12)
axes[0].set_ylabel('Normalized field B/B₀', fontsize=12)
axes[0].set_title('Magnetic Field Penetration into Superconductor', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(0, 300)

# Plot 2: Screening current
axes[1].plot(x_nm, j_MA / j_MA[0], 'g-', linewidth=2)
axes[1].axvline(x=lambda_nm, color='r', linestyle='--',
               label=f'x = λ = {lambda_nm:.0f} nm')
axes[1].fill_between(x_nm, 0, j_MA / j_MA[0], where=(x_nm <= lambda_nm),
                    alpha=0.3, color='yellow', label='Screening current region')
axes[1].set_xlabel('Distance from surface x (nm)', fontsize=12)
axes[1].set_ylabel('Normalized current density j/j₀', fontsize=12)
axes[1].set_title('Meissner Screening Current', fontsize=14)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim(0, 300)

plt.tight_layout()
plt.show()

# Print penetration depth info
print(f"\nPenetration depth λ = {lambda_nm:.1f} nm")
print(f"At x = λ, field decays to {np.exp(-1)*100:.1f}% of initial value")
print(f"At x = 3λ, field decays to {np.exp(-3)*100:.2f}% of initial value")

1.7 GL Parameter κ = λ/ξ

Definition and Significance

The Ginzburg-Landau parameter κ is defined as the ratio of the two characteristic lengths:

κ = λ / ξ

κ is a material-specific dimensionless parameter and one of the most important quantities determining superconductor properties.

Physical Meaning of κ

• ξ is the scale of order parameter variation (superconducting "stiffness")
• λ is the scale of magnetic field penetration (magnetic response)
• κ represents the competition between field penetration and order parameter suppression

1.8 Type I/II Classification (κ < 1/√2 vs κ > 1/√2)

Concept of Surface Energy

The boundary between superconducting and normal states has an associated energy cost. The sign of this surface energy determines the superconductor type.

Origin of Surface Energy

• Positive energy cost: Order parameter decreases to 0 over distance ξ (loss of condensation energy)
• Negative energy gain: Field penetrates over distance λ (reduction of magnetic energy)

Surface energy σns ∝ (ξ/λ - λ/ξ) ∝ (1/κ - κ)

Critical Value κc = 1/√2 ≈ 0.707

Type I Superconductors (κ < 1/√2)

• Positive surface energy minimizes boundaries between superconducting and normal regions
• Sharp transition to normal state at single critical field Hc
• Intermediate state (mixed but with macroscopic domain structure)

Type II Superconductors (κ > 1/√2)

• Negative surface energy favors boundary formation
• Two critical fields Hc1 and Hc2
• Mixed state (vortex state, Abrikosov vortex lattice) appears for Hc1 < H < Hc2
• Each vortex carries quantized flux Φ0 = h/(2e)

import numpy as np
import matplotlib.pyplot as plt

# κ values
kappa_values = [0.3, 0.707, 2.0, 10.0]
labels = ['Type I (κ=0.3)', 'Critical (κ=0.707)', 'Type II (κ=2)', 'Type II (κ=10)']
colors = ['blue', 'green', 'orange', 'red']

# Dimensionless distance from boundary
x = np.linspace(0, 5, 500)

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

for i, (kappa, label, color) in enumerate(zip(kappa_values, labels, colors)):
    ax = axes[i // 2, i % 2]

    # Order parameter profile (approximate)
    # ψ(x) ≈ ψ∞ tanh(x / (√2 ξ))
    psi_profile = np.tanh(x / np.sqrt(2))

    # Magnetic field profile (approximate)
    # B(x) ≈ B0 (1 - exp(-x/λ))
    # In units of ξ: x/λ = x/(κξ)
    B_profile = 1 - np.exp(-x / kappa)

    # Plot both profiles
    ax.plot(x, psi_profile, 'b-', linewidth=2, label='|ψ|/ψ∞ (order parameter)')
    ax.plot(x, B_profile, 'r-', linewidth=2, label='B/B₀ (magnetic field)')

    # Mark characteristic lengths
    ax.axvline(x=1, color='b', linestyle='--', alpha=0.5, label='ξ')
    ax.axvline(x=kappa, color='r', linestyle='--', alpha=0.5, label='λ = κξ')

    # Shading
    ax.fill_between(x, 0, psi_profile, alpha=0.2, color='blue')
    ax.fill_between(x, 0, B_profile, alpha=0.2, color='red')

    ax.set_xlabel('Distance (ξ units)', fontsize=11)
    ax.set_ylabel('Normalized value', fontsize=11)
    ax.set_title(label, fontsize=12, color=color)
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)
    ax.set_xlim(0, 5)
    ax.set_ylim(0, 1.1)

plt.tight_layout()
plt.show()

# Surface energy comparison
print("\nSurface energy sign:")
for kappa, label in zip(kappa_values, labels):
    sigma_sign = (1/kappa - kappa)
    if sigma_sign > 0:
        energy_type = "Positive (avoids boundaries)"
    elif sigma_sign < 0:
        energy_type = "Negative (favors boundaries)"
    else:
        energy_type = "Zero (critical)"
    print(f"{label}: σns ∝ {sigma_sign:.3f} → {energy_type}")

1.9 Surface Energy and Boundary Solutions

Details of Superconductor-Normal Boundary

Consider a planar boundary (x = 0) between a superconductor and normal metal. x < 0 is superconducting, x > 0 is normal.

Behavior at the Boundary

• For κ < 1/√2: ψ and B vary over different length scales (ξ < λ), resulting in net energy cost
• For κ > 1/√2: ψ and B vary over different length scales (ξ > λ), resulting in net energy gain

Vortex State (Abrikosov Lattice)

In Type II superconductors, for Hc1 < H < Hc2, flux penetrates as quantized vortices:

Φ0 = h / (2e) = 2.07 × 10⁻¹⁵ Wb (flux quantum)

Each vortex has:

• ψ = 0 at center (normal core, radius ≈ ξ)
• Magnetic field concentrated around vortex (penetration ≈ λ)
• Carries flux Φ0

import numpy as np
import matplotlib.pyplot as plt

# Simulate a single Abrikosov vortex
# Parameters
xi = 10  # Coherence length (nm)
lambda_L = 50  # Penetration depth (nm)
kappa = lambda_L / xi

# 2D grid
x = np.linspace(-100, 100, 200)
y = np.linspace(-100, 100, 200)
X, Y = np.meshgrid(x, y)
R = np.sqrt(X**2 + Y**2)

# Order parameter profile (approximate)
# |ψ(r)| ≈ ψ∞ tanh(r / (√2 ξ))
psi = np.tanh(R / (np.sqrt(2) * xi))

# Magnetic field profile (approximate)
# B(r) ≈ B0 K0(r/λ) (K0: modified Bessel function of second kind)
# Approximation: B(r) ≈ exp(-r/λ) for visualization
B = np.exp(-R / lambda_L)

# Plotting
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# Plot 1: Order parameter
im1 = axes[0].contourf(X, Y, psi, levels=20, cmap='Blues')
axes[0].contour(X, Y, psi, levels=10, colors='black', alpha=0.3, linewidths=0.5)
axes[0].set_xlabel('x (nm)', fontsize=11)
axes[0].set_ylabel('y (nm)', fontsize=11)
axes[0].set_title(f'Order parameter |ψ|/ψ∞\n(vortex core ~ ξ = {xi} nm)', fontsize=12)
axes[0].set_aspect('equal')
plt.colorbar(im1, ax=axes[0], label='|ψ|/ψ∞')

# Add circle at ξ
circle_xi = plt.Circle((0, 0), xi, fill=False, color='red',
                       linestyle='--', linewidth=2, label='r = ξ')
axes[0].add_patch(circle_xi)
axes[0].legend(fontsize=9)

# Plot 2: Magnetic field
im2 = axes[1].contourf(X, Y, B, levels=20, cmap='Reds')
axes[1].contour(X, Y, B, levels=10, colors='black', alpha=0.3, linewidths=0.5)
axes[1].set_xlabel('x (nm)', fontsize=11)
axes[1].set_ylabel('y (nm)', fontsize=11)
axes[1].set_title(f'Magnetic field B/B₀\n(penetration ~ λ = {lambda_L} nm)', fontsize=12)
axes[1].set_aspect('equal')
plt.colorbar(im2, ax=axes[1], label='B/B₀')

# Add circle at λ
circle_lambda = plt.Circle((0, 0), lambda_L, fill=False, color='blue',
                          linestyle='--', linewidth=2, label='r = λ')
axes[1].add_patch(circle_lambda)
axes[1].legend(fontsize=9)

# Plot 3: Radial profiles
r_profile = np.linspace(0, 100, 500)
psi_profile = np.tanh(r_profile / (np.sqrt(2) * xi))
B_profile = np.exp(-r_profile / lambda_L)

axes[2].plot(r_profile, psi_profile, 'b-', linewidth=2, label='|ψ|/ψ∞')
axes[2].plot(r_profile, B_profile, 'r-', linewidth=2, label='B/B₀')
axes[2].axvline(x=xi, color='b', linestyle='--', alpha=0.5, label=f'ξ = {xi} nm')
axes[2].axvline(x=lambda_L, color='r', linestyle='--', alpha=0.5,
               label=f'λ = {lambda_L} nm')
axes[2].set_xlabel('Radius r (nm)', fontsize=11)
axes[2].set_ylabel('Normalized value', fontsize=11)
axes[2].set_title(f'Radial profiles\n(κ = λ/ξ = {kappa:.1f})', fontsize=12)
axes[2].legend(fontsize=9)
axes[2].grid(True, alpha=0.3)
axes[2].set_xlim(0, 100)

plt.tight_layout()
plt.show()

# Print vortex info
Phi0 = 2.067833848e-15  # Wb
print(f"\nAbrikosov vortex structure (κ = {kappa:.1f}):")
print(f"Vortex core size (normal state): ~ξ = {xi} nm")
print(f"Magnetic field penetration: ~λ = {lambda_L} nm")
print(f"Flux quantum: Φ₀ = {Phi0:.3e} Wb")

1.10 Phase Diagram: Type I vs Type II

Complete Classification Based on κ

import numpy as np
import matplotlib.pyplot as plt

# Range of kappa values
kappa = np.logspace(-1, 2, 500)  # 0.1 to 100

# Critical kappa
kappa_c = 1 / np.sqrt(2)

# Surface energy (normalized)
sigma = 1/kappa - kappa

# Thermodynamic critical field (normalized to Hc)
# For Type I: H = Hc
# For Type II: Hc1 < H < Hc2
Hc1_ratio = np.log(kappa) / (np.sqrt(2) * kappa)  # Hc1/Hc (approximate)
Hc2_ratio = np.sqrt(2) * kappa  # Hc2/Hc

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

# Plot 1: Surface energy
axes[0, 0].plot(kappa, sigma, 'b-', linewidth=2)
axes[0, 0].axhline(y=0, color='k', linestyle='-', linewidth=1)
axes[0, 0].axvline(x=kappa_c, color='r', linestyle='--', linewidth=2,
                  label=f'κc = 1/√2 ≈ {kappa_c:.3f}')
axes[0, 0].fill_between(kappa, sigma, 0, where=(sigma > 0), alpha=0.3,
                       color='blue', label='Type I (σ > 0)')
axes[0, 0].fill_between(kappa, sigma, 0, where=(sigma < 0), alpha=0.3,
                       color='red', label='Type II (σ < 0)')
axes[0, 0].set_xlabel('GL parameter κ', fontsize=11)
axes[0, 0].set_ylabel('Surface energy σns (arb. units)', fontsize=11)
axes[0, 0].set_title('Surface Energy vs κ', fontsize=12)
axes[0, 0].set_xscale('log')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3, which='both')
axes[0, 0].set_xlim(0.1, 100)
axes[0, 0].set_ylim(-5, 5)

# Plot 2: Critical fields
axes[0, 1].plot(kappa[kappa > kappa_c], Hc1_ratio[kappa > kappa_c],
               'g-', linewidth=2, label='Hc1/Hc')
axes[0, 1].plot(kappa[kappa > kappa_c], Hc2_ratio[kappa > kappa_c],
               'r-', linewidth=2, label='Hc2/Hc')
axes[0, 1].axhline(y=1, color='b', linestyle='--', linewidth=1, label='Hc')
axes[0, 1].axvline(x=kappa_c, color='orange', linestyle='--', linewidth=2,
                  label=f'κc = {kappa_c:.3f}')
axes[0, 1].set_xlabel('GL parameter κ', fontsize=11)
axes[0, 1].set_ylabel('Normalized critical field H/Hc', fontsize=11)
axes[0, 1].set_title('Critical Fields vs κ (Type II only)', fontsize=12)
axes[0, 1].set_xscale('log')
axes[0, 1].set_yscale('log')
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3, which='both')
axes[0, 1].set_xlim(0.7, 100)
axes[0, 1].set_ylim(0.01, 100)

# Plot 3: Phase diagram (H vs κ)
fig3 = axes[1, 0]
kappa_type1 = np.linspace(0.1, kappa_c, 100)
kappa_type2 = np.linspace(kappa_c, 10, 100)

# Type I
fig3.fill_between([0, kappa_c], [0, 0], [2, 2], alpha=0.3, color='blue',
                  label='Type I SC region')
fig3.plot([kappa_c, kappa_c], [0, 2], 'k-', linewidth=2)

# Type II
Hc1_type2 = np.log(kappa_type2) / (np.sqrt(2) * kappa_type2)
Hc2_type2 = np.sqrt(2) * kappa_type2

fig3.fill_between(kappa_type2, 0, Hc1_type2, alpha=0.3, color='cyan',
                  label='Meissner state')
fig3.fill_between(kappa_type2, Hc1_type2, np.minimum(Hc2_type2, 2),
                  alpha=0.3, color='yellow', label='Mixed state (vortex)')
fig3.fill_between(kappa_type2, np.minimum(Hc2_type2, 2), 2,
                  alpha=0.3, color='lightgray', label='Normal state')

fig3.plot(kappa_type2, Hc1_type2, 'g-', linewidth=2, label='Hc1')
fig3.plot(kappa_type2, Hc2_type2, 'r-', linewidth=2, label='Hc2')
fig3.axvline(x=kappa_c, color='orange', linestyle='--', linewidth=2)

fig3.set_xlabel('GL parameter κ', fontsize=11)
fig3.set_ylabel('Magnetic field H/Hc', fontsize=11)
fig3.set_title('Superconducting Phase Diagram (κ-H)', fontsize=12)
fig3.legend(fontsize=8, loc='upper left')
fig3.grid(True, alpha=0.3)
fig3.set_xlim(0, 10)
fig3.set_ylim(0, 2)

# Plot 4: Material examples
fig4 = axes[1, 1]
materials = {
    'Al': 0.15, 'Pb': 0.4, 'Sn': 0.2, 'Nb': 0.8,
    'NbTi': 1.5, 'Nb3Sn': 2.5, 'YBCO': 50, 'Bi2212': 80
}
type1 = {k: v for k, v in materials.items() if v < kappa_c}
type2 = {k: v for k, v in materials.items() if v >= kappa_c}

y_pos_type1 = np.arange(len(type1))
y_pos_type2 = np.arange(len(type2)) + len(type1) + 0.5

fig4.barh(y_pos_type1, list(type1.values()), color='blue', alpha=0.7,
         label='Type I')
fig4.barh(y_pos_type2, list(type2.values()), color='red', alpha=0.7,
         label='Type II')

fig4.axvline(x=kappa_c, color='green', linestyle='--', linewidth=2,
            label=f'κc = {kappa_c:.3f}')

all_materials = list(type1.keys()) + list(type2.keys())
y_pos_all = list(y_pos_type1) + list(y_pos_type2)
fig4.set_yticks(y_pos_all)
fig4.set_yticklabels(all_materials)
fig4.set_xlabel('GL parameter κ', fontsize=11)
fig4.set_title('κ Values of Real Superconductors', fontsize=12)
fig4.set_xscale('log')
fig4.legend(fontsize=9)
fig4.grid(True, alpha=0.3, axis='x')
fig4.set_xlim(0.1, 100)

plt.tight_layout()
plt.show()

# Print classification
print("\nSuperconductor classification:")
print(f"Critical value: κc = 1/√2 = {kappa_c:.4f}\n")
print("Type I superconductors (κ < κc):")
for mat, k in sorted(type1.items(), key=lambda x: x[1]):
    print(f"  {mat}: κ = {k:.2f}")
print("\nType II superconductors (κ > κc):")
for mat, k in sorted(type2.items(), key=lambda x: x[1]):
    print(f"  {mat}: κ = {k:.1f}")

Summary

Exercises