理想固体の状態方程式 ― Grüneisen 定数と Mie–Grüneisen の状態方程式

固体を温めると膨らみます．あまりに当たり前ですが， ばねが理想的なばねである限り，固体は絶対に膨らみません．膨らむためには，原子の間の力が「伸びと縮みで対称でない」必要があります ―― 非調和性です．ここではその非調和性を Grüneisen 定数 という 1 つの数に押し込め，熱力学の関係式とフォノンの自由エネルギーから Mie–Grüneisen の状態方程式を導きます．そのうえで擬調和近似（QHA）を置くと，体積熱膨張係数が α_V = γC_V/(B_TV) という驚くほど簡単な形になります．最後に，この式のどの因子が実際に熱膨張を支配しているのか，そして温めると縮む物質（負熱膨張）がなぜ存在するのかを見ます．

1. 「膨らむ」を式にする

体積熱膨張係数は，圧力を一定に保ったまま温度を上げたときの体積の増え方です．

α_{V} \equiv \frac{1}{V} {(\frac{\partial V}{\partial T})}_{P}

(1)

ここで V は平均原子体積，つまり格子の体積を格子内の原子数で割ったものとします ―― こうしておくと，組成の違う物質どうしを同じ土俵で比べられます．

この定義のままでは扱いにくいので，三重積の関係 （ ${(\partial V / \partial T)}_{P} {(\partial T / \partial P)}_{V} {(\partial P / \partial V)}_{T} = - 1$ ）を使って書き直します．

α_{V} = - \frac{1}{V} {(\frac{\partial V}{\partial P})}_{T} {(\frac{\partial P}{\partial T})}_{V} = \frac{1}{B_{T}} {(\frac{\partial P}{\partial T})}_{V}

(2)

前半の $- \frac{1}{V} {(\partial V / \partial P)}_{T}$ は圧縮率 κ で，体積弾性率 B_T の逆数です．つまり（1）を知るには，体積を固定したまま温度を上げたとき，結晶の内圧がどれだけ上がるか ―― それだけを計算すればよい，ということになります．

なぜ内圧なのか温めると原子は激しく振動し，まわりを押し広げようとします．これが熱圧力です．実際に膨らむ量は「押す力（熱圧力の増え方）」÷「押し返す力（体積弾性率）」で決まります．硬い物質ほど膨らまないのは，この分母が大きいからです ―― 話の結末はもうここに書いてあります．

2. 理想的なばねは膨らまない

原子 1 対を取り出して，つり合いの位置からのずれを x とします．ポテンシャルが完全な放物線 $U (x) = \frac{1}{2} k x^{2}$ なら，Boltzmann 分布 $e^{- U / k_{B} T}$ は x について左右対称です．どんなに温度を上げても

⟨ x ⟩ = 0

―― 原子は激しく振動するだけで，平均の位置は 1 ミリも動きません．調和振動子をいくら並べても熱膨張は出てこないのです．

そこで 3 次の項を足します．伸びる側（x > 0）のほうが登りやすい ―― 逆に縮む側は急に立ち上がる，という非対称性です．ただし 3 次の項だけを足すと，x を大きくしたときポテンシャルが 下へ落ちていって底がなくなります（U → −∞）．それでは Boltzmann 因子 $e^{- U / k_{B} T}$ の積分（分配関数）が発散してしまい，平均そのものが定義できません．そこで底を作る 4 次の項を足しておきます．

U (x) = \frac{1}{2} k x^{2} - g x^{3} + λ x^{4}, ⟨ x ⟩ ≃ \frac{3 g k_{B} T}{k^{2}}

(3)

4 次の項は平均をずらしません λx⁴ は x について偶関数なので，最低次の摂動では

⟨ x ⟩

を動かしません．平均をずらすのは奇関数の 3 次の項だけで，そこから出るのが（3）の右の式です． 4 次の項は「底を作って，分配関数を収束させる」ためだけに入れてあります．図 1 では λ = g²/k と取っていて，このとき U′ = x(k − 3gx + 4λx²) の括弧が つねに正（判別式 9g² − 16λk = −7g² < 0）なので，谷は x = 0 の 1 つだけです．

平均の位置が温度に比例してずれます．これが熱膨張の最も素朴な姿です．図 1 で確かめてください ―― g = 0 にすると，温度をいくら上げても平均位置は動きません．ただし（3）は Boltzmann 分布で平均を取った古典近似なので，高温側でしか使えません ―― 量子的には k_BT のところに零点振動を含む振動子のエネルギーが入り， T → 0 では平均位置のずれが温度によらなくなります（7 節の α_V → 0 に対応します）．

図 1．非調和なポテンシャルの中では，平均の位置がずれます． 薄い色は Boltzmann 分布 $e^{- U / k_{B} T}$ ，縦の線がその平均 $⟨ x ⟩$ です． k = 15 eV/Å² に固定し，安定化の 4 次項は λ = g²/k としてあります（換算質量を 30 u とすれば，このばね定数は約 11 THz にあたります ―― 振動数はばね定数だけでは決まらず，質量が要ります）．ずれがとても小さいことにも注目してください ―― 熱膨張とは，もともとこの大きさの現象です．

温度 T = 600 K 非調和の強さ g = 10 eV/Å³

結晶では，非調和性は「振動数の体積依存性」として現れます 1 対の原子なら 3 次の項でよいのですが，結晶では原子が N 個つながっています．そこでふつうは見方を変えて，結晶の体積を変えるとフォノンの振動数がどう変わるかで非調和性を測ります．理想的なばねなら，体積を変えても振動数は変わりません ―― ばね定数はばねの持ち主の都合で決まっているからです．実際の結晶では，広げると結合がゆるんで振動数が下がります．

3. Grüneisen 定数

波数 q・バンド ν のフォノンについて，モード Grüneisen 定数を次で定義します．

γ_{q ν} \equiv \frac{- V}{ω_{q ν}} {(\frac{\partial ω_{q ν}}{\partial V})}_{T} = - \frac{\partial \ln ω_{q ν}}{\partial \ln V}

(4)

右の形が読みやすいでしょう ―― 体積を 1 % 増やしたとき，振動数が何 % 下がるかです．マイナス符号は「広げると下がる」を正にするためのものです．

理想的なばね（調和固体）なら振動数は体積によらないので γ = 0．熱膨張も 0．
ふつうの固体は γ ≒ 1 〜 3．E. Grüneisen が 1912 年に導入しました．
γ が負のモードもあります．広げると振動数が上がるモードで，これが負熱膨張の正体です（9 節）．

4. 結晶の Helmholtz 自由エネルギー

体積 V・温度 T の結晶の Helmholtz 自由エネルギーは， 0 K の静的なエネルギー E(V) と， フォノンのぶんの和で書けます．

量の規格をそろえておきます以下の E, F, U, C_V とモード和

\sum q, ν

は， すべて同じ規格で書きます ―― 結晶全体なら全部を全体で， 1 原子あたりなら全部を 1 原子あたりで（このとき和には

1 / N_{atom}

が付きます）．最後に出てくる α_V は示強量なので， C_V と V を同じ規格でそろえてさえいれば，どちらで書いても式（11）はそのままです． 1 節で V を平均原子体積と決めたので，このページではすべて 1 原子あたりで読んでください．

F (V, T) = E (V) + \sum_{q, ν} [\frac{ℏ ω_{q ν}}{2} + k_{B} T \ln (1 - e^{- ℏ ω_{q ν} / k_{B} T})]

(5)

第 1 項が零点振動，第 2 項が熱で励起されたぶんです．この式そのものは調和振動子の分配関数から出てきます ―― 調和振動子の導出と比熱の導出に書きました． 1 つのモードが持つ内部エネルギーは

U_{q ν} = (\frac{1}{e^{ℏ ω_{q ν} / k_{B} T} - 1} + \frac{1}{2}) ℏ ω_{q ν}

(6)

丸括弧の中は「そのモードに乗っているフォノンの数 + 1/2」です． 絶対零度でも 1/2 が残ります ―― 零点振動は消せません．

5. Mie–Grüneisen の状態方程式

圧力は自由エネルギーを体積で微分すれば出ます． $P = - {(\partial F / \partial V)}_{T}$ です．（5）の第 2 項は V をあらわに含まず，ω_qν を通してだけ体積に依るので，連鎖律で

\frac{\partial F_{ph}}{\partial V} = \sum_{q, ν} \frac{\partial F_{ph}}{\partial ω_{q ν}} \frac{\partial ω_{q ν}}{\partial V}, \frac{\partial F_{ph}}{\partial ω_{q ν}} = \frac{ℏ}{2} + \frac{ℏ}{e^{ℏ ω / k_{B} T} - 1} = \frac{U_{q ν}}{ω_{q ν}}

ここが気持ちのよいところです ―― 微分してみると，ちょうど（6）の U_qν/ω_qν が現れます．あとは（4）から $\partial ω_{q ν} / \partial V = - γ_{q ν} ω_{q ν} / V$ を入れるだけです．

P + \frac{d E}{d V} = \frac{1}{V} \sum_{q, ν} γ_{q ν} U_{q ν}

(7)

これが Mie–Grüneisen の状態方程式です（G. Mie, 1903／E. Grüneisen, 1912）．理想気体の PV = Nk_BT にあたるものが，固体ではこれです．

（7）の読み方左辺の

d E / d V

は結合の弾性だけから来る項，つまり 0 K の圧力です（E(V) は静的なエネルギーなので，温度を含みません）．右辺はフォノンが押し広げる圧力です． U_qν には零点の ħω/2 が入っているので，これは 0 K でも消えません ―― だから「熱圧力」ではなくフォノン圧力と呼ぶのが正確です（温度で上がったぶんだけを指したいなら P_vib(T) − P_vib(0) を取ります）． 1 つ 1 つのモードが「自分のエネルギー U_qν × 自分の γ_qν」だけ寄与します． γ_qν が正のモードは外へ押し，負のモードは内へ引きます．シミュレーターの①で見ているのは，この綱引きの結果です．

外圧をかけずに（P = 0）置いておくと，結晶は（7）の右辺と左辺がつり合う体積に落ち着きます． 0 K の弾性エネルギーを $E (V) = \frac{1}{2} B_{T} {(V - V_{0})}^{2} / V_{0}$ と放物線で近似すれば，膨らむ量そのものが書けます．

V - V_{0} = \frac{V_{0}}{B_{T} V} \sum_{q, ν} γ_{q ν} U_{q ν}

(8)

U_qν は必ず正なので， 膨らむか縮むかは γ_qν の符号だけで決まります ―― ただし決めるのは和のほうです．負の γ_qν を持つモードが 1 本あるだけでは足りず， $\sum γ_{q ν} U_{q ν} < 0$ が要ります．

6. 擬調和近似（QHA）とは何か

（7）はきれいですが，そのままでは使えません． ω_qν も γ_qν も，本当は温度と体積の両方に依るからです．そこで置く仮定が擬調和近似（quasi-harmonic approximation, QHA）です．

擬調和近似の中身は 2 行です

体積を止めておけば，フォノンの振動数は温度によらない． つまり各体積では，結晶はふつうの調和振動子の集まりとして扱う．
振動数の変化は，すべて体積変化に由来する． 温度が効くのは，フォノンの占有数を通してだけ．

「調和」ではなく「擬調和」と呼ぶのは，各体積では調和なのに，体積が変わると振動数が変わる（＝非調和性が残っている）からです．

実際の第一原理計算は，この仮定のもとで次の手順になります．

格子定数を少しずつ変えた何通りかの体積で，調和フォノンを計算する（ω_qν(V) が得られます）．
各体積・各温度について，（5）から F(V, T) を組み立てる．
温度ごとに F(V, T) を最小にする体積を探す．それが平衡体積 V(T)．

虚数のフォノンがあってはいけません振動数が虚数になる体積では，調和フォノンの自由エネルギーが実数として定義できません．遠くまで圧縮・膨張した試行体積の 1 点だけが不安定なら，そこを外して残りで組み立てられますが， 平衡体積を探すのに必要な体積範囲で虚数モードが出てしまうと，標準的な QHA はそのままでは使えません．論文で SiO₂ の 4 つの多形が解析から外れているのは，この理由です．

図 2 で，その「谷が動く」ようすを見てください． 0 K の谷（V₀）に対して，温度を上げると谷底がずれていきます． γ を負にすると，谷底は左へ動きます ―― 温めたのに縮む，これが負熱膨張です．

図 2．擬調和近似では，自由エネルギーの谷底が温度とともに動きます． Einstein 模型（1 原子あたり 3 モード，Θ_E = 500 K）で $ω (V) = ω_{0} {(V / V_{0})}^{- γ}$ とした，（5）そのものです．V₀ = 12 Å³/原子． 2 本の曲線は，それぞれ自分の谷底が 0 になるようにそろえて描いてあります ―― 見たいのは谷の深さではなく，谷底の位置が横にずれることだからです．谷底を数値で探して出した α_V と，式（11）から出した α_V が一致することも表示しています．なお「0 K の谷」は，零点振動まで含めた F(V, 0) の谷底です ―― 静的な E(V) の谷 V₀ = 12 Å³ とは， γ ≠ 0 なら零点圧力のぶんだけずれます．

温度 T = 600 K Grüneisen 定数 γ = 1.5 体積弾性率 B_T = 150 GPa

図 2 では，温めると B_T が少し増えます読み取り値を見ると，150 GPa と置いたのに 600 K では 157 GPa になっています．まず，つまみの 150 GPa は静的な放物線 E(V) の曲率パラメーターで，読み取りの B_T = V∂²F/∂V² は フォノン自由エネルギーの曲率まで含んだ量です ―― まず別の量です．そのうえで増える主な理由は，放物線では E″ が一定でも B = V E″ なので体積とともに増えてしまうことにあります．放物線には「押し縮めると硬くなる」という性質（

\partial B_{T} / \partial P \approx 4

）が入っていないので，膨らんでも柔らかくなりません．実際の固体は温めると柔らかくなります．シミュレーターのほうは， 595 K までは QHA 計算が出した B_T(T) をそのまま表示し，そこから上だけ Anderson–Grüneisen の関係

B_{T} = B_{0} {(V / V_{0})}^{- δ_{T}}

（δ_T = 5）で外挿しています．そちらでは温度とともに B_T が下がります ―― MgO なら 300 K の 150 GPa が 1000 K で 126 GPa まで落ちます（B_T(1000)/B_T(300) = 0.840．実験値は 0.84 です）．

7. 体積熱膨張係数が簡単な形になる

（2）に必要なのは ${(\partial P / \partial T)}_{V}$ でした．（7）を体積一定のまま温度で微分します． E(V) は温度を含まないので消えて，

\begin{aligned} {(\frac{\partial P}{\partial T})}_{V} & = \frac{1}{V} \sum_{q, ν} [{(\frac{\partial U_{q ν}}{\partial T})}_{V} γ_{q ν} + U_{q ν} {(\frac{\partial γ_{q ν}}{\partial T})}_{V}] \\ ≃ \sum_{q, ν} \frac{C_{V, q ν} γ_{q ν}}{V} = \frac{γ C_{V}}{V} \end{aligned}

(9)

2 行目で擬調和近似を使いました ―― 体積一定なら振動数は温度によらないので， γ_qν も温度によらず，第 2 項が落ちます．残った第 1 項の ${(\partial U_{q ν} / \partial T)}_{V}$ は，そのモードの定積比熱 C_V,qν にほかなりません．最後に，全体の比熱と，比熱で重みをつけた平均の Grüneisen 定数を

C_{V} = \sum_{q, ν} C_{V, q ν}, γ = \frac{\sum_{q, ν} C_{V, q ν} γ_{q ν}}{\sum_{q, ν} C_{V, q ν}}

(10)

と定義すれば，（2）と（9）から一行で出ます．

α_{V} = \frac{γ C_{V}}{B_{T} V}

(11)

固体の熱膨張は，たった 4 つの量で書けてしまいます．しかも C_V・B_T・V は必ず正なので， α_V の符号を決めているのは γ だけです．

絶対零度では膨らみません T → 0 で C_V → 0 なので（比熱の導出），（11）から α_V → 0 です．これは熱力学第三法則が要求することとぴったり合っています ―― 詳しくは調和振動子の導出の 5 節に書きました．逆に高温では C_V → 3k_B（Dulong–Petit）で頭打ちになります．ただし α_V が一定になるとまでは言えません ―― （11）には γ, B_T, V も入っていて，高温では B_T が下がり V が増え，固有の非調和性で γ も動くからです． これらの温度変化を無視できる範囲でだけ， α_V はほぼ一定に落ち着きます．

8. 熱膨張を支配しているのは何か

（11）の右辺には 4 つの量が並んでいます．どれがいちばん効いているのでしょうか．実際の物質でどれくらいばらつくかを並べると，答えははっきりしています．

因子	室温での値の幅	何倍ちがうか	役割
C_V	0.7 〜 3 k_B／原子	約 4 倍	高温では Dulong–Petit で頭打ち．ダイヤモンドのように Θ_D の高い物質は室温ではまだ頭打ちにならない
γ	−0.06 〜 3.3	符号が変わる正の 37 物質では 0.38 〜 3.34，約 9 倍	符号を決める．結晶構造と配位数で決まる
V	5.7 〜 71 Å³	12 倍	B_T と逆相関しているので単独では効きにくい
B_T	2.9 〜 438 GPa	150 倍	これが支配的

幅の値は，論文の Table S2／S3 にある立方晶 38 物質から取りました（シミュレーターの②で同じ表を見られます）． B_T だけが 2 桁以上ちがいます． だから α_V と B_T を散布図にすると， α_V ∝ 1/B_T の反比例曲線がきれいに現れます．正確に言えば，この 38 種類の立方晶で，正の熱膨張の「大きさ」のばらつきにいちばん強く効くのが B_T，ということです ―― 符号を決めるのはあくまで γ ですし，強誘電体・磁気転移・電子状態の転移のように，格子振動以外が熱膨張を支配する物質もあります．

ダイヤモンド（B_T = 438 GPa）の α_V は 3.3 ppm/K，カリウム（B_T = 2.9 GPa）は 245 ppm/K ―― 75 倍です． 硬いものは膨らまない．柔らかいものはよく膨らむ． 日常の感覚とも合いますし，式（11）の分母がそう言っています．

平均原子体積 V は，見かけほど効きません（11）だけ見ると V が大きいほど α_V は小さくなりそうです．ところが実際の物質では，V が大きい結晶ほど B_T が小さい という強い逆相関があります（原子どうしが遠いほど結合がゆるい）．分母の B_TV のうち B_T の変化のほうが激しいので，結果として α_V は V と弱い正の相関を示します．シミュレーターの②で横軸を「平均原子体積 V」に切り替えると確かめられます．

9. 温めると縮む物質 ― 負熱膨張

（11）で C_V, B_T, V はすべて正でした．したがって α_V < 0 になるには γ < 0 しかありません．これが負熱膨張（negative thermal expansion, NTE）の必要条件です．

γ_qν < 0 とは，（4）から 体積を広げると振動数が上がるモードのことです．直感に反しますが，起こります．そして負のモードが何本かあるだけでは足りません ―― 比熱で重みをつけた和 $\sum C_{V, q ν} γ_{q ν} < 0$ になって初めて，結晶全体が縮みます．重みは温度で変わるので，同じ物質でも負熱膨張になる温度範囲が限られるのはこのためです（Si が低温だけ縮むのがその例です）．

張力効果（tension effect） A–O–A のように 2 つの原子に挟まれた酸素が，結合を横切る向きに振れる場面を考えます．酸素が横に逃げると，A–O–A の距離は縮みます（弦をつまむと両端が引き寄せられるのと同じです）．ここで結晶を引き伸ばすと，弦が張られて横振動しにくくなる ―― 振動数が上がる．これがまさに γ_qν < 0 です．そして温度が上がるほど横振動が激しくなり，結晶は縮みます．

論文が 46 種の単体元素と 45 種の二元系酸化物を系統的に調べて示したのは，次の 5 点です．

大きな V をもつ結晶は，小さな B_T をもつ傾向がある．
γ の主な役割は α_V の符号を決めることであり，その値は元素そのものよりも結晶構造と配位数によって決まる．
横波音響（TA）フォノンや低振動数の横波フォノンが，負熱膨張に大きく寄与する．
負熱膨張材料は比較的大きな V をもち，結晶構造内に空隙がある．
負熱膨張材料に共通するのは，空隙に向かって振動するフォノンが存在することである．

Y. Mochizuki, H. Koiso, K. Nagamatsu, S. Bae, T. Isobe, A. Nakajima, J. Phys. Chem. C 128, 525–535 (2024) の結論（要約）．

つまり負熱膨張の設計指針は，「空隙があり，配位数が小さく，柔らかい」ということになります．空隙は振れる先を与え，小さい配位数は結合角を曲げやすくし，柔らかさ（小さい B_T）は（11）の分母を小さくして効果を増幅します．

物質	構造	γ	B_T (GPa)	α_V (ppm/K)	ひとこと
Cu₂O	赤銅鉱型 Pn3̄m	−0.061	117	−1.40	Cu の配位数が 2．O–Cu–O が横に振れる
ReO₃	ReO₃ 型 Pm3̄m	0.384	244	3.47	A サイトが空いたペロブスカイト．八面体が回る
Si	ダイヤモンド型	0.39	90.0	7.22	低温では α_V が負になる
MgO	岩塩型	1.63	150	35.9	密に詰まっていて空隙がない
K	体心立方	1.30	2.89	245	とにかく柔らかい

同じ組成でも構造が変われば γ が変わる，という点は特に大事です．論文では CaO の岩塩型と閃亜鉛鉱型を比べていて，岩塩型は正の γ_qν しか持たないのに，閃亜鉛鉱型は X 点と L 点に負の γ_qν を持ちます． γ は元素の性質ではなく，並べ方の性質なのです．

10. 擬調和近似はどこで破れるか

（9）の 2 行目で落とした ${(\partial γ_{q ν} / \partial T)}_{V}$ を，落とさずに残してみます． 熱 Grüneisen 定数 δ_qν と非調和定数 Δ_qν を

δ_{q ν} \equiv \frac{- T}{ω_{q ν}} {(\frac{\partial ω_{q ν}}{\partial T})}_{V}, Δ_{q ν} \equiv \frac{- V T}{ω_{q ν}} \frac{\partial^{2} ω_{q ν}}{\partial T \partial V}

(12)

と定義します．この δ_qν は，シミュレーターで 595 K より上の B_T の外挿に使った Anderson–Grüneisen パラメーター δ_T とは別物です ―― 同じ文字ですが，定義も意味も違います．

この 2 つを使うと，γ_qν の温度微分がきれいにまとまります．

{(\frac{\partial γ_{q ν}}{\partial T})}_{V} = \frac{1}{T} (γ_{q ν} δ_{q ν} + Δ_{q ν})

(13)

（13）の出どころ

γ_{q ν} = - (V / ω) ω_{V}

を温度で微分するだけです．

T (\partial γ / \partial T) = - V T ω_{V T} / ω + V T ω_{V} ω_{T} / ω^{2}

で，第 1 項がそのまま Δ_qν，第 2 項が

[- (V / ω) ω_{V}] [- (T / ω) ω_{T}] = γ_{q ν} δ_{q ν}

です．δ の前に γ が付くのがこの式の要点で，次の平均の定義にもそのまま効いてきます．

これを（9）の 1 行目に戻すと，体積熱膨張係数は次のように拡張されます．ここでモード和を U_qν で重みづけした平均を

δ \equiv \frac{\sum_{q, ν} U_{q ν} γ_{q ν} δ_{q ν}}{\sum_{q, ν} U_{q ν}}, Δ \equiv \frac{\sum_{q, ν} U_{q ν} Δ_{q ν}}{\sum_{q, ν} U_{q ν}}

(14)

と書きます ―― δ の中には γ_qν が入っています． $U = \sum U_{q ν}$ と置けば，

\begin{aligned} α_{V} & = \frac{1}{B_{T} V} \sum_{q, ν} [C_{V, q ν} γ_{q ν} + \frac{U_{q ν}}{T} (γ_{q ν} δ_{q ν} + Δ_{q ν})] \\ = \frac{γ C_{V}}{B_{T} V} + \frac{1}{B_{T} V} \frac{U}{T} (δ + Δ) \end{aligned}

(15)

第 1 項が QHA，第 2 項がそこからのはみ出しです．

高温極限の読み方に注意 1 本のフォノンモードについて，高温（ħω ≪ k_BT）では

平均占有数 $⟨ n ⟩ ≃ k_{B} T / ℏ ω$ ―― これは個数で，無次元です．
モードのエネルギー $U_{q ν} ≃ k_{B} T$
モードの比熱 $C_{V, q ν} ≃ k_{B}$

です．3k_B は 1 原子あたり 3 モードを足した値で， 1 モードの比熱ではありません．したがって（15）の第 2 項の重みは高温でどのモードも

U_{q ν} / T ≃ k_{B}

で同じになり，効くのは

γ_{q ν} δ_{q ν} + Δ_{q ν}

の大きさそのものです．

低い振動数のモードほど，温度を上げたときの $\partial ω / \partial T$ が相対的に大きくなります（δ_qν は ω で割った量です）．だから低振動数のモードをたくさん持つ柔らかい物質ほど，QHA から外れます．

この拡張式そのものの限界（15）は，ω に温度依存を持たせたうえで形式的に微分したものです．厳密にやるには，C_V,qν を「振動数を固定した調和比熱」と見るのか，「温度依存する振動数まで含めた全微分」と見るのかを決めなければならず，さらに調和自由エネルギーに ω(V,T) を代入しただけでは，エントロピーや自由エネルギーとの熱力学的な整合が保証されません．きちんと扱うには自己無撞着フォノン理論のような，温度依存フォノンと熱力学ポテンシャルが整合した枠組みが要ります．ここでは「QHA には決まった形の外れ方がある」という道しるべとして読んでください．

QHA が当てにならない物質岩塩型 NaBr のように非常に柔らかい物質では，QHA だけでは実測の α_V を再現できません． Pb や Bi のような重くて柔らかい元素を含む化合物，臭化物・ヨウ化物・セレン化物・テルル化物も同じ理由で怪しくなります．一方，この導出で扱っている単体元素と二元系酸化物の多くは NaBr ほど柔らかくないので，QHA でよく合います．ただし K（B_T = 2.9 GPa）をはじめ，Li・Na・Ba・Sr・Ca といったアルカリ金属・アルカリ土類金属は NaBr と同程度かそれ以上に柔らかく，すぐ上に挙げた Pb ともども，同じ注意が要ります．

11. まとめ

理想的なばね（調和固体）は，どんなに温めても膨らみません．熱膨張は非調和性そのものです．
結晶では非調和性を Grüneisen 定数 $γ_{q ν} = - \partial \ln ω / \partial \ln V$ にまとめます．
自由エネルギーを体積で微分すると Mie–Grüneisen の状態方程式（7）が出ます．右辺はフォノンが押し広げる圧力です ―― 零点振動のぶんが入るので，0 K でも消えません（5 節）．
擬調和近似（体積を止めれば振動数は温度によらない）を置くと， α_V = γC_V/(B_TV)（11）という 1 行になります．
4 つの因子のうち，物質によって桁で変わるのは B_T だけです． 熱膨張を支配しているのは硬さです．
C_V, B_T, V は正なので，温めて縮むには γ < 0 しかありません．空隙・小さい配位数・横波フォノンが，その条件を作ります．

①で揺れている大きさは，熱膨張とは別の量ですシミュレーターの①で原子が揺れる振幅は， Debye 模型から見積もった原子変位パラメーター（ADP） U_iso（1 方向あたりの平均二乗変位）の平方根です． U_iso には零点振動のぶんが残るので，0 K でも揺れは止まりません ―― このページで追ってきた「膨らむ」（平衡体積 V(T) が動くこと）とは別の量です．式そのものは比熱の導出の 8 節に書きました．なお①の箱の大きさ（格子定数）は QHA 計算の V(T) から作った計算値で，揺れ方だけが Debye 模型による模式表示です．ただし画面上の伸びは既定で 25 倍に誇張してあります（倍率は画面に「膨張 ×25」と出ていて，スライダーで変えられます）―― 実際の熱膨張は，2 節で見たとおりもっとずっと小さい量です．

実際の数値は理想固体の状態方程式シミュレーターで動かせます ―― ①で結晶を回しながら温度を上げ，②で 38 物質の α_V–B_T 相関を眺め，③で温度依存を追えます． ①〜③に出てくる V(T)・α_V(T)・ B_T(T)・C_V(T) は， 300 K の断面だけでなく 0 〜 595 K の QHA 計算の出力そのものです．データは Y. Mochizuki, H. Koiso, K. Nagamatsu, S. Bae, T. Isobe, A. Nakajima, “Thermal Properties of the Element and Binary Oxides toward Negative Thermal Expansion: A First-Principles Lattice-Dynamics Study,” J. Phys. Chem. C 128, 525–535 (2024)（オープンアクセス，CC BY 4.0）の Table S2／S3 と，その計算に使った phonopy-qha の出力によります．関連ページ：比熱の導出，調和振動子と第 2 量子化，エントロピーと自由エネルギー．