第3章：Josephson効果

学習目標

• Josephson接合の種類（トンネル接合、弱結合、点接触）を理解する
• DC Josephson効果の理論的導出と臨界電流を説明する
• AC Josephson効果の電圧-周波数関係を導出する
• RCSJモデル（抵抗・容量分流接合モデル）を理解し、数値計算できる
• Stewart-McCumberパラメータとダンピング特性を解析する
• Josephson接合のI-V特性とヒステリシスを予測する
• DC SQUIDとRF SQUIDの動作原理と感度を理解する
• Pythonで位相ダイナミクスとSQUID特性をシミュレーションできる

3.1 Josephson接合の種類

弱結合と接合の定義

2つの超伝導体が弱い結合によって接続されている構造をJosephson接合と呼びます。この結合が「弱い」ことにより、2つの超伝導体の位相差 \(\delta = \phi_1 - \phi_2\) が独立した変数となり、顕著な量子効果が現れます。

主要な接合タイプ

1. トンネル接合（Tunnel Junction）

薄い絶縁層（厚さ \(d \sim 1-2\) nm）を介した量子トンネル効果による結合です。

• 構造：S-I-S（超伝導体-絶縁体-超伝導体）
• 臨界電流密度：\(J_c \sim 10^2 - 10^4\) A/cm²
• 特徴：高品質で再現性が良い、集積回路に適する
• 材料例：Nb/Al₂O₃/Nb、Nb/AlN/Nb

2. 弱結合（Weak Link）

超伝導特性が弱められた薄い領域を介した結合です。

• 構造：S-N-S（超伝導体-常伝導金属-超伝導体）
• メカニズム：近接効果により常伝導金属に超伝導相関が誘起される
• 特徴：臨界電流が温度・磁場に敏感
• 材料例：Nb/Cu/Nb、Al/Au/Al

3. 点接触（Point Contact）

超伝導体間の微小な導電性接触です。

• 構造：ナノスケールの金属接触（直径 \(\sim 10-100\) nm）
• 臨界電流：接触面積に比例（\(I_c \sim 1-100\) μA）
• 特徴：作製が容易だが再現性が低い
• 用途：研究用、簡易SQUID

3.2 DC Josephson効果

結合超伝導体からの導出

2つの超伝導体 1, 2 の秩序パラメータを \(\psi_1 = |\psi_1| e^{i\phi_1}\)、\(\psi_2 = |\psi_2| e^{i\phi_2}\) とします。弱い結合により、時間依存Schrödinger方程式は次のように結合します：

ここで、\(K\) は結合エネルギー（トンネルマトリクス要素）です。

電流-位相関係の導出

クーパー対の粒子数 \(n_1, n_2\) の時間変化は：

電流は \(I = 2e \frac{dn}{dt}\) なので、位相差 \(\delta = \phi_2 - \phi_1\) を導入すると：

臨界電流の物理的意味

\(I_c\) は、接合に電圧をかけずに流せる最大の超電流です：

• \(|I| < I_c\)：電圧ゼロで電流が流れる（ゼロ電圧状態）
• \(|I| > I_c\)：接合に電圧が発生し、抵抗状態に転移

Ambegaokar-Baratoff関係式

トンネル接合の臨界電流は、常伝導抵抗 \(R_n\) とエネルギーギャップ \(\Delta\) に関係します：

T = 0 では：

import numpy as np
import matplotlib.pyplot as plt

# Physical constants
e = 1.602e-19  # Elementary charge (C)
kB = 1.381e-23  # Boltzmann constant (J/K)

# Superconductor parameters (Nb example)
Tc = 9.2  # Critical temperature (K)
Delta_0 = 1.76 * kB * Tc  # Energy gap at T=0

# Normal resistance of junction
Rn = 10  # Ohms

# Temperature range
T = np.linspace(0.1, Tc, 100)
Delta_T = Delta_0 * np.sqrt(1 - (T/Tc)**2)  # BCS temperature dependence

# Ambegaokar-Baratoff formula
Ic_Rn = (np.pi * Delta_T / (2*e)) * np.tanh(Delta_T / (2*kB*T))
Ic = Ic_Rn / Rn  # Critical current

# Plot 1: Ic(T)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

axes[0].plot(T, Ic*1e6, 'b-', linewidth=2)
axes[0].axvline(x=Tc, color='r', linestyle='--', label=f'Tc = {Tc} K')
axes[0].set_xlabel('温度 T (K)', fontsize=12)
axes[0].set_ylabel('臨界電流 Ic (μA)', fontsize=12)
axes[0].set_title('Josephson接合の臨界電流の温度依存性', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)

# Plot 2: Current-phase relation
delta = np.linspace(-np.pi, np.pi, 500)
Ic_example = 100e-6  # 100 μA

I = Ic_example * np.sin(delta)

axes[1].plot(delta, I*1e6, 'g-', linewidth=2)
axes[1].axhline(y=Ic_example*1e6, color='r', linestyle='--',
                label=f'Ic = {Ic_example*1e6:.0f} μA')
axes[1].axhline(y=-Ic_example*1e6, color='r', linestyle='--')
axes[1].axvline(x=0, color='k', linestyle='-', alpha=0.3)
axes[1].axhline(y=0, color='k', linestyle='-', alpha=0.3)
axes[1].set_xlabel('位相差 δ (rad)', fontsize=12)
axes[1].set_ylabel('電流 I (μA)', fontsize=12)
axes[1].set_title('DC Josephson関係：I = Ic sin(δ)', fontsize=14)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)

# Mark key points
axes[1].plot([np.pi/2, -np.pi/2], [Ic_example*1e6, -Ic_example*1e6],
             'ro', markersize=8, label='最大電流点')
axes[1].legend(fontsize=10)

plt.tight_layout()
plt.show()

print(f"\nNb Josephson接合の特性 (Rn = {Rn} Ω):")
print(f"  T = 4.2 K における Ic = {Ic[42]*1e6:.2f} μA")
print(f"  T = 0 K における IcRn積 = {Ic_Rn[0]*1e3:.2f} mV")

3.3 AC Josephson効果

電圧-周波数関係の導出

接合に電圧 \(V\) が印加されると、2つの超伝導体のエネルギー差が生じます：

位相差 \(\delta = \phi_2 - \phi_1\) の時間発展は：

Josephson周波数

周波数 \(f_J = \omega_J / (2\pi)\) は：

ここで、Josephson定数は：

数値例

• V = 1 μV → f_J ≈ 483.6 MHz
• V = 1 mV → f_J ≈ 483.6 GHz
• V = 10 mV → f_J ≈ 4.836 THz

交流超電流

一定電圧 \(V\) 印加時、位相は \(\delta(t) = \delta_0 + \omega_J t\) と時間発展し、電流は振動します：

Shapiroステップ

マイクロ波照射（周波数 \(f_{RF}\)）下では、I-V特性に一定電圧のShapiroステップが現れます：

これは、Josephson振動とマイクロ波の同期（位相ロック）により生じます。

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

# Josephson constant
h = 6.626e-34  # Planck constant (J·s)
e = 1.602e-19  # Elementary charge (C)
K_J = 2*e / h  # Josephson constant (Hz/V)

print(f"Josephson定数 KJ = {K_J*1e-9:.3f} GHz/mV")

# AC Josephson effect: phase evolution under constant voltage
def phase_evolution(t, y, V):
    """dy/dt = 2eV/hbar"""
    delta = y[0]
    omega_J = 2 * np.pi * K_J * V  # Angular frequency
    return [omega_J]

# Voltage values
V_values = [1e-6, 10e-6, 100e-6]  # 1, 10, 100 μV

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

# Plot AC Josephson oscillations
for i, V in enumerate(V_values):
    f_J = K_J * V  # Josephson frequency (Hz)

    # Time range (several periods)
    t_span = (0, 5/f_J)
    t_eval = np.linspace(0, 5/f_J, 1000)

    # Solve phase evolution
    sol = solve_ivp(phase_evolution, t_span, [0], t_eval=t_eval, args=(V,))
    delta = sol.y[0]
    I_norm = np.sin(delta)  # Normalized current I/Ic

    # Plot current oscillation
    ax = axes[i//2, i%2]
    ax.plot(t_eval*1e9, I_norm, 'b-', linewidth=1.5)
    ax.set_xlabel('時間 (ns)', fontsize=11)
    ax.set_ylabel('規格化電流 I/Ic', fontsize=11)
    ax.set_title(f'AC Josephson振動\nV = {V*1e6:.0f} μV, fJ = {f_J*1e-6:.1f} MHz',
                 fontsize=12)
    ax.grid(True, alpha=0.3)
    ax.axhline(y=0, color='k', linestyle='-', alpha=0.3)

# Plot 4: Shapiro steps
f_RF = 10e9  # RF frequency: 10 GHz
n_steps = np.arange(-5, 6)
V_shapiro = n_steps / K_J / f_RF  # Voltage of each step

ax = axes[1, 1]
# Simulated I-V with Shapiro steps
I_bias = np.linspace(-2, 2, 500)  # Normalized bias current

# Simplified I-V with steps (phenomenological)
V_IV = np.zeros_like(I_bias)
for n in n_steps:
    V_step = n / (K_J * f_RF)
    # Add step at each Shapiro voltage
    mask = (np.abs(I_bias) > np.abs(n)*0.2) & (np.abs(I_bias) < np.abs(n)*0.2 + 0.3)
    V_IV[mask] = V_step

ax.plot(I_bias, V_IV*1e6, 'r-', linewidth=2, label='マイクロ波照射下')
ax.axhline(y=0, color='k', linestyle='-', alpha=0.3)
ax.axvline(x=0, color='k', linestyle='-', alpha=0.3)

# Mark Shapiro steps
for n in [-3, -2, -1, 0, 1, 2, 3]:
    V_step = n / (K_J * f_RF)
    ax.axhline(y=V_step*1e6, color='orange', linestyle='--', alpha=0.5)
    ax.text(2.1, V_step*1e6, f'n={n}', fontsize=9)

ax.set_xlabel('規格化バイアス電流 I/Ic', fontsize=11)
ax.set_ylabel('電圧 V (μV)', fontsize=11)
ax.set_title(f'Shapiroステップ (fRF = {f_RF*1e-9:.0f} GHz)', fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim([-2.5, 2.5])

plt.tight_layout()
plt.show()

# Print Shapiro step voltages
print(f"\nShapiroステップ電圧 (fRF = {f_RF*1e-9} GHz):")
for n in range(-3, 4):
    V_n = n / (K_J * f_RF)
    print(f"  n = {n:2d}: V = {V_n*1e6:8.4f} μV")

3.4 RCSJモデル（抵抗・容量分流接合モデル）

等価回路モデル

Resistively and Capacitively Shunted Junction (RCSJ)モデルは、Josephson接合を以下の並列回路として記述します：

• Josephson素子：超電流 \(I_J = I_c \sin\delta\)
• 抵抗 \(R\)：準粒子電流 \(I_R = V/R\)
• 容量 \(C\)：変位電流 \(I_C = C dV/dt\)

RCSJの回路方程式

電流保存則より：

AC Josephson関係式 \(V = (\hbar/2e) d\delta/dt\) を代入すると：

無次元化

Josephsonプラズマ周波数 \(\omega_p\) と特性電圧 \(V_c = I_c R\) を導入：

無次元時間 \(\tau = \omega_p t\) とすると：

ここで、Stewart-McCumberパラメータは：

物理的アナロジー：洗濯板ポテンシャル

RCSJ方程式は、傾いた洗濯板上の粘性減衰を受ける粒子の運動方程式と等価です：

• 粒子の位置：位相 \(\delta\)
• ポテンシャル：\(U(\delta) = -E_J \cos\delta - (\hbar I/2e) \delta\)
• 減衰係数：\(\gamma = 1/\beta_c\)
• 質量：\(m \propto C\)

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

# RCSJ model simulation
def rcsj_ode(t, y, I_bias, beta_c):
    """
    RCSJ model ODEs in dimensionless form.
    y[0] = delta (phase)
    y[1] = d(delta)/dtau (phase velocity)
    """
    delta, delta_dot = y

    # d²δ/dτ² + (1/βc) dδ/dτ + sin(δ) = I/Ic
    delta_ddot = I_bias - np.sin(delta) - delta_dot / beta_c

    return [delta_dot, delta_ddot]

# Simulation parameters
beta_c_values = [0.5, 1.0, 5.0, 20.0]  # Overdamped to underdamped
I_bias_range = np.linspace(0, 2.5, 50)

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

for idx, beta_c in enumerate(beta_c_values):
    ax = axes[idx // 2, idx % 2]

    V_avg_list = []

    for I_bias in I_bias_range:
        # Time span (dimensionless)
        t_span = (0, 200)
        t_eval = np.linspace(0, 200, 2000)

        # Initial conditions
        y0 = [0, 0]  # delta=0, delta_dot=0

        # Solve RCSJ equations
        sol = solve_ivp(rcsj_ode, t_span, y0, t_eval=t_eval,
                        args=(I_bias, beta_c), method='RK45')

        # Average voltage (proportional to average phase velocity)
        # Skip transient (first 50 time units)
        transient_idx = np.where(sol.t > 50)[0][0]
        delta_dot_avg = np.mean(sol.y[1, transient_idx:])
        V_avg = delta_dot_avg  # V/Vc = 

        V_avg_list.append(V_avg)

    # Plot I-V curve
    ax.plot(I_bias_range, V_avg_list, 'b-', linewidth=2)
    ax.axhline(y=0, color='k', linestyle='-', alpha=0.3)
    ax.axvline(x=1, color='r', linestyle='--', alpha=0.5, label='I = Ic')

    ax.set_xlabel('バイアス電流 I/Ic', fontsize=11)
    ax.set_ylabel('平均電圧 V/Vc', fontsize=11)
    ax.set_title(f'RCSJ I-V特性 (βc = {beta_c})', fontsize=12)
    ax.legend(fontsize=10)
    ax.grid(True, alpha=0.3)
    ax.set_xlim([0, 2.5])

    # Indicate damping regime
    if beta_c < 1:
        regime = "オーバーダンプ"
    elif beta_c < 5:
        regime = "中間減衰"
    else:
        regime = "アンダーダンプ"
    ax.text(0.05, 0.95, regime, transform=ax.transAxes,
            fontsize=11, verticalalignment='top',
            bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

plt.tight_layout()
plt.show()

print("\nStewart-McCumberパラメータ βc の影響:")
print("βc < 1: オーバーダンプ → ヒステリシスなし、なめらかなI-V")
print("βc > 1: アンダーダンプ → ヒステリシスあり、急峻な転移")

3.5 I-V特性とヒステリシス

ゼロ電圧状態

バイアス電流 \(I < I_c\) では、位相が一定値に固定され、電圧はゼロです：

電圧状態

\(I > I_c\) では、位相が時間的に回転し、有限の時間平均電圧が発生します：

ヒステリシス（アンダーダンプ接合）

\(\beta_c > 1\) の場合、I-V特性にヒステリシスが現れます：

• 上昇時：\(I > I_c\) で急激に電圧状態へ転移
• 下降時：\(I < I_r\)（リトラップ電流、\(I_r < I_c\)）で再びゼロ電圧状態に戻る

リトラップ電流

アンダーダンプ極限（\(\beta_c \gg 1\)）では：

熱揺らぎと位相拡散

有限温度では、熱エネルギー \(k_B T\) により位相が確率的に揺らぎます。これにより：

• I-V特性が丸まる（sharpnessの低下）
• 臨界電流が減少
• ゼロ電圧状態の安定性が低下

熱揺らぎの強さは、Josephsonエネルギー \(E_J = \hbar I_c / 2e\) と比較されます：

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

# Hysteretic I-V characteristics for underdamped junction
def rcsj_hysteresis(beta_c):
    """
    Simulate I-V with hysteresis by sweeping current up and down.
    """
    # Current sweep: up and down
    I_up = np.linspace(0, 2.5, 100)
    I_down = np.linspace(2.5, 0, 100)

    V_up = []
    V_down = []

    # Initial state
    delta = 0
    delta_dot = 0

    # Sweep up
    for I_bias in I_up:
        t_span = (0, 100)
        t_eval = np.linspace(0, 100, 1000)

        sol = solve_ivp(lambda t, y: rcsj_ode(t, y, I_bias, beta_c),
                        t_span, [delta, delta_dot], t_eval=t_eval, method='RK45')

        # Update state
        delta = sol.y[0, -1] % (2*np.pi)
        delta_dot = sol.y[1, -1]

        # Average voltage
        V_avg = np.mean(sol.y[1, -200:])  # Last portion
        V_up.append(V_avg)

    # Sweep down (start from voltage state)
    for I_bias in I_down:
        t_span = (0, 100)
        t_eval = np.linspace(0, 100, 1000)

        sol = solve_ivp(lambda t, y: rcsj_ode(t, y, I_bias, beta_c),
                        t_span, [delta, delta_dot], t_eval=t_eval, method='RK45')

        delta = sol.y[0, -1] % (2*np.pi)
        delta_dot = sol.y[1, -1]

        V_avg = np.mean(sol.y[1, -200:])
        V_down.append(V_avg)

    return I_up, V_up, I_down, V_down

# Compare different beta_c values
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Low beta_c (no hysteresis)
beta_c_low = 0.5
I_up, V_up, I_down, V_down = rcsj_hysteresis(beta_c_low)

axes[0].plot(I_up, V_up, 'b-', linewidth=2, label='電流上昇')
axes[0].plot(I_down, V_down, 'r--', linewidth=2, label='電流下降')
axes[0].set_xlabel('バイアス電流 I/Ic', fontsize=12)
axes[0].set_ylabel('電圧 V/Vc', fontsize=12)
axes[0].set_title(f'オーバーダンプ接合 (βc = {beta_c_low})\nヒステリシスなし', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim([0, 2.5])
axes[0].set_ylim([0, 2])

# High beta_c (with hysteresis)
beta_c_high = 10.0
I_up, V_up, I_down, V_down = rcsj_hysteresis(beta_c_high)

axes[1].plot(I_up, V_up, 'b-', linewidth=2, label='電流上昇')
axes[1].plot(I_down, V_down, 'r--', linewidth=2, label='電流下降')
axes[1].axvline(x=1, color='orange', linestyle='--', alpha=0.5, label='Ic')

# Estimate retrap current
I_r_estimate = 4 / (np.pi * beta_c_high)
axes[1].axvline(x=I_r_estimate, color='purple', linestyle='--',
                alpha=0.5, label=f'Ir ≈ {I_r_estimate:.2f}')

axes[1].set_xlabel('バイアス電流 I/Ic', fontsize=12)
axes[1].set_ylabel('電圧 V/Vc', fontsize=12)
axes[1].set_title(f'アンダーダンプ接合 (βc = {beta_c_high})\nヒステリシスあり', fontsize=14)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim([0, 2.5])
axes[1].set_ylim([0, 2])

# Shade hysteresis region
axes[1].fill_betweenx([0, 2], I_r_estimate, 1, alpha=0.2, color='yellow',
                      label='ヒステリシス領域')

plt.tight_layout()
plt.show()

print(f"\nヒステリシス特性 (βc = {beta_c_high}):")
print(f"  臨界電流: Ic/Ic = 1.0")
print(f"  リトラップ電流: Ir/Ic ≈ {I_r_estimate:.3f}")
print(f"  ヒステリシス幅: ΔI/Ic ≈ {1 - I_r_estimate:.3f}")

3.6 DC SQUID（二接合SQUIDループ）

構成と動作原理

DC SQUID（Superconducting Quantum Interference Device）は、2つのJosephson接合を超伝導ループで結合した構成です。

構造

• インダクタンス \(L\) のループ
• 2つの同一接合（臨界電流 \(I_0\)）
• ループを貫く磁束 \(\Phi\)

臨界電流変調

2つの接合の位相差 \(\delta_1, \delta_2\) は、磁束量子化条件により関係付けられます：

総電流 \(I = I_0 \sin\delta_1 + I_0 \sin\delta_2\) を最大化すると：

臨界電流は磁束に対して周期 \(\Phi_0\) で振動します。

電圧変調

バイアス電流 \(I_b\) を一定に保つと、磁束変化に応じて電圧が変調されます：

磁束-電圧変換感度

最適バイアス点での変換係数は：

ノイズと感度

DC SQUIDの磁束ノイズは：

典型的な磁束感度：\(\Phi_n \sim 10^{-6} \Phi_0 / \sqrt{\text{Hz}}\)

import numpy as np
import matplotlib.pyplot as plt

# DC SQUID characteristics
Phi_0 = 2.067e-15  # Flux quantum (Wb)
I_0 = 10e-6  # Single junction critical current (A)
R = 10  # Shunt resistance (Ohm)

# Flux range
Phi = np.linspace(-2*Phi_0, 2*Phi_0, 500)

# Critical current modulation
I_c_squid = 2 * I_0 * np.abs(np.cos(np.pi * Phi / Phi_0))

# Bias currents
I_bias_values = [1.2*I_0, 1.5*I_0, 1.8*I_0]

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

# Plot 1: Critical current modulation
axes[0, 0].plot(Phi/Phi_0, I_c_squid*1e6, 'b-', linewidth=2)
axes[0, 0].set_xlabel('磁束 Φ/Φ₀', fontsize=12)
axes[0, 0].set_ylabel('臨界電流 Ic (μA)', fontsize=12)
axes[0, 0].set_title('DC SQUID臨界電流の磁束変調', fontsize=14)
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].axhline(y=2*I_0*1e6, color='r', linestyle='--',
                   alpha=0.5, label='最大値 2I₀')
axes[0, 0].axhline(y=0, color='k', linestyle='-', alpha=0.3)
axes[0, 0].legend(fontsize=10)

# Plot 2: Voltage modulation for different bias currents
for I_bias in I_bias_values:
    V_squid = np.zeros_like(Phi)
    for i, phi in enumerate(Phi):
        I_c = 2 * I_0 * np.abs(np.cos(np.pi * phi / Phi_0))
        if I_bias > I_c:
            V_squid[i] = R * np.sqrt(I_bias**2 - I_c**2)
        else:
            V_squid[i] = 0

    axes[0, 1].plot(Phi/Phi_0, V_squid*1e6, linewidth=2,
                    label=f'Ib = {I_bias/I_0:.1f}I₀')

axes[0, 1].set_xlabel('磁束 Φ/Φ₀', fontsize=12)
axes[0, 1].set_ylabel('電圧 V (μV)', fontsize=12)
axes[0, 1].set_title('DC SQUID電圧変調（異なるバイアス電流）', fontsize=14)
axes[0, 1].legend(fontsize=10)
axes[0, 1].grid(True, alpha=0.3)

# Plot 3: Transfer function V_Phi
# Choose optimal bias
I_bias_opt = 1.5 * I_0
Phi_fine = np.linspace(-Phi_0, Phi_0, 1000)
V_fine = np.zeros_like(Phi_fine)

for i, phi in enumerate(Phi_fine):
    I_c = 2 * I_0 * np.abs(np.cos(np.pi * phi / Phi_0))
    if I_bias_opt > I_c:
        V_fine[i] = R * np.sqrt(I_bias_opt**2 - I_c**2)

# Calculate transfer function (numerical derivative)
V_Phi = np.gradient(V_fine, Phi_fine)

axes[1, 0].plot(Phi_fine/Phi_0, V_Phi, 'g-', linewidth=2)
axes[1, 0].set_xlabel('磁束 Φ/Φ₀', fontsize=12)
axes[1, 0].set_ylabel('変換係数 ∂V/∂Φ (V/Wb)', fontsize=12)
axes[1, 0].set_title(f'磁束-電圧変換特性 (Ib = {I_bias_opt/I_0:.1f}I₀)', fontsize=14)
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].axhline(y=0, color='k', linestyle='-', alpha=0.3)

# Mark working point (maximum sensitivity)
max_idx = np.argmax(np.abs(V_Phi))
axes[1, 0].plot(Phi_fine[max_idx]/Phi_0, V_Phi[max_idx], 'ro',
                markersize=10, label='最大感度点')
axes[1, 0].legend(fontsize=10)

# Plot 4: SQUID applications
ax = axes[1, 1]
ax.axis('off')

applications_text = """
DC SQUIDの応用

【高感度磁気センサー】
• 磁束感度: ~10⁻⁶ Φ₀/√Hz
• 磁場感度: ~fT/√Hz レベル
• 生体磁気計測（MEG, MCG）
• 非破壊検査（NDE）

【量子情報】
• 超伝導量子ビット読み出し
• 量子状態測定
• 量子非破壊測定

【基礎物理】
• 微小磁気モーメント測定
• 重力波検出器の変位センサー
• 暗黒物質探索

【主要パラメータ】
• 最大感度: ∂V/∂Φ ~ R/Φ₀
• 最適動作点: Ib ≈ 1.5 I₀
• ダイナミックレンジ: ~Φ₀
• 周波数帯域: DC - GHz
"""

ax.text(0.1, 0.9, applications_text, transform=ax.transAxes,
        fontsize=11, verticalalignment='top', family='monospace',
        bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

plt.tight_layout()
plt.show()

# Print sensitivity
V_Phi_max = np.max(np.abs(V_Phi))
print(f"\nDC SQUID感度解析 (I₀ = {I_0*1e6:.0f} μA, R = {R} Ω):")
print(f"  最大変換係数: ∂V/∂Φ ≈ {V_Phi_max:.2e} V/Wb")
print(f"  理論限界: R/Φ₀ = {R/Phi_0:.2e} V/Wb")
print(f"  最適バイアス: Ib = {I_bias_opt*1e6:.1f} μA")

3.7 RF SQUID（単一接合SQUIDループ）

構成

RF SQUIDは、単一のJosephson接合を含む超伝導ループです。DC SQUIDより構造が単純ですが、感度はやや劣ります。

主要構成要素

• インダクタンス \(L\) のループ
• 単一接合（臨界電流 \(I_c\)）
• RFタンク回路（共振周波数 \(f_{RF} \sim 10-100\) MHz）

動作原理

RF SQUIDは、ループのインダクタンス変化を共振器のQ値変化として検出します。

磁束量子化条件

ループを貫く総磁束は：

ここで、\(I = I_c \sin\delta\) は循環電流です。量子化条件より：

ポテンシャルエネルギー

無次元パラメータ \(\beta_L = 2LI_c/\Phi_0\) を導入すると、系のエネルギーは：

βL の役割

• \(\beta_L < 1\)：単一最小値（モノステーブル）
• \(\beta_L > 1\)：二重最小値（バイステーブル）、ヒステリシスあり

タンク回路結合

RFタンク回路の共振周波数は、SQUIDのインダクタンスに依存します：

磁束変化により \(L_{\text{SQUID}}\) が変調され、共振周波数がシフトします。

import numpy as np
import matplotlib.pyplot as plt

# RF SQUID potential energy
def rf_squid_potential(delta, Phi_ext, beta_L):
    """
    Calculate RF SQUID potential energy (normalized).
    U(delta) ~ (Phi_ext/Phi_0 - delta/(2pi))^2 + beta_L(1 - cos(delta))
    """
    phi_norm = Phi_ext  # Already normalized to Phi_0
    U = (phi_norm - delta/(2*np.pi))**2 + beta_L * (1 - np.cos(delta))
    return U

# Phase range
delta = np.linspace(-np.pi, 3*np.pi, 500)

# Different beta_L values
beta_L_values = [0.5, 1.0, 2.0, 5.0]

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

for idx, beta_L in enumerate(beta_L_values):
    ax = axes[idx // 2, idx % 2]

    # Plot potential for different external flux
    Phi_ext_values = [0, 0.25, 0.5, 0.75]

    for Phi_ext in Phi_ext_values:
        U = rf_squid_potential(delta, Phi_ext, beta_L)
        ax.plot(delta/np.pi, U, linewidth=2, label=f'Φ/Φ₀ = {Phi_ext}')

    ax.set_xlabel('位相 δ/π', fontsize=11)
    ax.set_ylabel('ポテンシャル U (任意単位)', fontsize=11)
    ax.set_title(f'RF SQUIDポテンシャル (βL = {beta_L})', fontsize=12)
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)
    ax.set_xlim([-1, 3])

    # Indicate stability regime
    if beta_L < 1:
        regime = "モノステーブル"
        color = 'lightblue'
    else:
        regime = "バイステーブル"
        color = 'lightcoral'

    ax.text(0.05, 0.95, regime, transform=ax.transAxes,
            fontsize=11, verticalalignment='top',
            bbox=dict(boxstyle='round', facecolor=color, alpha=0.7))

plt.tight_layout()
plt.show()

# RF SQUID vs DC SQUID comparison
print("\nRF SQUID vs DC SQUID 比較:")
print("\n【RF SQUID】")
print("  利点: 構造が単純、1つの接合のみ")
print("  欠点: 感度がやや低い、RF回路が必要")
print("  典型的感度: ~10⁻⁴ - 10⁻⁵ Φ₀/√Hz")
print("  用途: 地磁気測定、岩石磁気、簡易センサー")

print("\n【DC SQUID】")
print("  利点: 高感度、DC動作可能")
print("  欠点: 2つの接合が必要、対称性が重要")
print("  典型的感度: ~10⁻⁶ Φ₀/√Hz")
print("  用途: 生体磁気、精密測定、量子情報")

まとめ

演習問題