音響ビリアル係数から密度ビリアル係数を求める公式 [1]○多田 康平 A) B)
A)

京都大学 大学院工学研究科、 B)京都大学 環境安全保健機構 低温物質管理部門ているが[5-7]、その数式は、一見使いやすいものとは

1. 序論
液化ガスの製造など高圧ガスを取り扱う場合には、用いる気体に対して、（理想気体としてではなく、）いえない。加えて、その他の分子については、そのような方法論自体、報告されていない。

実在気体としての性質がしばしば利用される。例え

そこで本研究では、γ0 が温度によらず定数である

ば、ヘリウム液化の最終段階で利用されるジュール・と仮定した条件の下で、式(3)を数学的に解き、βa にトムソン効果は、定義により理想気体では現れず、ついて適当な関数形を仮定し、さらに物理学的に適[2]実在気体でこそ現れる効果である 。このような背当な初期条件を課すことで、単原子分子のみならず景から、用いる気体の実在気体としてのふるまいをその他の分子にも適用できる、βa の温度依存性からよく理解しておくことは重要であると考えられる。

B の温度依存性を計算する公式を導く。最近の筆者

希薄な実在気体の記述には、物理量のビリアル展の論文 [1] では非直線多原子分子の例としてメタン開がしばしば用いられる。主なものに、圧縮因子 Z（CH4）についても議論しているが、本稿では、極低を気体のモル密度 ρ のべきで展開したもの：
𝑃𝑍== 1 + 𝐵𝜌 + 𝐶𝜌2 + ⋯𝜌𝑅𝑇

温寒剤として広く利用されている単原子分子のヘリウム 4（4He）と二原子分子の窒素（N2）に限定して、(1)

音速 u の二乗を圧力 P のべきで展開したもの：
𝛾0 𝑅𝑇
𝛽a
𝛾a 2
𝑢2 =
(1 +𝑃+𝑃 + ⋯)𝑀𝑅𝑇𝑅𝑇(2)がある。ここに、R は気体定数、T は熱力学温度、Mはモル質量、γ0 は理想気体における定圧熱容量と定積熱容量の比（比熱比）を表す。B, C, …, は密度ビリアル係数、βa, γa, …, は音響ビリアル係数と呼ばれ、これらはいずれも温度の関数である。比熱比 γ0 は、単原子分子では温度によらず一定の値 5/3 をとるが、それ以外の分子では、主に分子内振動や他の分子内自由度のために、温度の関数である。分子回転のみが励起され、他の自由度を無視してよいとすると、γ0 は、二原子分子および直線多原子分子では 7/5、非直線多原子分子では 4/3 の定数値をとる。
密度と音響のビリアル係数を結ぶ関係式は広く知られており[3, 4]、第二ビリアル係数については、以下の式で与えられる。
𝛽a = 2𝐵 + 2(𝛾0 − 1)𝑇d𝐵 (𝛾0 − 1)2 2 d2 𝐵+𝑇d𝑇
𝛾0d𝑇 2(3)B の温度依存性が分かれば βa の温度依存性は計算で得られた公式により βa の温度依存性から B の温度依存性をよく計算できることを確認する。
2. 公式の導出γ0 が定数であるとすれば、式(3)はオイラー型の微分方程式である。変数 T を t = lnT で置換すると、
d2 𝐵 𝛾0 + 1 d𝐵
2𝛾0𝛾0++
𝐵=
𝛽22(𝛾0 − 1)2 ad𝑡𝛾0 − 1 d𝑡 (𝛾0 − 1)(4)が得られ、この微分方程式の一般解として𝐵 = e𝑝𝑡 [𝐶1 sin(𝑞𝑡) + 𝐶2 cos(𝑞𝑡)]
𝛾0 e𝑝𝑡 sin(𝑞𝑡) 𝛽a cos(𝑞𝑡)+∫d𝑡𝑞(𝛾0 − 1)2e𝑝𝑡
𝛾0 e𝑝𝑡 cos(𝑞𝑡) 𝛽a sin(𝑞𝑡)−∫d𝑡𝑞(𝛾0 − 1)2e𝑝𝑡(5)が得られる。式(4)は 2 階の微分方程式なので、積分定数が 2 個（C1 および C2）存在する。これらはあとで初期条件を課して決定する。p と q はそれぞれ𝑝=−𝛾0 + 1,2(𝛾0 − 1)𝑞=√−𝛾02 + 6𝛾0 − 12(𝛾0 − 1)(6)で与えられ、γ0 が定数のとき p と q も定数である。
式(4)の一般解である式(5)は、βa の具体的な関数形きる一方で、βa の温度依存性が分かっていてもこのが与えられておらず、そのままでは利用することが式だけでは B の温度依存性は自明でない。単原子分難しい。そこで βa に適当な関数形を仮定することに子については式(3)の微分方程式が解かれ、βa の温度する。ビリアル係数の温度依存性は、例えば論文[8, 9]依存性から B の温度依存性を計算する数式が示されのように、しばしば熱力学温度のべき級数で表される。これを踏まえ、βa の温度依存性が
𝛽a = ∑ 𝛽𝑗 𝑇𝑗 = ∑ 𝛽𝑗 e𝑗𝑡𝑗𝑗とが極めて難しいからである。
いくつかの検討の結果、物理学的にも適当な初期条件として、「T = TB で B = 0」および「TB とは別の(7)で表される場合を考える。βj は、温度に依存しない、温度 T = TB/m で B = Bm」を選ぶことにした。m は 1T j の係数とする。この場合、式(5)は簡単になり、
𝐵 = e𝑝𝑡 [𝐶1 sin(𝑞𝑡) + 𝐶2 cos(𝑞𝑡)]
𝛾0 𝛽𝑗 e𝑗𝑡+∑
(𝛾0 − 1)2 𝑗 2 + (𝛾02 − 1)𝑗 + 2𝛾0

期条件から、C1 と C2 は以下のように与えられる。
𝑚 𝑝
𝑇B
𝐶1 = − {[𝐵𝑚 ( ) − 𝛺 ( )] cos(𝑞 ln 𝑇B )𝑇B
𝑚𝑇B+ 𝛺(𝑇B ) cos (𝑞 ln )}⁄sin(𝑞 ln 𝑚),𝑚𝑚 𝑝𝑇B
𝐶2 = {[𝐵𝑚 ( ) − 𝛺 ( )] sin(𝑞 ln 𝑇B )𝑇B𝑚𝑇B+ 𝛺(𝑇B ) sin (𝑞 ln )}⁄sin(𝑞 ln 𝑚)𝑚(10)ここに、𝑗(8)を得る。変数を T に戻せば、βa の温度依存性から Bの温度依存性を計算する公式として次式を得る。
𝐵 = 𝑇 𝑝 [𝐶1 sin(𝑞 ln 𝑇) + 𝐶2 cos(𝑞 ln 𝑇)]
𝛾0 𝛽𝑗 𝑇𝑗+∑
(𝛾0 − 1)2 𝑗 2 + (𝛾02 − 1)𝑗 + 2𝛾0とは明確に異なる任意の正の数とする。これらの初𝑗(9)式(9)の第 1 項および第 2 項は減衰振動のようにふるまい、その振幅は C1 と C2 で特徴づけられる。これ
らの項を以下では「減衰振動項」と呼ぶことにする。
これらは低温領域では B の温度依存性に影響するが、温度上昇につれて急激に減衰する。γ0 = 5/3 の場合 p= –2、γ0 = 7/5 の場合 p = –3、γ0 = 4/3 の場合 p = –7/2
だからである。高温では、B の温度依存性は式(9)の最後の項のみで与えられる。この最後の項を、以下「シグマ項」と呼ぶことにする。
C1 と C2 を与える初期条件を検討するため、典型的な B の関数形を思い起こすことにする。低温では主に分子間引力のために B は負の値をとる。温度上昇につれて、分子の運動エネルギーが増加し分子間力と競合し始める。また、実際の分子は排除体積をもち、高い温度ではその効果が大きくなる。そのため B は温度上昇とともに増加し、ボイル温度 TB でゼロになり、さらに正の値をとるようになる。実際の分子は完全な剛体ではなく、さらなる高温では実効的な排除体積は減少し始める。温度が無限大に近づくにつれ、排除体積はゼロに漸近し、それゆえ実在気体のふるまいは理想気体のそれに近づいていくと考えられる。したがって、B はある温度で極大をとった後、緩やかにゼロに漸近すると考えられる。
このような B の温度依存性を考慮すると、「T = TBで B = 0」および「極大をとる温度 Tmax で dB/dT = 0」が適当な 2 つの初期条件と思われるかもしれない。
しかしながら、後者の条件は、数学的には適当でも物理学的には適当でない。なぜならば、ほとんどの分子において Tmax が高温すぎるために、その温度に到達する前に分子が熱解離してしまうか、仮に分子が安定だとしても、実験的にその温度を達成するこ𝑗−𝑝

𝛺(𝑇B ) = ∑𝑗

𝛾0 𝛽𝑗 𝑇B,
22(𝛾0 − 1) 𝑗 + (𝛾02 − 1)𝑗 + 2𝛾0
𝑇

𝑗−𝑝𝛾0 𝛽𝑗 ( B)
𝑇B𝑚
𝛺( ) = ∑(𝛾0 − 1)2 𝑗 2 + (𝛾02 − 1)𝑗 + 2𝛾0𝑚𝑗(11)である。
3. ヘリウム 4（4He）
ヘリウム原子は 2 個しか電子をもたないため、近年、高精度な第一原理計算が可能となってきた。ビリアル係数においても理論計算の方が実験よりも高精度である。本稿では、最近発表された理論計算の論文[10]のビリアル係数のデータを利用した。文献[10]に報告されている 0.5 K から 1000 K までの 118 個のβa のデータを、熱力学温度のべき級数で、標準偏差0.29 cm3 mol-1 で、次のようにフィッティングした。
𝛽a𝑇 3−8=−4.1420198
×10()
cm3 mol−1K
𝑇 2+ 8.2295908 × 10−5 ( )K
𝑇− 5.9056759 × 10−2 ( )K𝑇 −1+ 34.541662 − 465.67255 ( )K−2𝑇− 91.990366 ( )K𝑇 −3− 4.7918563 ( )K𝑇 −4− 10.250695 ( )K𝑇 −5+ 2.3086230 ( )K(12)図 1 0.5 – 10 K での 4He の B の温度依存性。シン

図 2 10 – 1000 K での 4He の B の温度依存性。シンボルは文献値[10]を表す。点線、破線、実線は、式(13)ボルは文献値[10]を表す。点線、破線、実線は、式(13)における「減衰振動項」、「シグマ項」、それらの和、における「減衰振動項」、「シグマ項」、それらの和、をそれぞれ表す。

をそれぞれ表す。破線と実線は、ほぼ重なっている。

4He の場合には、γ0 は 5/3 の定数である。B の文献[10]（T = 0.5 – 10 K の範囲）および図 2（T = 10 – 1000 Kから、ボイル温度は TB = 23.3 K と見積もることの範囲）に示す。全体的に見て、式(13)は文献値[10]をができ、また、m = 2 を選ぶと、T = TB/2 = 11.65 K によく再現している。
「減衰振動項」を考慮しなければ、

値

3-1

おいて Bm=2 = –17.2 cm mol と見積もることができ

本公式による温度依存性すなわち「シグマ項」は、た。これらのパラメータを用いて、前章で求めた公1 K 以下の領域では、温度が下がるにつれて次第に式により、B の温度依存性は、文献値[10]から離れていき、0.5 K で両者の差は 122

𝐵𝑇 −27𝑇=() [31.63767 sin (√ ln ( ))3−1cm mol
K2Kcm3 mol-1 まで広がる。
「減衰振動項」と「シグマ項」との和が、数学的にも物理学的にも適当な、式(3)の微分方程式の解であるため、「減衰振動項」はこの差

7𝑇− 4.345371 cos (√ ln ( ))]2Kを補償し、この温度での両者の差を 5 cm3 mol-1 まで5𝑇 3+× [−4.1420198 × 10−8 ( ) ]38K

は、絶対値で 0.3 cm3 mol-1 より小さくなる。図 2 の

5𝑇 2× 8.2295908 × 10−5 ( )26K
3𝑇−2+× [−5.9056759 × 10 ( )]10K1+ × 34.541662
25𝑇 −1+ × [−465.67255 ( ) ]6Kの標準偏差で、B の温度依存性が再現されている。

抑えている。
「減衰振動項」の寄与は、10 K 以上で全温度領域では、「シグマ項」のみで 0.13 cm3 mol-1+15𝑇 −2+× [−91.990366 ( ) ]14K
5𝑇 −3+ × [−4.7918563 ( ) ]6K
1𝑇 −4+ × [−10.250695 ( ) ]2K
3𝑇 −5+× 2.3086230 ( )10K(13)と計算できた。B の計算値[10]と式(13)との比較を図 1

4. 窒素（N2）N2 については、80 K から 400 K の範囲で、βa の実験値が 16 点報告されている[11, 12]。これらの βa の値を、熱力学温度のべき級数で、0.19 cm3 mol-1 の標準偏差で、次のようにフィッティングした。
𝛽a𝑇 −14=85.00−
1.6265×10()
cm3 mol−1K
𝑇 −2− 8.078 × 105 ( )K𝑇 −3− 1.707091 × 107 ( )K(14)アメリカ国立標準技術研究所（NIST）の開発したREFPROP ソフトウェア[13, 14]を用いると、N2 の γ0 は、400 K 以下では 1.400 ± 0.003 の範囲にあると求めら5. 結論γ0 が温度によらない仮定の下で、βa と B を結ぶ微分方程式を解くことで、βa の温度依存性から B の温度依存性を計算する公式を導いた。極低温寒剤として広く利用されている単原子分子の 4He と二原子分子の N2 について、本公式により B の温度依存性をよく計算できることを確認した。このことから、本公式は、γ0 が近似的に温度によらない場合に、単原子分子のみならずその他の分子においても有効であることを示した。本公式で低温領域を計算する場合には「減衰振動項」の寄与を無視できないため、初図 3 80 – 400 K での N2 の B の温度依存性。シンボ期条件として 2 温度点での B の値が必要である。一ルは実験値[15-18]を表す。点線、破線、実線は、式(15)方、「減衰振動項」を無視できる程度の高温領域では、における「減衰振動項」、「シグマ項」、それらの和、そのような初期条件は必要なくなり、βa の温度依存をそれぞれ表す。破線と実線は、ほぼ重なっている。

性のみから B の温度依存性を容易に計算できる。

れた。したがって論文[11, 12]の温度範囲では近似的にγ0 は温度によらず 7/5 で一定とみなせる。B の実験
値[15-18]から、ボイル温度は TB = 327 K と見積もられ、また、m = 2 を選ぶと、T = TB/2 = 163.5 K においてBm=2 = –59 cm3 mol-1 と見積もられた。これらのパラメータにより、B の温度依存性は、
𝐵𝑇 −317𝑇=() [1.724942 × 104 sin (√ ln ( ))3−1cm molK2K+ 2.000459 × 104 cos (√17𝑇ln ( ))]2K1+ × 85.00
27𝑇 −1+× [−1.6265 × 104 ( ) ]10K+35𝑇 −2× [−8.078 × 105 ( ) ]38K

35𝑇 −3× [−1.707091 × 107 ( ) ]34K(15)[15-18]と計算できた。文献値と式(15)との比較を図 3 に+示す。図 3 のスケールでは両者はよく一致している。
式(15)は、文献値[15-18]を 0.82 cm3 mol-1 の標準偏差で再現した。この標準偏差は、式(15)のフィッティングの標準偏差および文献[11, 12, 15-18]の不確かさの和と同程度である。したがって、式(9)が単原子分子のみならず二原子分子にも有効であることを示した。このN2 の場合では、「減衰振動項」の寄与は 80 K 以上で0.05 cm3 mol-1 以下であり、図 3 の全温度領域で、「シグマ項」のみで B の温度依存性が再現されている。

参考文献[1] K. Tada, Int. J. Thermophys. 43, 64 (2022).
[2] K. Tada, Int. J. Refrig. 127, 157 (2021).
[3] K. A. Gillis and M. R. Moldover, Int. J. Thermophys.
17, 1305 (1996).
[4] 宝積勉、佐藤春樹、渡部康一 熱物性 11, 10(1997).
[5] L. W. Bruch, Phys. Rev. 178, 303 (1969).
[6] L. W. Bruch, Phys. Rev. A 2, 2167 (1970).
[7] M. E. Boyd and R. D. Mountain, Phys. Rev. A 2, 2164
(1970).
[8] F. C. Matacotta, G. T. McConville, P. P. M. Steur andM. Durieux, Metrologia 24, 61 (1987).
[9] G. T. McConville and J. J. Hurly, Metrologia 28, 375(1991).
[10] P. Czachorowski, M. Przybytek, M. Lesiuk, M.
Puchalski and B. Jeziorski, Phys. Rev. A 102, 042810(2020).
[11] M. B. Ewing and J. P. M. Trusler, Physica A 184, 415(1992).
[12] J. F. Estela-Uribe and J. P. M. Trusler, Int. J.
Thermophys. 21, 1033 (2000).
[13] E. W. Lemmon, I. H. Bell, M. L. Huber and M. O.
McLinden, NIST Standard reference database 23:
Reference fluid thermodynamic and transport propertiesREFPROP, Version 10.0, National Institute of Standards
and Technology, Standard reference data program,Gaithersburg, 2018.
[14] R. Span, E. W. Lemmon, R. T. Jacobsen, W. Wagnerand A. Yokozeki, J. Phys. Chem. Ref. Data 29, 1361
(2000).
[15] W. Duschek, R. Kleinrahm, W. Wagner and M.
Jaeschke, J. Chem. Thermodyn. 20, 1069 (1988).
[16] P. Nowak, R. Kleinrahm and W. Wagner, J. Chem.
Thermodyn. 29, 1137 (1997).
[17] M. O. McLinden and C. Lösch-Will, J. Chem.
Thermodyn. 39, 507 (2007).
[18] I. D. Mantilla, D. E. Cristancho, S. Ejaz, K. R. Hall,
M. Atilhan and G. A. Iglesias-Silva, J. Chem. Eng. Data55, 4227 (2010).