第1章: 強結合理論とEliashberg形式

import numpy as np
from scipy.integrate import quad
from scipy.interpolate import interp1d

class EliashbergSolver:
    """
    Eliashberg方程式ソルバー（虚軸形式）

    Parameters:
    -----------
    alpha2F : callable
        電子-フォノン結合スペクトル関数 α²F(ω)
    mu_star : float
        Coulombの擬ポテンシャル
    omega_c : float
        Coulomb相互作用のカットオフ周波数
    n_matsubara : int
        Matsubara周波数の数
    """

    def __init__(self, alpha2F, mu_star=0.13, omega_c=500, n_matsubara=1024):
        self.alpha2F = alpha2F
        self.mu_star = mu_star
        self.omega_c = omega_c
        self.n_matsubara = n_matsubara

        # 電子-フォノン結合定数を計算
        self.lambda_ep = self.compute_lambda()

    def compute_lambda(self):
        """電子-フォノン結合定数λの計算"""
        integrand = lambda omega: 2 * self.alpha2F(omega) / omega if omega > 0 else 0
        lambda_val, _ = quad(integrand, 0, 1000, limit=100)
        return lambda_val

    def kernel_lambda(self, nu, T):
        """
        電子-フォノンカーネル λ(ν)

        Parameters:
        -----------
        nu : float
            Matsubara周波数差
        T : float
            温度 [K]
        """
        def integrand(omega):
            if omega == 0:
                return 0
            return 2 * omega * self.alpha2F(omega) / (nu**2 + omega**2)

        result, _ = quad(integrand, 0, 1000, limit=100)
        return result

    def solve_at_temperature(self, T, max_iter=100, tol=1e-6):
        """
        指定温度でEliashberg方程式を解く

        Parameters:
        -----------
        T : float
            温度 [K]
        max_iter : int
            最大反復回数
        tol : float
            収束判定閾値

        Returns:
        --------
        Z : ndarray
            質量繰り込み関数
        phi : ndarray
            ギャップ関数
        omega_n : ndarray
            Matsubara周波数
        converged : bool
            収束したかどうか
        """
        # Matsubara周波数の設定
        pi_T = np.pi * T * 8.617e-5  # eV単位（kB = 8.617e-5 eV/K）
        n = np.arange(self.n_matsubara)
        omega_n = (2*n + 1) * pi_T

        # 初期値
        Z = np.ones(self.n_matsubara)
        phi = np.full(self.n_matsubara, 1.76 * T * 8.617e-5)  # BCS初期値

        # カーネル行列を事前計算（高速化）
        lambda_kernel = np.zeros((self.n_matsubara, self.n_matsubara))
        for i in range(self.n_matsubara):
            for j in range(self.n_matsubara):
                lambda_kernel[i, j] = self.kernel_lambda(
                    omega_n[i] - omega_n[j], T
                )

        # 反復解法
        for iteration in range(max_iter):
            Z_old = Z.copy()
            phi_old = phi.copy()

            # Zの更新
            for i in range(self.n_matsubara):
                sum_Z = 0
                for j in range(self.n_matsubara):
                    denom = np.sqrt(omega_n[j]**2 * Z[j]**2 + phi[j]**2)
                    sum_Z += lambda_kernel[i, j] * omega_n[j] * Z[j] / denom * np.sign(omega_n[j])

                Z[i] = 1 + pi_T / omega_n[i] * sum_Z

            # φの更新
            for i in range(self.n_matsubara):
                sum_phi = 0
                for j in range(self.n_matsubara):
                    # Coulomb項
                    mu_term = self.mu_star if abs(omega_n[j]) < self.omega_c else 0

                    denom = np.sqrt(omega_n[j]**2 * Z[j]**2 + phi[j]**2)
                    sum_phi += (lambda_kernel[i, j] - mu_term) * phi[j] / denom

                phi[i] = pi_T / Z[i] * sum_phi

            # 収束判定
            delta_Z = np.max(np.abs(Z - Z_old))
            delta_phi = np.max(np.abs(phi - phi_old))

            if delta_Z < tol and delta_phi < tol:
                return Z, phi, omega_n, True

        # 収束しなかった
        return Z, phi, omega_n, False

    def find_Tc(self, T_min=1.0, T_max=50.0, n_points=50):
        """
        転移温度T_cを二分法で探索

        Parameters:
        -----------
        T_min : float
            探索下限温度 [K]
        T_max : float
            探索上限温度 [K]
        n_points : int
            探索点数

        Returns:
        --------
        Tc : float
            転移温度 [K]
        """
        # 粗い探索
        T_array = np.linspace(T_min, T_max, n_points)
        phi0_array = []

        for T in T_array:
            _, phi, _, converged = self.solve_at_temperature(T)
            if converged:
                phi0_array.append(phi[0])
            else:
                phi0_array.append(np.nan)

        phi0_array = np.array(phi0_array)

        # φ(ω_0) ≈ 0 となる温度を探す
        valid = ~np.isnan(phi0_array)
        if not np.any(valid):
            raise ValueError("収束した解が見つかりませんでした")

        # 二分法
        T_low = T_min
        T_high = T_max

        while T_high - T_low > 0.01:  # 0.01 K精度
            T_mid = (T_low + T_high) / 2
            _, phi, _, converged = self.solve_at_temperature(T_mid)

            if not converged:
                T_high = T_mid
                continue

            if phi[0] > 1e-6:  # まだ有限のギャップ
                T_low = T_mid
            else:
                T_high = T_mid

        return (T_low + T_high) / 2


# 使用例: Pb型のα²F(ω)
def alpha2F_Pb(omega):
    """
    Pbのα²F(ω)のモデル（簡略化）
    3つのGaussianピークで近似
    """
    # パラメータ（単位: meV）
    peaks = [
        {'center': 4.0, 'width': 1.5, 'height': 0.30},  # TA
        {'center': 7.0, 'width': 1.5, 'height': 0.25},  # LA
        {'center': 9.0, 'width': 1.5, 'height': 0.20},  # LO
    ]

    result = 0
    for peak in peaks:
        result += peak['height'] * np.exp(
            -((omega - peak['center']) / peak['width'])**2
        )

    return result


# ソルバーの初期化と実行
print("=== Eliashberg方程式ソルバー ===\n")

solver = EliashbergSolver(
    alpha2F=alpha2F_Pb,
    mu_star=0.13,
    omega_c=500,  # meV
    n_matsubara=512
)

print(f"電子-フォノン結合定数 λ = {solver.lambda_ep:.3f}")

# T_cの探索
print("\nT_cを探索中...")
Tc = solver.find_Tc(T_min=1.0, T_max=15.0, n_points=30)
print(f"転移温度 T_c = {Tc:.2f} K")

# 低温での解を計算
print(f"\nT = {Tc/2:.2f} K でのギャップ構造を計算中...")
Z, phi, omega_n, converged = solver.solve_at_temperature(Tc/2)

if converged:
    print(f"収束成功")
    print(f"ゼロ周波数ギャップ φ(0) = {phi[0]*1000:.2f} meV")
    print(f"質量繰り込み Z(0) = {Z[0]:.3f}")
    print(f"ギャップ比 2Δ₀/k_B T_c = {2*phi[0]/(8.617e-5*Tc):.2f}")
else:
    print("収束失敗")


# 結果の可視化
import matplotlib.pyplot as plt

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

# (a) α²F(ω)
omega_plot = np.linspace(0, 15, 300)
alpha2F_plot = [alpha2F_Pb(w) for w in omega_plot]

axes[0, 0].plot(omega_plot, alpha2F_plot, 'b-', linewidth=2)
axes[0, 0].set_xlabel('ω (meV)', fontsize=12)
axes[0, 0].set_ylabel('α²F(ω)', fontsize=12)
axes[0, 0].set_title('(a) 電子-フォノン結合スペクトル', fontsize=13)
axes[0, 0].grid(True, alpha=0.3)

# (b) Z(ωn)
axes[0, 1].plot(omega_n[:50]*1000, Z[:50], 'ro-', markersize=4)
axes[0, 1].axhline(y=1, color='k', linestyle='--', alpha=0.5)
axes[0, 1].set_xlabel('ω_n (meV)', fontsize=12)
axes[0, 1].set_ylabel('Z(iω_n)', fontsize=12)
axes[0, 1].set_title(f'(b) 質量繰り込み関数 (T = {Tc/2:.1f} K)', fontsize=13)
axes[0, 1].grid(True, alpha=0.3)

# (c) φ(ωn)
axes[1, 0].plot(omega_n[:50]*1000, phi[:50]*1000, 'go-', markersize=4)
axes[1, 0].set_xlabel('ω_n (meV)', fontsize=12)
axes[1, 0].set_ylabel('φ(iω_n) (meV)', fontsize=12)
axes[1, 0].set_title(f'(c) ギャップ関数 (T = {Tc/2:.1f} K)', fontsize=13)
axes[1, 0].grid(True, alpha=0.3)

# (d) 温度依存性
T_range = np.linspace(0.3*Tc, 0.95*Tc, 15)
phi0_T = []

for T in T_range:
    _, phi_T, _, conv = solver.solve_at_temperature(T, max_iter=50)
    if conv:
        phi0_T.append(phi_T[0]*1000)
    else:
        phi0_T.append(np.nan)

# BCS理論との比較
T_BCS = np.linspace(0, Tc, 100)
Delta_BCS = 1.76 * 8.617e-5 * Tc * 1000 * np.sqrt(1 - (T_BCS/Tc)**2)

axes[1, 1].plot(T_range, phi0_T, 'ro-', markersize=6, label='Eliashberg')
axes[1, 1].plot(T_BCS, Delta_BCS, 'b--', linewidth=2, label='BCS')
axes[1, 1].set_xlabel('T (K)', fontsize=12)
axes[1, 1].set_ylabel('Δ(T) (meV)', fontsize=12)
axes[1, 1].set_title('(d) ギャップの温度依存性', fontsize=13)
axes[1, 1].legend(fontsize=11)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('eliashberg_solution.png', dpi=150, bbox_inches='tight')
print("\n図を 'eliashberg_solution.png' として保存しました")