Chapter 2: Type II Superconductivity and Vortex Physics

Learning Objectives

Understand why κ > 1/√2 leads to fundamentally different magnetic behavior
Derive the critical fields H_c1, H_c2, and their relationship to H_c
Master the concept of the mixed state (vortex state)
Derive flux quantization from single-valuedness of the wavefunction
Understand Abrikosov vortex lattice structure and energy minimization
Analyze vortex dynamics: Lorentz force, flux flow, and flux creep
Understand flux pinning mechanisms and critical current density
Apply the Bean critical state model
Implement numerical simulations of vortex physics

2.1 Type II Superconductors: The κ > 1/√2 Regime

Why κ Matters: Surface Energy Revisited

From Chapter 1, we learned that the GL parameter κ = λ/ξ determines the sign of the surface energy between superconducting and normal regions:

Surface Energy Scaling:

$$\sigma_{ns} \propto \xi(1 - \kappa\sqrt{2})$$

κ < 1/√2: σ_ns > 0 (positive surface energy) → Type I
κ > 1/√2: σ_ns < 0 (negative surface energy) → Type II

Physical Consequences of Negative Surface Energy

When σ_ns < 0, the system gains energy by creating superconducting-normal interfaces. This has profound implications:

Complete Meissner expulsion becomes unstable above a threshold field H_c1
Magnetic flux penetrates in quantized units (vortices)
Mixed state forms with coexisting superconducting and normal regions
Much higher critical fields H_c2 >> H_c1 become possible

Historical Context

Alexei Abrikosov predicted the vortex state in 1957 using GL theory, explaining Shubnikov's 1937 experimental observations of Type II behavior. Abrikosov received the Nobel Prize in Physics in 2003 for this work.

2.2 Critical Magnetic Fields

Three Critical Fields

Type II superconductors are characterized by three fundamental field scales:

Lower Critical Field H_c1:

$$H_{c1} = \frac{\Phi_0}{4\pi\lambda^2} \ln(\kappa)$$

where Φ₀ = h/2e is the flux quantum. Below H_c1, the Meissner state is stable.

Upper Critical Field H_c2:

$$H_{c2} = \frac{\Phi_0}{2\pi\xi^2}$$

Above H_c2, superconductivity is completely destroyed.

Thermodynamic Critical Field H_c:

$$\frac{B_c^2}{2\mu_0} = f_n - f_s$$

This represents the condensation energy density. It relates to H_c1 and H_c2 through:

$$H_{c1} \approx \frac{H_c}{\kappa} \ln(\kappa)$$

$$H_{c2} = \sqrt{2}\,\kappa\,H_c$$

Derivation of H_c1

The lower critical field corresponds to the field at which the first vortex becomes energetically favorable. Consider a single vortex in an infinite superconductor:

Energy of single vortex:

$$E_v \approx \frac{\Phi_0^2}{4\pi\mu_0\lambda^2} \ln\left(\frac{\lambda}{\xi}\right) = \frac{\Phi_0^2}{4\pi\mu_0\lambda^2} \ln(\kappa)$$

This energy must be compared to the Zeeman energy saved by admitting flux:

$$E_{Zeeman} = -\Phi_0 H_{ext}$$

Setting E_v + E_Zeeman = 0 gives:

$$H_{c1} = \frac{\Phi_0}{4\pi\lambda^2} \ln(\kappa)$$

Derivation of H_c2

The upper critical field is reached when the vortex cores overlap completely. At this point, the average order parameter goes to zero. From the first GL equation in a uniform field:

$$\alpha\psi + \beta|\psi|^2\psi + \frac{1}{2m^*}\left(-i\hbar\nabla - e^*\mathbf{A}\right)^2\psi = 0$$

For the lowest Landau level (at the onset of superconductivity), |ψ|² → 0, and:

$$\frac{1}{2m^*}\left(-i\hbar\nabla - e^*\mathbf{A}\right)^2\psi = -\alpha\psi$$

The eigenvalue of the kinetic operator is ℏe*B/m* = ℏe*μ₀H/m*, giving:

$$\alpha = -\frac{\hbar e^* H_{c2}}{m^*}$$

Using α = -ℏ²/(2m*ξ²) and e* = 2e:

$$H_{c2} = \frac{\hbar}{2e\xi^2} = \frac{\Phi_0}{2\pi\xi^2}$$

Temperature Dependence

All three critical fields vanish at T_c and grow as temperature decreases:

import numpy as np
import matplotlib.pyplot as plt

# Temperature dependence of critical fields
def critical_fields_vs_T(T, Tc, Hc0, kappa):
    """
    Calculate critical fields as function of temperature

    Parameters:
    T: temperature array
    Tc: critical temperature
    Hc0: thermodynamic critical field at T=0
    kappa: GL parameter

    Returns:
    Hc, Hc1, Hc2 at each temperature
    """
    t = T / Tc  # Reduced temperature

    # Near Tc: Hc ∝ (1 - t)
    # At T=0: Hc = Hc0
    # Approximate form: Hc(T) = Hc0 * (1 - t²)
    Hc = Hc0 * (1 - t**2)

    # H_c1 = (Hc / kappa) * ln(kappa)
    Hc1 = (Hc / kappa) * np.log(kappa)

    # H_c2 = sqrt(2) * kappa * Hc
    Hc2 = np.sqrt(2) * kappa * Hc

    return Hc, Hc1, Hc2

# Create figure with two subplots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# Temperature array
T = np.linspace(0, 1, 200)  # Reduced temperature T/Tc

# Parameters for different materials
materials = [
    ('Nb (κ=1.0)', 1.0, 0.2, 'blue'),
    ('NbTi (κ=50)', 50, 0.15, 'red'),
    ('Nb₃Sn (κ=22)', 22, 0.3, 'green'),
]

# Left plot: Critical fields vs temperature
for name, kappa, Hc0, color in materials:
    Hc, Hc1, Hc2 = critical_fields_vs_T(T, 1.0, Hc0, kappa)

    ax1.plot(T, Hc, '--', color=color, alpha=0.6, linewidth=1.5, label=f'{name}: Hc')
    ax1.plot(T, Hc1, '-', color=color, alpha=0.8, linewidth=2, label=f'{name}: Hc1')
    ax1.plot(T, Hc2, '-', color=color, linewidth=2.5, label=f'{name}: Hc2')

ax1.set_xlabel('Reduced temperature T/Tc', fontsize=12)
ax1.set_ylabel('Critical fields (normalized)', fontsize=12)
ax1.set_title('Temperature Dependence of Critical Fields', fontsize=14)
ax1.legend(fontsize=9, ncol=2, loc='upper right')
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 1)
ax1.set_ylim(0, None)

# Right plot: H-T phase diagram for Type II
T_phase = np.linspace(0, 1, 100)
kappa_demo = 10  # Representative Type II
Hc0_demo = 1.0

Hc_phase, Hc1_phase, Hc2_phase = critical_fields_vs_T(T_phase, 1.0, Hc0_demo, kappa_demo)

# Create field array
H_max = 15
H_phase = np.linspace(0, H_max, 300)
T_mesh, H_mesh = np.meshgrid(T_phase, H_phase)

# Calculate Hc1 and Hc2 for mesh
Hc_mesh = Hc0_demo * (1 - T_mesh**2)
Hc1_mesh = (Hc_mesh / kappa_demo) * np.log(kappa_demo)
Hc2_mesh = np.sqrt(2) * kappa_demo * Hc_mesh

# Create phase diagram
phase_diagram = np.zeros_like(H_mesh)
phase_diagram[H_mesh < Hc1_mesh] = 1  # Meissner
phase_diagram[(H_mesh >= Hc1_mesh) & (H_mesh < Hc2_mesh)] = 2  # Mixed
phase_diagram[H_mesh >= Hc2_mesh] = 3  # Normal

ax2.contourf(T_mesh, H_mesh, phase_diagram, levels=[0.5, 1.5, 2.5, 3.5],
             colors=['cyan', 'yellow', 'lightgray'], alpha=0.6)
ax2.plot(T_phase, Hc1_phase, 'b-', linewidth=3, label='Hc1(T)')
ax2.plot(T_phase, Hc2_phase, 'r-', linewidth=3, label='Hc2(T)')
ax2.plot(T_phase, Hc_phase, 'k--', linewidth=2, alpha=0.7, label='Hc(T)')

# Add labels for phases
ax2.text(0.5, 0.3, 'Meissner State\n(Complete flux expulsion)',
         ha='center', va='center', fontsize=11,
         bbox=dict(boxstyle='round', facecolor='cyan', alpha=0.8))
ax2.text(0.5, 4, 'Mixed State\n(Vortex lattice)',
         ha='center', va='center', fontsize=11,
         bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.8))
ax2.text(0.5, 10, 'Normal State',
         ha='center', va='center', fontsize=11,
         bbox=dict(boxstyle='round', facecolor='lightgray', alpha=0.8))

ax2.set_xlabel('Reduced temperature T/Tc', fontsize=12)
ax2.set_ylabel('Applied field H (normalized)', fontsize=12)
ax2.set_title(f'H-T Phase Diagram (κ = {kappa_demo})', fontsize=14)
ax2.legend(fontsize=11, loc='upper right')
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, 1)
ax2.set_ylim(0, H_max)

plt.tight_layout()
plt.savefig('critical_fields_temperature.png', dpi=150, bbox_inches='tight')
plt.show()

print("Critical Field Relationships:")
print(f"κ = {kappa_demo}")
print(f"At T = 0:")
print(f"  Hc2/Hc1 = {(np.sqrt(2) * kappa_demo**2 / np.log(kappa_demo)):.1f}")
print(f"  Hc2/Hc = {np.sqrt(2) * kappa_demo:.1f}")
print(f"  Hc1/Hc = {np.log(kappa_demo) / kappa_demo:.3f}")

2.3 Flux Quantization

Derivation from Single-Valuedness

One of the most remarkable properties of superconductors is flux quantization. The magnetic flux through any closed loop in a superconductor must be quantized in units of Φ₀.

Flux Quantum:

$$\Phi_0 = \frac{h}{2e} = 2.067 \times 10^{-15}\text{ Wb}$$

This fundamental constant arises from the quantum nature of Cooper pairs.

Rigorous Derivation

Consider a closed loop C deep inside the superconductor where the order parameter has uniform magnitude |ψ| = |ψ_∞| and only the phase varies:

$$\psi(\mathbf{r}) = |\psi_\infty| e^{i\theta(\mathbf{r})}$$

The supercurrent density from GL theory is:

$$\mathbf{J}_s = -\frac{e^*}{m^*}\text{Re}\left[\psi^*\left(-i\hbar\nabla - e^*\mathbf{A}\right)\psi\right]$$

Substituting ψ = |ψ|e^iθ:

$$\mathbf{J}_s = \frac{e^*|\psi|^2}{m^*}\left(\hbar\nabla\theta - e^*\mathbf{A}\right)$$

Deep inside the superconductor, J_s = 0, so:

$$\nabla\theta = \frac{e^*}{\hbar}\mathbf{A}$$

Integrating around closed loop C:

$$\oint_C \nabla\theta \cdot d\mathbf{l} = \frac{e^*}{\hbar}\oint_C \mathbf{A} \cdot d\mathbf{l} = \frac{e^*}{\hbar}\Phi$$

The left side must be 2πn (n = integer) due to single-valuedness of ψ:

$$2\pi n = \frac{e^*}{\hbar}\Phi = \frac{2e}{\hbar}\Phi$$

Therefore:

$$\Phi = n\frac{h}{2e} = n\Phi_0$$

Key Insight: The Factor of 2

The factor of 2 in Φ₀ = h/2e is experimental proof that superconducting charge carriers have charge 2e (Cooper pairs), not e (single electrons). This was confirmed by flux quantization experiments in 1961.

import numpy as np
import scipy.constants as const

# Physical constants
h = const.h          # Planck constant (J·s)
e = const.e          # Elementary charge (C)

# Flux quantum
Phi_0 = h / (2 * e)

print(f"Flux quantum Φ₀ = {Phi_0:.4e} Wb")
print(f"Φ₀ = {Phi_0*1e15:.4f} × 10⁻¹⁵ Wb")

# How many flux quanta in a typical magnetic field?
# Example: 1 Tesla field through a 1 mm² area
B = 1.0  # Tesla
A = 1e-6  # m² (1 mm²)
Phi = B * A
n = Phi / Phi_0

print(f"\nFlux through 1 mm² in 1 T field:")
print(f"Φ = {Phi:.4e} Wb")
print(f"Number of flux quanta n = {n:.2e}")

# Superconducting ring example
r_ring = 1e-3  # 1 mm radius
A_ring = np.pi * r_ring**2
B_weak = 0.01  # 10 mT
Phi_ring = B_weak * A_ring
n_ring = Phi_ring / Phi_0

print(f"\nSuperconducting ring (r = 1 mm) in 10 mT field:")
print(f"Flux through ring = {Phi_ring:.4e} Wb")
print(f"Number of flux quanta = {n_ring:.2f}")

2.4 Abrikosov Vortex Structure

Anatomy of a Single Vortex

An Abrikosov vortex (or flux vortex) is a topological defect carrying exactly one flux quantum Φ₀:

Vortex Structure:

Core region (r < ξ): Order parameter suppressed, |ψ(r)| ≈ 0 at r = 0
Magnetic field (r < λ): Field concentrated near core, B(r) ~ exp(-r/λ)
Phase circulation: θ(φ) = φ (winds by 2π around core)
Supercurrent: Circulates around core with J_s(r) ~ 1/r for ξ < r < λ

Detailed Vortex Profile

For an isolated vortex, the order parameter and magnetic field can be approximated as:

$$|\psi(r)| \approx |\psi_\infty| \tanh\left(\frac{r}{\sqrt{2}\xi}\right)$$

$$B(r) \approx \frac{\Phi_0}{2\pi\lambda^2} K_0\left(\frac{r}{\lambda}\right)$$

where K₀ is the modified Bessel function of the second kind.

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import k0, k1  # Modified Bessel functions

def vortex_profile(r, xi, lambda_L):
    """
    Calculate order parameter and field for isolated vortex

    Parameters:
    r: radial distance
    xi: coherence length
    lambda_L: penetration depth

    Returns:
    psi: normalized order parameter |ψ(r)|/|ψ∞|
    B: normalized magnetic field B(r)/(Φ₀/2πλ²)
    J: supercurrent density (normalized)
    """
    # Order parameter (tanh profile)
    psi = np.tanh(r / (np.sqrt(2) * xi))

    # Magnetic field (Bessel function, with cutoff for r→0)
    r_safe = np.where(r < 0.01 * lambda_L, 0.01 * lambda_L, r)
    B = k0(r_safe / lambda_L)

    # Supercurrent J ∝ dB/dr
    J = -k1(r_safe / lambda_L) / lambda_L

    return psi, B, J

# Physical parameters
xi = 10  # nm
kappa = 5
lambda_L = kappa * xi  # nm

# Radial coordinate
r = np.linspace(0, 10 * lambda_L, 1000)

psi, B, J = vortex_profile(r, xi, lambda_L)

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

# Plot 1: Order parameter
ax1 = axes[0, 0]
ax1.plot(r / xi, psi, 'b-', linewidth=2.5)
ax1.axvline(x=1, color='b', linestyle='--', alpha=0.5, label=f'r = ξ')
ax1.axvline(x=kappa, color='r', linestyle='--', alpha=0.5, label=f'r = λ = {kappa}ξ')
ax1.axhline(y=1/np.sqrt(2), color='gray', linestyle=':', alpha=0.5,
            label=f'|ψ| = 1/√2')
ax1.fill_between(r / xi, 0, psi, alpha=0.2, color='blue')
ax1.set_xlabel('Radial distance r/ξ', fontsize=12)
ax1.set_ylabel('|ψ(r)| / |ψ∞|', fontsize=12)
ax1.set_title('Order Parameter Profile', fontsize=14)
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 50)
ax1.set_ylim(0, 1.1)

# Annotate core region
ax1.axvspan(0, 1, alpha=0.15, color='red', label='Vortex core')
ax1.text(0.5, 0.5, 'Core\n(Normal)', ha='center', va='center',
         fontsize=10, bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

# Plot 2: Magnetic field
ax2 = axes[0, 1]
ax2.plot(r / xi, B / B[0], 'r-', linewidth=2.5)
ax2.axvline(x=1, color='b', linestyle='--', alpha=0.5, label='r = ξ')
ax2.axvline(x=kappa, color='r', linestyle='--', alpha=0.5, label='r = λ')
ax2.axhline(y=1/np.e, color='gray', linestyle=':', alpha=0.5,
            label='B = B₀/e')
ax2.fill_between(r / xi, 0, B / B[0], alpha=0.2, color='red')
ax2.set_xlabel('Radial distance r/ξ', fontsize=12)
ax2.set_ylabel('B(r) / B(0)', fontsize=12)
ax2.set_title('Magnetic Field Profile (Bessel K₀)', fontsize=14)
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, 50)
ax2.set_ylim(0, 1.1)
ax2.set_yscale('log')

# Plot 3: Supercurrent density
ax3 = axes[1, 0]
J_normalized = -J / np.max(np.abs(J))
ax3.plot(r / xi, J_normalized, 'g-', linewidth=2.5)
ax3.axvline(x=1, color='b', linestyle='--', alpha=0.5, label='r = ξ')
ax3.axvline(x=kappa, color='r', linestyle='--', alpha=0.5, label='r = λ')
ax3.fill_between(r / xi, 0, J_normalized, alpha=0.2, color='green')
ax3.set_xlabel('Radial distance r/ξ', fontsize=12)
ax3.set_ylabel('Jₛ(r) (normalized)', fontsize=12)
ax3.set_title('Supercurrent Density (Bessel K₁)', fontsize=14)
ax3.legend(fontsize=11)
ax3.grid(True, alpha=0.3)
ax3.set_xlim(0, 50)

# Add annotation for peak
peak_idx = np.argmax(np.abs(J))
ax3.plot(r[peak_idx] / xi, J_normalized[peak_idx], 'ro', markersize=8)
ax3.annotate(f'Peak at r ≈ {r[peak_idx]/xi:.1f}ξ',
             xy=(r[peak_idx]/xi, J_normalized[peak_idx]),
             xytext=(r[peak_idx]/xi + 10, J_normalized[peak_idx] - 0.2),
             arrowprops=dict(arrowstyle='->', color='red'),
             fontsize=10)

# Plot 4: 2D visualization
ax4 = axes[1, 1]
# Create 2D grid
x = np.linspace(-5*lambda_L, 5*lambda_L, 200)
y = np.linspace(-5*lambda_L, 5*lambda_L, 200)
X, Y = np.meshgrid(x, y)
R = np.sqrt(X**2 + Y**2)

# Calculate field
_, B_2d, _ = vortex_profile(R, xi, lambda_L)

# Plot
im = ax4.contourf(X/xi, Y/xi, B_2d/B_2d.max(), levels=20, cmap='hot')
ax4.contour(X/xi, Y/xi, B_2d/B_2d.max(), levels=[0.5, 0.1, 0.01],
            colors='white', alpha=0.5, linewidths=1)

# Add circles for ξ and λ
circle_xi = plt.Circle((0, 0), 1, fill=False, edgecolor='cyan',
                       linewidth=2, linestyle='--', label='r = ξ')
circle_lambda = plt.Circle((0, 0), kappa, fill=False, edgecolor='blue',
                           linewidth=2, linestyle='--', label='r = λ')
ax4.add_patch(circle_xi)
ax4.add_patch(circle_lambda)

ax4.set_xlabel('x/ξ', fontsize=12)
ax4.set_ylabel('y/ξ', fontsize=12)
ax4.set_title('2D Magnetic Field Map', fontsize=14)
ax4.axis('equal')
ax4.legend(fontsize=11, loc='upper right')
cbar = plt.colorbar(im, ax=ax4)
cbar.set_label('B(r) / B(0)', fontsize=11)

plt.tight_layout()
plt.savefig('vortex_structure.png', dpi=150, bbox_inches='tight')
plt.show()

# Calculate integrated flux
flux_numerical = 2 * np.pi * np.trapz(r * B, r)
print(f"Vortex core size: ξ = {xi} nm")
print(f"Field penetration depth: λ = {lambda_L} nm")
print(f"GL parameter: κ = λ/ξ = {kappa}")
print(f"\nIntegrated flux (numerical): {flux_numerical:.3f} × (Φ₀/2πλ²)")
print(f"Expected: 1.0 (one flux quantum)")
print(f"\nKey length scales:")
print(f"  - Order parameter varies over ~ξ = {xi} nm")
print(f"  - Magnetic field decays over ~λ = {lambda_L} nm")
print(f"  - Supercurrent peaks at r ≈ {r[peak_idx]:.1f} nm")

2.5 Hexagonal Vortex Lattice

Lattice Formation

When multiple vortices are present, they arrange themselves in a periodic lattice to minimize energy. Abrikosov showed that for κ >> 1, the hexagonal (triangular) lattice has the lowest energy:

Lattice Spacing:

$$a \approx 1.075 \sqrt{\frac{\Phi_0}{B}}$$

where B is the average magnetic field. As B increases (H → H_c2), the vortices pack closer together.

import numpy as np
import matplotlib.pyplot as plt

def hexagonal_lattice(B, Phi0, Lx, Ly):
    """
    Generate hexagonal vortex lattice positions

    Parameters:
    B: average magnetic field
    Phi0: flux quantum
    Lx, Ly: system size

    Returns:
    vortex_positions: array of (x, y) coordinates
    """
    # Lattice constant
    a = 1.075 * np.sqrt(Phi0 / B)

    # Basis vectors for hexagonal lattice
    a1 = np.array([a, 0])
    a2 = np.array([a/2, a*np.sqrt(3)/2])

    # Number of unit cells
    n1 = int(Lx / a) + 2
    n2 = int(Ly / (a*np.sqrt(3)/2)) + 2

    positions = []
    for i in range(-n1, n1):
        for j in range(-n2, n2):
            pos = i * a1 + j * a2
            if -Lx/2 < pos[0] < Lx/2 and -Ly/2 < pos[1] < Ly/2:
                positions.append(pos)

    return np.array(positions)

def calculate_field_map(vortex_positions, X, Y, xi, lambda_L):
    """
    Calculate total magnetic field from all vortices
    """
    from scipy.special import k0

    B_total = np.zeros_like(X)

    for vx, vy in vortex_positions:
        R = np.sqrt((X - vx)**2 + (Y - vy)**2)
        R_safe = np.where(R < 0.01 * lambda_L, 0.01 * lambda_L, R)
        B_total += k0(R_safe / lambda_L)

    return B_total

# Physical parameters
Phi0 = 2.067e-15  # Wb (flux quantum)
xi = 10  # nm
kappa = 5
lambda_L = kappa * xi
Lx, Ly = 500, 500  # System size in nm

# Different field values
B_values = [0.01, 0.05, 0.2]  # Normalized to Φ₀/λ²
field_labels = ['B ≈ 0.05 Hc2', 'B ≈ 0.3 Hc2', 'B ≈ 0.8 Hc2']

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

for ax, B_norm, label in zip(axes, B_values, field_labels):
    # Generate vortex lattice
    B_physical = B_norm * Phi0 / lambda_L**2  # Approximate scaling
    positions = hexagonal_lattice(B_physical, Phi0, Lx, Ly)

    # Plot vortices
    ax.scatter(positions[:, 0], positions[:, 1], c='red', s=100,
               marker='o', edgecolors='darkred', linewidths=2,
               label=f'{len(positions)} vortices', zorder=5)

    # Calculate lattice constant
    a = 1.075 * np.sqrt(Phi0 / B_physical)

    # Draw hexagonal unit cell for first few vortices
    if len(positions) > 0:
        # Find neighbors
        center = positions[len(positions)//2]
        distances = np.sqrt(np.sum((positions - center)**2, axis=1))
        neighbors_idx = np.where((distances > 0.1) & (distances < 1.5*a))[0]

        # Draw lines to nearest neighbors
        for idx in neighbors_idx[:6]:
            ax.plot([center[0], positions[idx, 0]],
                   [center[1], positions[idx, 1]],
                   'b--', alpha=0.5, linewidth=1.5)

    # Add circles showing vortex core sizes
    for vx, vy in positions[:min(20, len(positions))]:
        circle = plt.Circle((vx, vy), xi, fill=False, edgecolor='cyan',
                           linewidth=1, alpha=0.5)
        ax.add_patch(circle)

    ax.set_xlabel('x (nm)', fontsize=12)
    ax.set_ylabel('y (nm)', fontsize=12)
    ax.set_title(f'{label}\na ≈ {a:.1f} nm', fontsize=13)
    ax.axis('equal')
    ax.set_xlim(-Lx/2, Lx/2)
    ax.set_ylim(-Ly/2, Ly/2)
    ax.grid(True, alpha=0.3)
    ax.legend(fontsize=10, loc='upper right')

    # Add text showing spacing
    ax.text(0.05, 0.95, f'Spacing: {a:.1f} nm\nCores: ξ = {xi} nm',
            transform=ax.transAxes, va='top', ha='left',
            fontsize=10, bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

plt.tight_layout()
plt.savefig('vortex_lattice.png', dpi=150, bbox_inches='tight')
plt.show()

print("Vortex Lattice Properties:")
for B_norm, label in zip(B_values, field_labels):
    B_physical = B_norm * Phi0 / lambda_L**2
    a = 1.075 * np.sqrt(Phi0 / B_physical)
    n_vortices = (Lx * Ly) * B_physical / Phi0  # Approximate
    print(f"\n{label}:")
    print(f"  Lattice constant: a = {a:.1f} nm")
    print(f"  Vortex density: {B_physical/Phi0:.2e} vortices/nm²")
    print(f"  Core separation: a/ξ = {a/xi:.1f}")

2.6 Vortex Dynamics

Lorentz Force on Vortices

When a current density J flows through a superconductor containing vortices, each vortex experiences a Lorentz force:

$$\mathbf{F}_L = \mathbf{J} \times \mathbf{\Phi}_0$$

where Φ₀ is the flux vector (magnitude Φ₀, direction parallel to B).

Per unit length of vortex:

$$f_L = J \Phi_0$$

Flux Flow Resistance

When vortices move with velocity v, they generate an electric field:

$$\mathbf{E} = \mathbf{B} \times \mathbf{v}$$

This leads to energy dissipation (resistance) even in the "superconducting" state.

In steady state, the Lorentz force is balanced by a viscous drag force:

$$f_L = \eta v$$

where η is the viscous drag coefficient. This gives:

$$v = \frac{J\Phi_0}{\eta}$$

The resulting electric field is:

$$E = Bv = \frac{B\Phi_0}{\eta}J = \rho_{ff} J$$

where ρ_ff = BΦ₀/η is the flux flow resistivity.

Flux Creep and Thermal Activation

At finite temperature, thermal fluctuations can cause vortices to hop between pinning sites even when J < J_c. This is called flux creep.

Anderson-Kim Model:

$$v_{creep} = v_0 \exp\left(-\frac{U_0}{k_B T}\left[1 - \frac{J}{J_c}\right]\right)$$

where U₀ is the pinning barrier height.

2.7 Flux Pinning and Critical Current

Pinning Mechanisms

Vortex motion (and thus energy dissipation) can be prevented by flux pinning. Defects in the material create local variations in the free energy that trap vortices:

Point defects: Vacancies, interstitials, impurity atoms
Line defects: Dislocations, columnar defects from irradiation
Planar defects: Grain boundaries, stacking faults, twin boundaries
Volume defects: Precipitates, second-phase particles, voids

Pinning Energy:

A vortex core in a region with suppressed T_c or enhanced ξ gains energy:

$$\Delta E_{pin} \sim \xi^2 \Delta f$$

where Δf is the change in condensation energy density.

Critical Current Density J_c

The critical current density is the maximum current that can flow without vortex motion:

$$J_c = \frac{F_{pin}}{\Phi_0}$$

where F_pin is the maximum pinning force per unit length of vortex.

Bean Critical State Model

The Bean model provides a simple description of current distribution in pinned superconductors:

Bean Model Assumptions:

J = ±J_c wherever vortices are present
J = 0 in field-free regions
Abrupt current reversal at boundaries

For a slab of thickness 2a in perpendicular field:

$$B(x) = B_a - \mu_0 J_c |x|$$

for |x| < a, where B_a is the field at the surface.

import numpy as np
import matplotlib.pyplot as plt

def bean_model_slab(x, Ba, Jc, mu0=4*np.pi*1e-7):
    """
    Bean critical state model for infinite slab

    Parameters:
    x: position array (distance from center)
    Ba: applied field at surface
    Jc: critical current density
    mu0: permeability

    Returns:
    B: magnetic field profile
    J: current density profile
    """
    # Field penetration depth
    d_pen = Ba / (mu0 * Jc)

    # Field profile
    B = np.zeros_like(x)
    mask = np.abs(x) < d_pen
    B[mask] = Ba - mu0 * Jc * np.abs(x[mask])

    # Current density
    J = np.zeros_like(x)
    J[mask] = Jc * np.sign(x[mask])

    return B, J

# Parameters
a = 1.0  # Half-thickness (mm)
Jc = 1e9  # A/m² (typical for NbTi)
mu0 = 4*np.pi*1e-7

# Different applied fields
Ba_values = [0.1, 0.5, 1.0, 2.0]  # Tesla
colors = ['blue', 'green', 'orange', 'red']

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

# Loop over applied fields
for idx, (Ba, color) in enumerate(zip(Ba_values, colors)):
    row = idx // 2
    col = idx % 2
    ax = axes[row, col]

    x = np.linspace(-a, a, 1000)  # mm
    B, J = bean_model_slab(x, Ba, Jc, mu0)

    # Plot field profile
    ax_twin = ax.twinx()

    l1 = ax.plot(x, B, color=color, linewidth=2.5, label=f'B(x)')
    l2 = ax_twin.plot(x, J/1e9, color='darkred', linewidth=2.5,
                      linestyle='--', label='J(x)')

    # Penetration depth
    d_pen = Ba / (mu0 * Jc)

    # Mark penetration depth
    if d_pen < a:
        ax.axvline(x=-d_pen, color='gray', linestyle=':', alpha=0.7)
        ax.axvline(x=d_pen, color='gray', linestyle=':', alpha=0.7)
        ax.fill_between([-d_pen, d_pen], 0, [B[np.argmin(np.abs(x + d_pen))],
                        B[np.argmin(np.abs(x - d_pen))]], alpha=0.1, color=color)

        ax.text(0, Ba/2, f'Field-free\ncore\n(J = 0)',
                ha='center', va='center', fontsize=10,
                bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
    else:
        ax.fill_between(x, 0, B, alpha=0.1, color=color)
        ax.text(0, Ba/2, f'Full penetration\n(J = ±Jc)',
                ha='center', va='center', fontsize=10,
                bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

    ax.set_xlabel('Position x (mm)', fontsize=12)
    ax.set_ylabel('Magnetic field B (T)', fontsize=12, color=color)
    ax_twin.set_ylabel('Current density J (GA/m²)', fontsize=12, color='darkred')
    ax.tick_params(axis='y', labelcolor=color)
    ax_twin.tick_params(axis='y', labelcolor='darkred')

    ax.set_title(f'Applied field Ba = {Ba} T\nPenetration depth = {d_pen*1000:.2f} mm',
                fontsize=13)
    ax.grid(True, alpha=0.3)
    ax.set_xlim(-a, a)
    ax.set_ylim(0, Ba*1.1)
    ax_twin.set_ylim(-Jc/1e9*1.2, Jc/1e9*1.2)

    # Combined legend
    lns = l1 + l2
    labs = [l.get_label() for l in lns]
    ax.legend(lns, labs, loc='upper right', fontsize=10)

plt.tight_layout()
plt.savefig('bean_model.png', dpi=150, bbox_inches='tight')
plt.show()

print("Bean Critical State Model:")
print(f"Critical current density: Jc = {Jc:.2e} A/m²")
print("\nField penetration depths:")
for Ba in Ba_values:
    d_pen = Ba / (mu0 * Jc)
    status = "Partial" if d_pen < a else "Full"
    print(f"  Ba = {Ba} T: d = {d_pen*1000:.2f} mm ({status} penetration)")

Magnetization Curves

The Bean model predicts hysteretic magnetization curves M(H) with characteristic features:

import numpy as np
import matplotlib.pyplot as plt

def bean_magnetization_slab(Ha_array, Jc, a, mu0=4*np.pi*1e-7):
    """
    Calculate magnetization for Bean model in slab geometry

    Parameters:
    Ha_array: applied field values
    Jc: critical current density
    a: half-thickness of slab
    mu0: permeability

    Returns:
    M: magnetization (A/m)
    """
    M = np.zeros_like(Ha_array)

    for i, Ha in enumerate(Ha_array):
        # Convert H to B at surface
        Ba = mu0 * Ha

        # Penetration depth
        d_pen = min(Ba / (mu0 * Jc), a)

        # Average field inside
        if d_pen < a:
            # Partial penetration
            B_avg = Ba * d_pen / (2*a)
        else:
            # Full penetration
            B_avg = Ba - mu0 * Jc * a / 2

        # Magnetization M = B/μ₀ - H
        M[i] = B_avg / mu0 - Ha

    return M

# Parameters
a = 1e-3  # 1 mm half-thickness
Jc = 1e9  # A/m²
mu0 = 4*np.pi*1e-7

# Applied field cycle: 0 → +Hmax → -Hmax → +Hmax
Hmax = 2e6  # A/m (about 2.5 T)
H_up1 = np.linspace(0, Hmax, 200)
H_down = np.linspace(Hmax, -Hmax, 400)
H_up2 = np.linspace(-Hmax, Hmax, 400)

# Calculate magnetization for each branch
M_up1 = bean_magnetization_slab(H_up1, Jc, a, mu0)
M_down = bean_magnetization_slab(np.abs(H_down), Jc, a, mu0) * np.sign(H_down)
M_up2 = bean_magnetization_slab(np.abs(H_up2), Jc, a, mu0) * np.sign(H_up2)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# Plot 1: M-H curve (hysteresis loop)
ax1.plot(H_up1/1e6, M_up1/1e6, 'b-', linewidth=2.5, label='Virgin curve')
ax1.plot(H_down/1e6, M_down/1e6, 'r-', linewidth=2.5, label='Decreasing H')
ax1.plot(H_up2/1e6, M_up2/1e6, 'g-', linewidth=2.5, label='Increasing H')

# Mark key points
ax1.plot([0], [0], 'ko', markersize=8)
ax1.annotate('Start', xy=(0, 0), xytext=(0.2, -0.2),
            arrowprops=dict(arrowstyle='->', color='black'),
            fontsize=11)

# Calculate width of hysteresis
H_mid = 0
idx_mid_down = np.argmin(np.abs(H_down - H_mid))
idx_mid_up = np.argmin(np.abs(H_up2 - H_mid))
delta_M = M_down[idx_mid_down] - M_up2[idx_mid_up]

ax1.plot([H_mid/1e6, H_mid/1e6],
         [M_up2[idx_mid_up]/1e6, M_down[idx_mid_down]/1e6],
         'k--', linewidth=1.5, alpha=0.7)
ax1.text(H_mid/1e6 + 0.1, (M_down[idx_mid_down] + M_up2[idx_mid_up])/2/1e6,
         f'ΔM = {delta_M/1e6:.2f} MA/m',
         fontsize=10, bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.7))

ax1.set_xlabel('Applied field H (MA/m)', fontsize=12)
ax1.set_ylabel('Magnetization M (MA/m)', fontsize=12)
ax1.set_title('Bean Model: Magnetization Hysteresis', fontsize=14)
ax1.legend(fontsize=11, loc='lower right')
ax1.grid(True, alpha=0.3)
ax1.axhline(y=0, color='k', linewidth=0.5)
ax1.axvline(x=0, color='k', linewidth=0.5)

# Plot 2: Different Jc values
Jc_values = [1e8, 5e8, 1e9, 2e9]
colors_jc = ['blue', 'green', 'orange', 'red']

for Jc_val, color in zip(Jc_values, colors_jc):
    M_cycle = bean_magnetization_slab(H_up1, Jc_val, a, mu0)
    ax2.plot(H_up1/1e6, M_cycle/1e6, color=color, linewidth=2.5,
            label=f'Jc = {Jc_val:.1e} A/m²')

ax2.set_xlabel('Applied field H (MA/m)', fontsize=12)
ax2.set_ylabel('Magnetization M (MA/m)', fontsize=12)
ax2.set_title('Effect of Critical Current on Magnetization', fontsize=14)
ax2.legend(fontsize=11, loc='lower right')
ax2.grid(True, alpha=0.3)
ax2.axhline(y=0, color='k', linewidth=0.5)
ax2.axvline(x=0, color='k', linewidth=0.5)

plt.tight_layout()
plt.savefig('magnetization_hysteresis.png', dpi=150, bbox_inches='tight')
plt.show()

print("Magnetization Hysteresis Analysis:")
print(f"Sample thickness: 2a = {2*a*1000} mm")
print(f"Critical current density: Jc = {Jc:.2e} A/m²")
print(f"\nHysteresis width at H = 0:")
print(f"  ΔM = {delta_M/1e6:.2f} MA/m")
print(f"  ΔM ∝ Jc × a (Bean model)")
print(f"\nEnergy loss per cycle ∝ ∫ M dH (area of loop)")

Summary

Key Takeaways

Type II superconductors (κ > 1/√2) have negative surface energy, leading to spontaneous vortex formation
Three critical fields: H_c1 (vortex entry), H_c (thermodynamic), H_c2 (superconductivity destruction)
Mixed state exists between H_c1 and H_c2 with quantized flux vortices
Flux quantization Φ = nΦ₀ arises from single-valuedness of the order parameter wavefunction
Abrikosov vortices have normal cores (~ξ) and circulating supercurrents (~λ)
Hexagonal lattice minimizes energy for vortex arrays
Vortex motion under current creates resistance (flux flow) unless prevented by pinning
Critical current J_c is determined by pinning force from defects
Bean model describes current distribution and magnetization in hard superconductors

Practice Problems

Problem 1: Critical Field Calculation

For Nb₃Sn: ξ(0) = 3 nm, λ(0) = 65 nm, T_c = 18 K.

(a) Calculate κ, H_c1(0), and H_c2(0).

(b) If H_c(0) = 0.4 T, verify the relationships between the critical fields.

Problem 2: Flux Quantization Verification

A superconducting ring has inner diameter 1 mm and traps magnetic flux.

(a) If the measured flux is 3.1 × 10⁻¹⁵ Wb, how many flux quanta are trapped?

(b) Calculate the minimum detectable field change ΔB if one additional flux quantum enters.

Problem 3: Vortex Lattice Spacing

A Type II superconductor with λ = 100 nm is in a field B = 1 T.

(a) Calculate the vortex lattice constant a.

(b) How many vortices are present per μm²?

Problem 4: Bean Model Application

A superconducting slab (thickness 2 mm, J_c = 10⁹ A/m²) is in field B_a = 0.5 T.

(a) Calculate the penetration depth.

(b) Sketch the B(x) and J(x) profiles.

Problem 5: Critical Current from Pinning

A vortex experiences pinning force F_pin = 10⁻¹¹ N per unit length.

(a) Calculate the critical current density J_c.

(b) If the pinning centers have density n_p = 10²¹ m⁻³, what is the average pinning energy per defect?