学習目標
本章を学習することで、以下の知識とスキルを習得できます：
- Bogoliubov-de Gennes (BdG) 方程式を導出し、その物理的意味を理解する
- Andreev反射の機構と超伝導/常伝導界面での準粒子変換を説明できる
- 近接効果の物理を理解し、SNSジャンクションでの応用を学ぶ
- 超伝導量子ドットとYu-Shiba-Rusinov (YSR) 状態の理論を理解する
- トポロジカル超伝導体とMajorana束縛状態の基礎を学ぶ
- タイトバインディングBdG方程式を数値的に解くPythonコードを実装できる
- メゾスコピック超伝導の量子コンピューティングへの応用を理解する
1. メゾスコピック超伝導とは
メゾスコピック超伝導は、系のサイズが超伝導のコヒーレンス長$\xi$と同程度または小さくなる領域における超伝導現象を扱います。このスケールでは、量子力学的な波動関数の位相コヒーレンスが重要な役割を果たし、バルク超伝導とは異なる新奇な物理が現れます。
メゾスコピック領域の特徴
- サイズ効果: 系のサイズ $L \sim \xi$ で準粒子スペクトルが離散化
- 界面効果: 超伝導/常伝導界面でのAndreev反射と近接効果
- 位相コヒーレンス: 超伝導位相の空間変化が顕著に
- 量子閉じ込め: 量子ドットや1次元系での離散準位構造
1.1 特徴的な長さスケール
メゾスコピック超伝導を理解するには、以下の特徴的長さスケールの理解が不可欠です：
| 長さスケール | 定義 | 典型的な値 | 物理的意味 |
|---|---|---|---|
| コヒーレンス長 $\xi_0$ | $\hbar v_F / \pi \Delta_0$ | 10-1000 nm | Cooper対のサイズ |
| ロンドン侵入長 $\lambda_L$ | $\sqrt{m/\mu_0 n_s e^2}$ | 50-500 nm | 磁場の侵入深さ |
| 熱的長さ $\xi_T$ | $\hbar v_F / 2\pi k_B T$ | T依存 | 有限温度での準粒子コヒーレンス |
| Thouless長さ $L_T$ | $\sqrt{\hbar D / E_T}$ | 拡散依存 | エネルギー緩和長 |
2. Bogoliubov-de Gennes方程式
2.1 BdG方程式の導出
BCS平均場ハミルトニアンから出発して、空間変化する超伝導秩序を記述するBogoliubov-de Gennes (BdG) 方程式を導出します。
出発点は平均場BCSハミルトニアン：
$$ H = \sum_{\sigma} \int d\mathbf{r} \, \psi_{\sigma}^{\dagger}(\mathbf{r}) \left[ -\frac{\hbar^2 \nabla^2}{2m} + V(\mathbf{r}) - \mu \right] \psi_{\sigma}(\mathbf{r}) + \int d\mathbf{r} \left[ \Delta(\mathbf{r}) \psi_{\uparrow}^{\dagger}(\mathbf{r}) \psi_{\downarrow}^{\dagger}(\mathbf{r}) + \text{h.c.} \right] $$ここで $\Delta(\mathbf{r})$ は空間に依存する超伝導ギャップ関数です。準粒子演算子を導入：
$$ \psi_{\uparrow}(\mathbf{r}) = \sum_n \left[ u_n(\mathbf{r}) \gamma_n - v_n^*(\mathbf{r}) \gamma_n^{\dagger} \right] $$ $$ \psi_{\downarrow}(\mathbf{r}) = \sum_n \left[ v_n(\mathbf{r}) \gamma_n + u_n^*(\mathbf{r}) \gamma_n^{\dagger} \right] $$対角化により、BdG方程式が得られます：
Bogoliubov-de Gennes方程式
$$ \begin{pmatrix} H_0(\mathbf{r}) & \Delta(\mathbf{r}) \\ \Delta^*(\mathbf{r}) & -H_0^*(\mathbf{r}) \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} $$ここで $H_0(\mathbf{r}) = -\frac{\hbar^2 \nabla^2}{2m} + V(\mathbf{r}) - \mu$ は単一粒子ハミルトニアンです。
2.2 BdG方程式の性質
BdG方程式は以下の重要な性質を持ちます：
- 粒子-正孔対称性: $(u_n, v_n, E_n)$ が解ならば $(v_n^*, u_n^*, -E_n)$ も解
- 正規化条件: $\int d\mathbf{r} \, (|u_n(\mathbf{r})|^2 + |v_n(\mathbf{r})|^2) = 1$
- 自己無撞着性: ギャップ関数は $\Delta(\mathbf{r}) = g \sum_{E_n > 0} u_n(\mathbf{r}) v_n^*(\mathbf{r}) \tanh(E_n/2k_B T)$ で更新
2.3 1次元BdG方程式の例
1次元系での離散化BdG方程式は、実際の計算で広く用いられます。タイトバインディング近似では：
$$ H_{\text{BdG}} = \sum_i \begin{pmatrix} c_i^{\dagger} & c_i \end{pmatrix} \begin{pmatrix} \epsilon_i - \mu & \Delta_i \\ \Delta_i^* & -(\epsilon_i - \mu) \end{pmatrix} \begin{pmatrix} c_i \\ c_i^{\dagger} \end{pmatrix} + \sum_{\langle i,j \rangle} \begin{pmatrix} c_i^{\dagger} & c_i \end{pmatrix} \begin{pmatrix} -t & 0 \\ 0 & t \end{pmatrix} \begin{pmatrix} c_j \\ c_j^{\dagger} \end{pmatrix} $$Python実装例: 1次元BdG方程式のソルバー
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh
class BdGSolver1D:
"""1次元Bogoliubov-de Gennes方程式のソルバー"""
def __init__(self, N_sites, t, mu, Delta):
"""
Parameters:
-----------
N_sites : int
格子点の数
t : float
ホッピング積分 (eV)
mu : float
化学ポテンシャル (eV)
Delta : float or array
超伝導ギャップ (eV), スカラーまたはサイト依存配列
"""
self.N = N_sites
self.t = t
self.mu = mu
# ギャップ関数の設定（スカラーの場合は配列に変換）
if np.isscalar(Delta):
self.Delta = Delta * np.ones(N_sites)
else:
self.Delta = np.array(Delta)
def build_hamiltonian(self, potential=None):
"""BdGハミルトニアン行列を構築
Parameters:
-----------
potential : array, optional
サイトポテンシャル V(i)
Returns:
--------
H_BdG : ndarray
2N × 2N BdGハミルトニアン行列
"""
N = self.N
H_BdG = np.zeros((2*N, 2*N), dtype=complex)
# ポテンシャルの設定
if potential is None:
potential = np.zeros(N)
# 対角要素（サイトエネルギー）
for i in range(N):
# 粒子ブロック
H_BdG[i, i] = potential[i] - self.mu
# 正孔ブロック
H_BdG[N+i, N+i] = -(potential[i] - self.mu)
# ギャップ関数（粒子-正孔結合）
for i in range(N):
H_BdG[i, N+i] = self.Delta[i]
H_BdG[N+i, i] = np.conj(self.Delta[i])
# ホッピング項
for i in range(N-1):
# 粒子ブロック
H_BdG[i, i+1] = -self.t
H_BdG[i+1, i] = -self.t
# 正孔ブロック
H_BdG[N+i, N+i+1] = self.t
H_BdG[N+i+1, N+i] = self.t
return H_BdG
def solve(self, potential=None):
"""BdG方程式を対角化
Returns:
--------
energies : ndarray
準粒子エネルギー（昇順）
u_vectors : ndarray
粒子成分 u_n(i)
v_vectors : ndarray
正孔成分 v_n(i)
"""
H_BdG = self.build_hamiltonian(potential)
# エルミート行列の対角化
energies, eigenvectors = eigh(H_BdG)
N = self.N
u_vectors = eigenvectors[:N, :]
v_vectors = eigenvectors[N:, :]
return energies, u_vectors, v_vectors
def calculate_ldos(self, energies, u_vectors, v_vectors,
energy_range, eta=0.01):
"""局所状態密度 (LDOS) を計算
Parameters:
-----------
energies : ndarray
準粒子エネルギー
u_vectors, v_vectors : ndarray
Bogoliubov波動関数
energy_range : ndarray
LDOSを計算するエネルギー範囲
eta : float
ブロードニングパラメータ (eV)
Returns:
--------
ldos : ndarray (N_sites, len(energy_range))
各サイトのLDOS
"""
N = self.N
ldos = np.zeros((N, len(energy_range)))
for i in range(N):
for n in range(2*N):
# Lorentzブロードニング
ldos[i, :] += (np.abs(u_vectors[i, n])**2 +
np.abs(v_vectors[i, n])**2) * \
(eta / np.pi) / ((energy_range - energies[n])**2 + eta**2)
return ldos
# 使用例：一様超伝導体の計算
if __name__ == "__main__":
# パラメータ設定
N_sites = 100
t = 1.0 # eV
mu = 0.0 # eV（ハーフフィリング）
Delta = 0.2 # eV
# ソルバーの初期化
solver = BdGSolver1D(N_sites, t, mu, Delta)
# BdG方程式を解く
energies, u_vecs, v_vecs = solver.solve()
# 準粒子スペクトルのプロット
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
# エネルギースペクトル
ax1.plot(energies, 'o', markersize=3)
ax1.axhline(y=Delta, color='r', linestyle='--', label=f'Δ = {Delta} eV')
ax1.axhline(y=-Delta, color='r', linestyle='--')
ax1.set_xlabel('State index')
ax1.set_ylabel('Energy (eV)')
ax1.set_title('BdG Quasiparticle Spectrum')
ax1.legend()
ax1.grid(True, alpha=0.3)
# LDOS計算
E_range = np.linspace(-0.5, 0.5, 500)
ldos = solver.calculate_ldos(energies, u_vecs, v_vecs, E_range)
# 中心サイトのLDOS
center = N_sites // 2
ax2.plot(E_range, ldos[center, :])
ax2.axvline(x=Delta, color='r', linestyle='--', alpha=0.5)
ax2.axvline(x=-Delta, color='r', linestyle='--', alpha=0.5)
ax2.set_xlabel('Energy (eV)')
ax2.set_ylabel('LDOS (states/eV)')
ax2.set_title(f'Local Density of States (site {center})')
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('bdg_spectrum_ldos.png', dpi=300, bbox_inches='tight')
plt.show()
print(f"最低励起エネルギー: {np.min(np.abs(energies)):.4f} eV")
print(f"理論ギャップ: {Delta:.4f} eV")3. Andreev反射と近接効果
3.1 Andreev反射の機構
Andreev反射は、超伝導/常伝導 (SN) 界面における特有の準粒子変換過程です。エネルギーが超伝導ギャップより小さい電子が超伝導体に入射すると、電子は正孔に変換されて反射され、Cooper対が超伝導体に注入されます。
graph LR
A[常伝導金属
N] -->|電子 e⁻
E < Δ| B[SN界面] B -->|正孔 h⁺| A B -->|Cooper対| C[超伝導体
S] style A fill:#e8f4f8,stroke:#2980b9,stroke-width:2px style B fill:#ffeaa7,stroke:#fdcb6e,stroke-width:2px style C fill:#dfe6e9,stroke:#636e72,stroke-width:2px
N] -->|電子 e⁻
E < Δ| B[SN界面] B -->|正孔 h⁺| A B -->|Cooper対| C[超伝導体
S] style A fill:#e8f4f8,stroke:#2980b9,stroke-width:2px style B fill:#ffeaa7,stroke:#fdcb6e,stroke-width:2px style C fill:#dfe6e9,stroke:#636e72,stroke-width:2px
Andreev反射の境界条件
SN界面（$x=0$）での波動関数のマッチング条件：
$$ u_N(0) = u_S(0), \quad v_N(0) = v_S(0) $$ $$ \left. \frac{du_N}{dx} \right|_0 = \left. \frac{du_S}{dx} \right|_0, \quad \left. \frac{dv_N}{dx} \right|_0 = \left. \frac{dv_S}{dx} \right|_0 $$Andreev近似（$E \ll E_F$）では、入射電子と反射正孔のフェルミ波数がほぼ等しく：
$$ k_e^- \approx k_h^+ = k_F \left(1 - \frac{E}{E_F}\right) \approx k_F $$3.2 Andreev反射確率とBTK理論
Blonder-Tinkham-Klapwijk (BTK) 理論は、SN界面でのAndreev反射を定量的に記述します。界面散乱を特徴づける無次元パラメータ $Z$ を導入：
$$ Z = \frac{H_b}{k_F \hbar v_F} $$ここで $H_b$ は界面障壁の高さです。Andreev反射確率は：
$$ A(E) = \frac{1 + Z^2 - \left(\frac{E}{\Delta}\right)^2}{\left[1 + Z^2 - \left(\frac{E}{\Delta}\right)^2\right]^2 + \left(\frac{E}{\Delta}\right)^2} $$Python実装例: BTK理論によるAndreev反射確率
import numpy as np
import matplotlib.pyplot as plt
def andreev_reflection_probability(E, Delta, Z):
"""Andreev反射確率（BTK理論）
Parameters:
-----------
E : float or array
準粒子エネルギー
Delta : float
超伝導ギャップ
Z : float
界面障壁パラメータ
Returns:
--------
A : float or array
Andreev反射確率
"""
# エネルギーを正規化
e = E / Delta
# |E| < Delta の範囲
mask_gap = np.abs(e) < 1.0
A = np.zeros_like(e)
# ギャップ内のAndreev反射
numerator = 1 + Z**2 - e[mask_gap]**2
denominator = (1 + Z**2 - e[mask_gap]**2)**2 + e[mask_gap]**2
A[mask_gap] = numerator / denominator
return A
def normal_reflection_probability(E, Delta, Z):
"""正常反射確率（BTK理論）"""
e = E / Delta
mask_gap = np.abs(e) < 1.0
B = np.zeros_like(e)
numerator = Z**2 * e[mask_gap]**2
denominator = (1 + Z**2 - e[mask_gap]**2)**2 + e[mask_gap]**2
B[mask_gap] = numerator / denominator
return B
def conductance_ns_junction(E, Delta, Z, T=0):
"""NS接合のコンダクタンス（BTK理論）
Parameters:
-----------
E : float or array
印加電圧 eV (= 準粒子エネルギー)
Delta : float
超伝導ギャップ
Z : float
界面障壁パラメータ
T : float
温度 (K)
Returns:
--------
G_NS : float or array
正規化コンダクタンス G_NS / G_N
"""
A = andreev_reflection_probability(E, Delta, Z)
B = normal_reflection_probability(E, Delta, Z)
# T=0でのコンダクタンス
# Andreev反射は電荷2eを運ぶため係数2
G_NS = 1 + A - B
return G_NS
# プロット
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
Delta = 1.0 # 正規化単位
E = np.linspace(-2, 2, 500)
# (a) Andreev反射確率 vs Z
Z_values = [0, 0.5, 1.0, 2.0, 5.0]
ax = axes[0, 0]
for Z in Z_values:
A = andreev_reflection_probability(E, Delta, Z)
ax.plot(E/Delta, A, label=f'Z = {Z}')
ax.axvline(x=-1, color='k', linestyle='--', alpha=0.3)
ax.axvline(x=1, color='k', linestyle='--', alpha=0.3)
ax.set_xlabel('E / Δ')
ax.set_ylabel('Andreev Reflection Probability A(E)')
ax.set_title('(a) Andreev Reflection vs Barrier Strength')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(-2, 2)
ax.set_ylim(0, 1.05)
# (b) 正常反射確率 vs Z
ax = axes[0, 1]
for Z in Z_values:
B = normal_reflection_probability(E, Delta, Z)
ax.plot(E/Delta, B, label=f'Z = {Z}')
ax.axvline(x=-1, color='k', linestyle='--', alpha=0.3)
ax.axvline(x=1, color='k', linestyle='--', alpha=0.3)
ax.set_xlabel('E / Δ')
ax.set_ylabel('Normal Reflection Probability B(E)')
ax.set_title('(b) Normal Reflection vs Barrier Strength')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(-2, 2)
ax.set_ylim(0, 1.05)
# (c) コンダクタンス vs バイアス電圧
ax = axes[1, 0]
for Z in Z_values:
G = conductance_ns_junction(E, Delta, Z)
ax.plot(E/Delta, G, label=f'Z = {Z}')
ax.axvline(x=-1, color='k', linestyle='--', alpha=0.3)
ax.axvline(x=1, color='k', linestyle='--', alpha=0.3)
ax.axhline(y=1, color='k', linestyle=':', alpha=0.3)
ax.set_xlabel('eV / Δ')
ax.set_ylabel('G$_{NS}$ / G$_N$')
ax.set_title('(c) NS Junction Conductance')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(-2, 2)
# (d) ゼロバイアスコンダクタンス vs Z
Z_range = np.linspace(0, 10, 200)
G_zero = np.array([conductance_ns_junction(0, Delta, Z) for Z in Z_range])
ax = axes[1, 1]
ax.plot(Z_range, G_zero, 'b-', linewidth=2)
ax.axhline(y=2, color='r', linestyle='--', label='Perfect Andreev (Z=0)')
ax.axhline(y=1, color='g', linestyle='--', label='Normal transmission')
ax.set_xlabel('Barrier Strength Z')
ax.set_ylabel('G(V=0) / G$_N$')
ax.set_title('(d) Zero-Bias Conductance vs Barrier')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('andreev_reflection_btk.png', dpi=300, bbox_inches='tight')
plt.show()
print("BTK理論の主要な結果:")
print(f"Z = 0 (透明界面): G(V=0)/G_N = {conductance_ns_junction(0, Delta, 0):.2f}")
print(f"Z = 1 (中程度障壁): G(V=0)/G_N = {conductance_ns_junction(0, Delta, 1):.2f}")
print(f"Z → ∞ (トンネル極限): G(V=0)/G_N → 1")3.3 近接効果
近接効果は、超伝導体に接する常伝導金属中にCooper対が染み出し、常伝導体に超伝導相関が誘起される現象です。この効果は以下のUsadelまたはEilenberger方程式で記述されます。
拡散極限でのUsadel方程式
平均自由行程 $\ell \ll \xi$ の汚れた極限では：
$$ \hbar D \nabla \left( \hat{g} \nabla \hat{g} \right) + \left[ i\omega_n \hat{\tau}_3 + \hat{\Delta}, \hat{g} \right] = 0 $$ここで $\hat{g}$ はGreen関数、$D$ は拡散定数、$\hat{\tau}_3$ はPauli行列です。境界条件として、SN界面では：
$$ \xi_N \left. \frac{\partial \theta}{\partial x} \right|_{x=0} = -\frac{R_b}{\rho_N \xi_N} \sin\theta(0) $$ここで $\theta$ は超伝導位相、$R_b$ は界面抵抗、$\xi_N = \sqrt{\hbar D / 2\pi k_B T}$ は常伝導側のコヒーレンス長です。
近接効果の特徴的長さスケール
- 清浄極限 ($\ell \gg \xi$): 誘起ギャップは $\sim \hbar v_F / L$ の範囲で振動しながら減衰
- 拡散極限 ($\ell \ll \xi$): 誘起ギャップは $\xi_N = \sqrt{\hbar D / 2\pi k_B T}$ で単調減衰
- 温度依存性: $\xi_N \propto T^{-1/2}$ なので、低温ほど近接効果は長距離に及ぶ
4. SNSジャンクションとJosephson効果
4.1 SNSジャンクションの物理
超伝導体-常伝導体-超伝導体 (SNS) ジャンクションでは、2つの超伝導体間の位相差 $\phi$ に依存した超電流が流れます。常伝導領域の長さ $L$ とコヒーレンス長 $\xi$ の比によって、短接合 ($L \ll \xi$) と長接合 ($L \gg \xi$) の2つの極限があります。
短接合極限
$L \ll \xi_N$ では、電流-位相関係は正弦的です：
$$ I_c(\phi) = I_c \sin\phi $$臨界電流は温度に強く依存：
$$ I_c(T) = \frac{\pi \Delta(T)}{2 e R_N} \tanh\left(\frac{\Delta(T)}{2 k_B T}\right) $$長接合極限（拡散型）
$L \gg \xi_N$ では、電流-位相関係は非正弦的になります。Kulik-Omelyanchuk理論によると：
$$ I(\phi, T) = \frac{2\pi k_B T}{e R_N} \sum_{n=0}^{\infty} \frac{\sin\phi}{\sqrt{\omega_n^2 + \Delta^2 \sin^2(\phi/2)}} $$Python実装例: SNSジャンクションの電流-位相関係
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
class SNSJunction:
"""SNSジャンクションシミュレーター"""
def __init__(self, L, xi_N, R_N, Delta_0):
"""
Parameters:
-----------
L : float
常伝導領域の長さ (nm)
xi_N : float
常伝導コヒーレンス長 (nm)
R_N : float
常伝導抵抗 (Ω)
Delta_0 : float
超伝導ギャップ (meV)
"""
self.L = L
self.xi_N = xi_N
self.R_N = R_N
self.Delta_0 = Delta_0
self.ratio = L / xi_N
def critical_current_short(self, T, T_c):
"""短接合の臨界電流（解析式）
Parameters:
-----------
T : float
温度 (K)
T_c : float
臨界温度 (K)
"""
# BCS温度依存性
if T >= T_c:
return 0.0
Delta_T = self.Delta_0 * np.tanh(1.74 * np.sqrt(T_c/T - 1))
k_B = 8.617e-5 # eV/K
e = 1.602e-19 # C
I_c = (np.pi * Delta_T * 1e-3) / (2 * e * self.R_N) * \
np.tanh(Delta_T * 1e-3 / (2 * k_B * T))
return I_c
def current_phase_relation_short(self, phi, T, T_c):
"""短接合の電流-位相関係（正弦的）"""
I_c = self.critical_current_short(T, T_c)
return I_c * np.sin(phi)
def current_phase_relation_long(self, phi, T, T_c, n_matsubara=50):
"""長接合の電流-位相関係（Kulik-Omelyanchuk）
Parameters:
-----------
phi : float or array
位相差
T : float
温度 (K)
T_c : float
臨界温度 (K)
n_matsubara : int
Matsubara和の項数
"""
if T >= T_c:
return 0.0
k_B = 8.617e-5 # eV/K
e = 1.602e-19 # C
Delta_T = self.Delta_0 * 1e-3 * np.tanh(1.74 * np.sqrt(T_c/T - 1))
# Matsubara周波数
omega_n = np.pi * k_B * T * (2 * np.arange(n_matsubara) + 1)
# 位相依存性
phi_array = np.atleast_1d(phi)
I = np.zeros_like(phi_array)
for i, ph in enumerate(phi_array):
summand = np.sin(ph) / np.sqrt(omega_n**2 + Delta_T**2 * np.sin(ph/2)**2)
I[i] = (2 * np.pi * k_B * T / (e * self.R_N)) * np.sum(summand)
return I if phi_array.size > 1 else I[0]
def plot_cpr(self, T, T_c):
"""電流-位相関係のプロット"""
phi = np.linspace(0, 2*np.pi, 200)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# (a) 短接合と長接合の比較
ax = axes[0]
I_short = self.current_phase_relation_short(phi, T, T_c)
I_long = np.array([self.current_phase_relation_long(p, T, T_c) for p in phi])
ax.plot(phi, I_short * 1e9, 'b-', linewidth=2, label='Short junction (L << ξ)')
ax.plot(phi, I_long * 1e9, 'r--', linewidth=2, label='Long junction (L >> ξ)')
ax.set_xlabel('Phase difference φ (rad)')
ax.set_ylabel('Supercurrent I (nA)')
ax.set_title(f'(a) Current-Phase Relation (T = {T} K)')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xticks([0, np.pi/2, np.pi, 3*np.pi/2, 2*np.pi])
ax.set_xticklabels(['0', 'π/2', 'π', '3π/2', '2π'])
# (b) 温度依存性
ax = axes[1]
T_range = np.linspace(0.1, T_c*0.95, 20)
I_c_short = [self.critical_current_short(t, T_c) * 1e9 for t in T_range]
phi_max = np.pi/2 # 最大電流の位相
I_c_long = [self.current_phase_relation_long(phi_max, t, T_c) * 1e9
for t in T_range]
ax.plot(T_range, I_c_short, 'b-', linewidth=2, label='Short junction')
ax.plot(T_range, I_c_long, 'r--', linewidth=2, label='Long junction')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Critical Current (nA)')
ax.set_title('(b) Temperature Dependence of I$_c$')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
return fig
# 使用例
if __name__ == "__main__":
# パラメータ設定（典型的なアルミニウムSNSジャンクション）
L = 100 # nm
xi_N = 500 # nm (T依存、ここでは代表値)
R_N = 100 # Ω
Delta_0 = 0.18 # meV (Al)
T_c = 1.2 # K (Al)
T = 0.3 # K
# 短接合の例
print("=== 短接合 (L << ξ_N) ===")
sns_short = SNSJunction(L=50, xi_N=500, R_N=R_N, Delta_0=Delta_0)
I_c_short = sns_short.critical_current_short(T, T_c)
print(f"L/ξ_N = {sns_short.ratio:.3f}")
print(f"臨界電流 I_c = {I_c_short*1e9:.2f} nA")
fig1 = sns_short.plot_cpr(T, T_c)
plt.savefig('sns_short_junction.png', dpi=300, bbox_inches='tight')
# 長接合の例
print("\n=== 長接合 (L >> ξ_N) ===")
sns_long = SNSJunction(L=2000, xi_N=500, R_N=R_N, Delta_0=Delta_0)
print(f"L/ξ_N = {sns_long.ratio:.3f}")
fig2 = sns_long.plot_cpr(T, T_c)
plt.savefig('sns_long_junction.png', dpi=300, bbox_inches='tight')
plt.show()4.2 π接合とスピン三重項超伝導
通常のSNSジャンクションでは位相差0で電流が最大（0接合）ですが、特定の条件下では位相πで電流が最大となる「π接合」が実現します。π接合の実現方法：
- 強磁性層の挿入: SFS接合で交換相互作用により位相シフト
- スピン軌道相互作用: Rashba型SOIによる異常Josephson効果
- トポロジカル超伝導: Majorana束縛状態による4π周期性
5. 超伝導量子ドットとYSR状態
5.1 Anderson不純物モデルと超伝導
超伝導体中または超伝導体に結合した量子ドット（磁性不純物）は、Anderson不純物モデルで記述されます：
$$ H = H_{\text{BCS}} + \epsilon_d d^{\dagger} d + U n_{\uparrow} n_{\downarrow} + \sum_k V_k (c_k^{\dagger} d + d^{\dagger} c_k) $$ここで $\epsilon_d$ は量子ドットの準位、$U$ はCoulomb相互作用、$V_k$ はトンネル結合です。
5.2 Yu-Shiba-Rusinov (YSR) 状態
磁性不純物（スピン $S$）が超伝導体中に置かれると、超伝導ギャップ内に準粒子束縛状態が形成されます。これがYu-Shiba-Rusinov状態です。
YSR状態のエネルギー
Born近似では、YSR状態のエネルギーは：
$$ E_{\text{YSR}} = \pm \Delta \frac{1 - \alpha^2}{1 + \alpha^2} $$ここで $\alpha = \pi \nu_0 J S / 2$ は無次元化された交換相互作用強度、$\nu_0$ は常伝導状態密度、$J$ は交換結合定数です。
YSR状態の量子相転移
$\alpha$ を変化させると、$\alpha = 1$ で量子相転移が起こります：
- $\alpha < 1$ (弱結合): 一重項基底状態、正のエネルギーにYSR状態
- $\alpha > 1$ (強結合): 二重項基底状態、負のエネルギーにYSR状態
- $\alpha = 1$: 零エネルギー束縛状態（量子相転移点）
Python実装例: YSR状態の計算
import numpy as np
import matplotlib.pyplot as plt
def ysr_energy(alpha):
"""YSR状態のエネルギー（Δで規格化）
Parameters:
-----------
alpha : float or array
無次元交換相互作用強度 πν₀JS/2
Returns:
--------
E_YSR : float or array
YSR準位エネルギー（単位：Δ）
"""
alpha = np.atleast_1d(alpha)
E_YSR = np.zeros_like(alpha)
# 正のエネルギー分岐
mask_weak = alpha < 1
E_YSR[mask_weak] = (1 - alpha[mask_weak]**2) / (1 + alpha[mask_weak]**2)
# 負のエネルギー分岐
mask_strong = alpha >= 1
E_YSR[mask_strong] = -(1 - alpha[mask_strong]**2) / (1 + alpha[mask_strong]**2)
return E_YSR if alpha.size > 1 else E_YSR[0]
def ysr_wavefunction(r, xi, alpha, Delta):
"""YSR状態の空間波動関数
Parameters:
-----------
r : float or array
不純物からの距離
xi : float
コヒーレンス長
alpha : float
交換相互作用強度
Delta : float
超伝導ギャップ
Returns:
--------
psi : complex array
YSR波動関数
"""
E_YSR = ysr_energy(alpha) * Delta
kappa = np.sqrt(Delta**2 - E_YSR**2) / (const.hbar * const.v_F)
# 球対称1/r × exp(-κr)型の減衰
psi = np.exp(-kappa * r) / r
return psi
# YSRエネルギーのプロット
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# (a) YSRエネルギー vs 交換相互作用
alpha_range = np.linspace(0.1, 3, 300)
E_YSR = ysr_energy(alpha_range)
ax = axes[0, 0]
ax.plot(alpha_range, E_YSR, 'b-', linewidth=2)
ax.axhline(y=0, color='k', linestyle='--', alpha=0.5)
ax.axvline(x=1, color='r', linestyle='--', alpha=0.5, label='QPT at α=1')
ax.fill_between(alpha_range, -1, 1, alpha=0.1, color='gray', label='SC gap')
ax.set_xlabel('Exchange coupling α = πν₀JS/2')
ax.set_ylabel('E$_{YSR}$ / Δ')
ax.set_title('(a) YSR State Energy vs Coupling Strength')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_ylim(-1.2, 1.2)
# (b) 基底状態の相図
ax = axes[0, 1]
alpha_phase = np.linspace(0, 2, 500)
phase = np.where(alpha_phase < 1, 0, 1) # 0: singlet, 1: doublet
ax.fill_between(alpha_phase, 0, 1, where=(phase==0),
alpha=0.3, color='blue', label='Singlet (|BCS⟩)')
ax.fill_between(alpha_phase, 0, 1, where=(phase==1),
alpha=0.3, color='red', label='Doublet (|S=1/2⟩)')
ax.axvline(x=1, color='k', linestyle='--', linewidth=2, label='QPT')
ax.set_xlabel('Exchange coupling α')
ax.set_ylabel('Ground State')
ax.set_title('(b) Quantum Phase Diagram')
ax.set_yticks([0.25, 0.75])
ax.set_yticklabels(['Singlet', 'Doublet'])
ax.legend(loc='upper right')
ax.grid(True, alpha=0.3, axis='x')
# (c) YSR状態密度
energy = np.linspace(-1.5, 1.5, 500)
eta = 0.02 # ブロードニング
ax = axes[1, 0]
for alpha_val in [0.3, 0.7, 1.0, 1.5, 2.0]:
E_YSR_val = ysr_energy(alpha_val)
dos = (eta / np.pi) / ((energy - E_YSR_val)**2 + eta**2)
ax.plot(energy, dos, label=f'α = {alpha_val}')
ax.axvline(x=-1, color='k', linestyle=':', alpha=0.3)
ax.axvline(x=1, color='k', linestyle=':', alpha=0.3)
ax.axvline(x=0, color='r', linestyle='--', alpha=0.5)
ax.set_xlabel('Energy / Δ')
ax.set_ylabel('Density of States (arb. units)')
ax.set_title('(c) YSR State Density of States')
ax.legend()
ax.grid(True, alpha=0.3)
# (d) YSR状態数 vs 不純物濃度（簡単なモデル）
n_imp = np.logspace(-4, -1, 100) # 不純物濃度
# 平均準位間隔 δ ∼ 1/(ν₀ V)
# 多数不純物でYSRバンド形成
delta_E = 2 * np.abs(ysr_energy(0.5)) # YSR準位幅
n_states = n_imp / (1e-3) # 簡単化したモデル
ax = axes[1, 1]
ax.loglog(n_imp, n_states, 'b-', linewidth=2)
ax.set_xlabel('Impurity concentration n$_{imp}$')
ax.set_ylabel('Number of YSR states in gap')
ax.set_title('(d) YSR Band Formation (simplified)')
ax.grid(True, alpha=0.3, which='both')
plt.tight_layout()
plt.savefig('ysr_states_analysis.png', dpi=300, bbox_inches='tight')
plt.show()
# 重要な物理的結果の出力
print("=== YSR状態の重要な特性 ===")
print(f"\n弱結合 (α=0.5):")
print(f" E_YSR/Δ = +{ysr_energy(0.5):.3f} (正エネルギー)")
print(f" 基底状態: 一重項 |BCS⟩")
print(f"\n量子相転移点 (α=1.0):")
print(f" E_YSR/Δ = {ysr_energy(1.0):.3f} (零エネルギー)")
print(f" 量子相転移: 一重項 ↔ 二重項")
print(f"\n強結合 (α=1.5):")
print(f" E_YSR/Δ = {ysr_energy(1.5):.3f} (負エネルギー)")
print(f" 基底状態: 二重項 |S=1/2⟩")6. トポロジカル超伝導とMajorana束縛状態
6.1 トポロジカル超伝導の基礎
トポロジカル超伝導体は、バルクにギャップがあり、エッジや界面に保護された零エネルギー状態（Majorana束縛状態）を持つ超伝導体です。最も単純なモデルは1次元Kitaev鎖です。
Kitaev鎖モデル
スピンレスフェルミオンの1次元鎖：
$$ H = -\mu \sum_j c_j^{\dagger} c_j - t \sum_j (c_{j+1}^{\dagger} c_j + \text{h.c.}) + \Delta \sum_j (c_{j+1}^{\dagger} c_j^{\dagger} + \text{h.c.}) $$ここで $\mu$ は化学ポテンシャル、$t$ はホッピング、$\Delta$ はp波超伝導ギャップです。
Kitaev鎖の位相相転移
化学ポテンシャル $\mu$ を変化させると、系は位相的トリビアル相とトポロジカル相の間を転移します：
- トリビアル相: $|\mu| > 2t$、エッジ状態なし
- トポロジカル相: $|\mu| < 2t$、エッジにMajorana零モード
位相不変量（Pfaffian数）で分類：$\nu = \text{sgn}[\Delta/t] \in \{-1, +1\}$
6.2 Majorana束縛状態
トポロジカル超伝導体の端や磁束渦中心には、Majorana束縛状態（MBS）が出現します。Majorana演算子は以下の性質を持ちます：
$$ \gamma = \gamma^{\dagger}, \quad \{\gamma_i, \gamma_j\} = 2\delta_{ij} $$2つのMajorana演算子から1つの通常のフェルミオンが構成されます：
$$ c = \frac{\gamma_1 + i\gamma_2}{2}, \quad c^{\dagger} = \frac{\gamma_1 - i\gamma_2}{2} $$Python実装例: Kitaev鎖とMajorana零モード
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh
class KitaevChain:
"""1次元Kitaev鎖のシミュレーター"""
def __init__(self, N_sites, t, mu, Delta):
"""
Parameters:
-----------
N_sites : int
格子サイト数
t : float
ホッピング積分
mu : float
化学ポテンシャル
Delta : float
p波超伝導ギャップ
"""
self.N = N_sites
self.t = t
self.mu = mu
self.Delta = Delta
def build_hamiltonian(self, periodic=False):
"""Kitaev鎖のBdGハミルトニアンを構築
Parameters:
-----------
periodic : bool
周期境界条件の有無
Returns:
--------
H : ndarray (2N, 2N)
BdGハミルトニアン
"""
N = self.N
H = np.zeros((2*N, 2*N), dtype=complex)
# サイトエネルギー
for i in range(N):
H[i, i] = -self.mu
H[N+i, N+i] = self.mu
# ホッピングとp波ペアリング
for i in range(N-1):
# 粒子ホッピング
H[i, i+1] = -self.t
H[i+1, i] = -self.t
# 正孔ホッピング
H[N+i, N+i+1] = self.t
H[N+i+1, N+i] = self.t
# p波ペアリング（異常項）
H[i, N+i+1] = self.Delta
H[i+1, N+i] = self.Delta
H[N+i, i+1] = np.conj(self.Delta)
H[N+i+1, i] = np.conj(self.Delta)
# 周期境界条件
if periodic:
# N-1サイトと0サイトを結合
H[N-1, 0] = -self.t
H[0, N-1] = -self.t
H[2*N-1, N] = self.t
H[N, 2*N-1] = self.t
H[N-1, N] = self.Delta
H[0, 2*N-1] = self.Delta
H[2*N-1, 0] = np.conj(self.Delta)
H[N, N-1] = np.conj(self.Delta)
return H
def solve(self, periodic=False):
"""ハミルトニアンを対角化
Returns:
--------
energies : ndarray
準粒子エネルギー
wavefunctions : ndarray
Nambu基底での波動関数
"""
H = self.build_hamiltonian(periodic)
energies, wavefunctions = eigh(H)
return energies, wavefunctions
def check_topological_phase(self):
"""位相相を判定
Returns:
--------
phase : str
'Topological' または 'Trivial'
"""
# 簡単な判定: |μ| < 2t かつ Δ ≠ 0
if abs(self.mu) < 2 * abs(self.t) and abs(self.Delta) > 1e-10:
return 'Topological'
else:
return 'Trivial'
def calculate_majorana_polarization(self, energies, wavefunctions,
zero_threshold=1e-3):
"""Majorana零モードの空間分布を計算
Parameters:
-----------
energies : ndarray
準粒子エネルギー
wavefunctions : ndarray
波動関数
zero_threshold : float
零エネルギーとみなす閾値
Returns:
--------
majorana_weight : ndarray
各サイトでのMajorana重み
"""
N = self.N
zero_modes_indices = np.where(np.abs(energies) < zero_threshold)[0]
if len(zero_modes_indices) == 0:
return np.zeros(N)
# 零エネルギーモードの粒子成分
majorana_weight = np.zeros(N)
for idx in zero_modes_indices:
psi = wavefunctions[:N, idx] # 粒子成分
majorana_weight += np.abs(psi)**2
# 正規化
majorana_weight /= np.sum(majorana_weight)
return majorana_weight
# 位相図の作成
def plot_phase_diagram():
"""Kitaev鎖の位相図"""
mu_range = np.linspace(-3, 3, 100)
Delta_range = np.linspace(0, 2, 100)
t = 1.0 # 固定
phase_diagram = np.zeros((len(Delta_range), len(mu_range)))
for i, Delta in enumerate(Delta_range):
for j, mu in enumerate(mu_range):
# 簡易的位相判定: |μ| < 2t かつ Δ ≠ 0
if abs(mu) < 2*t and abs(Delta) > 1e-10:
phase_diagram[i, j] = 1 # Topological
else:
phase_diagram[i, j] = 0 # Trivial
fig, ax = plt.subplots(figsize=(10, 6))
im = ax.imshow(phase_diagram, extent=[mu_range[0], mu_range[-1],
Delta_range[0], Delta_range[-1]],
origin='lower', aspect='auto', cmap='RdBu_r', alpha=0.7)
# 相境界を強調
ax.axvline(x=-2*t, color='k', linestyle='--', linewidth=2, label='Phase boundary')
ax.axvline(x=2*t, color='k', linestyle='--', linewidth=2)
ax.set_xlabel('Chemical potential μ / t', fontsize=12)
ax.set_ylabel('Pairing Δ / t', fontsize=12)
ax.set_title('Kitaev Chain Phase Diagram', fontsize=14)
# カラーバー
cbar = plt.colorbar(im, ax=ax)
cbar.set_ticks([0, 1])
cbar.set_ticklabels(['Trivial', 'Topological'])
ax.text(-1, 1.5, 'Topological\n(MBS at edges)',
fontsize=11, ha='center', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
ax.text(2.5, 1.5, 'Trivial\n(No edge states)',
fontsize=11, ha='center', bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.5))
plt.tight_layout()
return fig
# メイン計算
if __name__ == "__main__":
# パラメータ
N_sites = 100
t = 1.0
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (a) トポロジカル相のスペクトル
kitaev_topo = KitaevChain(N_sites, t, mu=0.5, Delta=1.0)
energies_topo, wf_topo = kitaev_topo.solve(periodic=False)
ax = axes[0, 0]
ax.plot(energies_topo, 'o', markersize=3)
ax.axhline(y=0, color='r', linestyle='--', alpha=0.5)
ax.set_xlabel('State index')
ax.set_ylabel('Energy / t')
ax.set_title(f'(a) Topological Phase (μ={kitaev_topo.mu}t, Δ={kitaev_topo.Delta}t)')
ax.grid(True, alpha=0.3)
ax.set_ylim(-2.5, 2.5)
# (b) トリビアル相のスペクトル
kitaev_triv = KitaevChain(N_sites, t, mu=2.5, Delta=1.0)
energies_triv, wf_triv = kitaev_triv.solve(periodic=False)
ax = axes[0, 1]
ax.plot(energies_triv, 'o', markersize=3, color='orange')
ax.axhline(y=0, color='r', linestyle='--', alpha=0.5)
ax.set_xlabel('State index')
ax.set_ylabel('Energy / t')
ax.set_title(f'(b) Trivial Phase (μ={kitaev_triv.mu}t, Δ={kitaev_triv.Delta}t)')
ax.grid(True, alpha=0.3)
ax.set_ylim(-2.5, 2.5)
# (c) Majorana零モードの空間分布
majorana_weight = kitaev_topo.calculate_majorana_polarization(
energies_topo, wf_topo, zero_threshold=1e-2)
ax = axes[1, 0]
ax.plot(range(N_sites), majorana_weight, 'b-', linewidth=2)
ax.set_xlabel('Site index')
ax.set_ylabel('Majorana weight |ψ|²')
ax.set_title('(c) Majorana Zero Mode Spatial Distribution')
ax.grid(True, alpha=0.3)
ax.set_yscale('log')
# (d) ギャップクロージング（位相転移）
mu_scan = np.linspace(-3, 3, 200)
min_gap = []
for mu in mu_scan:
kitaev_scan = KitaevChain(N_sites, t, mu, Delta=1.0)
energies_scan, _ = kitaev_scan.solve(periodic=True) # PBC
# バルクギャップ（最小正エネルギー）
positive_E = energies_scan[energies_scan > 0]
if len(positive_E) > 0:
min_gap.append(np.min(positive_E))
else:
min_gap.append(0)
ax = axes[1, 1]
ax.plot(mu_scan, min_gap, 'g-', linewidth=2)
ax.axvline(x=-2*t, color='r', linestyle='--', label='Phase boundary')
ax.axvline(x=2*t, color='r', linestyle='--')
ax.fill_between(mu_scan, 0, 3,
where=(np.abs(mu_scan) < 2*t),
alpha=0.2, color='blue', label='Topological phase')
ax.set_xlabel('Chemical potential μ / t')
ax.set_ylabel('Bulk gap / t')
ax.set_title('(d) Gap Closing at Phase Transition')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 2.5)
plt.tight_layout()
plt.savefig('kitaev_chain_majorana.png', dpi=300, bbox_inches='tight')
# 位相図
fig_phase = plot_phase_diagram()
plt.savefig('kitaev_phase_diagram.png', dpi=300, bbox_inches='tight')
plt.show()
# 結果の要約
print("=== Kitaev鎖の計算結果 ===")
print(f"\nトポロジカル相 (μ={kitaev_topo.mu}t):")
print(f" 相: {kitaev_topo.check_topological_phase()}")
zero_E_topo = energies_topo[np.abs(energies_topo) < 1e-2]
print(f" 零エネルギーモード数: {len(zero_E_topo)}")
if len(zero_E_topo) > 0:
print(f" 零モードエネルギー: {zero_E_topo}")
print(f"\nトリビアル相 (μ={kitaev_triv.mu}t):")
print(f" 相: {kitaev_triv.check_topological_phase()}")
zero_E_triv = energies_triv[np.abs(energies_triv) < 1e-2]
print(f" 零エネルギーモード数: {len(zero_E_triv)}")6.3 半導体-超伝導体ハイブリッド系
実験的には、以下の系でMajorana束縛状態が追求されています：
- 半導体ナノワイヤー + s波超伝導: InSb/InAsナノワイヤーとAl超伝導体
- トポロジカル絶縁体 + 超伝導: Bi₂Se₃表面状態と超伝導近接効果
- 鉄系超伝導体: FeTe₀.₅Se₀.₅のトポロジカル表面状態
半導体ナノワイヤー系の有効ハミルトニアン：
$$ H = \frac{p^2}{2m^*} - \mu + \alpha (\sigma_x p_y - \sigma_y p_x) + \frac{1}{2} g \mu_B B \sigma_z + \Delta \sigma_y \tau_y $$ここで $\alpha$ はRashbaスピン軌道相互作用、$B$ はZeeman磁場、$\Delta$ は誘起された超伝導ギャップです。
7. BdG方程式の数値解法
7.1 実空間離散化とタイトバインディング
実空間でのBdG方程式の数値解法は、一般に以下のステップで行います：
- 空間離散化: 連続空間を格子点に離散化（格子定数 $a$）
- 差分近似: 微分を有限差分で近似 $\nabla^2 \to (V_{i+1} - 2V_i + V_{i-1})/a^2$
- 行列構築: $2N \times 2N$ エルミート行列（$N$ は格子点数）
- 対角化: エルミート行列の固有値問題を解く
- 自己無撞着: ギャップ関数を更新し収束まで反復
7.2 自己無撞着計算アルゴリズム
Python実装例: 自己無撞着BdG計算
import numpy as np
from scipy.linalg import eigh
import matplotlib.pyplot as plt
class SelfConsistentBdG:
"""自己無撞着BdG方程式ソルバー"""
def __init__(self, N_sites, t, U, T, mu=0):
"""
Parameters:
-----------
N_sites : int
格子点数
t : float
ホッピング積分 (eV)
U : float
引力相互作用強度 (eV, 負の値)
T : float
温度 (K)
mu : float
化学ポテンシャル (eV)
"""
self.N = N_sites
self.t = t
self.U = U # 引力: U < 0
self.T = T
self.mu = mu
self.k_B = 8.617e-5 # eV/K
# 初期ギャップ（小さなランダム値）
self.Delta = 0.1 * np.random.rand(N_sites)
def build_hamiltonian(self, Delta):
"""BdGハミルトニアンを構築"""
N = self.N
H = np.zeros((2*N, 2*N), dtype=complex)
# サイトエネルギー
for i in range(N):
H[i, i] = -self.mu
H[N+i, N+i] = self.mu
# ホッピング
for i in range(N-1):
H[i, i+1] = -self.t
H[i+1, i] = -self.t
H[N+i, N+i+1] = self.t
H[N+i+1, N+i] = self.t
# ギャップ関数
for i in range(N):
H[i, N+i] = Delta[i]
H[N+i, i] = np.conj(Delta[i])
return H
def fermi_function(self, E):
"""Fermi-Dirac分布関数"""
if self.T < 1e-6:
# T=0極限
return np.where(E < 0, 1.0, 0.0)
else:
beta = 1 / (self.k_B * self.T)
# 数値安定性のための処理
x = beta * E
return 1 / (1 + np.exp(np.clip(x, -500, 500)))
def update_gap(self, energies, u_vectors, v_vectors):
"""ギャップ関数を自己無撞着に更新
Parameters:
-----------
energies : ndarray
準粒子エネルギー
u_vectors, v_vectors : ndarray
Bogoliubov波動関数
Returns:
--------
Delta_new : ndarray
更新されたギャップ関数
"""
N = self.N
Delta_new = np.zeros(N, dtype=complex)
for i in range(N):
sum_val = 0
for n in range(2*N):
# 負エネルギー状態のみを足す（基底状態からの寄与）
if energies[n] < 0:
sum_val += u_vectors[i, n] * np.conj(v_vectors[i, n])
# 有限温度の場合
elif self.T > 1e-6:
f = 1 - self.fermi_function(energies[n])
sum_val += u_vectors[i, n] * np.conj(v_vectors[i, n]) * f
Delta_new[i] = -self.U * sum_val
return Delta_new
def solve_self_consistent(self, max_iter=100, tol=1e-6, mixing=0.5):
"""自己無撞着ループ
Parameters:
-----------
max_iter : int
最大反復回数
tol : float
収束判定の閾値
mixing : float
混合パラメータ（0-1、小さいほど安定）
Returns:
--------
converged : bool
収束したかどうか
Delta_history : list
各反復でのギャップ関数
"""
Delta_history = []
for iteration in range(max_iter):
# BdGハミルトニアンを解く
H = self.build_hamiltonian(self.Delta)
energies, eigvecs = eigh(H)
N = self.N
u_vecs = eigvecs[:N, :]
v_vecs = eigvecs[N:, :]
# ギャップ関数を更新
Delta_new = self.update_gap(energies, u_vecs, v_vecs)
# 混合（収束安定化）
Delta_mixed = mixing * Delta_new + (1 - mixing) * self.Delta
# 収束判定
diff = np.max(np.abs(Delta_mixed - self.Delta))
Delta_history.append(np.copy(self.Delta))
if diff < tol:
print(f"収束しました（{iteration+1}回の反復）")
self.Delta = Delta_mixed
return True, Delta_history
self.Delta = Delta_mixed
if (iteration + 1) % 10 == 0:
print(f"反復 {iteration+1}: 差分 = {diff:.6e}")
print(f"警告: {max_iter}回の反復で収束しませんでした")
return False, Delta_history
# 使用例
if __name__ == "__main__":
# パラメータ設定
N_sites = 50
t = 1.0 # eV
U = -1.5 # eV (引力)
T = 0.01 # K (ほぼゼロ温度)
# ソルバー初期化
solver = SelfConsistentBdG(N_sites, t, U, T, mu=0)
# 自己無撞着計算
print("自己無撞着BdG計算を開始...")
converged, Delta_history = solver.solve_self_consistent(
max_iter=100, tol=1e-6, mixing=0.3)
# 最終結果の計算
H_final = solver.build_hamiltonian(solver.Delta)
energies_final, eigvecs_final = eigh(H_final)
# 結果のプロット
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (a) ギャップ関数の空間分布
ax = axes[0, 0]
ax.plot(range(N_sites), np.abs(solver.Delta), 'b-', linewidth=2)
ax.set_xlabel('Site index')
ax.set_ylabel('|Δ(i)| (eV)')
ax.set_title('(a) Self-Consistent Gap Function')
ax.grid(True, alpha=0.3)
# (b) 準粒子スペクトル
ax = axes[0, 1]
ax.plot(energies_final, 'o', markersize=3)
ax.axhline(y=0, color='r', linestyle='--', alpha=0.5)
ax.set_xlabel('State index')
ax.set_ylabel('Energy (eV)')
ax.set_title('(b) Quasiparticle Spectrum')
ax.grid(True, alpha=0.3)
# (c) 収束過程
ax = axes[1, 0]
Delta_max_history = [np.max(np.abs(D)) for D in Delta_history]
ax.plot(Delta_max_history, 'g-', linewidth=2)
ax.set_xlabel('Iteration')
ax.set_ylabel('max|Δ| (eV)')
ax.set_title('(c) Convergence of Gap Function')
ax.grid(True, alpha=0.3)
# (d) 状態密度
ax = axes[1, 1]
E_range = np.linspace(-3, 3, 500)
dos = np.zeros_like(E_range)
eta = 0.05 # ブロードニング
for E_n in energies_final:
dos += (eta/np.pi) / ((E_range - E_n)**2 + eta**2)
ax.plot(E_range, dos, 'r-', linewidth=2)
ax.axvline(x=np.max(np.abs(solver.Delta)), color='b',
linestyle='--', label=f'Δ = {np.max(np.abs(solver.Delta)):.3f} eV')
ax.axvline(x=-np.max(np.abs(solver.Delta)), color='b', linestyle='--')
ax.set_xlabel('Energy (eV)')
ax.set_ylabel('Density of States')
ax.set_title('(d) Self-Consistent DOS')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('self_consistent_bdg.png', dpi=300, bbox_inches='tight')
plt.show()
print(f"\n最終結果:")
print(f"最大ギャップ: {np.max(np.abs(solver.Delta)):.4f} eV")
print(f"平均ギャップ: {np.mean(np.abs(solver.Delta)):.4f} eV")
print(f"最小準粒子エネルギー: {np.min(np.abs(energies_final)):.4f} eV")7.3 大規模系への拡張と並列計算
実際の研究では、3次元系や不純物を多数含む系など、より大規模な計算が必要です。効率化のテクニック：
- 疎行列手法: scipy.sparse を用いたメモリ効率化
- Krylov部分空間法: 全ての固有値ではなく、低エネルギー状態のみを求める
- 並列化: MPI/OpenMPによるマルチコア計算
- GPU加速: CuPy/PyTorchによる大規模行列計算
8. 量子コンピューティングへの応用
8.1 トポロジカル量子ビット
Majorana束縛状態は、フォールトトレラント量子計算のためのトポロジカル量子ビットの候補です。2つのMajorana演算子 $\gamma_1, \gamma_2$ から1つのフェルミオン準位が形成され、これが量子ビットとして機能します。
トポロジカル保護の原理
Majorana量子ビットの優位性：
- 非局所性: 量子情報が空間的に分離したMajorana対に符号化
- トポロジカル保護: 局所的な摂動に対して量子情報が安定
- 非アーベルエニオン統計: ブレイディング操作で量子ゲートを実装
- デコヒーレンス抑制: 環境との結合が指数関数的に抑制
Majoranaブレイディング
Majorana束縛状態の交換（ブレイディング）は非アーベル統計に従います。4つのMajorana演算子 $\gamma_1, \gamma_2, \gamma_3, \gamma_4$ があるとき、$\gamma_1$ と $\gamma_2$ を交換する操作は：
$$ \gamma_1 \to \gamma_2, \quad \gamma_2 \to -\gamma_1, \quad \gamma_3, \gamma_4 \to \text{不変} $$これはユニタリ演算子 $U_{12} = e^{\pi \gamma_1 \gamma_2 / 4}$ で表されます。
8.2 Microsoftのトポロジカル量子ビット
Microsoft Azure Quantumは、InAsナノワイヤー上のMajorana束縛状態を用いたトポロジカル量子ビットの開発を進めています。実験的に以下が観測されています：
- 零バイアスコンダクタンスピーク（Majoranaの証拠候補）
- 磁場・化学ポテンシャル依存性の位相図
- 非局所伝導度測定によるMajorana相関
8.3 課題と展望
トポロジカル量子コンピューティングの実現には、以下の技術的課題があります：
| 課題 | 現状 | 必要な進展 |
|---|---|---|
| Majorana検出の確実性 | 間接的証拠のみ | 確定的な同定手法の確立 |
| 量子ビットの読み出し | 基礎実験段階 | 高忠実度読み出しプロトコル |
| ブレイディング操作 | 理論提案段階 | T字ジャンクション等での実証 |
| スケーラビリティ | 数個のMajorana | 多数量子ビットアレイ |
まとめ
本章では、メゾスコピック超伝導の理論と応用を包括的に学習しました。主要なポイントを振り返ります：
本章の重要事項
- Bogoliubov-de Gennes方程式は、空間変化する超伝導を記述する基礎理論であり、粒子-正孔対称性を持つ
- Andreev反射は、SN界面での電子-正孔変換過程であり、ゼロバイアスコンダクタンスの倍化をもたらす
- 近接効果により、常伝導体中にCooper対相関が誘起され、SNSジャンクションでJosephson効果が発現
- Yu-Shiba-Rusinov状態は、磁性不純物が超伝導ギャップ内に作る束縛状態であり、量子相転移を示す
- Majorana束縛状態は、トポロジカル超伝導体のエッジに現れ、フォールトトレラント量子計算への応用が期待される
- 自己無撞着BdG計算は、実際の超伝導系の準粒子スペクトルとギャップ関数を第一原理的に求める手法
- トポロジカル量子ビットは、Majoranaの非局所性とトポロジカル保護により、デコヒーレンスに強い
次章への接続
第4章「超伝導量子デバイス」では、本章で学んだメゾスコピック物理を実際のデバイス応用に展開します。超伝導量子干渉計（SQUID）、超伝導量子ビット（transmon, fluxoniumなど）、超伝導検出器の動作原理と設計を詳しく学びます。
演習問題
演習1: BdG方程式の粒子-正孔対称性
問題: BdG方程式が粒子-正孔対称性を持つことを示し、この対称性がエネルギースペクトルに与える影響を説明せよ。
ヒント: 粒子-正孔演算子 $\mathcal{C}$ を定義し、BdGハミルトニアンとの交換関係を調べよ。
演習2: BTK理論の数値計算
問題: 本章のBTKコードを用いて、以下を計算せよ：
- 界面障壁 $Z$ を変化させたときのゼロバイアスコンダクタンスの変化
- $Z = 0$ （完全透明界面）での電流-電圧特性
- 有限温度効果を考慮した場合のコンダクタンスの温度依存性
演習3: SNSジャンクションの臨界電流
問題: 短接合SNSジャンクションの臨界電流 $I_c(T)$ が、温度 $T$ の関数としてどのように変化するか、数値的に計算し、BCS理論のギャップ温度依存性 $\Delta(T)$ との関係を議論せよ。
演習4: YSR状態の相転移
問題: YSR状態のエネルギーが $E_{\text{YSR}} = \Delta (1-\alpha^2)/(1+\alpha^2)$ で与えられることを用いて：
- $\alpha = 1$ で零エネルギー状態が実現することを確認せよ
- 基底状態の一重項-二重項転移が起こる物理的機構を説明せよ
- STMトンネル分光でYSR状態をどのように観測できるか議論せよ
演習5: Kitaev鎖の位相不変量（発展問題）
問題: Kitaev鎖のトポロジカル相とトリビアル相を区別する位相不変量（Pfaffian不変量）を計算し、エッジ状態の存在との関係を示せ。周期境界条件と開放境界条件での違いを議論せよ。
演習6: 自己無撞着BdG計算の収束（プログラミング課題）
問題: 本章の自己無撞着BdGコードを拡張し、以下を実装せよ：
- 不純物ポテンシャル $V_{\text{imp}}(i)$ を含む系の計算
- 2次元正方格子への拡張
- 収束性を改善するためのAnderson混合法の実装
演習7: Majoranaブレイディング（理論問題）
問題: 4つのMajorana束縛状態 $\gamma_1, \gamma_2, \gamma_3, \gamma_4$ があるとき：
- $\gamma_1$ と $\gamma_2$ のブレイディング演算子 $U_{12} = e^{\pi \gamma_1 \gamma_2 / 4}$ が、フェルミオン占有数演算子 $n = c^{\dagger} c$ （ただし $c = (\gamma_1 + i\gamma_2)/2$）にどう作用するか計算せよ
- 連続ブレイディング $U_{12} U_{23}$ と $U_{23} U_{12}$ が非可換であることを示し、非アーベル統計の物理的意味を説明せよ
参考文献
- de Gennes, P. G. "Superconductivity of Metals and Alloys" (1966) - BdG方程式の古典的教科書
- Tinkham, M. "Introduction to Superconductivity" 2nd ed. (2004) - メゾスコピック超伝導の包括的解説
- Blonder, G. E., Tinkham, M., & Klapwijk, T. M. "Transition from metallic to tunneling regimes in superconducting microconstrictions: Excess current, charge imbalance, and supercurrent conversion" Phys. Rev. B 25, 4515 (1982) - BTK理論の原著論文
- Buzdin, A. I. "Proximity effects in superconductor-ferromagnet heterostructures" Rev. Mod. Phys. 77, 935 (2005) - 近接効果の詳細レビュー
- Franke, K. J., Schulze, G., & Pascual, J. I. "Competition of Superconducting Phenomena and Kondo Screening at the Nanoscale" Science 332, 940 (2011) - YSR状態のSTM観測
- Kitaev, A. Yu. "Unpaired Majorana fermions in quantum wires" Physics-Uspekhi 44, 131 (2001) - Kitaev鎖の原著論文
- Alicea, J. "New directions in the pursuit of Majorana fermions in solid state systems" Rep. Prog. Phys. 75, 076501 (2012) - Majorana束縛状態の包括的レビュー
- Nayak, C., Simon, S. H., Stern, A., Freedman, M., & Das Sarma, S. "Non-Abelian anyons and topological quantum computation" Rev. Mod. Phys. 80, 1083 (2008) - トポロジカル量子計算の理論的基礎
- Lutchyn, R. M., Bakkers, E. P. A. M., Kouwenhoven, L. P., Krogstrup, P., Marcus, C. M., & Oreg, Y. "Majorana zero modes in superconductor–semiconductor heterostructures" Nat. Rev. Mater. 3, 52 (2018) - 半導体-超伝導ハイブリッド系の最新レビュー
- Microsoft Quantum: "Topological Quantum Computing" (2023) - トポロジカル量子ビットの実験的進展