第3章: 熱輸送計算

3.1 フォノンボルツマン輸送方程式

3.1.1 分布関数とボルツマン方程式

フォノンの非平衡状態を記述するために、フォノン分布関数 $n(\mathbf{q}, j, \mathbf{r}, t)$ を導入します。 これは、時刻 $t$ に位置 $\mathbf{r}$ において、波数 $\mathbf{q}$、ブランチ $j$ のフォノンの 平均占有数を表します。

3.1.2 定常状態と線形応答

温度勾配 $\nabla T$ が小さい線形応答領域では、分布関数を平衡分布からの 小さな偏差として展開できます：

ここで、$n_0 = [\exp(\hbar\omega/k_B T) - 1]^{-1}$ はBose-Einstein分布です。 定常状態（$\partial n/\partial t = 0$）かつ外力がない場合、方程式は：

3.1.3 フーリエの法則との関係

熱流束 $\mathbf{J}_Q$ は分布関数を用いて以下のように表されます：

線形応答理論により、$\mathbf{J}_Q = -\kappa \nabla T$ の関係から熱伝導率テンソル $\kappa$ が得られます。

3.1.4 散乱項の構造

散乱項は、様々な散乱メカニズムの寄与の和として表されます：

主要な散乱メカニズム：

• フォノン-フォノン散乱: 非調和効果による三フォノン・四フォノン過程
• 境界散乱: 試料の有限サイズ効果
• 欠陥散乱: 点欠陥、転位、粒界など
• 電子-フォノン散乱: 金属や高温での寄与

3.2 緩和時間近似と反復解法

3.2.1 緩和時間近似（RTA）

最も単純な近似は、散乱項を緩和時間 $\tau(\mathbf{q}, j)$ で特徴づける方法です：

3.2.2 反復的ボルツマン輸送方程式（Iterative BTE）

より厳密なアプローチは、BTE を直接解く反復法です。分布関数の偏差を 変分パラメータ $\mathbf{F}(\mathbf{q}, j)$ で表します：

$\mathbf{F}$ に対する自己無撞着方程式は：

ここで、$\Gamma$ は散乱行列です。この方程式を反復的に解くことで、 運動量保存則を正確に取り入れた熱伝導率が得られます。

3.2.3 Normal vs Umklapp 過程

反復解法では Normal 過程を適切に扱うことで、特に高純度結晶や 低温領域で RTA よりも正確な熱伝導率を予測できます。

3.3 第一原理からの三フォノン散乱率

3.3.1 三次の力定数

三フォノン散乱は、結晶ポテンシャルの三次項から生じます：

3.3.2 散乱振幅とFermiの黄金律

三次の力定数から、フォノン相互作用ハミルトニアンの行列要素が得られます：

Fermi の黄金律により、散乱率は：

3.3.3 選択則と保存則

3.3.4 有限変位法による三次力定数計算

実用的には、対称性を考慮した有限変位パターンを用いて三次力定数を計算します：

1. スーパーセルを構築（典型的には 2×2×2 から 4×4×4）
2. 原子ペアを変位させた配置を生成（対称性により独立な配置のみ）
3. 各配置で力を DFT 計算
4. 力の差分から三次力定数を数値的に求める

# 1. 変位パターン生成
phono3py --dim="2 2 2" -c POSCAR-unitcell

# 2. 各変位配置で DFT 計算（VASP例）
# POSCAR-{00001..00XXX} に対して自己無撞着計算

# 3. 力を収集して力定数を計算
phono3py --cf3 vasprun.xml-{00001..00XXX}

# 4. 散乱率と熱伝導率を計算
phono3py --br --mesh="20 20 20" --ts="300"

3.4 四フォノン散乱とその重要性

3.4.1 なぜ四フォノン過程が重要か

従来、三フォノン過程のみが考慮されてきましたが、近年の研究で 以下の系で四フォノン散乱が支配的になることが明らかになりました：

• 高温領域: フォノン占有数が大きくなると四フォノン過程の確率が増大
• 軽元素材料: Si, C, BN など高周波数フォノンを持つ系
• 弱い三フォノン散乱: 高対称性により三フォノン過程が抑制される系
• 二次元材料: グラフェン、h-BN など

3.4.2 四次の力定数と散乱率

四フォノン過程は四次の力定数 $\Phi^{(4)}$ から生じます：

四次力定数の計算は計算コストが高い（変位配置数が $O(N^3)$）ため、 機械学習ポテンシャルや対称性を最大限活用した効率的手法が開発されています。

3.4.3 実装と計算ツール

四フォノン散乱を扱えるソフトウェア：

• Phono3py (v2.0+): 四フォノン散乱に対応（FourPhonon オプション）
• ShengBTE: 一部の四フォノン過程をサポート
• almaBTE: 完全な四フォノン実装

# 四次力定数の計算（大きなスーパーセルが必要）
phono3py --dim="3 3 3" -c POSCAR --fc4

# 四フォノン散乱を含む熱伝導率計算
phono3py --br --mesh="24 24 24" --ts="100 200 300 500 1000" \
         --fc3 --fc4 --write-gamma

3.5 ShengBTE: 理論と実用

3.5.1 ShengBTE の概要

ShengBTE は、第一原理三次力定数から熱伝導率を計算するための 高性能Fortranコードです。Wu Li らによって開発され、広く使われています。

3.5.2 ShengBTE のワークフロー

3.5.3 CONTROL ファイルの設定

# ShengBTE CONTROL file for Silicon
&allocations
  nelements=1      # 元素数
  natoms=2         # 単位格子内の原子数
  ngrid(:)=20 20 20  # q点メッシュ
  norientations=0  # 異方性を考慮しない場合
&end

&crystal
  lfactor=0.5431   # 格子定数 (nm)
  lattvec(:,1)=0.0 0.5 0.5  # 格子ベクトル
  lattvec(:,2)=0.5 0.0 0.5
  lattvec(:,3)=0.5 0.5 0.0
  types(:)=1 1     # 原子タイプ
  elements='Si'    # 元素名
  positions(:,1)=0.00 0.00 0.00  # 原子位置
  positions(:,2)=0.25 0.25 0.25
  scell(:)=5 5 5   # スーパーセル（力定数計算用）
  masses(:)=28.0855  # 原子質量
&end

¶meters
  T=300.0          # 温度 (K)
  scalebroad=0.5   # ブロードニング係数
  isotopes=.TRUE.  # 同位体散乱を含む
  onlyharmonic=.FALSE.  # 非調和効果を含む
  convergence=1e-4 # 反復解法の収束条件
  maxiter=100      # 最大反復回数
&end

&flags
  espresso=.FALSE. # VASP を使用
  nonanalytic=.TRUE.  # LO-TO分裂を考慮
&end

3.5.4 結果の解釈

ShengBTE の主要な出力ファイル：

• BTE.kappa: 温度の関数としての熱伝導率テンソル
• BTE.w_final: 各モードの散乱率（逆緩和時間）
• BTE.kappa_tensor: 完全な3×3熱伝導率テンソル
• BTE.cumulative_kappa_*: 平均自由行程の関数としての累積熱伝導率

3.6 almaBTE と Phono3py

3.6.1 Phono3py の特徴

Phono3py は Phonopy の拡張として開発され、Python で記述された 使いやすいフォノン-フォノン相互作用計算ツールです。

3.6.2 Phono3py の基本的な使い方

# ステップ1: 変位ファイル生成
phono3py --dim="2 2 2" -c POSCAR-unitcell

# ステップ2: DFT計算（各 POSCAR-{00001..00XXX} で）
# （VASP INCAR設定: PREC=Accurate, ENCUT=高め, tight convergence）

# ステップ3: 力定数の計算
phono3py --cf3 vasprun.xml-{00001..00XXX}

# ステップ4: 熱伝導率計算
phono3py --mesh="19 19 19" --br --ts="200 300 500 800"

# ステップ5: RTA vs 反復解法の比較
phono3py --mesh="19 19 19" --br --ts="300" --write-gamma
phono3py --mesh="19 19 19" --br --ts="300" --write-gamma --full-pp

# ステップ6: 累積熱伝導率（MFP分布）
phono3py --mesh="19 19 19" --br --ts="300" --write-kappa

# ステップ7: 四フォノン散乱を含む計算
phono3py --dim="3 3 3" --fc4 -c POSCAR-unitcell  # 変位生成
# DFT計算...
phono3py --cf4 vasprun.xml-{00001..XXXXX}  # 四次力定数
phono3py --mesh="19 19 19" --br --ts="300 500 1000" --fc3 --fc4

3.6.3 almaBTE: 高度な解析ツール

almaBTE は、より詳細な熱輸送解析のための C++ ライブラリとツール群です。 特徴：

• サイズ効果: ナノ構造や薄膜の熱伝導率計算
• コヒーレント効果: 波動的輸送の取り扱い
• 時間依存輸送: 過渡的熱応答のシミュレーション
• 可視化ツール: 散乱率、MFP、寄与解析の高度な可視化

import almabte

# ShengBTEやPhono3pyの結果を読み込み
solver = almabte.Solver.from_shengbte("./")

# 詳細な平均自由行程分布を計算
mfp_spectrum = solver.compute_mfp_spectrum(T=300)

# サイズ効果の評価（Casimir限界）
kappa_L = solver.compute_kappa_vs_length(T=300, L_range=[10e-9, 10e-6])

# バリスティック-拡散的輸送の転移
regime_map = solver.analyze_transport_regime(T=300)

3.7 収束問題と解決策

3.7.1 主要な収束パラメータ

第一原理熱伝導率計算では、以下のパラメータに対する収束確認が必須です：

3.7.2 q点メッシュ収束テスト

import numpy as np
import subprocess

# テストするメッシュサイズ
mesh_sizes = [11, 15, 19, 23, 27, 31]
kappa_values = []

for mesh in mesh_sizes:
    # Phono3py実行
    cmd = f"phono3py --mesh='{mesh} {mesh} {mesh}' --br --ts=300"
    subprocess.run(cmd, shell=True, check=True)

    # 結果読み込み
    kappa = np.loadtxt("kappa-m{0:02d}{0:02d}{0:02d}.hdf5".format(mesh))
    kappa_avg = np.mean(kappa[0, 1:4])  # 300Kでの平均
    kappa_values.append(kappa_avg)

    print(f"Mesh {mesh}³: κ = {kappa_avg:.2f} W/mK")

# 収束判定
kappa_values = np.array(kappa_values)
convergence = np.abs(kappa_values[-1] - kappa_values[-2]) / kappa_values[-1] * 100
print(f"\n最終収束: {convergence:.2f}%")

if convergence < 5.0:
    print(f"✓ 収束達成（推奨メッシュ: {mesh_sizes[-2]}³）")
else:
    print(f"✗ さらに高密度メッシュが必要")

3.7.3 収束を妨げる一般的な問題

3.7.4 四面体法とスミアリング法

エネルギー保存デルタ関数 $\delta(\omega \pm \omega' \mp \omega'')$ の 評価方法が精度に大きく影響します：

• Gaussian スミアリング: 実装簡単、収束遅い
• 適応的ブロードニング: $\sigma(\mathbf{q}) = C \cdot \omega(\mathbf{q})$
• 四面体法: 最も正確だが実装複雑（ShengBTE の強み）

3.8 平均自由行程分布解析

3.8.1 平均自由行程の定義

各フォノンモード $(\mathbf{q}, j)$ の平均自由行程（MFP）は：

これは、フォノンが散乱されるまでに移動する平均距離を表します。 MFP 分布は材料のナノスケール熱輸送特性を理解する鍵となります。

3.8.2 累積熱伝導率

MFP が $\Lambda$ 以下のモードの寄与を積算した累積熱伝導率：

ここで、$\Theta$ はステップ関数です。$\kappa_\text{cum}(\Lambda)$ を プロットすることで、どの長さスケールのフォノンが熱伝導に寄与しているかが分かります。

3.8.3 サイズ効果の予測

試料サイズ $L$ が MFP と同程度になると、境界散乱が支配的になり 熱伝導率が低下します（Casimir 限界）：

これにより、ナノワイヤや薄膜の熱伝導率をバルク値から予測できます。

import numpy as np
import matplotlib.pyplot as plt
import h5py

# Phono3pyのHDF5出力を読み込み
with h5py.File("kappa-m191919.hdf5", "r") as f:
    temperatures = f["temperature"][:]
    gamma = f["gamma"][:]  # 散乱率
    group_velocity = f["group_velocity"][:]
    heat_capacity = f["heat_capacity"][:]
    qpoints = f["qpoint"][:]
    weights = f["weight"][:]

T_idx = 1  # 300 K を選択
v_norm = np.linalg.norm(group_velocity, axis=2)
mfp = v_norm / (gamma[:, :, T_idx] + 1e-12)  # MFP = v / Gamma

# 各モードの熱伝導率寄与
kappa_mode = heat_capacity[:, :, T_idx] * v_norm**2 / gamma[:, :, T_idx, None]
kappa_mode = np.nanmean(kappa_mode, axis=2) * weights[:, None]

# MFPをソートして累積分布計算
mfp_flat = mfp.flatten()
kappa_flat = kappa_mode.flatten()
sort_idx = np.argsort(mfp_flat)
mfp_sorted = mfp_flat[sort_idx]
kappa_cumulative = np.cumsum(kappa_flat[sort_idx])
kappa_cumulative /= kappa_cumulative[-1]  # 正規化

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

# 累積分布
ax1.semilogx(mfp_sorted * 1e9, kappa_cumulative * 100)
ax1.set_xlabel("Mean Free Path (nm)")
ax1.set_ylabel("Cumulative κ (%)")
ax1.set_title("Cumulative Thermal Conductivity")
ax1.grid(True, alpha=0.3)

# ヒストグラム（微分分布）
bins = np.logspace(-1, 4, 50)  # 0.1 nm to 10 μm
hist, bin_edges = np.histogram(mfp_sorted * 1e9, bins=bins, weights=kappa_flat[sort_idx])
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
ax2.semilogx(bin_centers, hist / np.sum(hist))
ax2.set_xlabel("Mean Free Path (nm)")
ax2.set_ylabel("κ Contribution (normalized)")
ax2.set_title("MFP Distribution")
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("mfp_distribution.png", dpi=300)
print(f"50% of κ from MFP < {mfp_sorted[np.argmin(np.abs(kappa_cumulative - 0.5))] * 1e9:.1f} nm")

3.9 コヒーレント vs インコヒーレント輸送

3.9.1 輸送レジームの分類

フォノン輸送は、長さスケールと温度に応じて異なるレジームを示します：

3.9.2 コヒーレント輸送の条件

フォノンが波としての性質を保つ（位相コヒーレンス）条件：

• 低温: フォノン-フォノン散乱が弱い（$k_B T \ll \hbar\omega$）
• 短距離: コヒーレンス長 $l_\phi$ 内での伝播
• 周期構造: 超格子、フォノニック結晶など
• クリーンな界面: 無秩序散乱が少ない

3.9.3 コヒーレント効果の例

3.9.4 AGF法: コヒーレント輸送の計算

原子グリーン関数（Atomistic Green's Function, AGF）法は、 波動輸送を直接扱う手法です：

ここで、$K$ は力定数行列、$\Sigma_{L/R}$ は左右のリードからの自己エネルギーです。 熱伝導率は透過係数 $T(\omega)$ から：

3.10 熱伝導率予測のための機械学習

3.10.1 なぜ機械学習か

第一原理BTE計算は高精度ですが、計算コストが高く（数日〜数週間/材料）、 大規模材料スクリーニングには不向きです。機械学習（ML）により：

• 高速予測（秒〜分/材料）
• 記述子ベースの物理的解釈
• 未知材料への外挿
• 逆設計（目標物性から構造を提案）

3.10.2 機械学習アプローチの分類

3.10.3 記述子ベース機械学習

材料の構造・組成情報から特徴量を抽出し、熱伝導率を予測：

import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from matminer.featurizers.composition import ElementProperty
from matminer.featurizers.structure import SiteStatsFingerprint

# データセット読み込み（例：Materials Projectから収集）
# compositions: 組成のリスト
# structures: pymatgen Structure objects
# kappa_values: 第一原理計算された熱伝導率

# 特徴量生成
featurizer = ElementProperty.from_preset("magpie")
X_comp = featurizer.featurize_many(compositions, ignore_errors=True)

struct_featurizer = SiteStatsFingerprint.from_preset("CoordinationNumber")
X_struct = [struct_featurizer.featurize(s) for s in structures]

# 特徴量結合
X = np.hstack([X_comp, X_struct])
y = np.array(kappa_values)

# 訓練・テスト分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# モデル訓練
model = RandomForestRegressor(n_estimators=200, max_depth=20, random_state=42)
model.fit(X_train, y_train)

# 評価
train_score = model.score(X_train, y_train)
test_score = model.score(X_test, y_test)
print(f"Train R²: {train_score:.3f}, Test R²: {test_score:.3f}")

# 特徴量重要度
feature_importance = model.feature_importances_
top_features = np.argsort(feature_importance)[-10:]
print("Top 10 important features:")
for idx in reversed(top_features):
    print(f"  {featurizer.feature_labels()[idx]}: {feature_importance[idx]:.3f}")

3.10.4 グラフニューラルネットワーク（GNN）

結晶構造を原子ノードと結合エッジからなるグラフとして表現し、 GNN で直接物性を予測する手法が注目されています。

3.10.5 機械学習力場 + BTE ハイブリッド

最先端のアプローチは、機械学習力場（MLIP）で力定数を高速計算し、 その後BTE解法で熱伝導率を求める二段階法です：

1. 少数の第一原理計算で MLIP を訓練（GAP, SNAP, NequIP など）
2. MLIP を用いて大きなスーパーセルでの力定数計算（DFT の 100-1000 倍高速）
3. 得られた力定数で通常の BTE 計算

3.10.6 ハイスループットスクリーニング

ML モデルを用いて、材料データベース（Materials Project, AFLOWなど）の 数十万材料の熱伝導率を一括予測し、有望候補を絞り込む戦略：

from pymatgen.ext.matproj import MPRester
import pandas as pd

# Materials Project から半導体を取得
mpr = MPRester("YOUR_API_KEY")
entries = mpr.query(
    criteria={"band_gap": {"$gt": 0.5, "$lt": 2.0}},  # 半導体
    properties=["material_id", "pretty_formula", "structure", "band_gap"]
)

# MLモデルで熱伝導率予測（上記で訓練済み）
structures = [e["structure"] for e in entries]
X_pred = generate_features(structures)  # 特徴量生成関数
kappa_pred = model.predict(X_pred)

# 結果を DataFrame に
df = pd.DataFrame({
    "material_id": [e["material_id"] for e in entries],
    "formula": [e["pretty_formula"] for e in entries],
    "band_gap": [e["band_gap"] for e in entries],
    "kappa_predicted": kappa_pred
})

# 低熱伝導率材料を抽出（熱電材料候補）
candidates = df[df["kappa_predicted"] < 2.0].sort_values("kappa_predicted")
print(f"Found {len(candidates)} low-κ candidates")
print(candidates.head(20))

# 上位候補を第一原理で検証
for idx, row in candidates.head(5).iterrows():
    print(f"\nValidating {row['formula']} (mp-{row['material_id']})")
    # ここで詳細なDFT+BTE計算を実行...

3.11 実践例：Si と PbTe の熱伝導率

3.11.1 ケーススタディ1: Silicon (Si)

シリコンは高熱伝導率材料の代表例で、熱伝導率計算の精度検証に 広く用いられています。

3.11.2 Si の解析スクリプト

import numpy as np
import h5py
import matplotlib.pyplot as plt

# Phono3py 結果の読み込み
def load_phono3py_data(filename):
    with h5py.File(filename, "r") as f:
        data = {
            "temperature": f["temperature"][:],
            "kappa": f["kappa"][:],  # shape: (T, 6)
            "gamma": f["gamma"][:],
            "group_velocity": f["group_velocity"][:],
            "frequencies": f["frequency"][:],
            "qpoints": f["qpoint"][:],
        }
    return data

# データ読み込み
data_rta = load_phono3py_data("kappa-m252525.hdf5")
data_full = load_phono3py_data("kappa-m252525-full.hdf5")

# 温度依存性プロット
T = data_rta["temperature"]
kappa_rta = np.mean(data_rta["kappa"][:, 1:4], axis=1)  # xx, yy, zz 平均
kappa_full = np.mean(data_full["kappa"][:, 1:4], axis=1)

plt.figure(figsize=(10, 6))
plt.loglog(T, kappa_rta, 'o-', label="RTA", linewidth=2)
plt.loglog(T, kappa_full, 's-', label="Iterative BTE", linewidth=2)

# 実験データ（文献値）
T_exp = np.array([100, 200, 300, 400, 500, 800, 1000])
kappa_exp = np.array([900, 260, 145, 95, 70, 40, 30])
plt.loglog(T_exp, kappa_exp, 'k^', label="Experiment", markersize=8)

plt.xlabel("Temperature (K)", fontsize=14)
plt.ylabel("Thermal Conductivity (W/mK)", fontsize=14)
plt.title("Silicon Thermal Conductivity", fontsize=16)
plt.legend(fontsize=12)
plt.grid(True, alpha=0.3)
plt.savefig("Si_kappa_vs_T.png", dpi=300)

# モード分解分析（300 K）
T_idx = np.argmin(np.abs(T - 300))
freq = data_rta["frequencies"]
gamma = data_rta["gamma"][:, :, T_idx]
v_g = np.linalg.norm(data_rta["group_velocity"], axis=2)

# 音響 vs 光学モード（6 branches: 3 acoustic + 3 optical）
freq_threshold = 8.0  # THz（LA/TA と LO/TO の境界）
acoustic_mask = freq < freq_threshold

gamma_ac = gamma[acoustic_mask]
gamma_op = gamma[~acoustic_mask]

print(f"300 K での散乱率統計:")
print(f"  音響モード: 平均 Γ = {np.mean(gamma_ac):.2e} THz")
print(f"  光学モード: 平均 Γ = {np.mean(gamma_op):.2e} THz")
print(f"  音響/光学 比: {np.mean(gamma_op) / np.mean(gamma_ac):.1f}")

3.11.3 ケーススタディ2: PbTe（熱電材料）

PbTe は代表的な熱電材料で、低熱伝導率と高い電気伝導率を併せ持ちます。

3.11.4 PbTe の計算実例

#!/bin/bash
# PbTe 熱伝導率計算の完全なワークフロー

# ステップ1: 構造最適化
cat > INCAR_relax << EOF
ISTART = 0
PREC = Accurate
ENCUT = 400
EDIFF = 1E-8
EDIFFG = -1E-6
NSW = 200
IBRION = 2
ISIF = 3
ISMEAR = 0
SIGMA = 0.05
LREAL = .FALSE.
EOF

# VASP 実行...

# ステップ2: Phono3py 変位生成
phono3py --dim="3 3 3" -c POSCAR-unitcell --pa="0 1/2 1/2 1/2 0 1/2 1/2 1/2 0"

# ステップ3: 力計算用 INCAR（タイトな収束条件）
cat > INCAR_forces << EOF
ISTART = 0
PREC = Accurate
ENCUT = 450
EDIFF = 1E-9
ISMEAR = 0
SIGMA = 0.02
LREAL = .FALSE.
ADDGRID = .TRUE.
EOF

# 各変位配置で力計算（並列化推奨）
for i in {00001..00XXX}; do
    mkdir disp-$i
    cp POSCAR-$i disp-$i/POSCAR
    cd disp-$i
    # VASP 実行
    cd ..
done

# ステップ4: 力定数計算
phono3py --cf3 disp-*/vasprun.xml

# ステップ5: 熱伝導率計算（RTA）
phono3py --mesh="21 21 21" --br --ts="300 400 500 600 700 800"

# ステップ6: 反復解法
phono3py --mesh="21 21 21" --br --ts="300 500 800" --full-pp --write-gamma

# ステップ7: MFP分布解析
phono3py --mesh="21 21 21" --br --ts="300" --write-kappa --gv-delta-q=1e-4

3.11.5 Si vs PbTe の比較