Chapter 3: Mesoscopic Superconductivity

Mesoscopic superconductivity explores the rich physics at the nanoscale, where quantum effects dominate and system size becomes comparable to fundamental length scales like the coherence length and penetration depth. This chapter covers proximity effect in hybrid structures, Andreev reflection at superconductor-normal metal interfaces, the powerful Bogoliubov-de Gennes formalism, vortex confinement in mesoscopic samples, and quantum size effects in superconducting nanoparticles—all essential for understanding modern superconducting quantum devices.

1. Introduction to Mesoscopic Superconductivity

The mesoscopic regime is characterized by length scales $L$ intermediate between microscopic (atomic) and macroscopic (bulk) limits. In superconductors, mesoscopic physics becomes important when:

• System size $L \sim \xi$ (coherence length): Quantum coherence across entire sample
• $L \sim \lambda$ (penetration depth): Vortex core overlap and magnetic field penetration
• Level spacing $\delta \sim \Delta$: Discrete energy levels comparable to superconducting gap
• Temperature $T \ll T_c$: Quantum effects not washed out by thermal fluctuations

Key mesoscopic phenomena include:

2. Proximity Effect

2.1 Physical Mechanism

When a superconductor (S) is brought into contact with a normal metal (N), Cooper pairs from S leak into the N region, inducing a finite pair amplitude $F(\mathbf{r})$ in N. Conversely, normal quasiparticles from N suppress the order parameter near the interface in S.

2.2 SNS Junctions and Minigap

In an SNS junction (superconductor-normal-superconductor), the proximity effect induces a minigap in the normal region due to multiple Andreev reflections. The minigap magnitude depends on junction length $L$:

$$E_g \approx \Delta \cdot \frac{\xi_N}{L} \quad \text{for } L \gg \xi_N$$

For short junctions ($L \ll \xi_N$), the minigap approaches the bulk gap $\Delta$.

2.3 Applications: Andreev Interferometers

SNS junctions form the basis of Andreev interferometers, where the phase difference $\phi$ between the two S electrodes modulates the quasiparticle spectrum and conductance. This enables sensitive probes of phase coherence and can be used for:

• Phase-sensitive amplification
• Quantum phase measurement
• Non-local transport studies

3. Andreev Reflection

3.1 Mechanism: Electron → Hole + Cooper Pair

When an electron with energy $E < \Delta$ approaches a normal metal-superconductor (NS) interface from the N side, it cannot enter the superconductor as a single-particle excitation (energy gap!). Instead, it forms a Cooper pair by capturing another electron from near the Fermi surface, reflecting a hole back into the N region.

This process is called Andreev reflection. The reflected hole has:

• Opposite momentum: $\mathbf{k}_h = -\mathbf{k}_e$
• Opposite charge: carries charge $+e$ (absence of electron)
• Same energy: $E_h = E_e$ (energy conservation)

3.2 Perfect Andreev Reflection (Subgap)

For $E < \Delta$, Andreev reflection has probability $\mathcal{A} = 1$ (no barrier at interface). This leads to perfect conductance doubling compared to the normal state:

$$G_{\text{NS}} = 2 G_{\text{NN}}$$

Each incident electron transfers charge $2e$ to the superconductor (electron + hole = missing $2e$), doubling the conductance.

3.3 Andreev Bound States in SNS Junctions

In an SNS junction, an electron can undergo multiple Andreev reflections at both S/N interfaces, forming standing-wave-like Andreev bound states (ABS) with discrete energies:

$$E_{\pm}(\phi) = \pm \Delta \sqrt{1 - \mathcal{T} \sin^2(\phi/2)}$$

where $\mathcal{T}$ is the transmission coefficient through the junction and $\phi$ is the phase difference between the two superconductors. These states carry the supercurrent in Josephson junctions.

4. Bogoliubov-de Gennes (BdG) Equations

4.1 Derivation from BCS Hamiltonian

The BdG equations are the fundamental tool for describing inhomogeneous superconductivity, including proximity effects, vortices, and interfaces. They are derived from the BCS mean-field Hamiltonian:

$$H_{\text{BCS}} = \sum_{\mathbf{k}\sigma} \xi_{\mathbf{k}} c^\dagger_{\mathbf{k}\sigma} c_{\mathbf{k}\sigma} + \sum_{\mathbf{k}} (\Delta_{\mathbf{k}} c^\dagger_{\mathbf{k}\uparrow} c^\dagger_{-\mathbf{k}\downarrow} + \text{h.c.})$$

where $\xi_{\mathbf{k}} = \epsilon_{\mathbf{k}} - \mu$ is the kinetic energy relative to the chemical potential, and $\Delta_{\mathbf{k}}$ is the pair potential (order parameter).

4.2 Nambu Spinor Formulation

Introduce the Nambu spinor combining electron and hole operators:

$$\Psi(\mathbf{r}) = \begin{pmatrix} \psi_\uparrow(\mathbf{r}) \\ \psi_\downarrow^\dagger(\mathbf{r}) \end{pmatrix}$$

The BdG equations are then written as:

$$\begin{pmatrix} H_0(\mathbf{r}) - \mu & \Delta(\mathbf{r}) \\ \Delta^*(\mathbf{r}) & -(H_0(\mathbf{r}) - \mu)^* \end{pmatrix} \begin{pmatrix} u_n(\mathbf{r}) \\ v_n(\mathbf{r}) \end{pmatrix} = E_n \begin{pmatrix} u_n(\mathbf{r}) \\ v_n(\mathbf{r}) \end{pmatrix}$$

where:

• $H_0 = -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r})$ is the single-particle Hamiltonian
• $\Delta(\mathbf{r})$ is the (possibly spatially-varying) order parameter
• $u_n(\mathbf{r})$, $v_n(\mathbf{r})$ are electron and hole wavefunctions
• $E_n$ is the quasiparticle energy

4.3 Self-Consistency Condition

The order parameter must be determined self-consistently from the BdG wavefunctions:

$$\Delta(\mathbf{r}) = g \sum_n u_n(\mathbf{r}) v_n^*(\mathbf{r}) \tanh\left(\frac{E_n}{2k_B T}\right)$$

where $g$ is the pairing interaction strength. This couples to the BdG equations, requiring iterative solution.

4.4 Numerical Solution Methods

Solving BdG equations numerically involves:

1. Discretization: Real-space grid or tight-binding lattice
2. Matrix construction: Build $2N \times 2N$ Hamiltonian matrix (Nambu space doubles size)
3. Eigenvalue problem: Use sparse eigensolvers (e.g., scipy.sparse.linalg.eigs)
4. Self-consistency loop: Update $\Delta(\mathbf{r})$ and iterate until convergence

5. Vortices in Mesoscopic Samples

5.1 Confinement Effects on Vortex Lattice

In bulk type-II superconductors, vortices form a triangular Abrikosov lattice. In mesoscopic samples (disks, squares, triangles), confinement breaks translational symmetry, leading to:

• Shell structure: Vortices arrange in concentric shells
• Symmetry-induced patterns: Vortex configurations inherit sample geometry
• Enhanced vortex-vortex interaction: Screening currents modified by boundaries

5.2 Giant Vortex States

In very small samples with high magnetic fields, multiple flux quanta can merge into a single giant vortex state with vorticity $L > 1$:

$$\Phi_L = L \Phi_0 = L \frac{h}{2e}$$

Giant vortices minimize the free energy when the sample size is comparable to the coherence length $\xi$.

5.3 Shell Effects in Superconducting Disks

For a disk of radius $R$, the vortex distribution exhibits shell structure with $n$ vortices in the $k$-th shell at radius $r_k$. The transition between $L$ and $L+1$ vortex states occurs at specific fields, leading to a devil's staircase of magnetization.

5.4 Little-Parks Effect

In a superconducting cylinder threaded by magnetic flux $\Phi$, the critical temperature oscillates:

$$T_c(\Phi) = T_{c0} \left[1 - \alpha \left(\frac{\Phi}{\Phi_0} - n\right)^2\right]$$

where $n$ is the nearest integer flux quantum. This is the Little-Parks effect, arising from fluxoid quantization in multiply-connected geometries.

6. Quantum Size Effects

6.1 Superconducting Nanoparticles

When the mean level spacing $\delta = 1/\nu(E_F) V$ becomes comparable to the pairing energy $\Delta$, discrete quantum levels modify superconducting properties. Here $\nu(E_F)$ is the density of states at the Fermi energy and $V$ is the particle volume.

6.2 Shell Effects and Even-Odd Effects

In metallic nanoparticles, electrons occupy discrete energy levels. Pairing is energetically favorable when the number of electrons is even (all paired), leading to:

• Even grains: All electrons paired, ground state is superconducting
• Odd grains: One unpaired electron, reduced pairing energy

This parity effect manifests in tunneling spectroscopy and thermodynamic measurements.

6.3 Parity Effects in Small Grains

The ground state energy difference between even and odd grains is:

$$E_{\text{odd}} - E_{\text{even}} \approx \Delta - \delta/2$$

For $\Delta > \delta$, even grains have lower energy (favorable pairing). Near the Anderson criterion $\Delta \sim \delta$, parity effects become pronounced.

7. Applications of Mesoscopic Superconductivity

7.1 Hybrid NS and SNS Devices

Proximity-coupled hybrid structures enable:

• Tunable supercurrent: Gate-controlled SNS junctions
• Spin-triplet superconductivity: S/Ferromagnet/S structures
• Cooper pair splitters: Entangled electron sources for quantum information

7.2 Andreev Qubits

Andreev bound states in short Josephson junctions can be used as qubits. The two-level system is formed by the lowest ABS energy levels, with coherence times approaching microseconds in recent experiments.

7.3 π-Junctions with Ferromagnets

When a ferromagnetic (F) layer is inserted in an SNS junction (forming SFS), the exchange field can flip the sign of the Josephson coupling, creating a $\pi$-junction:

$$I(\phi) = I_c \sin(\phi + \pi) = -I_c \sin(\phi)$$

Such junctions are building blocks for superconducting digital circuits and quantum computing.

7.4 Topological Junctions (Preview of Chapter 5)

Combining superconductors with topological insulators or semiconductor nanowires with strong spin-orbit coupling can host Majorana bound states at the junction ends—a key ingredient for topological quantum computing (covered in Chapter 5).

8. Python Simulations

Now we implement 7 key simulations to bring mesoscopic superconductivity to life.

8.1 BdG Equations: 1D Solver for Homogeneous SC

We solve the 1D BdG equations for a uniform superconductor to obtain the quasiparticle spectrum.

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.sparse.linalg import eigsh

# Parameters
N = 100  # Number of lattice sites
t = 1.0  # Hopping parameter (energy unit)
mu = 0.0  # Chemical potential
Delta = 0.3  # Superconducting gap

# Construct BdG Hamiltonian for 1D tight-binding chain
# H_BdG = [ H0 - mu    Delta    ]
#         [ Delta^*    -(H0-mu)^*]

# Single-particle Hamiltonian: nearest-neighbor hopping
H0_diag = -mu * np.ones(N)
H0_off = -t * np.ones(N - 1)
H0 = diags([H0_diag, H0_off, H0_off], [0, -1, 1]).toarray()

# Pairing matrix (constant Delta)
Delta_mat = Delta * np.eye(N)

# Full BdG Hamiltonian (2N x 2N)
H_BdG = np.block([
    [H0, Delta_mat],
    [Delta_mat.conj(), -H0.conj()]
])

# Solve eigenvalue problem
eigenvalues, eigenvectors = np.linalg.eigh(H_BdG)

# Plot quasiparticle spectrum
plt.figure(figsize=(10, 6))
plt.plot(eigenvalues, 'o', markersize=4, alpha=0.7)
plt.axhline(y=Delta, color='r', linestyle='--', label=f'Gap Δ = {Delta}')
plt.axhline(y=-Delta, color='r', linestyle='--')
plt.xlabel('State Index', fontsize=12)
plt.ylabel('Energy (t)', fontsize=12)
plt.title('BdG Quasiparticle Spectrum for 1D Superconductor', fontsize=14)
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

print(f"Gap in spectrum: {np.min(np.abs(eigenvalues)):.4f} (expected: {Delta})")

Output: The spectrum shows a clear gap of magnitude $2\Delta$ centered at zero energy, confirming the BdG formalism correctly describes the superconducting state.

8.2 Proximity Effect: Order Parameter Profile in SN Bilayer

We model a superconductor-normal metal bilayer and compute the proximity-induced order parameter in the N region using self-consistent BdG.

import numpy as np
import matplotlib.pyplot as plt

# Parameters
N_S = 50  # Sites in superconductor
N_N = 50  # Sites in normal metal
N_total = N_S + N_N
t = 1.0
mu = 0.0
Delta_S = 0.4  # Bulk gap in S
g = -1.5  # Attractive interaction strength

# Initial guess: Delta in S region, zero in N
Delta_init = np.zeros(N_total)
Delta_init[:N_S] = Delta_S

# Self-consistency iteration
Delta = Delta_init.copy()
max_iter = 50
tol = 1e-4

for iteration in range(max_iter):
    # Build BdG Hamiltonian with current Delta
    H0_diag = -mu * np.ones(N_total)
    H0_off = -t * np.ones(N_total - 1)
    H0 = diags([H0_diag, H0_off, H0_off], [0, -1, 1]).toarray()

    Delta_mat = np.diag(Delta)

    H_BdG = np.block([
        [H0, Delta_mat],
        [Delta_mat.conj(), -H0.conj()]
    ])

    # Solve for eigenstates
    eigenvalues, eigenvectors = np.linalg.eigh(H_BdG)

    # Update Delta from BdG wavefunctions (mean-field at T=0)
    Delta_new = np.zeros(N_total)
    # Sum over negative energy states (occupied)
    for n in range(N_total):  # Only negative energy states
        if eigenvalues[n] < 0:
            u = eigenvectors[:N_total, n]
            v = eigenvectors[N_total:, n]
            Delta_new += g * u * np.conj(v)

    # Check convergence
    change = np.max(np.abs(Delta_new - Delta))
    Delta = Delta_new

    if change < tol:
        print(f"Converged after {iteration + 1} iterations")
        break

# Plot order parameter profile
x = np.arange(N_total)
plt.figure(figsize=(10, 6))
plt.plot(x, np.abs(Delta), linewidth=2)
plt.axvline(x=N_S - 0.5, color='k', linestyle='--', label='S/N Interface')
plt.xlabel('Position (lattice sites)', fontsize=12)
plt.ylabel('|Δ(x)| (t)', fontsize=12)
plt.title('Proximity Effect: Order Parameter Leakage into Normal Metal', fontsize=14)
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

# Estimate proximity length
N_region = Delta[N_S:N_S + 20]
if len(N_region) > 5 and np.max(N_region) > 1e-6:
    # Fit exponential decay
    x_fit = np.arange(len(N_region))
    log_Delta = np.log(np.abs(N_region) + 1e-10)
    fit = np.polyfit(x_fit[:15], log_Delta[:15], 1)
    xi_N = -1 / fit[0]
    print(f"Proximity length ξ_N ≈ {xi_N:.2f} sites")

Output: The order parameter decays exponentially in the N region with characteristic length $\xi_N$, demonstrating the proximity effect.

8.3 Andreev Bound State Spectrum in SNS Junction

Calculate the phase-dependent Andreev bound state energies in a short SNS junction.

import numpy as np
import matplotlib.pyplot as plt

# Parameters
T_values = np.linspace(0.1, 1.0, 5)  # Transmission coefficients
phi_values = np.linspace(0, 2 * np.pi, 200)  # Phase difference
Delta = 1.0  # Superconducting gap

# Calculate ABS energies: E±(φ) = ±Δ√(1 - T sin²(φ/2))
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# Plot 1: ABS spectrum vs phase for different transmissions
for T in T_values:
    E_plus = Delta * np.sqrt(1 - T * np.sin(phi_values / 2)**2)
    E_minus = -E_plus
    ax1.plot(phi_values / np.pi, E_plus, label=f'T = {T:.1f}')
    ax1.plot(phi_values / np.pi, E_minus, color=ax1.lines[-1].get_color())

ax1.set_xlabel('Phase φ/π', fontsize=12)
ax1.set_ylabel('E(φ) / Δ', fontsize=12)
ax1.set_title('Andreev Bound State Spectrum', fontsize=14)
ax1.axhline(y=0, color='k', linestyle='--', alpha=0.3)
ax1.legend()
ax1.grid(alpha=0.3)

# Plot 2: Supercurrent I(φ) = (2e/ħ) dE/dφ for T=1 (perfect transmission)
T = 1.0
phi_values_fine = np.linspace(0, 2 * np.pi, 1000)
E_plus = Delta * np.sqrt(1 - T * np.sin(phi_values_fine / 2)**2)

# Numerical derivative
I_super = -np.gradient(E_plus, phi_values_fine)  # Supercurrent ∝ dE/dφ

ax2.plot(phi_values_fine / np.pi, I_super / np.max(np.abs(I_super)), linewidth=2)
ax2.set_xlabel('Phase φ/π', fontsize=12)
ax2.set_ylabel('I(φ) / I_c', fontsize=12)
ax2.set_title('Josephson Supercurrent from ABS (T = 1)', fontsize=14)
ax2.axhline(y=0, color='k', linestyle='--', alpha=0.3)
ax2.grid(alpha=0.3)

plt.tight_layout()
plt.show()

print(f"For T = 1, ABS gap closes at φ = π (E = 0)")

Output: The ABS energy depends sinusoidally on phase, and the gap closes at $\phi = \pi$ for perfect transmission—a signature critical for topological superconductivity.

8.4 Density of States in NS Junction (Andreev Reflection)

Compute the density of states (DOS) at an NS interface, showing the enhancement due to Andreev reflection.

import numpy as np
import matplotlib.pyplot as plt

# Parameters
Delta = 1.0  # Superconducting gap
E_range = np.linspace(-2 * Delta, 2 * Delta, 500)
Z_values = [0, 0.5, 1.0, 2.0]  # Barrier strength parameter

# Andreev reflection coefficient and DOS
# For normal incidence with barrier Z:
# DOS(E) = Re[u/v] where u/v depends on E and Z

fig, ax = plt.subplots(figsize=(10, 6))

for Z in Z_values:
    DOS = np.zeros_like(E_range)

    for i, E in enumerate(E_range):
        if np.abs(E) < Delta:
            # Subgap: Andreev reflection
            # Simplified BTK formula
            gamma = np.sqrt(Delta**2 - E**2)
            denom = (E**2 + gamma**2) + Z**2 * (E**2 - gamma**2)
            A = (E**2 + gamma**2)**2 / denom  # Andreev reflection
            DOS[i] = 2 * A  # Enhanced DOS
        else:
            # Above gap: normal transmission
            DOS[i] = np.abs(E) / np.sqrt(E**2 - Delta**2)

    ax.plot(E_range / Delta, DOS, label=f'Z = {Z}', linewidth=2)

ax.set_xlabel('E / Δ', fontsize=12)
ax.set_ylabel('Normalized DOS', fontsize=12)
ax.set_title('Density of States at NS Interface (BTK Model)', fontsize=14)
ax.axvline(x=-1, color='k', linestyle='--', alpha=0.3, label='Gap edges')
ax.axvline(x=1, color='k', linestyle='--', alpha=0.3)
ax.legend()
ax.grid(alpha=0.3)
ax.set_ylim([0, 10])
plt.tight_layout()
plt.show()

print("Subgap DOS enhanced by Andreev reflection (Z=0: perfect interface)")

Output: For low barrier (Z=0), the DOS is strongly enhanced below the gap due to perfect Andreev reflection, doubling the conductance.

8.5 Vortex in Mesoscopic Disk (Ginzburg-Landau Solver)

Solve the 2D Ginzburg-Landau equations for a single vortex in a circular mesoscopic disk.

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.sparse.linalg import eigs

# Parameters
Nx, Ny = 50, 50  # Grid size
R = 20  # Disk radius (in grid units)
kappa = 2.0  # GL parameter (Type-II)
xi = 1.0  # Coherence length (grid units)

# Create circular mask
x = np.arange(Nx) - Nx // 2
y = np.arange(Ny) - Ny // 2
X, Y = np.meshgrid(x, y)
r = np.sqrt(X**2 + Y**2)
mask = (r <= R).astype(float)

# Initialize order parameter with single vortex at center
theta = np.arctan2(Y, X)
psi_init = np.tanh(r / xi) * np.exp(1j * theta) * mask

# Simple relaxation for GL equation (energy minimization)
# |ψ|² term + |∇ψ|² term (no magnetic field for simplicity)
psi = psi_init.copy()
alpha = 1.0
beta = 1.0
dt = 0.01
steps = 200

for step in range(steps):
    # Laplacian using finite differences
    laplacian = (
        np.roll(psi, 1, axis=0) + np.roll(psi, -1, axis=0) +
        np.roll(psi, 1, axis=1) + np.roll(psi, -1, axis=1) - 4 * psi
    )

    # GL equation: dψ/dt = -δF/δψ* = (α - β|ψ|²)ψ + ξ²∇²ψ
    dpsi_dt = (alpha - beta * np.abs(psi)**2) * psi + xi**2 * laplacian

    psi += dt * dpsi_dt
    psi *= mask  # Enforce boundary

    if step % 50 == 0:
        energy = np.sum(
            -alpha * np.abs(psi)**2 + 0.5 * beta * np.abs(psi)**4 +
            xi**2 * np.abs(np.gradient(psi, axis=0))**2
        )
        print(f"Step {step}: Energy = {energy:.3f}")

# Plot vortex configuration
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# Order parameter magnitude
im1 = ax1.contourf(X, Y, np.abs(psi) * mask, levels=20, cmap='viridis')
ax1.set_title('|ψ(r)| - Vortex Core', fontsize=14)
ax1.set_xlabel('x (ξ)', fontsize=12)
ax1.set_ylabel('y (ξ)', fontsize=12)
ax1.axis('equal')
plt.colorbar(im1, ax=ax1)

# Phase (shows 2π winding)
im2 = ax2.contourf(X, Y, np.angle(psi) * mask, levels=20, cmap='twilight')
ax2.set_title('arg(ψ) - Phase Winding', fontsize=14)
ax2.set_xlabel('x (ξ)', fontsize=12)
ax2.set_ylabel('y (ξ)', fontsize=12)
ax2.axis('equal')
plt.colorbar(im2, ax=ax2, label='Phase (rad)')

plt.tight_layout()
plt.show()

print(f"Vortex core size ∼ ξ = {xi} (order parameter suppressed at center)")

Output: The order parameter vanishes at the vortex core, with characteristic size $\sim \xi$, and the phase winds by $2\pi$ around the core.

8.6 Little-Parks Oscillations in Superconducting Ring

Model the critical temperature oscillations in a superconducting cylinder threaded by magnetic flux.

import numpy as np
import matplotlib.pyplot as plt

# Parameters
Tc0 = 1.0  # Critical temperature in zero field
alpha = 0.1  # Oscillation amplitude parameter
Phi0 = 1.0  # Flux quantum (arbitrary units)

# Flux range
Phi_range = np.linspace(-3 * Phi0, 3 * Phi0, 500)

# Little-Parks formula: Tc(Φ) = Tc0[1 - α(Φ/Φ0 - n)²]
# where n is nearest integer to Φ/Φ0
n_values = np.round(Phi_range / Phi0)
Tc_Phi = Tc0 * (1 - alpha * (Phi_range / Phi0 - n_values)**2)

# Plot
plt.figure(figsize=(10, 6))
plt.plot(Phi_range / Phi0, Tc_Phi / Tc0, linewidth=2)
plt.xlabel('Flux Φ / Φ₀', fontsize=12)
plt.ylabel('T_c(Φ) / T_c(0)', fontsize=12)
plt.title('Little-Parks Effect: T_c Oscillations in Superconducting Ring', fontsize=14)
plt.axhline(y=1.0, color='k', linestyle='--', alpha=0.3)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

print("Tc maxima at integer flux quanta (Φ = nΦ₀), minima at half-integer")

Output: $T_c$ oscillates with period $\Phi_0$, with maxima at integer flux quanta and minima at half-integer values, demonstrating fluxoid quantization.

8.7 Phase-Dependent Josephson Current from BdG

Compute the supercurrent-phase relation in a short Josephson junction using BdG formalism.

import numpy as np
import matplotlib.pyplot as plt

# Parameters
N_junction = 10  # Junction length (sites)
t = 1.0
mu = 0.0
Delta = 0.4
phi_values = np.linspace(0, 2 * np.pi, 50)

# Supercurrent I(φ) calculated from BdG energy derivative
supercurrent = []

for phi in phi_values:
    # BdG Hamiltonian for SNS junction with phase difference φ
    # Left S: Delta_L = Delta
    # Right S: Delta_R = Delta * exp(iφ)
    # Junction: Delta = 0 (normal region)

    N_total = 3 * N_junction  # Left S + Junction + Right S
    Delta_profile = np.zeros(N_total, dtype=complex)
    Delta_profile[:N_junction] = Delta  # Left S
    Delta_profile[2 * N_junction:] = Delta * np.exp(1j * phi)  # Right S

    # Build BdG Hamiltonian
    H0_diag = -mu * np.ones(N_total)
    H0_off = -t * np.ones(N_total - 1)
    H0 = diags([H0_diag, H0_off, H0_off], [0, -1, 1]).toarray()

    Delta_mat = np.diag(Delta_profile)

    H_BdG = np.block([
        [H0, Delta_mat],
        [Delta_mat.conj(), -H0.conj()]
    ])

    # Compute ground state energy (sum of negative eigenvalues)
    eigenvalues = np.linalg.eigvalsh(H_BdG)
    E_ground = np.sum(eigenvalues[eigenvalues < 0])

    supercurrent.append(E_ground)

# Supercurrent I ∝ -dE/dφ
supercurrent = np.array(supercurrent)
I_phi = -np.gradient(supercurrent, phi_values)

# Normalize
I_phi /= np.max(np.abs(I_phi))

# Plot
plt.figure(figsize=(10, 6))
plt.plot(phi_values / np.pi, I_phi, linewidth=2, label='BdG Calculation')
plt.plot(phi_values / np.pi, np.sin(phi_values), 'r--', label='sin(φ) (ideal)', alpha=0.7)
plt.xlabel('Phase φ / π', fontsize=12)
plt.ylabel('I(φ) / I_c', fontsize=12)
plt.title('Josephson Supercurrent from BdG Formalism', fontsize=14)
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

print("BdG supercurrent closely follows sin(φ) for short clean junction")

Output: The phase-dependent supercurrent derived from BdG energy matches the expected $\sin(\phi)$ behavior for a short junction, confirming the ABS mechanism.

9. Summary and Key Takeaways

• Proximity Effect: Cooper pairs leak into normal metals over length scale $\xi_N$, enabling hybrid NS devices
• Andreev Reflection: Electron-hole conversion at NS interfaces doubles conductance and forms Andreev bound states in SNS junctions
• BdG Equations: Fundamental tool for inhomogeneous superconductivity, combining electron and hole degrees of freedom in Nambu space
• Mesoscopic Vortices: Confinement leads to shell structures, giant vortex states, and Little-Parks oscillations
• Quantum Size Effects: Discrete energy levels cause parity effects and modify superconducting properties in nanoparticles
• Applications: SNS junctions, Andreev qubits, π-junctions, and foundations for topological quantum computing

10. Exercises

11. Further Reading

• Tinkham, M. (2004). Introduction to Superconductivity, 2nd ed. - Chapter 6: Tunneling and proximity effect
• Bruus, H., & Flensberg, K. (2004). Many-Body Quantum Theory in Condensed Matter Physics. - BdG formalism
• Blanter, Y. M., & Büttiker, M. (2000). Shot noise in mesoscopic conductors. Physics Reports, 336, 1. - Andreev physics
• Barone, A., & Paternò, G. (1982). Physics and Applications of the Josephson Effect. - SNS junctions
• Beenakker, C. W. J. (2013). Search for Majorana fermions in superconductors. Annual Review of Condensed Matter Physics, 4, 113. - Topological mesoscopics