シリコンの負熱膨張とグリューナイゼン定数の第一原理計算

熱膨張は、精密光学・半導体デバイス・構造材料の設計を左右する基本的な材料物性です。物質の熱膨張は結晶格子の非調和性に由来し、原子間ポテンシャルを 2 次で打ち切る調和近似の範囲では記述できません（調和振動子は温度によらず平均原子間距離が一定のため、熱膨張はゼロになります）。特にシリコン（Si）は、室温付近では正の熱膨張を示す一方、約 18–125 K の低温域で線膨張係数が負になる「負熱膨張（NTE: Negative Thermal Expansion）」を示すことで知られます [1]。これは横音響（TA）フォノンの異常な体積依存性に起因する現象です [2, 3]。本事例では、第一原理計算ソフトウェア Advance/PHASE を力・エネルギーの評価エンジンとし、フォノン解析パッケージ ALAMODE と連携させて、Si のモードグリューナイゼン定数 $\gamma_j(\mathbf{q})$ を準調和近似（3 体積・中心差分）で求めました。そこから熱力学的グリューナイゼン定数 $\bar\gamma(T)$ と線膨張係数 $\alpha(T)$ を算出し、Si の負熱膨張を第一原理から再現します。

Keywords: 第一原理計算 (DFT), グリューナイゼン定数, 熱膨張係数, 負熱膨張 (NTE), 準調和近似 (QHA), フォノン, ALAMODE, シリコン

1. 計算方法：準調和近似によるグリューナイゼン定数と熱膨張

熱膨張は、フォノン振動数 $\omega_j(\mathbf{q})$ が結晶体積 $V$ に依存すること（＝非調和性の一次の効果）から生じます [2, 3]。この体積依存性の無次元指標がモードグリューナイゼン定数 $\gamma$ であり、体積 $V$ を微小に変えた 3 点のフォノン計算から中心差分で評価します。

$$\gamma_j(\mathbf{q}) = -\dfrac{\partial \ln \omega_j(\mathbf{q})}{\partial \ln V} \;\approx\; -\dfrac{\ln \omega_j(\mathbf{q};V_+) - \ln \omega_j(\mathbf{q};V_-)}{\ln V_+ - \ln V_-}$$

グリューナイゼン定数から、各モードの比熱 $c_i(T)$ を重みとしてマクロな熱物性へ接続します。

$$\bar\gamma(T) = \dfrac{\sum_i \gamma_i\, c_i(T)}{\sum_i c_i(T)}, \qquad c_i(T) = k_B \dfrac{x^2 e^{x}}{(e^{x}-1)^2},\ \ x=\dfrac{\hbar\omega_i}{k_B T}, \qquad \alpha(T) = \dfrac{\bar\gamma(T)\, C_V(T)}{3 B_0 V_0}$$

ここで $\alpha(T)$ は線膨張係数、$C_V(T)$ は定積比熱です。体積弾性率 $B_0$ と平衡体積 $V_0$ は、エネルギー–体積（E–V）曲線の状態方程式フィットから決定します。式が示すとおり、$\gamma_i$ が負のモードは熱膨張に負の寄与をします。

全体のワークフローを図1に示します。Advance/PHASE を「力・エネルギーの評価」エンジン、ALAMODE（alm, anphon） [4] を「力定数フィッティング＋フォノン物性」エンジンとして連携させ、状態方程式の決定からグリューナイゼン定数・熱膨張係数の算出までを自動化しています。なお Si は非極性結晶（Born 有効電荷がゼロ）であるため、LO-TO 非解析補正は不要です。

図1.　Si のモードグリューナイゼン定数・負熱膨張の算出ワークフロー: Advance/PHASE × ALAMODE。①状態方程式: 慣用セル（8 原子）で E–V 曲線 9 点を SCF 計算し、Birch–Murnaghan フィットで $a_0, B_0, V_0$ を決定。②フォノン計算: 平衡格子定数まわりの 3 体積で 2×2×2 スーパーセル（64 原子）を生成し、alm suggest の変位構造に対して Advance/PHASE で原子間力を評価、alm optimize で調和力定数を、anphon でフォノン $\omega_j(\mathbf{q})$ を体積ごとに求める。③解析: 3 体積の中心差分で $\gamma_j(\mathbf{q})$ を、モード比熱を重みに $\bar\gamma(T), \alpha(T)$ を導出。

2. 計算モデルと計算条件

状態方程式（E–V 曲線）の計算には Si ダイヤモンド構造の慣用立方セル（8 原子）を、フォノン計算には慣用セルを 2×2×2 に拡張した 64 原子スーパーセルを用いました。力の評価はいずれも構造緩和を伴わない 1 ショット SCF で行います。数値微分でグリューナイゼン定数を求めるため、3 体積で計算条件を完全に統一し、SCF 収束を通常のフォノン計算より厳しくとることが精度の要になります（振動数の「差」を扱うため）。主な計算条件を表1に示します。

表1.　計算条件の概要

項目	設定
対象	Si（ダイヤモンド構造, 空間群 \(Fd\bar{3}m\)）
状態方程式用セル	慣用立方セル 8 原子（格子スケール 9 点: 0.980–1.020）
フォノン用スーパーセル	2×2×2 慣用セル（64 原子）
フォノン計算体積	3 体積（格子スケール 0.995 / 1.000 / 1.005）
交換相関汎関数	LDA（Perdew-Wang）
カットオフ（波動関数 / 電荷密度）	25 / 100 Rydberg
k 点メッシュ（EOS 8 原子 / フォノン 64 原子）	8×8×8 / 2×2×2
変位量（有限変位法・調和 NORDER=1）	0.01 Å（3 体積で共通）
SCF 収束閾値	1×10−10 hartree（数値微分のため厳格化）
擬ポテンシャル	ノルム保存型
グリューナイゼン用 q メッシュ	16×16×16

交換相関汎関数には、Si の熱物性・フォノンで定量的に定評があり先行 LDA 計算とも比較しやすい LDA を採用しました。

3. 計算結果と考察

3.1 状態方程式と平衡格子定数

E–V 曲線 9 点の Birch–Murnaghan フィットから、平衡格子定数 $a_0 = 5.380$ Å、体積弾性率 $B_0 = 96.5$ GPa、その圧力微分 $B_0' = 4.24$ を得ました（図2）。実験値（$a_0 = 5.431$ Å [5], $B_0 = 97.9$ GPa [6]）に対して $a_0$ は $-0.94\%$、$B_0$ は $-1.43\%$ と、いずれも 1% 台の高い一致を示します。この $B_0, V_0$ が後段の $\alpha(T)$ の分母を規定します。

図2.　Si のエネルギー–体積（E–V）曲線と Birch–Murnaghan（3 次）フィット。赤点が Advance/PHASE（LDA）による 9 体積の全エネルギー、実線がフィット。平衡点で $a_0=5.3798$ Å、$B_0=96.5$ GPa、$B_0'=4.24$。緑破線は実験格子定数 5.431 Å。

3.2 フォノン分散の体積依存性

3 体積のフォノン分散を重ねて図3に示します。注目すべきは横音響（TA）枝の振る舞いです。通常のモードでは圧縮（体積減少）で振動数が上がりますが、Si の TA 枝（最低 2 分岐, ゾーン境界付近で 100–150 cm⁻¹）は逆に、圧縮した格子（スケール 0.995, 青破線）の振動数が膨張した格子（1.005, 赤点線）より低くなっています。すなわち $\partial\omega_{\rm TA}/\partial V > 0$ であり、これが $\gamma_{\rm TA} = -\partial\ln\omega/\partial\ln V \le 0$（負のグリューナイゼン定数）の微視的な起源です。一方、光学枝・LA 枝は圧縮で振動数が上がる通常の応答を示します。この TA 枝の異常な軟化こそが、Si の負熱膨張の出発点です。

図3.　3 体積（格子スケール 0.995 / 1.000 / 1.005）における Si のフォノン分散。TA 枝（最低 2 分岐, ゾーン境界付近で 100–150 cm⁻¹）では圧縮側（青破線）が膨張側（赤点線）より低振動数となり、体積依存が通常と逆転しています。光学枝は圧縮で高振動数化する通常応答。

3.3 モードグリューナイゼン定数の分散

中心差分で得たモードグリューナイゼン定数 $\gamma_j(\mathbf{q})$ の分散を図4に示します。TA 枝はブリルアンゾーンの広い範囲で負となり、X 点で $\gamma_{\rm TA}\approx-2.0$、L 点で $\approx-2.3$ に達します。対照的に、LA 枝・光学枝は全域で正の値をとります（光学枝はおよそ 0.8–1.55、LA 枝は L 点の 0.42 などやや小さめの正値）。高対称点での値を実験値 [7, 8] とともに表2にまとめます。$\Gamma$ 点 LTO の $\gamma=0.91$ は、圧力下ラマン散乱による実験値 0.98 ± 0.06 [7] と 7% 程度の差で整合します。

図4.　Si のモードグリューナイゼン定数 $\gamma_j(\mathbf{q})$ の分散（中心差分）。青の TA 枝が広範囲で負値（$\gamma\le 0$ 領域を淡色で強調）をとり、これが低温の負熱膨張をもたらします。$\Gamma$ 点近傍の音響枝は $\omega\to0$ で $\gamma$ が数値的に発散するため除外しています。

表2.　高対称点におけるフォノン振動数とモードグリューナイゼン定数（計算 vs 実験）

モード	\(\omega\) 計算 [cm−1]	\(\omega\) 実験 [cm−1]	\(\gamma\) 計算	\(\gamma\) 実験
LTO (Γ)	514.1	520	0.91	0.98 ± 0.06
TA (X)	140.1	150	−2.03	−1.40 ± 0.30
LA/LO (X)	409.9	410	0.82	—
TO (X)	461.3	463	1.48	—
TA (L)	107.7	114	−2.33	—
TO (L)	489.6	487	1.20	—

表2の実験値の出典: フォノン振動数は $\Gamma$ 点がラマン散乱 [7]、ゾーン境界（X, L 点）が中性子非弾性散乱 [8]。モードグリューナイゼン定数は圧力下ラマン測定 [7] による。

フォノン振動数自体はいずれも実験値の数 % 以内に収まっています。$\gamma_{\rm TA}(X)=-2.03$ は実験値（−1.4 [7]）より絶対値が大きめですが、これは LDA 準調和近似に共通する傾向です。実際、先行研究におけるLDA線形応答計算 [3] においても、$\gamma_{\rm TA}(X)$ は実験格子定数を用いた場合で $-1.78$、LDA平衡格子定数を用いた場合で $-2.30$ と、いずれも実験値より負の大きな値を示しています。本計算の結果である $-2.03$ は、この範囲内に収まっています。重要な点は、TA 枝の $\gamma$ が明確に負であり、その符号と分散の形が堅牢に再現されていることで、後述の負熱膨張はこの符号に支配されます。

3.4 熱膨張係数と負熱膨張

$\gamma_j(\mathbf{q})$ をモード比熱で重み付けして得た熱力学的グリューナイゼン定数 $\bar\gamma(T)$、比熱 $C_V(T)$、線膨張係数 $\alpha(T)$ を図5に示します。$\bar\gamma(T)$ は低温で負に落ち込み、約 150 K で符号を変えて 300 K で 0.37 となり、高温では 0.50 程度（1000 K）まで単調に増加します。比熱は高温で Dulong–Petit 極限（$3R$）に正しく収束します。線膨張係数 $\alpha(T)$ は、12–149 K の温度域で負となり（最小 $-0.92\times10^{-6}$ K⁻¹ at 82 K）、Si の負熱膨張を第一原理から再現しました。実験 [1] の負熱膨張域は 18–125 K、最小は $-0.48\times10^{-6}$ K⁻¹ at 76 K です。負熱膨張の符号反転と極小の位置はよく対応する一方、本計算は谷の深さを約 1.9 倍過大評価し、負の温度域も両端でやや広く出ています。これは $\gamma_{\rm TA}$ の絶対値の過大評価と整合する、LDA 準調和近似に共通の傾向です。300 K での $\alpha=2.14\times10^{-6}$ K⁻¹ は実験値 $2.57\times10^{-6}$（Okada–Tokumaru 1984 [5]、式を 300 K で評価）に対して約 17% の過小評価です。なお、近年の非弾性中性子散乱と第一原理計算の比較から、準調和近似は個々のフォノン振動数の温度シフトを正しく再現できない（多くのモードで符号が逆になる）にもかかわらず、モード間の寄与の大幅な相殺の結果として熱膨張自体はほぼ正しく予測することが指摘されています [9]。さらなる定量精度が必要な場合は、格子非調和性や核量子効果を陽に扱う枠組みが選択肢になります。

図5.　モードグリューナイゼン定数から導いた Si の熱物性。(a) 熱力学的グリューナイゼン定数 $\bar\gamma(T)$（低温で負, 星は実験の $\alpha, B_T, C_V$ から導かれる $\bar\gamma(300\,{\rm K})\approx0.45$）、(b) 定積比熱 $C_V(T)$（破線は Dulong–Petit $3R$）、(c) 線膨張係数 $\alpha(T)$（計算 vs Okada–Tokumaru 実験式 [5]・Lyon らの低温実測 [1]）、(d) 負熱膨張域の拡大（NTE 窓 12–149 K, 四角は実験値 [1]）。

4. まとめ

第一原理計算ソフトウェア Advance/PHASE をフォノン解析パッケージ ALAMODE と連携させ、3 体積・中心差分の準調和近似によってシリコンのモードグリューナイゼン定数 $\gamma_j(\mathbf{q})$ を算出しました。TA 枝が負のグリューナイゼン定数をもつことを確認し、そこから導いた線膨張係数 $\alpha(T)$ は 12–149 K の領域で負となるなど、Si の負熱膨張を第一原理計算から再現しました。調和近似では得られない熱膨張という非調和物性を、状態方程式からグリューナイゼン解析まで一貫した自動ワークフローで評価できることを示しました。第一原理フォノン計算は、熱膨張整合設計や低膨張材料の探索など、幅広い温度域にわたる材料設計を原子レベルで支える基盤となります。

本解析の詳細や、研究への適用可能性に関するご相談はこちら

お問い合わせ

参考文献

1. K. G. Lyon, G. L. Salinger, C. A. Swenson, and G. K. White, "Linear thermal expansion measurements on silicon from 6 to 340 K", J. Appl. Phys. 48, 865 (1977).
2. S. Biernacki and M. Scheffler, "Negative Thermal Expansion of Diamond and Zinc-Blende Semiconductors", Phys. Rev. Lett. 63, 290 (1989).
3. G.-M. Rignanese, J.-P. Michenaud, and X. Gonze, "Ab initio study of the volume dependence of dynamical and thermodynamical properties of silicon", Phys. Rev. B 53, 4488 (1996).
4. T. Tadano, Y. Gohda, and S. Tsuneyuki, "Anharmonic force constants extracted from first-principles molecular dynamics: applications to heat transfer simulations", J. Phys.: Condens. Matter 26, 225402 (2014).
5. Y. Okada and Y. Tokumaru, "Precise determination of lattice parameter and thermal expansion coefficient of silicon between 300 and 1500 K", J. Appl. Phys. 56, 314 (1984).
6. H. J. McSkimin, "Measurement of Elastic Constants at Low Temperatures by Means of Ultrasonic Waves — Data for Silicon and Germanium Single Crystals, and for Fused Silica", J. Appl. Phys. 24, 988 (1953).
7. B. A. Weinstein and G. J. Piermarini, "Raman scattering and phonon dispersion in Si and GaP at very high pressure", Phys. Rev. B 12, 1172 (1975).
8. G. Nilsson and G. Nelin, "Study of the Homology between Silicon and Germanium by Thermal-Neutron Spectrometry", Phys. Rev. B 6, 3777 (1972).
9. D. S. Kim, O. Hellman, J. Herriman, H. L. Smith, J. Y. Y. Lin, N. Shulumba, J. L. Niedziela, C. W. Li, D. L. Abernathy, and B. Fultz, "Nuclear quantum effect with pure anharmonicity and the anomalous thermal expansion of silicon", Proc. Natl. Acad. Sci. USA 115, 1992 (2018).

関連ページ

• 第一原理計算ソフトウェア Advance/PHASE
• 解析分野：ナノ・バイオ
• 産業分野：材料・化学