第5章：フォノン分光法

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
M = 1.0  # 原子質量（任意単位）
K = 1.0  # ばね定数（任意単位）
a = 1.0  # 格子定数（任意単位）

# 波数ベクトル
q = np.linspace(-np.pi/a, np.pi/a, 200)

# 分散関係
omega = 2 * np.sqrt(K/M) * np.abs(np.sin(q*a/2))

# プロット
plt.figure(figsize=(10, 6))
plt.plot(q*a/np.pi, omega, 'b-', linewidth=2)
plt.xlabel('q (π/a)', fontsize=14)
plt.ylabel('ω (√K/M)', fontsize=14)
plt.title('1次元単原子鎖のフォノン分散', fontsize=16)
plt.grid(True, alpha=0.3)
plt.axvline(x=0, color='k', linestyle='--', alpha=0.3)
plt.tight_layout()
plt.show()

print(f"ゾーン中心での群速度: {2*np.sqrt(K/M)*a/2:.3f} √(K/M)·a")
print(f"ゾーン境界での周波数: {2*np.sqrt(K/M):.3f} √(K/M)")

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
M1 = 1.0   # 原子1の質量
M2 = 2.0   # 原子2の質量
K = 1.0    # ばね定数
a = 2.0    # 格子定数（2原子含む）

# 波数ベクトル
q = np.linspace(-np.pi/a, np.pi/a, 200)

# 分散関係の計算
def phonon_dispersion_diatomic(q, M1, M2, K, a):
    """2原子鎖のフォノン分散を計算"""
    M_sum = M1 + M2
    M_prod = M1 * M2

    # ω²の2つの解（光学・音響ブランチ）
    term1 = K * M_sum / M_prod
    term2 = K * np.sqrt(M_sum**2 - 4*M_prod*np.sin(q*a/2)**2) / M_prod

    omega_optical = np.sqrt(term1 + term2)
    omega_acoustic = np.sqrt(term1 - term2)

    return omega_acoustic, omega_optical

omega_ac, omega_op = phonon_dispersion_diatomic(q, M1, M2, K, a)

# プロット
plt.figure(figsize=(10, 6))
plt.plot(q*a/np.pi, omega_ac, 'b-', linewidth=2, label='音響ブランチ')
plt.plot(q*a/np.pi, omega_op, 'r-', linewidth=2, label='光学ブランチ')
plt.xlabel('q (π/a)', fontsize=14)
plt.ylabel('ω (√K/M)', fontsize=14)
plt.title('1次元2原子鎖のフォノン分散', fontsize=16)
plt.legend(fontsize=12)
plt.grid(True, alpha=0.3)
plt.axvline(x=0, color='k', linestyle='--', alpha=0.3)
plt.tight_layout()
plt.show()

# ゾーン中心の周波数
omega_op_center = np.sqrt(2*K*(M1+M2)/(M1*M2))
omega_ac_center = 0

print(f"光学フォノン（q=0）: {omega_op_center:.3f} √(K/M)")
print(f"音響フォノン（q=0）: {omega_ac_center:.3f}")
print(f"ゾーン境界ギャップ: {omega_op[-1] - omega_ac[-1]:.3f}")

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

def calculate_vdos(q, omega, omega_bins):
    """
    フォノン分散から振動状態密度を計算

    Parameters:
    -----------
    q : array
        波数ベクトル
    omega : array or tuple
        角周波数（単一配列または複数ブランチのタプル）
    omega_bins : array
        DOSを計算する周波数ビン

    Returns:
    --------
    dos : array
        振動状態密度
    """
    dos = np.zeros_like(omega_bins)

    # 複数ブランチの場合
    if isinstance(omega, tuple):
        omega_list = omega
    else:
        omega_list = [omega]

    for omega_branch in omega_list:
        # ヒストグラムでDOSを計算
        hist, _ = np.histogram(omega_branch, bins=omega_bins)
        dos += hist

    # 群速度による補正（より正確なDOS）
    # DOS(ω) ∝ 1/|dω/dq|

    return dos / (omega_bins[1] - omega_bins[0])

# 2原子鎖の例
q = np.linspace(-np.pi/a, np.pi/a, 1000)
omega_ac, omega_op = phonon_dispersion_diatomic(q, M1, M2, K, a)

# DOSの計算
omega_bins = np.linspace(0, 2.0, 200)
dos = calculate_vdos(q, (omega_ac, omega_op), omega_bins)

# プロット
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 分散関係
ax1.plot(q*a/np.pi, omega_ac, 'b-', linewidth=2, label='音響')
ax1.plot(q*a/np.pi, omega_op, 'r-', linewidth=2, label='光学')
ax1.set_xlabel('q (π/a)', fontsize=12)
ax1.set_ylabel('ω (√K/M)', fontsize=12)
ax1.set_title('フォノン分散', fontsize=14)
ax1.legend()
ax1.grid(True, alpha=0.3)

# VDOS
ax2.plot(omega_bins, dos, 'k-', linewidth=2)
ax2.fill_between(omega_bins, dos, alpha=0.3)
ax2.set_xlabel('ω (√K/M)', fontsize=12)
ax2.set_ylabel('状態密度（任意単位）', fontsize=12)
ax2.set_title('振動状態密度', fontsize=14)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def lorentzian(omega, omega0, gamma, intensity=1.0):
    """
    ローレンツ関数（フォノンピークの線形）

    Parameters:
    -----------
    omega : array
        周波数
    omega0 : float
        ピーク中心周波数
    gamma : float
        半値全幅（FWHM）
    intensity : float
        ピーク強度
    """
    return intensity * (gamma/2)**2 / ((omega - omega0)**2 + (gamma/2)**2)

def bose_einstein(omega, T):
    """
    Bose-Einstein分布関数

    Parameters:
    -----------
    omega : array
        周波数
    T : float
        温度（K）
    """
    hbar = 1.054571817e-34  # J·s
    kB = 1.380649e-23       # J/K

    # cm⁻¹からJへの変換
    omega_J = omega * 100 * 3e8 * hbar * 2 * np.pi

    if T == 0:
        return np.zeros_like(omega)
    else:
        x = omega_J / (kB * T)
        return 1.0 / (np.exp(x) - 1)

def raman_spectrum(omega, peaks, T=300):
    """
    ラマンスペクトルをシミュレート

    Parameters:
    -----------
    omega : array
        ラマンシフト（cm⁻¹）
    peaks : list of tuples
        (周波数, 線幅, 強度)のリスト
    T : float
        温度（K）
    """
    spectrum_stokes = np.zeros_like(omega)
    spectrum_antistokes = np.zeros_like(omega)

    for omega0, gamma, intensity in peaks:
        # Stokes散乱
        n = bose_einstein(np.array([omega0]), T)[0]
        stokes_intensity = intensity * (n + 1)
        spectrum_stokes += lorentzian(omega, omega0, gamma, stokes_intensity)

        # Anti-Stokes散乱
        antistokes_intensity = intensity * n
        spectrum_antistokes += lorentzian(-omega, omega0, gamma, antistokes_intensity)

    return spectrum_stokes + spectrum_antistokes

# シリコンのラマンスペクトルをシミュレート
omega = np.linspace(-600, 600, 2000)

# Si: ゾーン中心の3重縮退F₂gモード
peaks_si = [
    (520, 3, 1.0),  # (周波数 cm⁻¹, FWHM cm⁻¹, 相対強度)
]

# 異なる温度でのスペクトル
temperatures = [77, 300, 500]
colors = ['blue', 'green', 'red']

plt.figure(figsize=(12, 6))

for T, color in zip(temperatures, colors):
    spectrum = raman_spectrum(omega, peaks_si, T)
    plt.plot(omega, spectrum, color=color, linewidth=2, label=f'T = {T} K')

plt.xlabel('ラマンシフト (cm⁻¹)', fontsize=14)
plt.ylabel('強度（任意単位）', fontsize=14)
plt.title('シリコンのラマンスペクトル（温度依存性）', fontsize=16)
plt.legend(fontsize=12)
plt.grid(True, alpha=0.3)
plt.axvline(x=0, color='k', linestyle='--', alpha=0.5, label='レイリー散乱')
plt.xlim(-600, 600)
plt.tight_layout()
plt.show()

# 温度によるピークシフトと線幅の効果を解析
print("\n温度依存性の分析：")
for T in temperatures:
    n = bose_einstein(np.array([520]), T)[0]
    ratio = n / (n + 1) if n > 0 else 0
    print(f"T = {T} K: n = {n:.3f}, Anti-Stokes/Stokes = {ratio:.3f}")

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def fit_raman_peak(omega, intensity):
    """
    ラマンピークをローレンツ関数でフィッティング

    Returns:
    --------
    fit_params : dict
        フィッティングパラメータ（中心周波数、線幅、強度）
    """
    # 初期推定値
    omega0_guess = omega[np.argmax(intensity)]
    gamma_guess = 5.0
    I_guess = np.max(intensity)
    baseline_guess = np.min(intensity)

    # フィッティング関数（ベースライン含む）
    def fit_func(x, omega0, gamma, I, baseline):
        return lorentzian(x, omega0, gamma, I) + baseline

    # フィッティング実行
    popt, pcov = curve_fit(
        fit_func,
        omega,
        intensity,
        p0=[omega0_guess, gamma_guess, I_guess, baseline_guess]
    )

    # 標準偏差
    perr = np.sqrt(np.diag(pcov))

    fit_params = {
        'omega0': popt[0],
        'gamma': popt[1],
        'intensity': popt[2],
        'baseline': popt[3],
        'omega0_err': perr[0],
        'gamma_err': perr[1]
    }

    return fit_params, fit_func(omega, *popt)

# 疑似実験データの生成（ノイズ付き）
omega_exp = np.linspace(500, 540, 100)
true_params = (520, 3, 100, 5)  # 真のパラメータ
signal = lorentzian(omega_exp, *true_params[:3]) + true_params[3]
noise = np.random.normal(0, 2, len(omega_exp))
intensity_exp = signal + noise

# フィッティング
fit_params, fitted_curve = fit_raman_peak(omega_exp, intensity_exp)

# 結果のプロット
plt.figure(figsize=(10, 6))
plt.plot(omega_exp, intensity_exp, 'o', label='実験データ', alpha=0.6)
plt.plot(omega_exp, fitted_curve, 'r-', linewidth=2, label='フィッティング')
plt.xlabel('ラマンシフト (cm⁻¹)', fontsize=14)
plt.ylabel('強度（任意単位）', fontsize=14)
plt.title('ラマンピークのフィッティング', fontsize=16)
plt.legend(fontsize=12)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# 結果の表示
print("フィッティング結果：")
print(f"ピーク位置: {fit_params['omega0']:.2f} ± {fit_params['omega0_err']:.2f} cm⁻¹")
print(f"線幅（FWHM）: {fit_params['gamma']:.2f} ± {fit_params['gamma_err']:.2f} cm⁻¹")
print(f"強度: {fit_params['intensity']:.2f}")
print(f"ベースライン: {fit_params['baseline']:.2f}")