第2章: フォノン相互作用

1. 調和近似を超えて

1.1 調和近似の限界

第1章で学んだ調和近似では、格子ポテンシャルエネルギーを変位の二次までで展開しました。 この近似下では、フォノンモードは互いに独立であり、エネルギーは完全に保存され、 フォノンの寿命は無限大となります。しかし、この描像には重大な物理的問題があります：

• 熱膨張の欠如: 調和近似では熱膨張係数がゼロになる
• 無限の熱伝導率: フォノン散乱がないため熱抵抗が生じない
• 温度独立性: フォノン周波数が温度に依存しない
• 無限の寿命: 実験で観測される分光線幅を説明できない

1.2 非調和ポテンシャル展開

結晶ポテンシャルエネルギーをより完全に展開すると：

ここで、$u_{i\alpha}$は原子$i$の方向$\alpha$における変位、$\Phi^{(n)}$は$n$次力定数です。

1.3 三次および四次力定数

力定数は、平衡位置近傍でのポテンシャルの導関数として定義されます：

graph LR
    A[調和近似] -->|三次項| B[三フォノン過程]
    A -->|四次項| C[四フォノン過程]
    B --> D[フォノン寿命]
    B --> E[熱膨張]
    C --> F[追加の散乱]
    D --> G[有限熱伝導率]
    E --> G
    F --> G
                

2. フォノン-フォノン散乱理論

2.1 フェルミの黄金則

非調和ハミルトニアンを摂動として扱い、時間依存摂動論を適用すると、 初期状態$|i\rangle$から終状態$|f\rangle$への遷移率は、フェルミの黄金則で与えられます：

ここで、$H_{\text{anh}}$は非調和ハミルトニアン、$\delta(E_f - E_i)$はエネルギー保存を表すデルタ関数です。

2.2 三フォノン過程

三次非調和項$\Phi^{(3)}$による主要な散乱過程は、三つのフォノンモードが関与します：

崩壊過程（Decay Process）

波数$\mathbf{q}_0$のフォノンが、二つのフォノン（$\mathbf{q}_1, \mathbf{q}_2$）に分裂します。

合体過程（Coalescence Process）

二つのフォノンが合体して一つのフォノンになります。

2.3 散乱率の導出

フォノンモード$(\mathbf{q}, j)$の三フォノン散乱率は、第二量子化された非調和ハミルトニアンから導かれます：

ここで、

• $n_{\mathbf{q}j} = [\exp(\hbar\omega_{\mathbf{q}j}/k_BT) - 1]^{-1}$: ボース・アインシュタイン分布
• $V_3^{\pm}$: 三次相互作用行列要素
• $\mathbf{G}$: 逆格子ベクトル
• $N$: 単位胞の数

3. 通常過程とウムクラップ過程

3.1 保存則

フォノン散乱過程は、次の二つの保存則に従います：

結晶運動量は逆格子ベクトル$\mathbf{G}$の不定性を持つことに注意が必要です。

3.2 通常過程（Normal Process, N-process）

$\mathbf{G} = 0$の場合、つまり：

この過程では、全フォノン運動量が保存されます。N過程は：

• フォノンのエネルギー分布を再配分する
• 局所平衡状態（ドリフトした分布）を確立する
• 熱抵抗を直接生じない（運動量が保存されるため）

3.3 ウムクラップ過程（Umklapp Process, U-process）

$\mathbf{G} \neq 0$の場合、つまり散乱後の波数ベクトルの和が第一ブリルアン域外に出る場合：

ウムクラップ過程は：

• 全運動量を保存しない（結晶格子が逆向きの運動量を受け取る）
• 熱抵抗の主要な原因となる
• 高温で重要になる（大きな波数ベクトルを持つフォノンが必要）
• 温度とともに指数関数的に増加する（$\propto e^{-\Theta_U/T}$）

graph TD
    A[フォノン散乱] --> B{G = 0?}
    B -->|Yes| C[通常過程 N]
    B -->|No| D[ウムクラップ過程 U]
    C --> E[運動量保存]
    D --> F[運動量非保存]
    E --> G[熱抵抗なし]
    F --> H[熱抵抗あり]
    C --> I[低温で支配的]
    D --> J[高温で重要]
                

3.4 温度依存性

低温（$T \ll \Theta_D$、デバイ温度）では：

• 励起されるフォノンは主に小さな$\mathbf{q}$に集中
• $\mathbf{q} + \mathbf{q}' \lesssim 2q_{\max} < |\mathbf{G}|_{\min}$となりU過程は稀
• N過程が支配的で、熱伝導率は境界散乱で制限される

高温（$T \gtrsim \Theta_D$）では：

• 大きな$\mathbf{q}$のフォノンも励起される
• U過程が頻繁に起こり、熱抵抗を支配する
• 熱伝導率は$\kappa \propto 1/T$の温度依存性を示す

4. フォノン寿命と線幅

4.1 不確定性原理との関連

フォノンの寿命$\tau$が有限であることは、エネルギー準位に不確定性を導入します。 ハイゼンベルクの不確定性原理より：

この不確定性は、分光スペクトルにおいて有限の線幅として現れます。

4.2 分光線幅

フォノンモード$(\mathbf{q}, j)$のスペクトル線形は、ローレンツ型となります：

ここで、半値全幅（FWHM）$\Gamma$は寿命と以下の関係があります：

4.3 寿命の温度依存性

フォノン寿命は温度とともに変化します。三フォノン散乱による寿命は、 ボース分布関数$n(\omega, T)$を通じて温度依存性を持ちます：

高温極限（$k_BT \gg \hbar\omega$）では、$n(\omega, T) \approx k_BT/\hbar\omega$となり：

つまり、温度上昇とともに線幅が線形に広がります。

5. 熱抵抗との関連

5.1 マティーセンの法則

異なる散乱機構は独立であると仮定すると、全散乱率は各機構の散乱率の和となります（マティーセンの法則）：

ここで、

• $\tau_{\text{ph-ph}}^{-1}$: フォノン-フォノン散乱（N過程とU過程）
• $\tau_{\text{boundary}}^{-1}$: 境界散乱（$\propto L^{-1}$、$L$は試料サイズ）
• $\tau_{\text{defect}}^{-1}$: 点欠陥、不純物散乱
• $\tau_{\text{electron}}^{-1}$: 電子-フォノン散乱

5.2 熱伝導率への寄与

ボルツマン輸送方程式（詳細は第3章）から、熱伝導率は：

ここで、$C_{\mathbf{q}j}$は比熱、$v_{\mathbf{q}j}$は群速度です。

5.3 温度依存性のパターン

6. Pythonによるフォノン散乱率計算

6.1 簡略化されたモデル

実際の散乱率計算は非常に複雑ですが、ここでは教育的な目的で簡略化されたモデルを実装します。 デバイモデルを仮定し、単一の音響ブランチを考えます。

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

# 物理定数
hbar = 1.054571817e-34  # J·s
kB = 1.380649e-23      # J/K

class PhononScatteringModel:
    """
    簡略化されたフォノン散乱率計算モデル
    デバイモデルを仮定し、三フォノン過程を考慮
    """

    def __init__(self, theta_D, mass, a, gamma):
        """
        パラメータ:
            theta_D: デバイ温度 [K]
            mass: 平均原子質量 [kg]
            a: 格子定数 [m]
            gamma: グリュナイゼンパラメータ (非調和性の尺度)
        """
        self.theta_D = theta_D
        self.mass = mass
        self.a = a
        self.gamma = gamma

        # デバイ周波数
        self.omega_D = kB * theta_D / hbar

        # 平均音速（簡略化）
        self.v_s = self.omega_D * a / (2 * np.pi)

    def bose_einstein(self, omega, T):
        """ボース・アインシュタイン分布"""
        x = hbar * omega / (kB * T)
        if x < 1e-3:  # 高温極限
            return kB * T / (hbar * omega)
        else:
            return 1.0 / (np.exp(x) - 1.0)

    def umklapp_rate(self, omega, T):
        """
        ウムクラップ散乱率（簡略化されたモデル）

        Klemens の近似式:
        τ_U^(-1) ∝ ω² T exp(-θ_D / 3T)
        """
        # ウムクラップ過程の活性化温度
        theta_U = self.theta_D / 3.0

        # 散乱率の前因子（現象論的）
        A_U = (self.gamma**2 * self.omega_D) / (self.mass * self.v_s**2 * self.a**3)

        # 温度と周波数依存性
        if T < 0.1:  # 数値安定性のため
            return 0.0
        else:
            rate = A_U * (omega / self.omega_D)**2 * T * np.exp(-theta_U / T)
            return rate

    def normal_rate(self, omega, T):
        """
        通常散乱率（簡略化されたモデル）

        N過程は熱抵抗に直接寄与しないが、
        平衡化に重要な役割を果たす
        """
        A_N = (self.gamma**2 * self.omega_D) / (self.mass * self.v_s**2 * self.a**3)

        # N過程はU過程よりも活性化障壁が低い
        rate = A_N * (omega / self.omega_D)**2 * T
        return rate

    def boundary_rate(self, L):
        """
        境界散乱率
        Casimir限界: τ_B^(-1) = v_s / L
        """
        return self.v_s / L

    def total_rate(self, omega, T, L=np.inf):
        """
        全散乱率（マティーセンの法則）
        """
        rate_U = self.umklapp_rate(omega, T)
        rate_N = self.normal_rate(omega, T)
        rate_B = self.boundary_rate(L) if L < np.inf else 0.0

        return rate_U + rate_N + rate_B

    def thermal_conductivity(self, T, L=np.inf, n_points=1000):
        """
        熱伝導率の計算（簡略化されたモデル）

        κ = (1/3) ∫ C(ω) v² τ(ω) g(ω) dω
        """
        # デバイ状態密度を仮定
        omega_max = self.omega_D
        omegas = np.linspace(0.01 * omega_max, omega_max, n_points)

        kappa = 0.0
        for omega in omegas:
            # 比熱（フォノンモードあたり）
            n = self.bose_einstein(omega, T)
            x = hbar * omega / (kB * T)
            if x < 1e-3:
                C_omega = kB
            else:
                C_omega = kB * x**2 * np.exp(x) / (np.exp(x) - 1)**2

            # 状態密度（デバイモデル）
            g_omega = 3 * omega**2 / omega_max**3

            # 緩和時間
            tau = 1.0 / self.total_rate(omega, T, L)

            # 寄与
            kappa += C_omega * self.v_s**2 * tau * g_omega

        # 数値積分（台形則）
        kappa *= (omega_max / n_points)

        return kappa / 3.0  # 因子1/3


# シリコンのパラメータ例
model_Si = PhononScatteringModel(
    theta_D=645,           # K
    mass=28.085 * 1.66054e-27,  # Si原子質量 [kg]
    a=5.43e-10,            # 格子定数 [m]
    gamma=1.0              # グリュナイゼンパラメータ
)

# 散乱率の温度依存性
temperatures = np.linspace(10, 600, 100)
omega_test = 0.5 * model_Si.omega_D  # テスト周波数

rates_U = [model_Si.umklapp_rate(omega_test, T) for T in temperatures]
rates_N = [model_Si.normal_rate(omega_test, T) for T in temperatures]

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

plt.subplot(1, 2, 1)
plt.semilogy(temperatures, rates_U, 'r-', label='U-process', linewidth=2)
plt.semilogy(temperatures, rates_N, 'b--', label='N-process', linewidth=2)
plt.xlabel('Temperature (K)', fontsize=12)
plt.ylabel('Scattering Rate (s⁻¹)', fontsize=12)
plt.title('Phonon-Phonon Scattering Rates', fontsize=14)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)

# 熱伝導率の温度依存性
kappas_bulk = [model_Si.thermal_conductivity(T, L=np.inf) for T in temperatures]
kappas_1mm = [model_Si.thermal_conductivity(T, L=1e-3) for T in temperatures]

plt.subplot(1, 2, 2)
plt.plot(temperatures, kappas_bulk, 'k-', label='Bulk (L=∞)', linewidth=2)
plt.plot(temperatures, kappas_1mm, 'g--', label='L=1mm', linewidth=2)
plt.xlabel('Temperature (K)', fontsize=12)
plt.ylabel('Thermal Conductivity (W/m·K)', fontsize=12)
plt.title('Temperature Dependence of κ', fontsize=14)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('phonon_scattering_rates.png', dpi=150, bbox_inches='tight')
plt.show()

# フォノン寿命の周波数依存性
omegas_scan = np.linspace(0.1 * model_Si.omega_D, 0.9 * model_Si.omega_D, 50)
T_room = 300  # K

lifetimes_300K = [1.0 / model_Si.total_rate(omega, T_room) for omega in omegas_scan]
linewidths_300K = [2 * hbar / tau for tau in lifetimes_300K]  # Γ = 2ℏ/τ

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

plt.subplot(1, 2, 1)
plt.semilogy(omegas_scan / model_Si.omega_D, np.array(lifetimes_300K) * 1e12, 'b-', linewidth=2)
plt.xlabel('ω/ω_D', fontsize=12)
plt.ylabel('Lifetime (ps)', fontsize=12)
plt.title(f'Phonon Lifetime at T = {T_room} K', fontsize=14)
plt.grid(True, alpha=0.3)

plt.subplot(1, 2, 2)
plt.plot(omegas_scan / model_Si.omega_D, np.array(linewidths_300K) * 1e12 / hbar, 'r-', linewidth=2)
plt.xlabel('ω/ω_D', fontsize=12)
plt.ylabel('Linewidth Γ/ℏ (THz)', fontsize=12)
plt.title('Spectral Linewidth', fontsize=14)
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('phonon_lifetime_linewidth.png', dpi=150, bbox_inches='tight')
plt.show()

print("計算完了！")
                

6.2 計算結果の解釈

7. 非調和性の実験的兆候

7.1 熱膨張

非調和性の最も直接的な証拠の一つは熱膨張です。体積熱膨張係数は：

ここで$\gamma$はグリュナイゼンパラメータで、非調和性の尺度です。調和系では$\gamma = 0$となり、熱膨張は起こりません。

7.2 フォノン周波数のシフト

温度上昇に伴い、フォノン周波数は変化します。これは主に：

• 体積効果: 熱膨張による格子定数の変化
• 純粋な非調和効果: フォノン-フォノン相互作用による繰り込み

周波数の温度依存性は：

7.3 分光線幅の温度依存性

ラマン散乱や中性子散乱スペクトルにおいて、線幅$\Gamma$は温度とともに増加します。 典型的な依存性は：

ここで、第二項は三フォノン過程、第三項は四フォノン過程に対応します。

7.4 熱伝導率の振る舞い

熱伝導率の温度依存性のパターン自体が、非調和性の証拠となります：

• 高温での$\kappa \propto 1/T$依存性はU過程の存在を示す
• 中間温度でのピークは、比熱増加とU過程増加の競合を反映
• 低温でのサイズ依存性は、内的散乱（U過程）が抑制されていることを示す

graph TD
    A[非調和性の実験的兆候] --> B[熱膨張]
    A --> C[周波数シフト]
    A --> D[線幅の温度依存性]
    A --> E[熱伝導率パターン]

    B --> F[α_V ≠ 0]
    C --> G[ω T依存性]
    D --> H[Γ T増加]
    E --> I[κ~1/T 高温]

    F --> J[グリュナイゼンパラメータ]
    G --> J
    H --> K[三・四フォノン過程]
    I --> K
                

7.5 測定技術