第3章：フォノン状態密度

1. 状態密度とは何か？

1.1 物理的イメージ

第2章で学んだフォノン分散関係 \(\omega(\mathbf{k})\) は、各波数ベクトル \(\mathbf{k}\) に対して どのような振動数のフォノンが存在するかを教えてくれます。しかし、材料の熱的性質を理解するには、 「特定の振動数 \(\omega\) を持つフォノンモードがどれだけ存在するか」という情報が重要です。

1.2 フォノン状態密度の定義

フォノン状態密度（Density of States, DOS）\(g(\omega)\) は、単位周波数あたりに 存在するフォノンモードの数として定義されます：

ここで：

• \(\mathbf{k}\)：波数ベクトル
• \(s\)：フォノンブランチのインデックス（音響・光学モード）
• \(\omega_s(\mathbf{k})\)：ブランチ \(s\)、波数 \(\mathbf{k}\) のフォノン振動数
• \(\delta(x)\)：ディラックのデルタ関数

実際の計算では、ブリルアンゾーン内の \(\mathbf{k}\) 点を離散的にサンプリングし、 積分形式に書き換えます：

1.3 DOSと分散関係の関係

状態密度は分散関係の「逆関数」的な情報を含んでいます。分散が平坦な領域 （\(\nabla_\mathbf{k} \omega \approx 0\)）では、多くの \(\mathbf{k}\) 点が 同じ振動数 \(\omega\) に対応するため、DOSが大きくなります。

1.4 規格化条件

全フォノンモード数は \(3N\)（\(N\) は単位胞あたりの原子数）に等しいため、 DOSは次の規格化条件を満たします：

2. Van Hove特異点

2.1 特異点の起源

Van Hove特異点は、分散関係が極値（最大値または最小値）を持つ点で発生します。 これらの点では、\(\nabla_\mathbf{k} \omega_s(\mathbf{k}) = 0\) となり、 状態密度が特異的に増大します。

2.2 3次元結晶における特異点の分類

3次元結晶では、以下の4種類のVan Hove特異点が存在します：

2.3 物理的重要性

Van Hove特異点は以下の現象に重要な役割を果たします：

• 比熱の温度依存性：特異点での状態密度の増大が比熱に寄与
• 中性子散乱実験：散乱強度がDOSに比例するため、特異点が顕著に観測される
• フォノン-フォノン散乱：特異点近傍でフォノン寿命が変化

3. Debyeモデル

3.1 モデルの仮定

Debyeモデルは、1912年にPeter Debyeによって提案された、 フォノンDOSを近似するための簡単なモデルです。以下の仮定に基づいています：

1. 線形分散近似：すべてのフォノンモードが音速 \(v\) で線形に分散する $$\omega = v|\mathbf{k}|$$
2. 等方的媒質：音速はすべての方向で同じ
3. 単一音速：縦波と横波の音速を平均化（または3つの音響ブランチを考慮）
4. カットオフ周波数：全モード数を \(3N\) に保つため、最大周波数 \(\omega_D\)（Debye周波数）を導入

3.2 Debye状態密度の導出

線形分散 \(\omega = v k\) を仮定すると、3次元での状態密度は：

\(k = \omega/v\) より \(dk/d\omega = 1/v\) なので：

ここで、規格化条件から Debye周波数 \(\omega_D\) を決定します：

これより：

3.3 Debye温度

Debye周波数に対応する温度をDebye温度 \(\Theta_D\) と呼びます：

Debye温度は材料の固有の特性であり、以下の意味を持ちます：

• \(T \ll \Theta_D\)：低温極限（量子効果が顕著）
• \(T \gg \Theta_D\)：高温極限（古典的挙動）
• 比熱の温度依存性を決定する重要なパラメータ

3.4 Debyeモデルの物理的解釈

Debyeモデルは、結晶を等方的弾性体として扱います。 この近似は以下の場合に良く機能します：

• 低周波数（長波長）フォノン：分散が線形に近い
• 高温での比熱：細かいDOS構造が平均化される
• 単純な結晶構造（単原子基底）

一方、以下の場合には不適切です：

• 光学フォノンが重要な多原子系
• Van Hove特異点の影響が大きい場合
• 強い異方性を持つ結晶

4. Einsteinモデル

4.1 モデルの仮定

Einsteinモデルは、1907年にAlbert Einsteinによって提案された、 さらに簡単な近似です。このモデルでは：

1. 単一周波数仮定：すべてのフォノンモードが同じ周波数 \(\omega_E\) を持つ
2. 独立調和振動子：各原子が独立に振動する

この仮定の下では、状態密度は以下のように書けます：

4.2 Einstein温度

Einstein周波数に対応する温度をEinstein温度 \(\Theta_E\) と呼びます：

4.3 物理的解釈と適用範囲

Einsteinモデルは、以下の状況で有用です：

• 光学フォノン：分散が平坦（\(\omega \approx \text{const}\)）なモードの記述
• 高温極限：古典的振動子の挙動
• 局在振動：不純物や欠陥周辺の局在モード

しかし、低温での比熱を正しく記述できないという重大な欠点があります。 実験的には、低温で \(C_V \propto T^3\) となる（Debyeの \(T^3\) 則）のに対し、 Einsteinモデルは指数関数的に減衰する予測を与えます。

5. DebyeモデルとEinsteinモデルの比較

5.1 定性的比較

5.2 視覚的比較

5.3 実材料への適用：混合モデル

実際の材料では、DebyeモデルとEinsteinモデルを組み合わせることで、 より良い近似が得られます：

ここで、\(g_{\text{Debye}}(\omega)\) は音響フォノンを記述し、 デルタ関数の項は光学フォノンを表します。

6. 実際のDOSと近似モデルの比較

6.1 単原子結晶（例：FCC金属）

面心立方（FCC）構造を持つアルミニウムのような単原子結晶では、 実際のDOSは以下の特徴を示します：

• 低周波数領域：Debye近似 \(g(\omega) \propto \omega^2\) に従う
• 中間領域：Van Hove特異点による鋭いピーク
• 高周波数端：分散が平坦になる領域で状態密度が増大

6.2 二原子結晶（例：NaCl型イオン結晶）

NaCl型構造を持つイオン結晶では：

• 音響フォノン：低周波数に3つのブランチ（1縦波、2横波）
• ギャップ：音響モードと光学モードの間に周波数ギャップ
• 光学フォノン：高周波数にピーク（Einstein的）

6.3 実DOSの計算方法

第一原理計算では、以下の手順でDOSを計算します：

1. 密度汎関数摂動論（DFPT）でフォノン分散 \(\omega_s(\mathbf{k})\) を計算
2. ブリルアンゾーン内で \(\mathbf{k}\) 点を密にサンプリング
3. 各 \(\mathbf{k}\) 点でのフォノン周波数をヒストグラム化
4. 適切な平滑化（Gaussian smearing など）を適用

7. PythonでフォノンDOSを計算する

7.1 Debye DOSのプロット

まず、Debyeモデルの状態密度を実装してみましょう。

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
N = 1000  # 原子数
omega_D = 20.0  # Debye周波数（THz）

# 周波数範囲
omega = np.linspace(0, omega_D * 1.5, 500)

# Debye DOS
def debye_dos(omega, omega_D, N):
    """
    Debyeモデルの状態密度

    Parameters:
    -----------
    omega : array
        周波数配列（THz）
    omega_D : float
        Debye周波数（THz）
    N : int
        原子数

    Returns:
    --------
    g : array
        状態密度
    """
    g = np.zeros_like(omega)
    mask = omega <= omega_D
    g[mask] = (9 * N / omega_D**3) * omega[mask]**2
    return g

# DOSを計算
g_debye = debye_dos(omega, omega_D, N)

# プロット
plt.figure(figsize=(10, 6))
plt.plot(omega, g_debye, linewidth=2, label='Debye model')
plt.axvline(omega_D, color='red', linestyle='--',
            label=f'$\\omega_D$ = {omega_D} THz')
plt.xlabel('Frequency ω (THz)', fontsize=12)
plt.ylabel('DOS g(ω) (states/THz)', fontsize=12)
plt.title('Debye Phonon Density of States', fontsize=14)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('debye_dos.png', dpi=300, bbox_inches='tight')
plt.show()

# 規格化チェック
total_modes = np.trapz(g_debye, omega)
print(f"Total modes: {total_modes:.1f} (expected: {3*N})")

7.2 EinsteinモデルとDebyeモデルの比較

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
N = 1000
omega_D = 20.0  # Debye周波数
omega_E = 15.0  # Einstein周波数

omega = np.linspace(0, 30, 1000)

# Debye DOS
def debye_dos(omega, omega_D, N):
    g = np.zeros_like(omega)
    mask = omega <= omega_D
    g[mask] = (9 * N / omega_D**3) * omega[mask]**2
    return g

# Einstein DOS（Gaussianで近似）
def einstein_dos(omega, omega_E, N, sigma=0.5):
    """
    EinsteinモデルのDOS（Gaussian近似）

    Parameters:
    -----------
    sigma : float
        Gaussianの幅（THz）
    """
    return 3 * N / (sigma * np.sqrt(2 * np.pi)) * \
           np.exp(-0.5 * ((omega - omega_E) / sigma)**2)

# 計算
g_debye = debye_dos(omega, omega_D, N)
g_einstein = einstein_dos(omega, omega_E, N)

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

# 左：個別プロット
axes[0].plot(omega, g_debye, linewidth=2, label='Debye model')
axes[0].plot(omega, g_einstein, linewidth=2, label='Einstein model')
axes[0].axvline(omega_D, color='blue', linestyle='--', alpha=0.5)
axes[0].axvline(omega_E, color='orange', linestyle='--', alpha=0.5)
axes[0].set_xlabel('Frequency ω (THz)', fontsize=12)
axes[0].set_ylabel('DOS g(ω)', fontsize=12)
axes[0].set_title('Debye vs Einstein Models', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)

# 右：積分DOS（累積モード数）
cumulative_debye = np.array([np.trapz(g_debye[:i+1], omega[:i+1])
                               for i in range(len(omega))])
cumulative_einstein = np.array([np.trapz(g_einstein[:i+1], omega[:i+1])
                                 for i in range(len(omega))])

axes[1].plot(omega, cumulative_debye, linewidth=2, label='Debye')
axes[1].plot(omega, cumulative_einstein, linewidth=2, label='Einstein')
axes[1].axhline(3*N, color='red', linestyle='--',
                label=f'Total modes = {3*N}')
axes[1].set_xlabel('Frequency ω (THz)', fontsize=12)
axes[1].set_ylabel('Cumulative modes', fontsize=12)
axes[1].set_title('Integrated DOS', fontsize=14)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('debye_vs_einstein.png', dpi=300, bbox_inches='tight')
plt.show()

7.3 実材料のDOS（単原子鎖モデル）

第2章で学んだ単原子鎖の分散関係から、実際のDOSを計算してみましょう。

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
M = 1.0      # 質量（任意単位）
K = 10.0     # バネ定数（任意単位）
a = 1.0      # 格子定数（任意単位）
N_atoms = 100  # 原子数

# 単原子鎖の分散関係
def monatomic_dispersion(k, K, M, a):
    """単原子鎖のフォノン分散"""
    return 2 * np.sqrt(K / M) * np.abs(np.sin(k * a / 2))

# ブリルアンゾーン内でk点をサンプリング
k_points = np.linspace(-np.pi/a, np.pi/a, 1000)
omega_k = monatomic_dispersion(k_points, K, M, a)

# DOSをヒストグラムで計算
omega_max = 2 * np.sqrt(K / M)
omega_bins = np.linspace(0, omega_max * 1.1, 200)
dos_hist, bin_edges = np.histogram(omega_k, bins=omega_bins, density=True)
omega_centers = 0.5 * (bin_edges[1:] + bin_edges[:-1])

# 規格化（全モード数 = N_atoms）
dos_normalized = dos_hist * N_atoms / np.trapz(dos_hist, omega_centers)

# 解析的DOS（1D単原子鎖）
def analytical_dos_1d(omega, omega_max, N):
    """
    1次元単原子鎖の解析的DOS
    g(ω) = N / (π√(ω_max² - ω²))
    """
    g = np.zeros_like(omega)
    mask = (omega > 0) & (omega < omega_max)
    g[mask] = N / (np.pi * np.sqrt(omega_max**2 - omega[mask]**2))
    return g

omega_analytical = np.linspace(0.01, omega_max * 0.99, 500)
dos_analytical = analytical_dos_1d(omega_analytical, omega_max, N_atoms)

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

# 左：分散関係
axes[0].plot(k_points * a / np.pi, omega_k, linewidth=2)
axes[0].set_xlabel('Wave vector ka/π', fontsize=12)
axes[0].set_ylabel('Frequency ω', fontsize=12)
axes[0].set_title('Phonon Dispersion (Monatomic Chain)', fontsize=14)
axes[0].grid(True, alpha=0.3)
axes[0].axhline(omega_max, color='red', linestyle='--', alpha=0.5,
                label=f'$\\omega_{{max}}$ = {omega_max:.2f}')
axes[0].legend()

# 右：DOS
axes[1].plot(omega_centers, dos_normalized, 'o-',
             label='Numerical (histogram)', markersize=3)
axes[1].plot(omega_analytical, dos_analytical, linewidth=2,
             label='Analytical', color='red', alpha=0.7)
axes[1].set_xlabel('Frequency ω', fontsize=12)
axes[1].set_ylabel('DOS g(ω)', fontsize=12)
axes[1].set_title('Phonon Density of States', fontsize=14)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)

# Van Hove特異点（ω = ω_maxで発散）を強調
axes[1].axvline(omega_max, color='green', linestyle='--', alpha=0.5,
                label='Van Hove singularity')

plt.tight_layout()
plt.savefig('monatomic_chain_dos.png', dpi=300, bbox_inches='tight')
plt.show()

print(f"Maximum frequency: {omega_max:.4f}")
print(f"Total modes (numerical): {np.trapz(dos_normalized, omega_centers):.1f}")

8. DOSの実験的測定方法

8.1 非弾性中性子散乱（INS）

中性子散乱は、フォノンDOSを直接測定する最も強力な実験手法です。 中性子が試料中のフォノンと相互作用し、エネルギーと運動量を交換します。

8.2 ラマン分光

ラマン分光では、光がフォノンと相互作用してエネルギーを交換します。 特に光学フォノンの測定に適しています。

8.3 比熱測定からの推定

比熱の温度依存性を測定し、逆問題を解くことでDOSを推定できます。 特に低温での \(C_V(T)\) データから、低周波数DOSの情報が得られます。

8.4 測定方法の比較

まとめ

この章では、フォノン状態密度（DOS）の概念と計算方法を学びました。 重要なポイントをまとめます：