Chapter 3: Josephson Effects

The Josephson effects, discovered by Brian Josephson in 1962 (Nobel Prize 1973), describe quantum tunneling of Cooper pairs through a weak link between two superconductors. This phenomenon leads to remarkable effects: supercurrent flow without voltage (DC Josephson effect) and oscillating supercurrent under applied voltage (AC Josephson effect). In this chapter, we rigorously derive these effects, explore the RCSJ model for realistic junctions, analyze phase dynamics, and investigate SQUID devices for ultra-sensitive magnetometry.

3.1 Josephson Junction Types and Structure

What is a Josephson Junction?

A Josephson junction (JJ) consists of two superconducting electrodes separated by a weak link through which Cooper pairs can tunnel quantum mechanically. The weak link can take several forms:

Tunnel Junction (SIS)

The most common type consists of two superconductors (S) separated by a thin insulating barrier (I):

• Structure: Superconductor / Insulator / Superconductor (SIS)
• Barrier thickness: 1-2 nm (typically Al₂O₃ or MgO)
• Tunneling mechanism: Quantum tunneling of Cooper pairs through potential barrier
• Critical current: Exponentially depends on barrier thickness

Weak Link Junction

A narrow superconducting constriction with width comparable to coherence length $\xi$:

• Cooper pairs flow through the constriction with reduced order parameter
• Critical current depends on constriction geometry
• Used in variable-thickness bridge junctions

SNS Junction

Normal metal (N) sandwiched between two superconductors:

• Cooper pairs diffuse through the normal metal via proximity effect
• Length scale: normal metal coherence length $\xi_N = \sqrt{\hbar D / k_B T}$
• Used in hybrid devices and quantum computing

3.2 DC Josephson Effect: Supercurrent Without Voltage

The Josephson Current-Phase Relation

The fundamental result of Josephson's theory is that a supercurrent flows through the junction even in the absence of voltage, with magnitude determined by the phase difference $\phi$ between the two superconductors:

where:

• $I$ is the supercurrent through the junction
• $I_c$ is the critical current (maximum supercurrent)
• $\phi = \phi_2 - \phi_1$ is the phase difference across the junction

Derivation from Coupled Ginzburg-Landau Equations

Consider two superconductors (labeled 1 and 2) with order parameters:

The weak coupling allows Cooper pairs to tunnel between the superconductors. The tunneling Hamiltonian couples the two order parameters:

where $T$ is the tunneling matrix element.

The time evolution of the order parameters is governed by the time-dependent Ginzburg-Landau equation. For the current density through the junction:

Integrating over the junction area $A$:

Critical Current Ic

For a tunnel junction with barrier thickness $d$ and area $A$, the Ambegaokar-Baratoff formula gives:

where:

• $\Delta(T)$ is the superconducting gap at temperature $T$
• $R_N$ is the normal-state resistance of the junction
• $e$ is the elementary charge

At $T = 0$:

Energy of the Junction

The energy stored in the junction due to the phase difference is:

where the Josephson coupling energy is:

Here $\Phi_0 = h/2e$ is the flux quantum.

Key Insight: The DC Josephson effect allows dissipationless current flow up to $I_c$. This is the basis for superconducting quantum interference devices (SQUIDs), Josephson voltage standards, and superconducting qubits.

import numpy as np
import matplotlib.pyplot as plt
import scipy.constants as const

# Physical constants
e = const.e          # Elementary charge (C)
h = const.h          # Planck constant (J·s)
hbar = const.hbar    # Reduced Planck constant (J·s)
Phi_0 = h / (2*e)    # Flux quantum (Wb)

# Junction parameters
I_c = 10e-6  # Critical current (A) - 10 μA typical for small junction
R_N = 10     # Normal resistance (Ω)

# Calculate Josephson energy
E_J = (Phi_0 * I_c) / (2*np.pi)

# Phase difference array
phi = np.linspace(-2*np.pi, 2*np.pi, 500)

# DC Josephson current
I = I_c * np.sin(phi)

# Josephson energy
U = -E_J * np.cos(phi)

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

# Current-phase relation
ax1.plot(phi, I*1e6, 'b-', linewidth=2.5, label='I(φ) = Ic sin(φ)')
ax1.axhline(I_c*1e6, color='red', linestyle='--', linewidth=2, alpha=0.7, label=f'Ic = {I_c*1e6:.0f} μA')
ax1.axhline(-I_c*1e6, color='red', linestyle='--', linewidth=2, alpha=0.7)
ax1.axhline(0, color='black', linestyle='-', linewidth=0.8, alpha=0.5)
ax1.axvline(0, color='black', linestyle='-', linewidth=0.8, alpha=0.5)

ax1.set_xlabel('Phase difference φ (rad)', fontsize=13, fontweight='bold')
ax1.set_ylabel('Current I (μA)', fontsize=13, fontweight='bold')
ax1.set_title('DC Josephson Effect: Current-Phase Relation', fontsize=14, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.legend(fontsize=11)
ax1.set_xticks([-2*np.pi, -np.pi, 0, np.pi, 2*np.pi])
ax1.set_xticklabels(['-2π', '-π', '0', 'π', '2π'])

# Energy-phase relation
ax2.plot(phi, U/const.e*1e6, 'r-', linewidth=2.5, label='U(φ) = -EJ cos(φ)')
ax2.axhline(0, color='black', linestyle='-', linewidth=0.8, alpha=0.5)
ax2.axvline(0, color='black', linestyle='-', linewidth=0.8, alpha=0.5)

# Mark minima and maxima
minima_phi = np.array([0])
maxima_phi = np.array([-np.pi, np.pi])
ax2.plot(minima_phi, -E_J/const.e*1e6*np.ones_like(minima_phi), 'go',
         markersize=10, label='Energy minimum (stable)')
ax2.plot(maxima_phi, E_J/const.e*1e6*np.ones_like(maxima_phi), 'ro',
         markersize=10, label='Energy maximum (unstable)')

ax2.set_xlabel('Phase difference φ (rad)', fontsize=13, fontweight='bold')
ax2.set_ylabel('Energy U (μeV)', fontsize=13, fontweight='bold')
ax2.set_title('Josephson Energy Landscape', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.legend(fontsize=11)
ax2.set_xticks([-2*np.pi, -np.pi, 0, np.pi, 2*np.pi])
ax2.set_xticklabels(['-2π', '-π', '0', 'π', '2π'])

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

print(f"Junction Parameters:")
print(f"  Critical current Ic = {I_c*1e6:.2f} μA")
print(f"  Normal resistance RN = {R_N:.1f} Ω")
print(f"  Josephson energy EJ = {E_J/const.e*1e6:.3f} μeV")
print(f"  IcRN product = {I_c*R_N*1e6:.2f} μV")
print(f"  Characteristic frequency = EJ/ℏ = {E_J/hbar*1e-9:.2f} GHz")

3.3 AC Josephson Effect: Voltage-to-Frequency Conversion

The Voltage-Frequency Relation

When a DC voltage $V$ is applied across the junction, the phase difference evolves in time according to:

This is the second Josephson equation, describing the AC Josephson effect.

Derivation from Energy Considerations

The Josephson energy difference between the two sides is related to the work done by the voltage:

The energy is also related to the phase:

However, from quantum mechanics, the phase evolution relates to energy as:

Josephson Frequency

Integrating the voltage-phase relation with constant voltage $V$:

where the Josephson frequency is:

Numerically:

Oscillating Supercurrent

Substituting the time-dependent phase into the DC Josephson relation:

The current oscillates at the Josephson frequency. For typical voltages ($V \sim 1$ mV), $f_J \sim 500$ GHz (microwave frequency).

Metrological Application: The AC Josephson effect provides the most accurate voltage standard. The frequency $f_J$ can be measured with atomic clock precision, allowing voltage to be defined through the relation $V = (h/2e) f_J$. This defines the volt in the SI system.

Shapiro Steps

When the junction is irradiated with microwave radiation at frequency $f_{RF}$, the I-V curve develops Shapiro steps (constant-voltage plateaus) at voltages:

These steps occur when the Josephson frequency synchronizes with the applied RF frequency: $f_J = n f_{RF}$.

import numpy as np
import matplotlib.pyplot as plt
import scipy.constants as const

# Physical constants
e = const.e
h = const.h
hbar = const.hbar

# Voltage-to-frequency conversion factor
K_J = 2*e / h  # Josephson constant (Hz/V)

# Voltage range
V = np.linspace(0, 2e-3, 100)  # 0 to 2 mV

# Josephson frequency
f_J = K_J * V

# Example: Fixed voltage for time-dependent current
V_fixed = 1e-3  # 1 mV
omega_J = 2*e*V_fixed / hbar
f_J_fixed = omega_J / (2*np.pi)

I_c = 10e-6  # 10 μA
phi_0 = 0

# Time array
t = np.linspace(0, 10e-12, 1000)  # 10 ps

# Oscillating current
phi_t = phi_0 + omega_J * t
I_t = I_c * np.sin(phi_t)

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

# Voltage-frequency relation
ax1.plot(V*1e3, f_J*1e-9, 'b-', linewidth=2.5)
ax1.set_xlabel('Voltage V (mV)', fontsize=13, fontweight='bold')
ax1.set_ylabel('Josephson Frequency fJ (GHz)', fontsize=13, fontweight='bold')
ax1.set_title('AC Josephson Effect: Voltage-Frequency Relation', fontsize=14, fontweight='bold')
ax1.grid(True, alpha=0.3)

# Add annotation for 483.6 GHz/mV slope
ax1.text(0.5, 400, f'fJ = {K_J*1e-9:.1f} GHz/mV × V\n= 483.6 GHz/mV × V',
         fontsize=11, bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))

# Oscillating current
ax2.plot(t*1e12, I_t*1e6, 'r-', linewidth=2)
ax2.set_xlabel('Time t (ps)', fontsize=13, fontweight='bold')
ax2.set_ylabel('Current I (μA)', fontsize=13, fontweight='bold')
ax2.set_title(f'Oscillating Supercurrent at V = {V_fixed*1e3:.1f} mV\nfJ = {f_J_fixed*1e-9:.1f} GHz',
             fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

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

print(f"AC Josephson Effect:")
print(f"  Josephson constant KJ = 2e/h = {K_J*1e-9:.2f} GHz/mV")
print(f"  For V = 1 mV:")
print(f"    Josephson frequency fJ = {K_J*1e-3*1e-9:.2f} GHz")
print(f"    Period T = {1/(K_J*1e-3)*1e12:.3f} ps")
print(f"\n  Voltage standard:")
print(f"    1 μV corresponds to fJ = {K_J*1e-6*1e-6:.3f} MHz")

3.4 RCSJ Model: Realistic Junction Dynamics

The Resistively and Capacitively Shunted Junction Model

Real Josephson junctions have not only the ideal Josephson element but also:

• Parallel resistance $R$ (due to quasiparticle tunneling)
• Parallel capacitance $C$ (due to junction geometry)

The RCSJ model represents the junction as a parallel combination of these elements:

Circuit Equation Derivation

Applying Kirchhoff's current law at the junction node:

Using the AC Josephson relation $V = (\hbar/2e)(d\phi/dt)$:

Defining characteristic time and damping:

• Plasma frequency: $\omega_p = \sqrt{2eI_c / (\hbar C)} = \sqrt{2\pi I_c / (\Phi_0 C)}$
• Quality factor: $Q = \omega_p R C$
• Stewart-McCumber parameter: $\beta_c = (2\pi/\Phi_0) I_c R^2 C = Q^2$

The normalized RCSJ equation becomes:

Or in dimensionless form with $\tau = \omega_p t$:

Washboard Potential Analogy

The RCSJ equation can be interpreted as a particle moving in a tilted washboard potential:

The tilt is proportional to the applied current $I$. For small $I < I_c$, the particle is trapped in a potential well (zero-voltage state). For $I > I_c$, the particle rolls down the potential (voltage state).

Stewart-McCumber Parameter βc

The parameter $\beta_c$ determines the junction dynamics:

• Underdamped ($\beta_c > 1$): Low resistance, high capacitance. Exhibits hysteresis in I-V curve.
• Overdamped ($\beta_c < 1$): High resistance, low capacitance. No hysteresis, smooth I-V curve.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import scipy.constants as const

# Physical constants
e = const.e
h = const.h
hbar = const.hbar
Phi_0 = h / (2*e)

# Junction parameters
I_c = 10e-6    # Critical current (A)
R = 10         # Resistance (Ω)
C = 1e-12      # Capacitance (F) - 1 pF

# Calculate characteristic parameters
omega_p = np.sqrt(2*np.pi*I_c / (Phi_0*C))
Q = omega_p * R * C
beta_c = (2*np.pi / Phi_0) * I_c * R**2 * C

print(f"Junction Parameters:")
print(f"  Ic = {I_c*1e6:.1f} μA")
print(f"  R = {R:.1f} Ω")
print(f"  C = {C*1e12:.1f} pF")
print(f"  Plasma frequency ωp = {omega_p*1e-9:.2f} GHz")
print(f"  Quality factor Q = {Q:.2f}")
print(f"  Stewart-McCumber parameter βc = {beta_c:.2f}")
print(f"  Junction type: {'Underdamped (hysteretic)' if beta_c > 1 else 'Overdamped (non-hysteretic)'}")

# RCSJ differential equation
def rcsj_equation(t, y, I_applied):
    """
    y[0] = phi (phase)
    y[1] = dphi/dt (phase velocity)
    """
    phi, dphi_dt = y

    # d²φ/dt² = -ωp²sin(φ) - (1/Q)·(dφ/dt) + (ωp²/Ic)·I
    d2phi_dt2 = -omega_p**2 * np.sin(phi) - (omega_p/Q) * dphi_dt + (omega_p**2 / I_c) * I_applied

    return [dphi_dt, d2phi_dt2]

# Time array
t_max = 50e-9  # 50 ns
t_eval = np.linspace(0, t_max, 5000)

# Applied current (step function)
I_applied = 1.2 * I_c  # 120% of critical current

# Initial conditions
y0 = [0.0, 0.0]  # phi=0, dphi/dt=0

# Solve RCSJ equation
sol = solve_ivp(rcsj_equation, [0, t_max], y0, args=(I_applied,),
                t_eval=t_eval, method='RK45', rtol=1e-8, atol=1e-10)

# Extract results
phi = sol.y[0]
dphi_dt = sol.y[1]

# Calculate voltage from phase velocity
V = (hbar / (2*e)) * dphi_dt

# Plotting
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(12, 12))

# Phase evolution
ax1.plot(sol.t*1e9, phi, 'b-', linewidth=2)
ax1.set_xlabel('Time (ns)', fontsize=12, fontweight='bold')
ax1.set_ylabel('Phase φ (rad)', fontsize=12, fontweight='bold')
ax1.set_title(f'Phase Dynamics (I = {I_applied/I_c:.2f}Ic, βc = {beta_c:.2f})',
             fontsize=13, fontweight='bold')
ax1.grid(True, alpha=0.3)

# Phase velocity
ax2.plot(sol.t*1e9, dphi_dt*1e-9, 'r-', linewidth=2)
ax2.set_xlabel('Time (ns)', fontsize=12, fontweight='bold')
ax2.set_ylabel('dφ/dt (GHz)', fontsize=12, fontweight='bold')
ax2.set_title('Phase Velocity', fontsize=13, fontweight='bold')
ax2.grid(True, alpha=0.3)

# Voltage
ax3.plot(sol.t*1e9, V*1e6, 'g-', linewidth=2)
ax3.axhline(I_applied*R*1e6, color='red', linestyle='--', linewidth=2,
           label=f'V = IR = {I_applied*R*1e6:.2f} μV')
ax3.set_xlabel('Time (ns)', fontsize=12, fontweight='bold')
ax3.set_ylabel('Voltage V (μV)', fontsize=12, fontweight='bold')
ax3.set_title('Junction Voltage', fontsize=13, fontweight='bold')
ax3.grid(True, alpha=0.3)
ax3.legend(fontsize=10)

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

# Time-averaged voltage
V_avg = np.mean(V[-1000:])  # Average over last part
print(f"\nTime-averaged voltage V = {V_avg*1e6:.3f} μV")
print(f"Ohmic voltage IR = {I_applied*R*1e6:.3f} μV")

3.5 I-V Characteristics: Zero-Voltage and Voltage States

Zero-Voltage State

For $I < I_c$, there exists a static solution with constant phase $\phi_0$ satisfying:

The voltage is zero: $V = (\hbar/2e)(d\phi/dt) = 0$. This is the superconducting state.

Voltage State

For $I > I_c$, no static solution exists. The phase continuously increases, and the time-averaged voltage is:

This is the voltage state with dissipation.

Hysteresis in Underdamped Junctions

For $\beta_c > 1$, the I-V curve exhibits hysteresis:

• Increasing current: Junction switches to voltage state at $I = I_c$ (switch current)
• Decreasing current: Junction returns to zero-voltage state at $I = I_r < I_c$ (return current)

The return current is approximately:

Thermal Fluctuations and Phase Diffusion

At finite temperature, thermal fluctuations cause the phase to diffuse even for $I < I_c$. The characteristic energy for thermal escape from the zero-voltage state is:

The switching rate from thermal activation is:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import scipy.constants as const

# Physical constants
e = const.e
h = const.h
hbar = const.hbar
Phi_0 = h / (2*e)

# Junction parameters
I_c = 10e-6
R = 10

def calculate_iv_curve(C, I_range, label, color):
    """Calculate I-V curve for given capacitance"""
    omega_p = np.sqrt(2*np.pi*I_c / (Phi_0*C))
    Q = omega_p * R * C
    beta_c = (2*np.pi / Phi_0) * I_c * R**2 * C

    def rcsj_equation(t, y, I_applied):
        phi, dphi_dt = y
        d2phi_dt2 = -omega_p**2 * np.sin(phi) - (omega_p/Q) * dphi_dt + (omega_p**2 / I_c) * I_applied
        return [dphi_dt, d2phi_dt2]

    voltages = []

    for I_applied in I_range:
        t_max = 100e-9
        t_eval = np.linspace(0, t_max, 1000)
        y0 = [0.0, 0.0]

        sol = solve_ivp(rcsj_equation, [0, t_max], y0, args=(I_applied,),
                       t_eval=t_eval, method='RK45', rtol=1e-8)

        dphi_dt = sol.y[1]
        V = (hbar / (2*e)) * dphi_dt
        V_avg = np.mean(V[-200:])  # Time-average over last part
        voltages.append(V_avg)

    return np.array(voltages), beta_c

# Current range for sweep
I_range_up = np.linspace(0, 2*I_c, 50)

# Different capacitances (different βc values)
capacitances = [0.1e-12, 0.5e-12, 2e-12]  # pF
labels = []
colors = ['blue', 'green', 'red']

fig, ax = plt.subplots(figsize=(12, 8))

for C, color in zip(capacitances, colors):
    V, beta_c = calculate_iv_curve(C, I_range_up, f'C={C*1e12:.1f}pF', color)
    ax.plot(I_range_up*1e6, V*1e6, color=color, linewidth=2.5,
           label=f'βc = {beta_c:.2f} ({"underdamped" if beta_c > 1 else "overdamped"})')

# Theoretical zero-voltage line
ax.axvline(I_c*1e6, color='black', linestyle='--', linewidth=2,
          label=f'Ic = {I_c*1e6:.0f} μA')

ax.set_xlabel('Current I (μA)', fontsize=13, fontweight='bold')
ax.set_ylabel('Voltage V (μV)', fontsize=13, fontweight='bold')
ax.set_title('I-V Characteristics: Effect of Stewart-McCumber Parameter βc',
            fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
ax.legend(fontsize=11, loc='upper left')
ax.set_xlim(0, 20)
ax.set_ylim(0, None)

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

print("I-V Characteristics calculated for different βc values")
print("Underdamped junctions (βc > 1) would show hysteresis in full bidirectional sweep")

3.6 Shapiro Steps: Microwave-Induced Phenomena

Phase Locking to External RF

When the junction is irradiated with microwave radiation at frequency $f_{RF}$, the voltage across the junction can lock to discrete values:

These Shapiro steps appear as horizontal plateaus in the I-V curve.

Modified RCSJ Equation with RF Drive

The current becomes:

The RCSJ equation with RF drive is:

where $i_{DC} = I_{DC}/I_c$, $i_{RF} = I_{RF}/I_c$, and $\Omega = \omega_{RF}/\omega_p$.

Physical Mechanism

Shapiro steps occur when the Josephson oscillation frequency synchronizes with harmonics of the RF frequency:

This phase-locking mechanism is analogous to mode-locking in lasers.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import scipy.constants as const

# Physical constants
e = const.e
h = const.h
hbar = const.hbar
Phi_0 = h / (2*e)

# Junction parameters
I_c = 10e-6
R = 10
C = 0.5e-12

omega_p = np.sqrt(2*np.pi*I_c / (Phi_0*C))
Q = omega_p * R * C
beta_c = (2*np.pi / Phi_0) * I_c * R**2 * C

# RF parameters
f_RF = 10e9  # 10 GHz
omega_RF = 2*np.pi*f_RF
I_RF = 0.3 * I_c  # RF amplitude

print(f"Shapiro Step Simulation:")
print(f"  RF frequency fRF = {f_RF*1e-9:.1f} GHz")
print(f"  Expected step voltage Vn = n × {(h*f_RF)/(2*e)*1e6:.2f} μV")

# RCSJ with RF drive
def rcsj_rf_equation(t, y, I_DC):
    phi, dphi_dt = y
    I_applied = I_DC + I_RF * np.sin(omega_RF * t)
    d2phi_dt2 = -omega_p**2 * np.sin(phi) - (omega_p/Q) * dphi_dt + (omega_p**2 / I_c) * I_applied
    return [dphi_dt, d2phi_dt2]

# DC current sweep
I_DC_range = np.linspace(0, 2*I_c, 60)
voltages = []

for I_DC in I_DC_range:
    t_max = 200e-9
    t_eval = np.linspace(0, t_max, 2000)
    y0 = [0.0, 0.0]

    sol = solve_ivp(rcsj_rf_equation, [0, t_max], y0, args=(I_DC,),
                   t_eval=t_eval, method='RK45', rtol=1e-8)

    dphi_dt = sol.y[1]
    V = (hbar / (2*e)) * dphi_dt
    V_avg = np.mean(V[-500:])
    voltages.append(V_avg)

voltages = np.array(voltages)

# Calculate theoretical step voltages
step_voltage = (h * f_RF) / (2*e)
step_numbers = np.arange(-2, 5)
step_voltages = step_numbers * step_voltage

# Plotting
fig, ax = plt.subplots(figsize=(12, 8))

# I-V curve with Shapiro steps
ax.plot(I_DC_range*1e6, voltages*1e6, 'b-', linewidth=2.5, label='I-V with RF drive')

# Mark theoretical step positions
for n, V_step in zip(step_numbers, step_voltages):
    if V_step > 0 and V_step < max(voltages):
        ax.axhline(V_step*1e6, color='red', linestyle='--', linewidth=1.5, alpha=0.6)
        ax.text(1, V_step*1e6 + 0.5, f'n={n}', fontsize=10, color='red')

ax.set_xlabel('DC Current IDC (μA)', fontsize=13, fontweight='bold')
ax.set_ylabel('Time-Averaged Voltage V (μV)', fontsize=13, fontweight='bold')
ax.set_title(f'Shapiro Steps under Microwave Irradiation\nfRF = {f_RF*1e-9:.0f} GHz, IRF = {I_RF/I_c:.2f}Ic',
            fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
ax.legend(fontsize=11)

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

print(f"\nShapiro step positions (n × {step_voltage*1e6:.2f} μV):")
for n in range(0, 4):
    print(f"  n = {n}: V = {n*step_voltage*1e6:.2f} μV")

3.7 DC SQUID: Two-Junction Quantum Interferometer

DC SQUID Structure

A DC SQUID (Superconducting Quantum Interference Device) consists of a superconducting loop interrupted by two Josephson junctions in parallel.

Flux Quantization in SQUID Loop

The total flux through the loop is:

where $I_s$ is the screening current and $L$ is the loop inductance.

Flux quantization requires:

The phase drops across the two junctions satisfy:

Critical Current Modulation

For a symmetric SQUID (identical junctions with $I_{c1} = I_{c2} = I_c$), the total critical current is:

The critical current oscillates periodically with applied flux, with period $\Phi_0$.

Voltage Modulation

When biased with constant current $I_b > 2I_c$, the time-averaged voltage also oscillates:

SQUID Sensitivity

The flux sensitivity (minimum detectable flux) is limited by thermal and quantum noise:

State-of-the-art DC SQUIDs achieve flux sensitivities of $\sim 10^{-6} \Phi_0/\sqrt{\text{Hz}}$, making them the most sensitive magnetic field detectors.

import numpy as np
import matplotlib.pyplot as plt
import scipy.constants as const

# Physical constants
e = const.e
h = const.h
Phi_0 = h / (2*e)

# SQUID parameters
I_c = 5e-6     # Critical current per junction (A)
R = 10         # Resistance per junction (Ω)
L = 100e-12    # Loop inductance (H) - 100 pH

# Screening parameter
beta_L = 2*np.pi * L * I_c / Phi_0

print(f"DC SQUID Parameters:")
print(f"  Critical current per junction Ic = {I_c*1e6:.1f} μA")
print(f"  Resistance per junction R = {R:.1f} Ω")
print(f"  Loop inductance L = {L*1e12:.0f} pH")
print(f"  Screening parameter βL = {beta_L:.2f}")
print(f"  Flux quantum Φ₀ = {Phi_0*1e15:.3f} fWb")

# Applied flux range
Phi_ext = np.linspace(-2*Phi_0, 2*Phi_0, 500)
phi_ext_normalized = Phi_ext / Phi_0

# Critical current modulation (symmetric SQUID)
I_c_tot = 2 * I_c * np.abs(np.cos(np.pi * phi_ext_normalized))

# Bias current
I_b = 1.5 * 2 * I_c  # 150% of max critical current

# Voltage modulation (simplified model)
# V = R * sqrt(Ib² - Ic_tot²) for Ib > Ic_tot
V = np.zeros_like(I_c_tot)
for i, Ic_t in enumerate(I_c_tot):
    if I_b > Ic_t:
        V[i] = R * np.sqrt(I_b**2 - Ic_t**2)

# Flux-to-voltage transfer function (sensitivity)
dV_dPhi = np.gradient(V, Phi_ext)

# Plotting
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(12, 14))

# Critical current modulation
ax1.plot(phi_ext_normalized, I_c_tot*1e6, 'b-', linewidth=2.5)
ax1.axhline(2*I_c*1e6, color='red', linestyle='--', linewidth=1.5,
           label=f'Max Ic = {2*I_c*1e6:.1f} μA')
ax1.axhline(0, color='black', linestyle='-', linewidth=0.8, alpha=0.5)
ax1.set_xlabel('Applied Flux Φext / Φ₀', fontsize=13, fontweight='bold')
ax1.set_ylabel('Critical Current Ic(Φ) (μA)', fontsize=13, fontweight='bold')
ax1.set_title('DC SQUID Critical Current Modulation', fontsize=14, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.legend(fontsize=11)

# Voltage modulation
ax2.plot(phi_ext_normalized, V*1e6, 'r-', linewidth=2.5)
ax2.set_xlabel('Applied Flux Φext / Φ₀', fontsize=13, fontweight='bold')
ax2.set_ylabel('Voltage V(Φ) (μV)', fontsize=13, fontweight='bold')
ax2.set_title(f'DC SQUID Voltage Modulation (Ib = {I_b*1e6:.1f} μA)',
             fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

# Flux-to-voltage transfer (sensitivity)
ax3.plot(phi_ext_normalized, dV_dPhi*1e-9, 'g-', linewidth=2.5)
ax3.axhline(0, color='black', linestyle='-', linewidth=0.8, alpha=0.5)
ax3.set_xlabel('Applied Flux Φext / Φ₀', fontsize=13, fontweight='bold')
ax3.set_ylabel('dV/dΦ (V/Φ₀)', fontsize=13, fontweight='bold')
ax3.set_title('Flux-to-Voltage Transfer Function (Sensitivity)', fontsize=14, fontweight='bold')
ax3.grid(True, alpha=0.3)

# Mark maximum sensitivity points
max_sensitivity = np.max(np.abs(dV_dPhi))
max_idx = np.argmax(np.abs(dV_dPhi))
ax3.plot(phi_ext_normalized[max_idx], dV_dPhi[max_idx]*1e-9, 'ro', markersize=10,
        label=f'Max |dV/dΦ| = {max_sensitivity*1e-9:.2e} V/Φ₀')
ax3.legend(fontsize=11)

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

print(f"\nSQUID Sensitivity:")
print(f"  Maximum |dV/dΦ| = {max_sensitivity*1e-9:.2e} V/Φ₀")
print(f"  At flux values: Φ/Φ₀ = ±0.5, ±1.5, ...")

3.8 RF SQUID: Single-Junction Interferometer

RF SQUID Structure

An RF SQUID consists of a superconducting loop with a single Josephson junction, inductively coupled to an external RF tank circuit.

Operating Principle

The RF tank circuit at frequency $f_{RF}$ (typically 20-30 MHz) drives an oscillating current through the SQUID loop. The energy dissipated in the junction depends on the applied flux through:

The junction switches between superconducting and resistive states, modulating the tank circuit Q-factor and thus the RF amplitude or phase.

Flux Quantization and Potential Wells

The effective potential energy is:

This creates a tilted washboard potential. The number of flux quanta trapped in the loop changes as $\Phi_{ext}$ varies, causing stepwise changes in the tank circuit response.

Comparison: DC vs RF SQUID

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

# Physical constants
e = const.e
h = const.h
Phi_0 = h / (2*e)

# RF SQUID parameters
I_c = 10e-6    # Critical current (A)
L = 500e-12    # Loop inductance (H) - 500 pH

E_J = (Phi_0 * I_c) / (2*np.pi)
beta_L = 2*np.pi * L * I_c / Phi_0

print(f"RF SQUID Parameters:")
print(f"  Critical current Ic = {I_c*1e6:.1f} μA")
print(f"  Loop inductance L = {L*1e12:.0f} pH")
print(f"  Josephson energy EJ = {E_J/const.e*1e6:.2f} μeV")
print(f"  Screening parameter βL = {beta_L:.2f}")

# Potential energy function
def potential_energy(phi, Phi_ext):
    """
    U(φ, Φext) = (1/2L)(Φ₀φ/2π - Φext)² - EJ cos(φ)
    """
    flux_term = (Phi_0 * phi / (2*np.pi) - Phi_ext)**2 / (2*L)
    josephson_term = -E_J * np.cos(phi)
    return flux_term + josephson_term

# Phase and external flux arrays
phi = np.linspace(-2*np.pi, 4*np.pi, 300)
Phi_ext_normalized = np.linspace(-0.5, 1.5, 200)
Phi_ext = Phi_ext_normalized * Phi_0

# Create meshgrid
PHI, PHIEXT = np.meshgrid(phi, Phi_ext)

# Calculate potential energy landscape
U = np.zeros_like(PHI)
for i in range(len(Phi_ext_normalized)):
    for j in range(len(phi)):
        U[i, j] = potential_energy(PHI[i, j], PHIEXT[i, j])

# Normalize energy to Josephson energy
U_normalized = U / E_J

# Plotting
fig = plt.figure(figsize=(16, 6))

# 3D surface plot
ax1 = fig.add_subplot(121, projection='3d')
surf = ax1.plot_surface(PHI, PHIEXT/Phi_0, U_normalized, cmap='viridis',
                        alpha=0.9, linewidth=0, antialiased=True)

ax1.set_xlabel('Phase φ (rad)', fontsize=11, fontweight='bold')
ax1.set_ylabel('Φext / Φ₀', fontsize=11, fontweight='bold')
ax1.set_zlabel('Energy U / EJ', fontsize=11, fontweight='bold')
ax1.set_title('RF SQUID Potential Energy Landscape', fontsize=13, fontweight='bold')
fig.colorbar(surf, ax=ax1, shrink=0.5, aspect=10)

# 2D contour plot
ax2 = fig.add_subplot(122)
contour = ax2.contourf(PHI, PHIEXT/Phi_0, U_normalized, levels=30, cmap='viridis')
ax2.contour(PHI, PHIEXT/Phi_0, U_normalized, levels=15, colors='black',
           linewidths=0.5, alpha=0.3)

ax2.set_xlabel('Phase φ (rad)', fontsize=12, fontweight='bold')
ax2.set_ylabel('Φext / Φ₀', fontsize=12, fontweight='bold')
ax2.set_title('RF SQUID Potential Energy Contours', fontsize=13, fontweight='bold')
fig.colorbar(contour, ax=ax2, label='Energy U / EJ')

# Mark flux quantum transitions
for n in range(-1, 3):
    ax2.axhline(n, color='red', linestyle='--', linewidth=1.5, alpha=0.6)
    ax2.text(phi[-1] + 0.2, n, f'n={n}', fontsize=10, color='red')

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

print(f"\nFlux quantization states visible as potential wells at integer Φ/Φ₀")
print(f"βL = {beta_L:.2f} determines well depth and barrier height")

3.9 SQUID Noise and Sensitivity Analysis

Noise Sources in SQUIDs

The ultimate sensitivity of SQUIDs is limited by various noise sources:

• Thermal (Johnson) noise: From resistive elements at finite temperature
• Shot noise: Quantum fluctuations in current
• 1/f noise: Low-frequency noise from junction critical current fluctuations
• Environmental noise: Magnetic field fluctuations, vibrations

Flux Noise Spectral Density

The flux noise power spectral density is:

The energy resolution is:

Quantum Limit

The quantum limit for SQUID energy resolution is:

where $\omega_p$ is the plasma frequency. For typical SQUIDs, this is $\sim 10^{-32}$ J, far below the thermal noise limit at mK temperatures.

Practical Sensitivity

State-of-the-art DC SQUIDs achieve:

• Flux noise: $\sim 10^{-6} \Phi_0 / \sqrt{\text{Hz}}$
• Magnetic field sensitivity: $\sim 1$ fT$/\sqrt{\text{Hz}}$ (femtotesla)
• Energy resolution: $\sim 10^{-31}$ J (near quantum limit at mK)

import numpy as np
import matplotlib.pyplot as plt
import scipy.constants as const

# Physical constants
k_B = const.k
e = const.e
h = const.h
hbar = const.hbar
Phi_0 = h / (2*e)

# SQUID parameters
I_c = 5e-6      # Critical current (A)
R = 10          # Resistance (Ω)
L = 100e-12     # Inductance (H)
T = 4.2         # Temperature (K) - liquid He

# Calculate transfer function sensitivity (simplified)
# Maximum dV/dΦ occurs near Φ = (n+1/2)Φ₀
# From DC SQUID theory: max(dV/dΦ) ~ R·Ic / Φ₀
dV_dPhi_max = (R * I_c) / Phi_0

# Thermal flux noise spectral density
S_Phi_thermal = (k_B * T * L) / (dV_dPhi_max**2)

# Flux noise (white noise approximation)
flux_noise_density = np.sqrt(S_Phi_thermal) / Phi_0  # in units of Φ₀/√Hz

# Energy resolution
epsilon = S_Phi_thermal / (2*L)

# Frequency range for noise analysis
freq = np.logspace(-2, 6, 100)  # 0.01 Hz to 1 MHz

# 1/f noise model (simplified)
S_Phi_1f = S_Phi_thermal * (1 + 1/(freq/10))  # 1/f knee at 10 Hz

# Total noise
S_Phi_total = S_Phi_thermal + S_Phi_1f

print(f"SQUID Noise and Sensitivity Analysis:")
print(f"  Temperature T = {T:.1f} K")
print(f"  Loop inductance L = {L*1e12:.0f} pH")
print(f"  Transfer function dV/dΦ = {dV_dPhi_max:.2e} V/Wb")
print(f"\nNoise Performance:")
print(f"  Thermal flux noise = {flux_noise_density:.2e} Φ₀/√Hz")
print(f"  Flux noise in SI = {np.sqrt(S_Phi_thermal)*1e15:.2f} fWb/√Hz")
print(f"  Energy resolution ε = {epsilon/const.e*1e6:.2e} μeV")
print(f"  Energy resolution ε = {epsilon:.2e} J")

# Magnetic field sensitivity
# B = Φ / A_eff, assuming effective area A_eff ~ 1 mm²
A_eff = 1e-6  # m²
B_noise = np.sqrt(S_Phi_thermal) / A_eff
print(f"\nMagnetic Field Sensitivity (assuming A = 1 mm²):")
print(f"  Field noise = {B_noise*1e15:.2f} fT/√Hz")

# Quantum limit estimate
C_junction = 1e-12  # Typical junction capacitance (F)
omega_p = np.sqrt(2*np.pi*I_c / (Phi_0*C_junction))
epsilon_quantum = hbar * omega_p

print(f"\nQuantum Limit:")
print(f"  Plasma frequency ωp = {omega_p*1e-9:.2f} GHz")
print(f"  Quantum energy ℏωp = {epsilon_quantum/const.e*1e6:.2e} μeV")
print(f"  Ratio ε/ℏωp = {epsilon/epsilon_quantum:.1f}")

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

# Flux noise spectral density
ax1.loglog(freq, np.sqrt(S_Phi_total)/Phi_0, 'b-', linewidth=2.5, label='Total noise')
ax1.loglog(freq, np.sqrt(S_Phi_thermal)/Phi_0 * np.ones_like(freq), 'r--',
          linewidth=2, label='White (thermal) noise floor')
ax1.set_xlabel('Frequency (Hz)', fontsize=13, fontweight='bold')
ax1.set_ylabel('Flux Noise (Φ₀/√Hz)', fontsize=13, fontweight='bold')
ax1.set_title('SQUID Flux Noise Spectral Density', fontsize=14, fontweight='bold')
ax1.grid(True, alpha=0.3, which='both')
ax1.legend(fontsize=11)

# Energy resolution vs inductance
L_range = np.logspace(-13, -10, 100)  # 0.1 pH to 10 nH
epsilon_vs_L = (k_B * T) / (2 * L_range * dV_dPhi_max**2 / L**2)

ax2.loglog(L_range*1e12, epsilon_vs_L/const.e*1e6, 'g-', linewidth=2.5)
ax2.axhline(epsilon_quantum/const.e*1e6, color='red', linestyle='--', linewidth=2,
           label=f'Quantum limit = {epsilon_quantum/const.e*1e6:.2e} μeV')
ax2.axvline(L*1e12, color='blue', linestyle='--', linewidth=2, alpha=0.7,
           label=f'Current L = {L*1e12:.0f} pH')
ax2.set_xlabel('Loop Inductance L (pH)', fontsize=13, fontweight='bold')
ax2.set_ylabel('Energy Resolution ε (μeV)', fontsize=13, fontweight='bold')
ax2.set_title('Energy Resolution vs Loop Inductance', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3, which='both')
ax2.legend(fontsize=11)

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

print(f"\nOptimization: Smaller L improves energy resolution but reduces flux coupling")

Summary

In this chapter, we explored the Josephson effects and SQUID devices in depth:

• Introduced Josephson junction types: tunnel junctions, weak links, point contacts, and SNS junctions
• Derived the DC Josephson effect current-phase relation $I = I_c \sin(\phi)$ from coupled GL equations
• Derived the AC Josephson effect voltage-frequency relation $d\phi/dt = 2eV/\hbar$ and calculated Josephson frequency $f_J = 483.6$ GHz/mV
• Explored Shapiro steps under microwave irradiation and their metrological applications
• Derived the RCSJ model circuit equations including resistance and capacitance
• Analyzed the washboard potential and Stewart-McCumber parameter $\beta_c$ determining junction behavior
• Investigated I-V characteristics showing zero-voltage state, voltage state, and hysteresis
• Studied thermal fluctuations and phase diffusion effects
• Analyzed DC SQUID critical current and voltage modulation with applied flux
• Introduced RF SQUID single-junction interferometer and compared with DC SQUID
• Examined SQUID sensitivity including noise sources, flux noise spectral density, and quantum limits
• Implemented comprehensive Python simulations for all phenomena using scipy.integrate.solve_ivp

The Josephson effects form the foundation of modern quantum electronics, enabling ultra-sensitive magnetometry (SQUIDs), voltage standards, superconducting qubits for quantum computing, and rapid single-flux-quantum (RSFQ) logic circuits. The next chapter will explore BCS theory in depth to understand the microscopic origin of superconductivity.

Exercises