第1章：Ginzburg-Landau理論

学習目標

現象論的理論としてのGL理論の基本的な考え方を理解する
秩序パラメータψの物理的意味を説明する
GL自由エネルギー汎関数を導出し、その各項の意味を理解する
第一GL方程式と第二GL方程式を導出する
コヒーレンス長ξと侵入長λの定義と物理的解釈を理解する
GLパラメータκ = λ/ξと第一種・第二種超伝導体の分類を説明する
表面エネルギーの符号が超伝導体のタイプを決定することを理解する
Pythonを使ってGL理論の予測を可視化する

1.1 現象論的理論の導入

GL理論の歴史的背景

1950年、ソビエトの物理学者Vitaly GinzburgとLev Landauは、超伝導の現象論的理論を提案しました。この理論は、BCS理論（1957年）が確立される前に開発されたため、微視的な機構を説明するのではなく、超伝導の巨視的な性質を記述することを目的としていました。

現象論的理論とは？

現象論的理論は、観測可能な量の間の関係を記述しますが、必ずしも微視的なメカニズムを説明するわけではありません。GL理論は、超伝導状態を秩序パラメータによって特徴づけ、その空間的および磁場的な振る舞いを予測します。

Landauの相転移理論

GL理論は、Landauの一般的な二次相転移理論に基づいています。この理論では、相転移は秩序パラメータによって特徴づけられます：

通常状態（T > Tc）: 秩序パラメータ = 0
超伝導状態（T < Tc）: 秩序パラメータ ≠ 0

臨界温度Tc近傍では、秩序パラメータは小さく、自由エネルギーをそのべき級数で展開できます。

1.2 秩序パラメータの概念

複素スカラー場ψ(r)

GL理論では、超伝導状態は複素秩序パラメータψ(r)によって記述されます：

ψ(r) = |ψ(r)| exp(iθ(r))

ここで：

|ψ(r)|²は超伝導電子対（クーパー対）の密度nsに比例します
θ(r)は位相で、超伝導電流の流れと関係します
rは空間座標です

BCS理論との関係

後にBCS理論が確立されると、ψは超伝導ギャップ関数Δに関連していることが明らかになりました。しかし、GL理論は微視的な詳細を知らなくても、多くの重要な予測を行うことができます。

規格化

秩序パラメータは通常、バルクの平衡値で規格化されます：

|ψ∞|² = ns0

ここで、ns0は均一な超伝導体における超伝導電子対の密度です。

1.3 Ginzburg-Landau自由エネルギー汎関数

自由エネルギーの構築

GL自由エネルギー汎関数は、通常状態からの超伝導自由エネルギー密度の超過分として書かれます：

Fs - Fn = ∫ dV [ α|ψ|² + (β/2)|ψ|⁴ + (1/2m*)|(−iℏ∇ − q*A)ψ|² + B²/(2μ0) ]

各項の物理的意味

項	意味	備考
`α\|ψ\|²`	秩序パラメータの二次項	α(T) = α0(T - Tc)、α0 > 0
`(β/2)\|ψ\|⁴`	秩序パラメータの四次項	β > 0（安定性のため）
`(1/2m)\|(−iℏ∇ − qA)ψ\|²`	運動エネルギーと電磁エネルギー	mはクーパー対の有効質量、q = 2eは電荷
`B²/(2μ0)`	磁場のエネルギー	B = ∇ × A

温度依存性

係数αは温度に依存し、臨界温度付近で線形に変化します：

α(T) = α0(T - Tc)

T > Tc: α > 0 → ψ = 0 が最小（通常状態）
T < Tc: α < 0 → ψ ≠ 0 が最小（超伝導状態）

Pythonによる自由エネルギーの可視化

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

# Parameters
alpha_0 = 1.0
beta = 1.0
Tc = 10.0  # Critical temperature (arbitrary units)

# Temperature range
temperatures = [12, 10, 8, 5]  # T > Tc, T = Tc, T < Tc

# Order parameter range
psi = np.linspace(-2, 2, 400)

# Plot free energy landscape
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
axes = axes.flatten()

for i, T in enumerate(temperatures):
    alpha = alpha_0 * (T - Tc)

    # Free energy density (excluding field terms)
    f = alpha * psi**2 + (beta/2) * psi**4

    axes[i].plot(psi, f, 'b-', linewidth=2)
    axes[i].axhline(y=0, color='k', linestyle='--', alpha=0.3)
    axes[i].axvline(x=0, color='k', linestyle='--', alpha=0.3)
    axes[i].set_xlabel('秩序パラメータ ψ', fontsize=11)
    axes[i].set_ylabel('自由エネルギー密度 f', fontsize=11)
    axes[i].set_title(f'T = {T} K (α = {alpha:.1f})', fontsize=12)
    axes[i].grid(True, alpha=0.3)

    # Mark minima
    if T > Tc:
        axes[i].plot(0, 0, 'ro', markersize=10, label='最小値（通常状態）')
    elif T == Tc:
        axes[i].plot(0, 0, 'yo', markersize=10, label='最小値（臨界点）')
    else:
        psi_eq = np.sqrt(-alpha / beta)
        f_min = alpha * psi_eq**2 + (beta/2) * psi_eq**4
        axes[i].plot([psi_eq, -psi_eq], [f_min, f_min], 'ro', markersize=10,
                    label='最小値（超伝導状態）')
        axes[i].plot(0, 0, 'go', markersize=8, alpha=0.5, label='局所最大値')

    axes[i].legend(fontsize=9)
    axes[i].set_ylim(-1, 3)

plt.tight_layout()
plt.show()

# Print equilibrium values
print("\n平衡状態の秩序パラメータ:")
for T in temperatures:
    alpha = alpha_0 * (T - Tc)
    if T <= Tc:
        psi_eq = np.sqrt(-alpha / beta)
        print(f"T = {T} K: |ψ| = {psi_eq:.3f}")
    else:
        print(f"T = {T} K: |ψ| = 0 (通常状態)")

コード出力の解釈

このコードは、異なる温度での自由エネルギーランドスケープを示します。T > Tcではψ = 0が唯一の最小値、T < Tcではψ ≠ 0に2つの等価な最小値が現れます（対称性の破れ）。

1.4 第一・第二GL方程式

変分原理

平衡状態では、自由エネルギーが最小化されます。変分原理を適用すると、2つのGL方程式が得られます：

第一GL方程式（秩序パラメータ方程式）

αψ + β|ψ|²ψ + (1/2m*)(−iℏ∇ − q*A)²ψ = 0

これは、ψに関する変分δF/δψ* = 0から得られます。

第二GL方程式（電流密度方程式）

j = (q*/2m*i) [ψ*(−iℏ∇ − q*A)ψ − ψ(−iℏ∇ − q*A)*ψ*]
  = (q*ℏ/m*) |ψ|² [∇θ − (q*/ℏ)A]

これは超伝導電流密度を与え、ベクトルポテンシャルAに関する変分から得られます。

第二GL方程式の解釈

電流密度は2つの寄与を持ちます：

位相項：∇θは超流体速度に対応
ベクトルポテンシャル項：磁場との相互作用を記述

London方程式との関係

均一な秩序パラメータ（|ψ| = const）を仮定すると、第二GL方程式はLondon方程式に帰着します：

j = −(ns q*²/m*) A
∇ × j = −(ns q*²/m*) B

1.5 コヒーレンス長（ξ）

定義と物理的意味

GL コヒーレンス長ξ(T)は、秩序パラメータが空間的に変化できる特徴的な長さスケールを表します：

ξ²(T) = ℏ² / (2m*|α(T)|)
      = ℏ² / (2m*α0(Tc - T))

物理的解釈

ξは、秩序パラメータが通常状態の値（ψ = 0）から超伝導状態の平衡値に変化する距離を特徴づけます
ξは温度が上昇するにつれて増加し、Tcで発散します
典型的な値：ξ(0) ≈ 1-100 nm（材料による）

境界近傍の秩序パラメータプロファイル

超伝導体と真空の境界を考えます。一次元で磁場がない場合、第一GL方程式は：

−(ℏ²/2m*) d²ψ/dx² + αψ + β|ψ|²ψ = 0

境界条件：x → ∞でψ → ψ∞、x = 0でdψ/dx = 0（境界条件）

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

# Parameters
hbar = 1.0545718e-34  # Reduced Planck constant (J·s)
m_star = 2 * 9.10938e-31  # Cooper pair mass (2 * electron mass)
alpha_0 = 1e-10  # Arbitrary units
Tc = 10.0
T = 5.0
alpha = alpha_0 * (T - Tc)
beta = 1e-11

# Coherence length
xi = np.sqrt(hbar**2 / (2 * m_star * abs(alpha)))

# Equilibrium order parameter
psi_infty = np.sqrt(-alpha / beta)

# Dimensionless GL equation
# Let ψ = ψ∞ · f(x/ξ), where x/ξ = s
# Then: d²f/ds² = f - f³

def gl_1d(y, s):
    """
    y[0] = f (dimensionless order parameter)
    y[1] = df/ds
    """
    f, df = y
    d2f = f - f**3  # Dimensionless GL equation
    return [df, d2f]

# Solve from bulk (s = large) back to boundary (s = 0)
s_span = np.linspace(10, 0, 1000)
y0 = [1.0, 0.0]  # Initial conditions: f = 1, df/ds = 0 at s = large

solution = odeint(gl_1d, y0, s_span)
f_profile = solution[:, 0]

# Reverse to go from boundary to bulk
s_plot = s_span[::-1]
f_plot = f_profile[::-1]

# Convert to real coordinates
x_plot = s_plot * xi * 1e9  # Convert to nm
psi_plot = f_plot * psi_infty

# Plotting
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Plot 1: Dimensionless profile
axes[0].plot(s_plot, f_plot, 'b-', linewidth=2)
axes[0].axhline(y=1, color='r', linestyle='--', label='f = 1 (バルク値)')
axes[0].axhline(y=0, color='k', linestyle='--', alpha=0.3)
axes[0].set_xlabel('x/ξ (無次元距離)', fontsize=12)
axes[0].set_ylabel('f = ψ/ψ∞ (無次元秩序パラメータ)', fontsize=12)
axes[0].set_title('境界近傍の秩序パラメータプロファイル', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(0, 5)

# Shading to show healing length
axes[0].fill_between(s_plot, 0, f_plot, where=(s_plot <= 2),
                     alpha=0.3, color='cyan', label='ξ領域')

# Plot 2: Temperature dependence of ξ
T_range = np.linspace(0.1, Tc*0.99, 100)
alpha_T = alpha_0 * (T_range - Tc)
xi_T = np.sqrt(hbar**2 / (2 * m_star * abs(alpha_T))) * 1e9  # nm

axes[1].plot(T_range, xi_T, 'r-', linewidth=2)
axes[1].axvline(x=Tc, color='b', linestyle='--', label=f'Tc = {Tc} K')
axes[1].set_xlabel('温度 T (K)', fontsize=12)
axes[1].set_ylabel('コヒーレンス長 ξ(T) (nm)', fontsize=12)
axes[1].set_title('温度依存するコヒーレンス長', fontsize=14)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
axes[1].set_ylim(0, 100)

plt.tight_layout()
plt.show()

print(f"\nコヒーレンス長 ξ(T={T} K) = {xi*1e9:.2f} nm")
print(f"平衡秩序パラメータ ψ∞ = {psi_infty:.3e}")

1.6 侵入長（λ）

定義と物理的意味

GL 侵入長λ(T)は、磁場が超伝導体の表面から侵入できる特徴的な深さを表します：

λ²(T) = m* β / (μ0 q*² |ψ∞|²)
      = m* / (μ0 q*² ns)

ここで、ns = |ψ∞|² / βは超伝導電子対密度です。

物理的解釈

λは、表面から内部への磁場の指数関数的減衰長さです
λはLondon侵入長と同等です
典型的な値：λ(0) ≈ 20-500 nm（材料による）

Meissner効果の数学的記述

半無限超伝導体（x > 0）に平行磁場B(0) = B0を印加すると、内部の磁場は：

B(x) = B0 exp(−x/λ)

import numpy as np
import matplotlib.pyplot as plt

# Parameters
lambda_L = 50e-9  # London penetration depth (m)
B0 = 0.1  # Applied field (T)

# Distance from surface
x = np.linspace(0, 300e-9, 500)  # 0 to 300 nm

# Magnetic field profile
B = B0 * np.exp(-x / lambda_L)

# Screening current density (from Ampere's law)
# ∇ × B = μ0 j  →  dB/dx = −μ0 j_y
mu0 = 4 * np.pi * 1e-7  # Vacuum permeability
j = (B0 / (mu0 * lambda_L)) * np.exp(-x / lambda_L)

# Convert to more convenient units
x_nm = x * 1e9  # nm
lambda_nm = lambda_L * 1e9  # nm
j_MA = j * 1e-6  # MA/m²

# Plotting
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Plot 1: Magnetic field penetration
axes[0].plot(x_nm, B / B0, 'b-', linewidth=2)
axes[0].axhline(y=np.exp(-1), color='r', linestyle='--',
               label=f'B/B₀ = 1/e (x = λ = {lambda_nm:.0f} nm)')
axes[0].axvline(x=lambda_nm, color='r', linestyle='--', alpha=0.5)
axes[0].fill_between(x_nm, 0, B / B0, where=(x_nm <= lambda_nm),
                    alpha=0.3, color='cyan', label='侵入領域')
axes[0].set_xlabel('表面からの距離 x (nm)', fontsize=12)
axes[0].set_ylabel('規格化磁場 B/B₀', fontsize=12)
axes[0].set_title('超伝導体への磁場侵入', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(0, 300)

# Plot 2: Screening current
axes[1].plot(x_nm, j_MA / j_MA[0], 'g-', linewidth=2)
axes[1].axvline(x=lambda_nm, color='r', linestyle='--',
               label=f'x = λ = {lambda_nm:.0f} nm')
axes[1].fill_between(x_nm, 0, j_MA / j_MA[0], where=(x_nm <= lambda_nm),
                    alpha=0.3, color='yellow', label='遮蔽電流領域')
axes[1].set_xlabel('表面からの距離 x (nm)', fontsize=12)
axes[1].set_ylabel('規格化電流密度 j/j₀', fontsize=12)
axes[1].set_title('Meissner遮蔽電流', fontsize=14)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim(0, 300)

plt.tight_layout()
plt.show()

# Print penetration depth info
print(f"\n侵入長 λ = {lambda_nm:.1f} nm")
print(f"x = λにおいて、磁場は初期値の {np.exp(-1)*100:.1f}% に減衰")
print(f"x = 3λにおいて、磁場は初期値の {np.exp(-3)*100:.2f}% に減衰")

1.7 GLパラメータ κ = λ/ξ

定義と意義

Ginzburg-Landauパラメータκは、2つの特徴的長さの比として定義されます：

κ = λ / ξ

κは材料固有の無次元パラメータで、超伝導体の性質を決定する最も重要な量の一つです。

κの物理的意味

ξは秩序パラメータの変化のスケール（超伝導の「硬さ」）
λは磁場侵入のスケール（磁気応答）
κは、磁場侵入と秩序パラメータ抑制の競合を表します

κの値の範囲

実際の超伝導体におけるκの値：

純粋金属（Al, Pb, Sn）: κ ≈ 0.1 - 0.5
合金（NbTi, Nb3Sn）: κ ≈ 1 - 100
高温超伝導体（YBCO, BSCCO）: κ ≈ 50 - 100

1.8 第一種・第二種の分類（κ < 1/√2 vs κ > 1/√2）

表面エネルギーの概念

超伝導体と通常状態の境界には、エネルギーコストが関連付けられます。この表面エネルギーの符号が、超伝導体のタイプを決定します。

表面エネルギーの起源

正のエネルギーコスト：秩序パラメータがξの距離で0に減少する（凝縮エネルギーの損失）
負のエネルギー利得：磁場がλの距離で侵入する（磁気エネルギーの削減）

表面エネルギー σns ∝ (ξ/λ - λ/ξ) ∝ (1/κ - κ)

臨界値 κc = 1/√2 ≈ 0.707

条件	表面エネルギー	超伝導体タイプ	振る舞い
κ < 1/√2	σns > 0（正）	第一種	境界形成を避ける、急激な転移
κ = 1/√2	σns = 0	臨界点	境界エネルギーなし
κ > 1/√2	σns < 0（負）	第二種	境界形成を好む、混合状態

第一種超伝導体（κ < 1/√2）

表面エネルギーが正のため、超伝導領域と通常領域の境界を最小化
単一の臨界磁場Hcで急激に通常状態に転移
中間状態（混合はするが、巨視的なドメイン構造）

第二種超伝導体（κ > 1/√2）

表面エネルギーが負のため、境界形成が有利
2つの臨界磁場Hc1とHc2
Hc1 < H < Hc2で混合状態（渦状態、Abrikosov渦格子）が出現
各渦は量子化された磁束Φ0 = h/(2e)を運ぶ

import numpy as np
import matplotlib.pyplot as plt

# κ values
kappa_values = [0.3, 0.707, 2.0, 10.0]
labels = ['第一種 (κ=0.3)', '臨界 (κ=0.707)', '第二種 (κ=2)', '第二種 (κ=10)']
colors = ['blue', 'green', 'orange', 'red']

# Dimensionless distance from boundary
x = np.linspace(0, 5, 500)

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

for i, (kappa, label, color) in enumerate(zip(kappa_values, labels, colors)):
    ax = axes[i // 2, i % 2]

    # Order parameter profile (approximate)
    # ψ(x) ≈ ψ∞ tanh(x / (√2 ξ))
    psi_profile = np.tanh(x / np.sqrt(2))

    # Magnetic field profile (approximate)
    # B(x) ≈ B0 (1 - exp(-x/λ))
    # In units of ξ: x/λ = x/(κξ)
    B_profile = 1 - np.exp(-x / kappa)

    # Plot both profiles
    ax.plot(x, psi_profile, 'b-', linewidth=2, label='|ψ|/ψ∞ (秩序パラメータ)')
    ax.plot(x, B_profile, 'r-', linewidth=2, label='B/B₀ (磁場)')

    # Mark characteristic lengths
    ax.axvline(x=1, color='b', linestyle='--', alpha=0.5, label='ξ')
    ax.axvline(x=kappa, color='r', linestyle='--', alpha=0.5, label='λ = κξ')

    # Shading
    ax.fill_between(x, 0, psi_profile, alpha=0.2, color='blue')
    ax.fill_between(x, 0, B_profile, alpha=0.2, color='red')

    ax.set_xlabel('距離 (ξ単位)', fontsize=11)
    ax.set_ylabel('規格化された値', fontsize=11)
    ax.set_title(label, fontsize=12, color=color)
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)
    ax.set_xlim(0, 5)
    ax.set_ylim(0, 1.1)

plt.tight_layout()
plt.show()

# Surface energy comparison
print("\n表面エネルギーの符号:")
for kappa, label in zip(kappa_values, labels):
    sigma_sign = (1/kappa - kappa)
    if sigma_sign > 0:
        energy_type = "正（境界を避ける）"
    elif sigma_sign < 0:
        energy_type = "負（境界を好む）"
    else:
        energy_type = "ゼロ（臨界）"
    print(f"{label}: σns ∝ {sigma_sign:.3f} → {energy_type}")

1.9 表面エネルギーと境界解

超伝導-通常境界の詳細

超伝導体と通常金属の平面境界（x = 0）を考えます。x < 0が超伝導、x > 0が通常状態です。

境界での振る舞い

κ < 1/√2の場合：ψとBは異なる長さスケール（ξ < λ）で変化し、エネルギーの正味のコストが生じます
κ > 1/√2の場合：ψとBは異なる長さスケール（ξ > λ）で変化し、エネルギーの正味の利得が生じます

渦状態（Abrikosov格子）

第二種超伝導体では、Hc1 < H < Hc2の範囲で、磁束が量子化された渦として侵入します：

Φ0 = h / (2e) = 2.07 × 10⁻¹⁵ Wb（磁束量子）

各渦は：

中心でψ = 0（通常状態の芯、半径 ≈ ξ）
磁場が渦周りに集中（侵入深さ ≈ λ）
Φ0の磁束を運ぶ

import numpy as np
import matplotlib.pyplot as plt

# Simulate a single Abrikosov vortex
# Parameters
xi = 10  # Coherence length (nm)
lambda_L = 50  # Penetration depth (nm)
kappa = lambda_L / xi

# 2D grid
x = np.linspace(-100, 100, 200)
y = np.linspace(-100, 100, 200)
X, Y = np.meshgrid(x, y)
R = np.sqrt(X**2 + Y**2)

# Order parameter profile (approximate)
# |ψ(r)| ≈ ψ∞ tanh(r / (√2 ξ))
psi = np.tanh(R / (np.sqrt(2) * xi))

# Magnetic field profile (approximate)
# B(r) ≈ B0 K0(r/λ) (K0: modified Bessel function of second kind)
# Approximation: B(r) ≈ exp(-r/λ) for visualization
B = np.exp(-R / lambda_L)

# Plotting
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# Plot 1: Order parameter
im1 = axes[0].contourf(X, Y, psi, levels=20, cmap='Blues')
axes[0].contour(X, Y, psi, levels=10, colors='black', alpha=0.3, linewidths=0.5)
axes[0].set_xlabel('x (nm)', fontsize=11)
axes[0].set_ylabel('y (nm)', fontsize=11)
axes[0].set_title(f'秩序パラメータ |ψ|/ψ∞\n(渦芯 ～ ξ = {xi} nm)', fontsize=12)
axes[0].set_aspect('equal')
plt.colorbar(im1, ax=axes[0], label='|ψ|/ψ∞')

# Add circle at ξ
circle_xi = plt.Circle((0, 0), xi, fill=False, color='red',
                       linestyle='--', linewidth=2, label='r = ξ')
axes[0].add_patch(circle_xi)
axes[0].legend(fontsize=9)

# Plot 2: Magnetic field
im2 = axes[1].contourf(X, Y, B, levels=20, cmap='Reds')
axes[1].contour(X, Y, B, levels=10, colors='black', alpha=0.3, linewidths=0.5)
axes[1].set_xlabel('x (nm)', fontsize=11)
axes[1].set_ylabel('y (nm)', fontsize=11)
axes[1].set_title(f'磁場 B/B₀\n(侵入深さ ～ λ = {lambda_L} nm)', fontsize=12)
axes[1].set_aspect('equal')
plt.colorbar(im2, ax=axes[1], label='B/B₀')

# Add circle at λ
circle_lambda = plt.Circle((0, 0), lambda_L, fill=False, color='blue',
                          linestyle='--', linewidth=2, label='r = λ')
axes[1].add_patch(circle_lambda)
axes[1].legend(fontsize=9)

# Plot 3: Radial profiles
r_profile = np.linspace(0, 100, 500)
psi_profile = np.tanh(r_profile / (np.sqrt(2) * xi))
B_profile = np.exp(-r_profile / lambda_L)

axes[2].plot(r_profile, psi_profile, 'b-', linewidth=2, label='|ψ|/ψ∞')
axes[2].plot(r_profile, B_profile, 'r-', linewidth=2, label='B/B₀')
axes[2].axvline(x=xi, color='b', linestyle='--', alpha=0.5, label=f'ξ = {xi} nm')
axes[2].axvline(x=lambda_L, color='r', linestyle='--', alpha=0.5,
               label=f'λ = {lambda_L} nm')
axes[2].set_xlabel('半径 r (nm)', fontsize=11)
axes[2].set_ylabel('規格化された値', fontsize=11)
axes[2].set_title(f'径方向プロファイル\n(κ = λ/ξ = {kappa:.1f})', fontsize=12)
axes[2].legend(fontsize=9)
axes[2].grid(True, alpha=0.3)
axes[2].set_xlim(0, 100)

plt.tight_layout()
plt.show()

# Print vortex info
Phi0 = 2.067833848e-15  # Wb
print(f"\nAbrikosov渦の構造（κ = {kappa:.1f}）:")
print(f"渦芯サイズ（通常状態）: ～ξ = {xi} nm")
print(f"磁場侵入深さ: ～λ = {lambda_L} nm")
print(f"磁束量子: Φ₀ = {Phi0:.3e} Wb")

1.10 相図：第一種 vs 第二種

κに基づく完全な分類

import numpy as np
import matplotlib.pyplot as plt

# Range of kappa values
kappa = np.logspace(-1, 2, 500)  # 0.1 to 100

# Critical kappa
kappa_c = 1 / np.sqrt(2)

# Surface energy (normalized)
sigma = 1/kappa - kappa

# Thermodynamic critical field (normalized to Hc)
# For Type I: H = Hc
# For Type II: Hc1 < H < Hc2
Hc1_ratio = np.log(kappa) / (np.sqrt(2) * kappa)  # Hc1/Hc (approximate)
Hc2_ratio = np.sqrt(2) * kappa  # Hc2/Hc

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

# Plot 1: Surface energy
axes[0, 0].plot(kappa, sigma, 'b-', linewidth=2)
axes[0, 0].axhline(y=0, color='k', linestyle='-', linewidth=1)
axes[0, 0].axvline(x=kappa_c, color='r', linestyle='--', linewidth=2,
                  label=f'κc = 1/√2 ≈ {kappa_c:.3f}')
axes[0, 0].fill_between(kappa, sigma, 0, where=(sigma > 0), alpha=0.3,
                       color='blue', label='第一種（σ > 0）')
axes[0, 0].fill_between(kappa, sigma, 0, where=(sigma < 0), alpha=0.3,
                       color='red', label='第二種（σ < 0）')
axes[0, 0].set_xlabel('GLパラメータ κ', fontsize=11)
axes[0, 0].set_ylabel('表面エネルギー σns (任意単位)', fontsize=11)
axes[0, 0].set_title('表面エネルギーとκ', fontsize=12)
axes[0, 0].set_xscale('log')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3, which='both')
axes[0, 0].set_xlim(0.1, 100)
axes[0, 0].set_ylim(-5, 5)

# Plot 2: Critical fields
axes[0, 1].plot(kappa[kappa > kappa_c], Hc1_ratio[kappa > kappa_c],
               'g-', linewidth=2, label='Hc1/Hc')
axes[0, 1].plot(kappa[kappa > kappa_c], Hc2_ratio[kappa > kappa_c],
               'r-', linewidth=2, label='Hc2/Hc')
axes[0, 1].axhline(y=1, color='b', linestyle='--', linewidth=1, label='Hc')
axes[0, 1].axvline(x=kappa_c, color='orange', linestyle='--', linewidth=2,
                  label=f'κc = {kappa_c:.3f}')
axes[0, 1].set_xlabel('GLパラメータ κ', fontsize=11)
axes[0, 1].set_ylabel('規格化臨界磁場 H/Hc', fontsize=11)
axes[0, 1].set_title('臨界磁場とκ（第二種のみ）', fontsize=12)
axes[0, 1].set_xscale('log')
axes[0, 1].set_yscale('log')
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3, which='both')
axes[0, 1].set_xlim(0.7, 100)
axes[0, 1].set_ylim(0.01, 100)

# Plot 3: Phase diagram (H vs κ)
fig3 = axes[1, 0]
kappa_type1 = np.linspace(0.1, kappa_c, 100)
kappa_type2 = np.linspace(kappa_c, 10, 100)

# Type I
fig3.fill_between([0, kappa_c], [0, 0], [2, 2], alpha=0.3, color='blue',
                  label='第一種超伝導領域')
fig3.plot([kappa_c, kappa_c], [0, 2], 'k-', linewidth=2)

# Type II
Hc1_type2 = np.log(kappa_type2) / (np.sqrt(2) * kappa_type2)
Hc2_type2 = np.sqrt(2) * kappa_type2

fig3.fill_between(kappa_type2, 0, Hc1_type2, alpha=0.3, color='cyan',
                  label='Meissner状態')
fig3.fill_between(kappa_type2, Hc1_type2, np.minimum(Hc2_type2, 2),
                  alpha=0.3, color='yellow', label='混合状態（渦状態）')
fig3.fill_between(kappa_type2, np.minimum(Hc2_type2, 2), 2,
                  alpha=0.3, color='lightgray', label='通常状態')

fig3.plot(kappa_type2, Hc1_type2, 'g-', linewidth=2, label='Hc1')
fig3.plot(kappa_type2, Hc2_type2, 'r-', linewidth=2, label='Hc2')
fig3.axvline(x=kappa_c, color='orange', linestyle='--', linewidth=2)

fig3.set_xlabel('GLパラメータ κ', fontsize=11)
fig3.set_ylabel('磁場 H/Hc', fontsize=11)
fig3.set_title('超伝導相図（κ-H）', fontsize=12)
fig3.legend(fontsize=8, loc='upper left')
fig3.grid(True, alpha=0.3)
fig3.set_xlim(0, 10)
fig3.set_ylim(0, 2)

# Plot 4: Material examples
fig4 = axes[1, 1]
materials = {
    'Al': 0.15, 'Pb': 0.4, 'Sn': 0.2, 'Nb': 0.8,
    'NbTi': 1.5, 'Nb3Sn': 2.5, 'YBCO': 50, 'Bi2212': 80
}
type1 = {k: v for k, v in materials.items() if v < kappa_c}
type2 = {k: v for k, v in materials.items() if v >= kappa_c}

y_pos_type1 = np.arange(len(type1))
y_pos_type2 = np.arange(len(type2)) + len(type1) + 0.5

fig4.barh(y_pos_type1, list(type1.values()), color='blue', alpha=0.7,
         label='第一種')
fig4.barh(y_pos_type2, list(type2.values()), color='red', alpha=0.7,
         label='第二種')

fig4.axvline(x=kappa_c, color='green', linestyle='--', linewidth=2,
            label=f'κc = {kappa_c:.3f}')

all_materials = list(type1.keys()) + list(type2.keys())
y_pos_all = list(y_pos_type1) + list(y_pos_type2)
fig4.set_yticks(y_pos_all)
fig4.set_yticklabels(all_materials)
fig4.set_xlabel('GLパラメータ κ', fontsize=11)
fig4.set_title('実際の超伝導体のκ値', fontsize=12)
fig4.set_xscale('log')
fig4.legend(fontsize=9)
fig4.grid(True, alpha=0.3, axis='x')
fig4.set_xlim(0.1, 100)

plt.tight_layout()
plt.show()

# Print classification
print("\n超伝導体の分類:")
print(f"臨界値: κc = 1/√2 = {kappa_c:.4f}\n")
print("第一種超伝導体 (κ < κc):")
for mat, k in sorted(type1.items(), key=lambda x: x[1]):
    print(f"  {mat}: κ = {k:.2f}")
print("\n第二種超伝導体 (κ > κc):")
for mat, k in sorted(type2.items(), key=lambda x: x[1]):
    print(f"  {mat}: κ = {k:.1f}")

まとめ

重要なポイント

GL理論は複素秩序パラメータψを用いた超伝導の現象論的記述である
GL自由エネルギー汎関数は秩序パラメータ、勾配項、磁場項を含む
第一GL方程式は秩序パラメータの空間分布を、第二GL方程式は超伝導電流を記述する
コヒーレンス長ξは秩序パラメータ変化のスケール、侵入長λは磁場侵入のスケールを表す
GLパラメータκ = λ/ξは材料固有の無次元量で、超伝導体のタイプを決定する
κ < 1/√2: 第一種（表面エネルギー正、急激な転移）
κ > 1/√2: 第二種（表面エネルギー負、混合状態、渦格子）
第二種超伝導体はHc1とHc2の間でAbrikosov渦状態を形成し、高磁場応用に不可欠

練習問題

問題1

GL自由エネルギー汎関数から、磁場がない場合の均一超伝導体における平衡秩序パラメータ|ψ∞|を求めてください。また、T < Tcにおける凝縮エネルギー密度（通常状態と超伝導状態の自由エネルギー差）を計算してください。

問題2

ある超伝導体がT = 0でξ(0) = 5 nm、λ(0) = 100 nmを持つとします。GLパラメータκを計算し、この材料が第一種か第二種かを判定してください。また、この材料がHc = 0.05 Tを持つ場合、Hc1とHc2（近似式を使用）を推定してください。

問題3

第二種超伝導体における単一のAbrikosov渦を考えます。渦の中心（r = 0）と、r >> λでの秩序パラメータ|ψ|および磁場Bの振る舞いを説明してください。また、なぜ各渦が磁束量子Φ0 = h/(2e)を運ぶのかを説明してください。

問題4（発展）

表面エネルギーσns ∝ (1/κ - κ)の導出を次元解析に基づいて説明してください。ヒント：エネルギー密度のスケールはα²/β、長さスケールはξとλです。