Superconducting Flux Qubit Josephson Junction Persistence Current Squid Quantum Coherence

# Superconducting Flux Qubits: Three-Junction SQUID Hamiltonians, Persistent Current States, and Tunneling Splitting

## 1. Introduction to Superconducting Flux Qubits

Flux qubits represent an alternative to charge and phase qubits, encoding quantum information in the persistent current state of a superconducting loop threaded by external magnetic flux. Unlike charge qubits (vulnerable to charge noise) and phase qubits (requiring large E_J/E_C ratios), flux qubits operate at an optimal "sweet spot" where first-order flux noise sensitivity vanishes (dE/dPhi = 0), enabling T₂ coherence times exceeding 50 microseconds—competitive with transmon qubits. The qubit consists of a superconducting loop interrupted by one to three Josephson junctions, creating a double-well potential with two degenerate or near-degenerate energy minima corresponding to clockwise and counterclockwise persistent currents. Tunneling through the barrier between these wells creates an energy gap proportional to the tunneling matrix element, defining the qubit transition frequency. This section develops the complete flux qubit Hamiltonian, derives the tunneling splitting, and explores the role of flux bias and junction asymmetry in optimizing coherence.

## 2. Superconducting Persistent Current and Flux Quantization

A superconducting loop threaded by magnetic flux Phi exhibits persistent current—a circulating current I_p that persists indefinitely (at zero temperature) without external driving. The persistent current is determined by flux quantization: the magnetic flux through a superconducting loop must be a multiple of the flux quantum:

where n is an integer, Phi_0 = h/(2e) approx 2.07 times 10^{-15} Wb is the flux quantum, and Phi_{text{ext}} is the applied external flux. The effective flux seen by the loop is the fractional part:

The persistent current magnitude is I_p = I_c |sin(phi/2)|, where I_c is the critical current of the junction(s). At phi = pi (half-flux-quantum bias), the persistent current is maximum: I_p = I_c. At phi = 0 or 2pi, the persistent current vanishes.

## 3. Single-Junction Flux Qubit and Effective Potential

For a simple superconducting loop interrupted by a single Josephson junction, the phase difference varphi across the junction relates to the persistent current via:

The Hamiltonian in terms of the gauge-invariant phase difference is:

where the charging energy term 4E_C n^2 is negligible for large loops (shunt capacitance dominated by stray capacitance, so E_C ll E_J), and the persistent-current term couples the phase to applied flux. The effective single-particle potential is:

This creates a tilted washboard: the tilt direction reverses at phi = pi (half-flux). At phi = pi, the potential has a deep symmetric double-well structure—ideal for qubits.

## 4. Three-Junction SQUID Loop (Flux Qubit)

Real flux qubits employ three-junction superconducting quantum interference device (3-junction SQUID) architectures, where three Josephson junctions are arranged in a triangular loop. Two junctions (with critical currents I_{c1} and I_{c2}, typically equal: I_{c1} = I_{c2} equiv I_c) form the "qubit" arm, while a third junction (with critical current I_{c3} = alpha I_c, where 0 < alpha < 1) serves as an asymmetry element. The effective potential energy is:

where f = Phi_{text{ext}}/Phi_0 is the normalized external flux (dimensionless, 0 to 1) and the phase variables varphi_1, varphi_2 satisfy the flux quantization constraint. The total flux through the loop is:

Near the half-flux-quantum point (f = 1/2), the potential develops a symmetric double-well structure, with minima at varphi_1 approx 0 (clockwise persistent current) and varphi_1 approx pi (counterclockwise persistent current).

## 5. Persistent Current States and Double-Well Quantization

The two potential minima at varphi_1 = 0 and varphi_1 = pi correspond to persistent current states |↺⟩ (clockwise) and |↻⟩ (counterclockwise) with opposite circulation directions. The persistent current magnitudes are:

At the minima (varphi_1 = 0, pi), the persistent current is either +I_p or -I_p, with magnitude typically I_p approx 0.1 text{ to } 1 microampere in practical devices. Quantum mechanically, the ground state and first excited state of this double-well potential are symmetric and antisymmetric combinations:

The energy difference (qubit transition frequency) arises from tunneling through the central barrier.

## 6. Tunneling Splitting and Flux-Tunable Frequency

The qubit Hamiltonian near the half-flux-quantum sweet spot can be projected to the two-level subspace spanned by |↺⟩ and |↻⟩:

where:
- varepsilon(Phi_{text{ext}}) = 2 I_p (Phi_{text{ext}} - frac{1}{2}Phi_0) is the flux-bias energy (tilt)
- Δ is the tunneling gap—the energy splitting at the symmetric point (varepsilon = 0)

The tunneling gap depends exponentially on the potential barrier height:

where S is the action integral through the classically forbidden region and A is a prefactor related to the attempt frequency. Typically, Delta/(2pi) approx 10-100 MHz for well-designed flux qubits.

At the optimal point where dvarepsilon/dPhi = 0 (half-flux quantum), the qubit is first-order insensitive to flux fluctuations. The transition frequency is simply:

at the sweet spot. Away from the sweet spot, omega_q increases linearly with flux bias: omega_q approx 2I_p DeltaPhi for small deviations DeltaPhi.

## 7. Flux Dephasing and Second-Order Noise Sensitivity

Although the flux qubit is first-order protected against flux noise (domega/dPhi = 0 at the sweet spot), second-order sensitivity remains:

Flux noise power spectral density S_Phi(f) (typically S_Phi^{1/f} sim 10^{-5} to 10^{-4} Phi_0^2/Hz) leads to dephasing:

At the sweet spot, the first-order term vanishes, reducing 1/T_2^* by a factor of (Delta/2I_pPhi_0)^2 approx 0.01, yielding T_2^* approx 10-100 microseconds compared to microseconds for charge qubits.

## 8. Junction Asymmetry and Design Optimization

The asymmetry parameter α (ratio of third junction critical current to pair critical current) profoundly affects qubit properties:

• Symmetric loop (alpha = 1): Potential has deeper wells but smaller tunneling gap Δ
• Asymmetric loop (alpha < 1): Reduced well depth but larger Δ, enabling faster gates

The optimal alpha approx 0.6-0.8 balances strong flux tunability with reasonable tunneling gap. Additionally, junction sizing asymmetry (slight differences in I_{c1} and I_{c2}) can create intentional approx 10% anticrossings, suppressing resonant two-photon coupling and improving isolation from unwanted transitions.

## 9. Readout via RF SQUID Coupling

Readout of flux qubits employs rf SQUID magnetometry: a separate rf-SQUID (superconducting quantum interference device operating at microwave frequencies) is capacitively coupled to the flux qubit. The rf SQUID's resonant frequency shifts based on the flux qubit's state-dependent persistent current:

where the sign and magnitude depend on the qubit state. Microwave measurement of the rf SQUID resonance frequency yields the qubit state with fidelities exceeding 99%.

Alternative readout: persistent current readout via dispersive interaction with a superconducting cavity. The qubit's state-dependent current shifts the cavity resonance frequency by chi = beta I_p / 2, enabling single-shot readout via cavity reflectance measurement.

## 10. Two-Qubit Gates and Scalability

Flux qubits can be coupled via:

1. Direct inductive coupling: Overlap magnetic fields from two qubit loops
2. Cavity-mediated coupling: Both qubits couple to a shared superconducting cavity
3. Tunable SQUID coupling: Bias lines enable dynamic activation/deactivation

The tunable-SQUID architecture has enabled >99% two-qubit gate fidelities and demonstrations of 2-3 qubit entanglement on single chips. Scalability to >10 qubits has been demonstrated; larger arrays face challenges in control electronics and crosstalk isolation.

## 11. Python Numerical Solver: Potential Energy and Tunneling Dynamics

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh
from scipy.integrate import quad
from scipy.optimize import fminbound

def three_junction_potential(varphi1, varphi2, flux_bias, alpha=0.7):
    """
    Three-junction SQUID potential energy.
    
    U(φ₁, φ₂) = -E_J [2 + α - cos(φ₁) - cos(φ₂) - α cos(2πf - φ₁ - φ₂)]
    
    Normalized to E_J units.
    
    Args:
        varphi1, varphi2: phase variables (radians)
        flux_bias: normalized external flux f = Φ_ext/Φ_0 (0 to 1)
        alpha: asymmetry parameter
    
    Returns:
        U: potential energy in units of E_J
    """
    phase_constraint = 2 * np.pi * flux_bias - varphi1 - varphi2
    
    U = -(2 + alpha - np.cos(varphi1) - np.cos(varphi2) - alpha * np.cos(phase_constraint))
    
    return U

def effective_qubit_potential_1d(varphi, flux_bias, alpha=0.7):
    """
    Effective 1D potential for flux qubit (varphi ≡ φ₁, φ₂ = π - φ₁ + δ).
    
    Simplified to single effective phase variable.
    """
    varphi2_eff = np.pi - varphi  # Approximation near symmetric point
    
    U = three_junction_potential(varphi, varphi2_eff, flux_bias, alpha)
    
    return U

def double_well_analysis(flux_bias, alpha=0.7, n_varphi=200):
    """
    Analyze double-well potential: extract well positions, heights, barrier.
    
    Returns:
        (varphi_min_left, varphi_min_right, barrier_height, well_separation)
    """
    varphi_range = np.linspace(-np.pi, np.pi, n_varphi)
    U_vals = [effective_qubit_potential_1d(v, flux_bias, alpha) for v in varphi_range]
    
    # Find minima (approximate)
    min_indices = np.where(np.diff(np.sign(np.diff(U_vals))) == -2)[0] + 1
    
    if len(min_indices) >= 2:
        varphi_min_left = varphi_range[min_indices[0]]
        varphi_min_right = varphi_range[min_indices[1]]
        U_min = min(U_vals[min_indices[0]], U_vals[min_indices[1]])
        
        # Barrier: maximum between wells
        barrier_region = np.argwhere((varphi_range > varphi_min_left) & (varphi_range < varphi_min_right))
        if len(barrier_region) > 0:
            barrier_idx = min_indices[0] + np.argmax([U_vals[i] for i in barrier_region.flatten()])
            U_barrier = U_vals[barrier_idx]
        else:
            U_barrier = max(U_vals)
        
        barrier_height = U_barrier - U_min
        well_separation = varphi_min_right - varphi_min_left
        
        return varphi_min_left, varphi_min_right, barrier_height, well_separation
    
    return None, None, None, None

def tunneling_gap_estimate(barrier_height, well_separation, prefactor=1.0):
    """
    Estimate tunneling gap via WKB approximation.
    
    Δ ~ A * exp(-S/ℏ) where S ≈ (1/π) * barrier_height * well_separation
    
    In normalized units (ℏ = 1, mass = m*).
    """
    if barrier_height is None or barrier_height < 0.01:
        return 0.0
    
    # Simplified WKB action
    action = barrier_height * well_separation / np.pi
    
    tunneling_gap = prefactor * np.exp(-action)
    
    return tunneling_gap

def qubit_energy_levels(flux_bias, alpha=0.7, n_levels=4):
    """
    Diagonalize the effective double-well Hamiltonian.
    
    Returns: array of energy eigenvalues
    """
    # Discretize 1D potential
    n_points = 100
    varphi_range = np.linspace(-np.pi, np.pi, n_points)
    dvarphi = varphi_range[1] - varphi_range[0]
    
    # Kinetic energy matrix (finite difference Laplacian)
    diag_KE = np.ones(n_points) * 2.0 / (dvarphi**2)
    off_diag_KE = -1.0 / (dvarphi**2) * np.ones(n_points - 1)
    T = np.diag(diag_KE) + np.diag(off_diag_KE, 1) + np.diag(off_diag_KE, -1)
    
    # Potential energy matrix (diagonal)
    V = np.diag([effective_qubit_potential_1d(v, flux_bias, alpha) for v in varphi_range])
    
    # Full Hamiltonian
    H = T + V
    
    # Diagonalize
    evals, evecs = eigh(H)
    
    return evals[:n_levels]

def flux_dephasing_time(tunneling_gap, persistent_current, flux_noise_spectral_density):
    """
    Estimate T2* from second-order flux noise sensitivity.
    
    T2* ~ 1 / [(d²ω/dΦ²)² * S_Φ(f)]
    
    d²ω/dΦ² ≈ 2I_p/Δ at sweet spot
    """
    if tunneling_gap < 1e-6:
        return 1e-6  # Min value
    
    d2omega_dPhi2 = 2 * persistent_current / tunneling_gap
    
    T2_star = 1.0 / ((d2omega_dPhi2)**2 * flux_noise_spectral_density)
    
    return np.clip(T2_star, 1e-6, 1e0)  # Clip to physical range (μs to s)

# Main execution
if __name__ == "__main__":
    print("="*70)
    print("SUPERCONDUCTING FLUX QUBITS: 3-JUNCTION SQUID")
    print("="*70)
    
    # Parameters
    alpha = 0.7  # Asymmetry parameter
    flux_bias_sweet = 0.5  # Half-flux quantum (sweet spot)
    flux_bias_scan = np.linspace(0.4, 0.6, 21)
    
    print(f"\nQubit Design Parameters:")
    print(f"  Junction asymmetry α = {alpha:.2f}")
    print(f"  Optimal bias point f = {flux_bias_sweet:.2f} (half-flux)")
    
    # Analysis at sweet spot
    vmin_L, vmin_R, barrier, well_sep = double_well_analysis(flux_bias_sweet, alpha)
    Delta = tunneling_gap_estimate(barrier, well_sep, prefactor=1.0)
    
    print(f"\nDouble-Well Analysis at Sweet Spot:")
    print(f"  Left minimum phase: {vmin_L:.3f} rad")
    print(f"  Right minimum phase: {vmin_R:.3f} rad")
    print(f"  Potential barrier height: {barrier:.3f} E_J")
    print(f"  Well separation: {well_sep:.3f} rad")
    print(f"  Tunneling gap Δ/(2π) ≈ {Delta/(2*np.pi)*1e3:.2f} MHz (normalized)")
    
    # Energy levels
    E_levels = qubit_energy_levels(flux_bias_sweet, alpha, n_levels=4)
    qubit_freq = E_levels[1] - E_levels[0]
    
    print(f"\nEnergy Spectrum at Sweet Spot:")
    for i, E in enumerate(E_levels[:4]):
        print(f"  E_{i} = {E:.4f} E_J")
    print(f"  Qubit transition E_01 = {qubit_freq:.4f} E_J")
    
    # Flux tunability
    print(f"\nFlux Tuning Performance:")
    tunneling_gaps = []
    qubit_freqs = []
    for f in flux_bias_scan:
        vmin_L, vmin_R, barrier, well_sep = double_well_analysis(f, alpha)
        Delta_f = tunneling_gap_estimate(barrier, well_sep, prefactor=1.0)
        tunneling_gaps.append(Delta_f)
        E_lev = qubit_energy_levels(f, alpha, n_levels=2)
        qubit_freqs.append(E_lev[1] - E_lev[0])
    
    max_tuning_range = max(qubit_freqs) - min(qubit_freqs)
    print(f"  Qubit frequency tuning range: {min(qubit_freqs):.4f} to {max(qubit_freqs):.4f} E_J")
    print(f"  Total tuning range: {max_tuning_range:.4f} E_J")
    
    # Dephasing
    I_p = 0.5  # Persistent current in units of I_c
    S_Phi = 1e-3  # Flux noise power spectral density (normalized)
    T2_star = flux_dephasing_time(Delta, I_p, S_Phi)
    
    print(f"\nCoherence Properties (at sweet spot):")
    print(f"  Persistent current I_p = {I_p} I_c")
    print(f"  Flux noise S_Φ = {S_Phi:.1e} Φ_0²/Hz")
    print(f"  Dephasing time T₂* ≈ {T2_star*1e6:.1f} μs")
    
    # Plotting
    fig, axes = plt.subplots(2, 2, figsize=(13, 10))
    
    # Plot 1: 2D Potential energy surface
    ax = axes[0, 0]
    n_grid = 100
    varphi1_grid = np.linspace(-np.pi, np.pi, n_grid)
    varphi2_grid = np.linspace(-np.pi, np.pi, n_grid)
    V_grid = np.zeros((n_grid, n_grid))
    
    for i, v1 in enumerate(varphi1_grid):
        for j, v2 in enumerate(varphi2_grid):
            V_grid[i, j] = three_junction_potential(v1, v2, flux_bias_sweet, alpha)
    
    im = ax.contourf(varphi1_grid, varphi2_grid, V_grid, levels=30, cmap='viridis')
    ax.set_xlabel('$\\varphi_1$ (rad)', fontsize=11)
    ax.set_ylabel('$\\varphi_2$ (rad)', fontsize=11)
    ax.set_title(f'3-Junction SQUID Potential at f = {flux_bias_sweet}', fontsize=12, fontweight='bold')
    plt.colorbar(im, ax=ax, label='$U / E_J$')
    
    # Plot 2: 1D Effective potential
    ax = axes[0, 1]
    varphi_1d = np.linspace(-np.pi, np.pi, 300)
    U_1d = [effective_qubit_potential_1d(v, flux_bias_sweet, alpha) for v in varphi_1d]
    
    ax.plot(varphi_1d, U_1d, 'b-', linewidth=2.5)
    ax.fill_between(varphi_1d, min(U_1d), U_1d, alpha=0.2)
    ax.axhline(y=min(U_1d), color='r', linestyle='--', alpha=0.5, label='Ground state')
    ax.set_xlabel('$\\varphi$ (rad)', fontsize=11)
    ax.set_ylabel('$U(\\varphi) / E_J$', fontsize=11)
    ax.set_title('Effective 1D Double-Well Potential (Sweet Spot)', fontsize=12, fontweight='bold')
    ax.set_xlim([-np.pi, np.pi])
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # Plot 3: Qubit frequency vs flux bias
    ax = axes[1, 0]
    ax.plot(flux_bias_scan, np.array(qubit_freqs)/max(qubit_freqs), 'go-', linewidth=2, markersize=6)
    ax.axvline(x=0.5, color='r', linestyle='--', alpha=0.7, label='Sweet spot')
    ax.set_xlabel('Flux Bias $f = \\Phi_{ext}/\\Phi_0$', fontsize=11)
    ax.set_ylabel('Normalized Qubit Frequency', fontsize=11)
    ax.set_title('Flux Tunability of Qubit Frequency', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # Plot 4: Tunneling gap and dephasing
    ax = axes[1, 1]
    T2_stars = [flux_dephasing_time(tg, I_p, S_Phi) for tg in tunneling_gaps]
    
    ax2 = ax.twinx()
    line1 = ax.semilogy(flux_bias_scan, tunneling_gaps, 'b-o', linewidth=2, markersize=6, label='Tunneling gap Δ')
    line2 = ax2.semilogy(flux_bias_scan, np.array(T2_stars)*1e6, 'r-s', linewidth=2, markersize=6, label='T_2*')
    
    ax.axvline(x=0.5, color='gray', linestyle='--', alpha=0.7)
    ax.set_xlabel('Flux Bias $f$', fontsize=11)
    ax.set_ylabel('Tunneling Gap Δ', color='b', fontsize=11)
    ax2.set_ylabel('Dephasing Time T_2* (μs)', color='r', fontsize=11)
    ax.set_title('Tunneling Gap and Coherence vs Flux', fontsize=12, fontweight='bold')
    ax.tick_params(axis='y', labelcolor='b')
    ax2.tick_params(axis='y', labelcolor='r')
    ax.grid(True, alpha=0.3)
    
    lines = line1 + line2
    labels = [l.get_label() for l in lines]
    ax.legend(lines, labels, loc='upper left', fontsize=9)
    
    plt.tight_layout()
    plt.savefig('flux_qubit_analysis.png', dpi=150, bbox_inches='tight')
    print(f"
[Plot saved as 'flux_qubit_analysis.png']")
    
    print("
" + "="*70)
    print("FLUX QUBIT ANALYSIS COMPLETE")
    print("="*70)

## References

1. Devoret, M. H., & Schoelkopf, R. J. (2013). "Superconducting circuits for quantum information." *Reviews of Modern Physics*, 85(3), 1143.

2. Krantz, P., et al. (2019). "A quantum engineer's guide to superconducting qubits." *Applied Physics Reviews*, 6(2), 021318.

3. Koch, J., et al. (2007). "Charge-insensitive qubit design derived from the Cooper pair box." *Physical Review A*, 76(4), 042319.

4. Siddiqi, I., et al. (2004). "Tunable Superconducting Qubits with Flux-Independent Coherence." *Physical Review Letters*, 93(20), 207002.

5. Clerk, A. A., et al. (2010). "Introduction to Quantum Noise, Measurement, and Amplification." *Reviews of Modern Physics*, 82(2), 1155.