強い斥力ポテンシャルが作用する 多粒子系の気体−固体相

強い斥力ポテンシャルが作用する
多粒子系の気体−固体相転移
に対する量子効果の研究
名古屋大学大学院 理学研究科 物性理論研究室
山下 耕平
概要
一般に与えられた系の相図に対する量子効果はきわめて重要である。本論文では、多粒子
系の気体―固体相転移に対する量子効果を明らかにすることを目的とする。
ヘリウムは、常圧下で絶対零度まで固化せずに液相にとどまる。これは相互作用が弱い
うえに、原子質量が小さく零点運動が大きいために、原子の局在が妨げられるためと考え
られている。つまり、量子効果によって固化が抑制されている。一方、古典および量子剛
体球に対する数値計算では、量子性が強いほど固化しやすいという結果が得られている。
本研究では、この一見すると相矛盾する結果をはじめとして、多粒子系の気体‐固体相転
移に対する量子効果の包括的理解を目指す。
このために、本研究では原子間距離が近づくに従い冪的に発散する斥力相互作用を及
ぼしあう粒子系を考える。粒子としては、区別できる粒子およびボソンの場合を考える。
引力成分は存在しないので液相への凝縮の可能性を考える必要はなく、広いパラメタ―範
囲で気体‐固体相転移を研究することができる。また、この系は、極限として剛体球系も
含む。この系に対して、拡散モンテカルロ法および経路積分モンテカルロ法を適用し、基
底状態および有限温度における相図を確定する。その際、粒子の質量を変化させることに
よって零点運動の大きさを変化させ、これが相図に与える影響を明らかにする。
数値計算に先立ち、まず、この系におけるスケーリング則を調べた。その結果、零点
運動には相反する２つの効果があることが明らかとなった、ひとつは有効密度を小さくす
る効果であり、これは固化を妨げる。もうひとつは、有効温度を減少させる効果であり、
これは固化を促進する。剛体球系においては後者のみが存在する。一般には、両者の競合
により、固化に対する零点エネルギーの効果が決まる。
数値計算の結果、温度ゼロでは零点運動による有効密度減少の効果のみが見られ、量
子効果によって固化は抑制されることが確かめられた。一方、熱エネルギーが零点エネル
ギーより大きい有限温度領域で零点エネルギーを大きくすると、有効温度減少の効果が勝
り、固化が促進されることが明らかとなった。つまり、熱エネルギーが支配的な古典領域
から零点エネルギーが支配的な量子領域へクロスオーバーする中間の過程において、零点
運動による固化の促進が観測されることが示された。
さらに、圧力‐温度相図において、融解圧曲線に弱い極小が現われることを見出し
た。これは固体相が気相に比べて高エントロピーの状態にあることを示唆する。この現象
を理解するために、固体相におけるずり弾性を調べ、相境界近傍でこれがソフト化するこ
とを見出した。このため、固体相での横波フォノンの速度が減少し、エントロピーの相対
的な増加の原因となっていると考えられる。
目次
第 1 章 はじめに
1
1.1
研究の背景 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.2
目的 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
第 2 章 相互作用と相平衡
6
2.1
相互作用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.2
相平衡 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
第 3 章 スケーリング則
10
第 4 章 温度ゼロにおける気体-固体相転移
13
4.1
手法：拡散モンテカルロ法
. . . . . . . . . . . . . . . . . . . . . . . . .
13
4.2
結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
4.3
考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
第 5 章 有限温度における気体-固体相転移
17
5.1
手法：経路積分モンテカルロ法 . . . . . . . . . . . . . . . . . . . . . . .
17
5.2
Helmholz の自由エネルギーの計算方法 . . . . . . . . . . . . . . . . . . .
18
5.3
synthetic data 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
22
5.4
区別出来る粒子の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
5.5
ボース粒子の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
第 6 章 まとめ
34
付録 A Monte Carlo integration
35
付録 B Monte Carlo simulations
37
付録 C Metropolis algorithm
42
付録 D Variational Monte Carlo
45
付録 E
Diffusion Monte Carlo
46
付録 F
Path Integral Monte Carlo
48
第1章
第1章
1.1
はじめに
はじめに
研究の背景
一般的に希ガスの単原子分子などでは、質量が軽く量子性が強い物質ほど零点運動が激
しくなり、固体になりにくいと考えられている（表 1）
。常圧で温度を下げてゆくと、通常
では固体への相転移が起きるのに対し、液体ヘリウムでは絶対零度まで液体のままにとど
まり、加圧下でなければ固化しない。一方でヘリウムは零度に近い温度では、他の物質で
は見られない超流動性を示す液体へと相転移する。これらの振る舞いは、ヘリウム原子の
質量が小さいことにより、通常は無視出来るような零点運動の寄与が重要になったためと
理解されている。この振る舞いについて議論するために、Lennard-Jones ポテンシャルが
働く系を考えてみる。
U (r) = 4
[( )
σ 12
r
−
( σ )6 ]
(1.1)
r
このポテンシャルは、r = σ で U (r) = 0 となり、原子同士がこの距離より近づくと Pauli
の原理によって強い斥力が働く。反対に距離が r > σ になると、van der Waals 引力
(−(r/σ)6 ) の寄与が大きくなり、r0 = 21/6 σ で極小値をとって、それより離れるにした
がって引力は弱まる。表 2 に示すように、 と σ の値は原子により異なっている。そこ
で、エネルギーを 、距離を σ でスケールしてみると、異なる原子の相互作用はお互いに
似た形になり、共通な関数に支配される系となる。
U (r) = U ∗ (r/σ)
(1.2)
このときの熱力学的性質では de Boer の対応原理が成立すると言われており、例えば温度
T と体積 V 、圧力 P を σ 、 で換算した無次元量 (reduced unit)
v∗ =
V
,
N σ3
T∗ =
kB T
,
p∗ =
P σ3
(1.3)
を用いて、Helmholtz の自由エネルギーは以下のように表される。
f ∗ = F/N = F ∗ (v ∗ , T ∗ )
(1.4)
対応原理ではすべての物質に対して F ∗ が共通な関数になり、U ∗ のみによって決まる。
これはネオンやアルゴンなどの重い希ガス原子に対してはよい精度で成立しているが、ヘ
リウム原子に対しては成り立っていない。これはヘリウム原子については零点エネルギー
が無視出来ない寄与をおよぼすためである。そこで零点運動と引力相互作用の比を表す、
量子パラメター η = ~2 /mσ 2 という量を導入する。量子パラメターが大きな物質ほど、
零点運動が激しくなり量子効果が顕著になると期待される。量子パラメターを導入したこ
とにより Helmholtz の自由エネルギーは
F ∗ = F ∗ (v ∗ , T ∗ , η)
1
(1.5)
第1章
表1
融点
表2
はじめに
希ガス単原子分子の融点
He
Ne
Ar
Kr
Xe
14.9K (1050atm)
24.56K
83.80K
115.79K
161.4K
さまざまな物質の量子パラメター η と Lennard-Jones パラメター
m (amu)
(K)
σ (A)
η
3
He
3.016
10.22
2.556
0.2409
4
He
4.003
10.22
2.556
0.1815
H2
2.016
37.0
2.92
0.0763
D2
4.028
37.0
2.92
0.0382
Ne
20.18
35.6
2.74
0.0085
Ar
39.95
120.0
3.41
0.00088
となり、一般に η 依存性を持ち、F ∗ は共通な関数ではなくなる。
Nosanow らは、変分モンテカルロ法 (Variational Monte Carlo method) を用いて、温
度ゼロにおいて量子パラメターを変化させたときの、固体−液体相転移を調べた 1) 。図 1
はボース粒子である希ガス原子の計算結果である。体積 (reduced volume) が大きく希薄
な領域に液体相、体積が小さい領域に固体相が存在し、量子パラメターが大きくなるにつ
れて、凝固体積 (VL∗ ) が減少している。つまり、量子パラメターが大きな物質ほど高密度
図1
ボース統計に従う粒子の固体−液体相転移：V ∗ , P ∗ は式 (1.3) と同じ定義。VL∗
は凝固曲線、VS∗ は融解曲線である。
（Nosanow ら 1) の論文より）
2
第1章
はじめに
にしなければ固化せず、量子性によって固化しにくくなっている。この結果は、一般に考
えられている「質量が軽く量子性が強い物質ほど、零点運動が激しくなり固化しにくく
なる」という事と定性的に一致している。なお小池らによって、より精度の良い方法（グ
リーン関数モンテカルロ法）で η ≤ 0.3 まで調べられ、同様な結果が得られている 2) 。
ここまでは Lennard-Jones ポテンシャルという、斥力と引力をもつ相互作用が働く系
について考えて来たが、粒子を剛体球として扱った場合について見ていこう。剛体球の計
算は古くから盛んにおこなわれている。その中でも、Alder らが分子動力学法 (Molecular
Dynamics method) を用いて、剛体球の系が圧縮していくと固体になるという結果を示し
たのが有名である 3) 。気体のとき、剛体球は図 2 のスナップショットに見られるように、
空間内を自由に動き回っているが、圧縮され固体となると図 3 のように格子点まわりを動
くだけになる。この相転移は Alder 転移と呼ばれており、凝固密度は 0.94/σ 3 と見積もら
れ、最密充填密度の約 0.67 倍に相当する 4) 。ここで σ は剛体球の直径のことである。後
図2
図3
気体状態のスナップショット
（Alder ら
3)
固体状態のスナップショット
（Alder ら 3) の論文より）
の論文より）
の研究で、剛体球だけでなく、冪乗の斥力ポテンシャルのみが働く系においても固化が起
こることが分かっており、固体−液体相転移の本質は斥力相互作用にあることが示唆され
非常に興味深い。（一方で気体−液体相転移の本質は引力相互作用にあることが分かって
いる。）14 年後、Kalos たちは温度ゼロにおける手法であるグリーン関数モンテカルロ法
(Green Function Monte Carlo method) を用いて、剛体球が量子的な枠組みのなかでも
気体−固体相転移を示すのか検証をおこなった 5, 6) 。その結果、量子剛体球も固化するこ
とが分かり、その凝固密度は 0.25/σ 3 で、最密充填密度に対し約 0.18 倍と見積もられた。
Alder らの結果と Kalos らの結果を比較すると、量子剛体球のほうが低密度で固化してい
るため、先ほど述べた一般的に言われている事とは反対の量子性によって（量子的な取り
扱いをすることによって）固化が助けられる事を示している。この量子性によって固化が
助けられる振る舞いは、古典系と温度ゼロにおける量子系の比較だけでなく、有限温度に
おける量子剛体球の気体−固体相転移で顕著な振る舞いとして報告されている。Runge
3
第1章
はじめに
0.8
transition density ρ σ
3
Solid
0.6
0.4
0
図4
Gas
1
de Broglie wavelength λ(T) / σ
2
量子剛体球の有限温度における気体−固体相転移：赤線と○シンボルは凝固曲
線、青線と△シンボルは融解曲線を表す。
ら 7) と Sese ら 8) は、有限温度における剛体球の気体−固体相転移を、経路積分モンテカ
ルロ法 (Path Integral Monte Carlo method) 用いて計算をおこなった。次章にて詳しく
述べるが、剛体球における特徴的な物理量はド・ブロイ波長 λ(T ) =
√
2π~2 /mkB T であ
り、様々なド・ブロイ波長の場合について転移密度をプロットしたのが図 4 である。ド・
ブロイ波長は質量が軽く量子性が強い物質ほど長くなる一方で、温度を下げても波長は長
くなるため、質量を小さくすることと、温度を下げることは同じ意味を持つ。彼らの結果
によると、ド・ブロイ波長が長くなると転移密度が小さくなっている。つまり量子性が強
い物質ほど低密度で固化し、固体相が量子性（零点運動）によって助けられて安定化して
いるという事を示している。以上のように、剛体球の計算においては一般に考えられてい
る事とは逆に、量子性によって固化しやすくなることが分かっている。
量子性によって固体相が安定化している同様の振る舞いが、２次元グラファイト上に吸
着された水素分子について観測され 9) 、また理論的 10, 11) にも示唆されている。図 5 はグ
ラファイト基盤上に吸着された重水素分子と水素分子に対して、比熱測定実験などをおこ
√
なって得られた状態相図である。この図において、 3 整合固体相について着目してみる
と、水素は重水素よりも質量が軽く零点運動が激しい（量子性が強い）にも関わらず、高
温まで整合相が安定である事が分かる。この
√
3 整合相については理論的な研究もおこな
われている。彼らは Nosanow らと同様な量子パラメターを導入し、経路積分モンテカル
ロ法を用いて計算をおこなった。その結果、量子パラメターが大きく、零点運動が激しい
物質ほど融解温度が高くるという、実験と定性的に一致した結果を得ている（図 6）
。２次
元面というバルクとは次元が異なる場合においても、量子性によって固体相が安定化して
いるというのは興味深い。
この他にも、零点運動の効果によって固体相が安定化していると考えられる振る舞い
が、さまざまな系において報告されている 12–14) 。
4
第1章
はじめに
図 6 ２次元面内に移動が制限された粒子
の融解温度：横軸の q は量子パラメター η
と同様な定義であり、σ はハードコア直径
図 5 グラファイト基盤上に吸着されたは
を表し、○は σ = 1.05、●は σ = 1.10、■
重水素 D2 の相図 (a)、水素 H2 の相図 (b)
は σ = 1.15 を表す。
（Freimuth ら
1.2
9)
（Hirashima ら 11) の論文より）
の論文より）
目的
ここまで見てきたように、Lennard-Jones ポテンシャルを用いた温度ゼロの計算では、
量子性によって固化しにくくなっている一方で、剛体球の計算では量子性が強くなると固
体相が安定化するという、相反する結果が先行研究によって得られている。先行研究の問
題点を考えると、(1)Nosanow らの計算は現実に近いポテンシャルを使っているにも関わ
らず温度ゼロのみの計算であり、またあまり精度のよくない手法を用いている、(2) 剛体
球系の計算は、有限温度についても精度のよい手法で計算されているが、粒子を剛体球と
して扱っているため、Lennard-Jones ポテンシャルなどを用いた計算と大きく異なる結果
を示す可能性がある。
そこで本研究では、冪乗ポテンシャル (式 (2.2)) が粒子間に働く系において、拡散モン
テカルロ法と経路積分モンテカルロ法を用いて、温度ゼロと有限温度のそれぞれの場合
に、量子パラメター η を指標として量子性の強さによって気体−固体相転移がどのように
変化するのかを調べる事を第一の目的とする。なお、この系は極限として剛体球系を含ん
でいる。また量子力学においては、粒子は区別出来ないものであり、ボース統計かフェル
ミ統計に従うはずである。しかし、シミュレーションの実行のし易さから、多くの研究で
粒子は区別出来るものとして扱われており、本研究においても区別出来る粒子に対しての
計算をメインにおこなう。そして区別出来る粒子で気体−固体相転移の大局的な振る舞い
を捉えてから、量子統計性（ボース統計性）を加味した計算を数点おこない、粒子の統計
性の違いが気体−固体相転移にどのような変化をもたらすのかを調べる事を第二の目的と
する。
5
第2章
第2章
相互作用と相平衡
相互作用と相平衡
2.1
相互作用
2.1.1
Inverse-Power-Law ポテンシャル
様々な物質の相互作用として、斥力と引力の両方を持った Lennard-Jones(LJ) ポテン
シャルが広く使われている。
[(
ULJ (rij ) = 4
σ
rij
(
)12
−
σ
rij
)6 ]
(2.1)
ここで rij は、i 番目の粒子と j 番目の粒子の間の距離を表す。LJ ポテンシャルは、物質
の定性的な振る舞いを議論をするのに非常に有効なポテンシャルであるが、引力相互作用
を持っているため液体−気体相転移が起こる場合がある。本研究では、物質の質量を変え
ながら広い密度領域で計算を行なうため、液体相が存在していると話が複雑になってし
まう。そのため、LJ ポテンシャルから引力の項を無くした、斥力相互作用のみのポテン
シャルを考えることにし、広いパラメター領域で固体−気体相転移に集中して調べる。
UIPL =
∑
(
4
i<j
σ
rij
)n
(2.2)
このポテンシャルは、Inverse-Power-Law ポテンシャルと呼ばれ、古くからこの相互作用
が働く系の振る舞いが調べられて来た 21) 。冪乗の指数 n が小さいときは長距離相互作用
であり、大きいときは短距離相互作用となる。特に n → ∞ の極限では剛体球の系を再現
する。本研究では、主に n = 12 の短距離相互作用について調べるが、温度ゼロにおいて
は n をさらに大きくした場合の振る舞いについても調べた。
2.1.2
相互作用のカットオフ
一般的に、粒子間の距離で決まるポテンシャルを用いてシミュレーションを行なう場
合、どこまでの距離の粒子を計算にいれるのか問題になる。なぜならば、N 個の粒子を考
えた場合、厳密に足し合わせるべき組み合わせの数は N C2 であり、粒子数が多い場合は
膨大な計算時間が必要となってしまうからである。そこでシミュレーション上で、r < rc
の距離までの粒子間相互作用はきちんと計算し、r > rc の距離の相互作用は働かないもの
とする、距離のカットオフを導入する。


u(rij ) =
(
4
0
σ
rij
)12
(rij < rc )
(2.3)
(rij > rc )
本研究では粒子間の相互作用として、LJ ポテンシャルよりも早く減衰する短距離の IPL
ポテンシャルを用いているが、計算精度向上のためカットオフを入れたことによる補正を
6
第2章
相互作用と相平衡
加える。距離 rc より遠方の粒子は均一に分布していると考えて、1 粒子あたりのポテン
シャルエネルギーの補正 Uc は２重数えを考慮し、
1
Uc =
2
∫
∞
rc
∫ ∞
u(r) ρ(r) dr3
(2.4)
)12
σ
4
4πr2 dr
rij
rc
( )9
8
σ
3
= ρσ π
9
rc
ρ
=
2
(
(2.5)
(2.6)
となる。また、圧力に対しても同様な補正が加わる。圧力は virial 定理
1
P =
3V
*
∑
du(rrij )
3K −
rij
drij
i<j
+
(2.7)
から求めることができる。ここで、V は系の体積であり、K は全体の運動エネルギーのを
表す。このときの圧力の補正は
]
∫ [
1 N ∞ dv
Pc = −
r
ρ(r)dr3
3V 2 rc
dr
[
∫
( σ )13 ]
ρ2 ∞
4 −12
4πr2 dr
=−
6 rc
r
( )9
32 2 3
σ
=
ρ σ π
9
rc
(2.8)
(2.9)
(2.10)
となる。カットオフを導入して計算された、ポテンシャルエネルギーと圧力のシミュレー
ション結果に、以上の補正を加えて計算結果の評価をおこなう。
7
第2章
2.2
相平衡
2.2.1
相境界の求め方
相互作用と相平衡
ここでは、気体相 (A) と固体相 (B) の相境界を求める方法について述べる。２つの相
が平衡である条件は、
TA = TB ,
PA = PB ,
µA = µB
(2.11)
である。化学ポテンシャルを圧力と温度の関数として表せば、両相の等しい温度と圧力を
T, P と書いて、方程式
µA (T, P ) = µB (T, P )
(2.12)
が得られる。この式から、平衡にある相の圧力と温度は、その一方を他方の関数として表
すことが出来ることを示している。式 (2.12) を温度で微分すると、
∂µA
∂µA dP
∂µB
∂µB dP
+
=
+
∂T
∂P dT
∂T
∂P dT
(2.13)
が得られ、さらに (∂µ/∂T )P = −s と (∂µ/∂P )T = v の関係式を使って、
dP
sA − sB
=
dT
vA − vB
(2.14)
と書き直される。ここで s は１粒子あたりのエントロピー、v は１粒子あたりの体積のこ
とである。これはクラジウス・クラペイロンの式と呼ばれ、平衡にある２相の圧力が温度
変化にともなってどう変わるかを決定する。
化学ポテンシャルは１粒子あたりの Gibbs の自由エネルギー G であるから、化学ポテ
ンシャルに対する条件は、
GA = GB
(2.15)
と表す事もできる。また Gibbs の自由エネルギー G は Helmholz の自由エネルギー F を
使って G = F + P V と書けるので、
FA + P VA = FB + P VB
(2.16)
と書き換えられる。これにより、
P =−
FA − FB
VA − VB
(2.17)
という関係式が得られる。これは、それぞれの相の Helmholz の自由エネルギーを体積に
対して求め、その共通接線の傾きが共存圧力 P となり、接点 VA と VB は相転移が始まる
（終わる）体積となることを意味する。また VA と VB の間は、状態 A と状態 B が混ざっ
た共存領域となっている。
なお本研究では、固体相と気体相の相転移を考えており、共存領域の固体側の相転移点
を融解体積 vml 、また共存領域の気体相側の相転移点を凝固体積 vfr と簡単のために呼ぶ
事にする。
8
第2章
2.2.2
相互作用と相平衡
粒子の初期状態と平衡状態
ここでは、シミュレーションに用いる粒子の初期状態と平衡状態について述べる。本研
究では固体−気体の相転移を調べるため、固体・気体それぞれの状態をシミュレーション
で再現しなければならない。本来であれば、粒子の初期状態はランダム構造からスタート
して、平衡状態を求めるべきである。しかし、ランダム構造から始めると、予期しない準
安定状態に落ち着いたりするなどして、シミュレーション時間が非常に長くなってしま
う。そのため、今回は固体相のシミュレーションを行なう場合には、FCC の構造を持っ
た状態から始め、また気体相については変分モンテカルロ法を用いて、あらかじめランダ
ム構造から気体の状態の初期粒子配置を作ることにした。これらの状態からシミュレー
ションを始めることで、少ない計算時間で平衡状態に落ち着かせることが出来る。
また平衡状態において、系が固体であるか気体であるか判別出来なければならない。そ
のため、状態の指標として静的構造因子を計算することによって、状態を判別する。FCC
の固体状態の場合、静的構造因子は波数
√
8 3π
k1 =
,
L
k2 =
16π
,
L
k3 = · · ·
(2.18)
において、デルタ関数的なピークを持つ。ここで L は系の長さ L = V 1/3 である。図 7
は実際のシミュレーションで得られた静的構造因子をプロットしたものであり、FCC 構
造を示す波数 k1 と k2 において大きな値を持つ事が確かめられた。また図 8 に示すよう
に、気体状態においては静的構造因子に固体のような鋭いピークは現れず、波数に対して
なだらかなカーブを描く。このように静的構造因子を計算することで、系がどのような状
態にあるのかを判別した後に、固体・気体の状態のサンプルとして採用した。
100
S(k)
S(k)
1
0
0
2
4
6
図7
2
4
6
k
k
T ∗ = 1, η = 0.1815, ρ∗ = 0.6 の固
図8
T ∗ = 1, η = 0.1815, ρ∗ = 0.33 の気
体状態における静的構造因子
体状態における静的構造因子
9
第3章
第3章
スケーリング則
スケーリング則
この章では粒子間に
U (r) = 4
( σ )n
r
(3.1)
の相互作用が働く系のスケーリング則を考えてゆく。この系の特徴的な長さスケールとし
て、零点エネルギーとポテンシャルエネルギーが同程度になる距離 r0 が考えられる。
( σ )n
~2
=
mr02
r
つまり、
−1
r0 = η n−2 σ
(3.2)
(3.3)
が距離のよいスケーリングとなる。このとき系のハミルトニアンは、
( )n
N
~2 ∑ 2 ∑
σ
H=−
∇i +
4
2m i=1
rij
i<j


( )n
N
∑
∑
n
1
1
˜ 2i +

= η n−2 −
∇
4
2 i=1
r̃
ij
i<j
n
= η n−2 H̃
(3.4)
(3.5)
(3.6)
となり、H̃ は η に明らかには依存しない距離 r0 でスケールされたハミルトニアンとなる。
n
また、くくりだされた η n−2 は系のエネルギースケール 0 と見なせ、温度 kB T と１粒子
あたりの体積 v = V /N は、
n kB T
kB T
= η − n−2
0
3
v
v
ṽ = 03 = η n−2 3
r
σ
T̃ =
(3.7)
(3.8)
とスケールされる。以上のことより、系の１粒子あたりの自由エネルギー f は、
f (N, V, T ) = 0 f˜(η − n−2
n
3
kB T n−2
v
,η
)
σ3
(3.9)
もしくは、
3
n kB T
v
f (N, V, T ) = kB T f˜(η − n−2
, η n−2 3 )
σ
(3.10)
と書く事ができる。これは IPL ポテンシャルが働く系においては、独立な変数が有効温
度 T̃ と有効体積 ṽ だけであることを示している。
なおスケーリング関数 f˜ は粒子の統計性によって異なることに注意してほしい。カノ
ニカルアンサンブルにおいて、系の自由エネルギーは分配関数を使って、F (N, V, T ) =
−kB T ln Z と表される。この分配関数 Z は粒子の統計性を反映しており、密度行列
ρ(R, R0 ; β) を使って書くと、
Z = Tr ρ(R, R0 ; β)
10
(3.11)
第3章
スケーリング則
{
R|e−βH |R
(区別出来る粒子)
ρ(R, R; β) = 1 ∑ −βH
|PR (ボース粒子)
P R|e
N!
(3.12)
となる。スケーリング関数は、この分配関数に関係しているため、考えている系の統計性
によって関数形が異なってくる。
有効温度と有効体積によって、系の固体−気体相転移が決まる事が分かったが、質量が
軽くなり量子パラメターが大きくなった場合に、２つの異なった量子効果が相転移に効い
てくることに注目したい。一つ目は、量子パラメター η が大きくなると有効体積が大きく
なる効果である。これは、零点運動が激しくなると、粒子同士が近づきやすくなり、粒子
の直径が有効的に小さくなるためであり、固体相を不安定化させる方向に働く。もう一つ
は、量子パラメター η が大きくなると有効温度が小さくなる効果である。こちらは反対に
固体相を安定化させる方向に働く。このことから、量子性が強くなった際に「有効体積の
増加効果」が支配的な領域では、固体相は不安定になり、また「有効温度の減少効果」が
支配的な領域では固体相が安定となる振る舞いをすることが示唆される。
ここまでは IPL ポテンシャルについて議論して来たが、剛体球についても考えみる。
IPL ポテンシャルでは、冪乗の指数 n を無限大に持っていけば、剛体球の系を再現する事
が出来る。n → ∞ とした場合、有効体積は η に依存しない定数となり、有効温度は
lim η − n−2
n
n→∞
kB T
mkB T σ 2
σ 2
kB T
=
=
∝(
)
2
η
~
λ(T )
(3.13)
と書き直せる。ここで λ(T ) はド・ブロイ波長を表す。つまり、剛体球の系においては、
量子性が強くなった場合には有効温度の減少効果しか現れず、固体相は零点運動によって
安定化することが期待される。実際、第 1 章で紹介した先行研究の計算結果では、ド・ブ
ロイ波長が長くなると、転移密度は小さくなり、固体相が安定化している。
最後に IPL ポテンシャルの古典極限について考えてみる。古典極限は η → 0 とすれば
良いのだが、独立変数を有効温度 T̃ と有効体積 ṽ の２変数としている場合、両方に η が
含まれているため、そのまま極限を考えるにはふさわしくない。そこで、スケーリング変
数としては２つのものを考えればよいので、η 依存のない変数を作り出すことにする。
(
T̃
3/n
ṽ =
kB T
)3/n
v
σ3
(3.14)
の関係が成り立つため、T̃ と ṽ を与える代わりに、T̃ 3/n ṽ と T̃ （もしくは ṽ ）をスケーリ
ング変数として与える。このときの元のスケーリングは、
f = 0 fˆ(T̃ , T̃ 3/n ṽ)
(3.15)
f = kB T fˆ(T̃ , T̃ 3/n ṽ)
(3.16)
または
と書き直され、二つ目の式において古典極限 η → 0 を考えると、
f = kB T fˆcl (T̃ 3/n ṽ)
11
(3.17)
第3章
スケーリング則
というスケーリング則が得られる。また圧力に関しては、自由エネルギーを体積で微分す
れば求まるので、
(
P =
kB T
)3/n
kB T cl 3/n
P̂ (T̃ ṽ)
σ3
(3.18)
というスケーリングが成り立つ。IPL ポテンシャル系の古典極限では、スケーリング変数
が T̃ 3/n ṽ だけであるため、固体−気体の相転移は T̃ 3/n ṽ がある値になったときに起こる
と期待される。つまり、
(
kB T
)3/n
v
=C
σ3
(3.19)
で転移体積 vtran は決まり、n = 12 のときは
vtran ∝ (kB T )−1/4
(3.20)
という温度依存性を持ち、同様に共存圧力 Pcoex に対しても、
Pcoex ∝ (kB T )5/4
となる。この温度依存性は、Hoover ら
21)
(3.21)
によって導かれた古典系の結果と一致して
いる。
以上のスケーリング則の議論をまとめると、有限温度の IPL ポテンシャルが働く系に
おいて、固体−気体相転移を特徴づける物理量は有効体積と有効温度であり、質量が軽く
なり量子性が強くなると、有効体積が増加して固体相を不安定化し、一方で有効温度は減
少し固体相を安定化させる。また温度ゼロでは、有効体積の増加の効果のみしかないた
め、量子性が強くなると固体相が不安定化させるようになる。
12
第4章
第4章
4.1
温度ゼロにおける気体-固体相転移
温度ゼロにおける気体-固体相転移
手法：拡散モンテカルロ法
温度ゼロにおいては、拡散モンテカルロ法 (DMC) を用いたシミュレーションを行なっ
ているが、ここではシミュレーションを実行する上で重要な試行関数について述べるだけ
にとどめる。手法の詳細については付録 E や、Reynolds ら 15) の論文を参照していただ
きたい。
拡散モンテカルロ法では、シミュレーションが上手くいくかは試行関数の選び方にか
かっている。試行関数は効率的にシミュレーションを行なうのに必要で、真の波動関数に
近いものを用意しなければならない。本研究では、気体相の試行関数には、Jastrow 関数
ΦJ (R) =
N
∏
(4.1)
f (rij )
i<j
を用い、固体相の試行関数には Nosanow-Jastrow 関数
ΦNJ (R) =
N
∏
N
∏
f (rij )
i<j
(4.2)
g(riI )
i,I=1
を用いた。ここで、i 番目の粒子と j 番目の粒子の距離を rij = |ri − rj | とし、riI は i
番目の粒子と結晶の I 番目の格子点との距離を表す。また f (r) = exp(−(a/r)(n−2)/2 )、
g(r) = exp(−0.5Ar2 ) であり、a と A はメインシミュレーションを行う前に VMC によっ
て最適化をおこなう。これらの関数は、ヘリウム４などの研究で広く使われており、基底
状態の真の波動関数に近いと考えられている。
また、気体相の試行関数 ΦJ は粒子の入れ替えに対して対称であるため、ボース粒子の
波動関数を表すが、固体相の試行関数 ΦNJ は対称になっておらず、ボース粒子の波動関
数になっていないことに注意してほしい。固体相における対称化された試行関数は幾つか
提案されているが、今回は Cazorla ら 16) によって提案された、
ΦSNJ (R) =
N
∏
f (rij )
i<j
N
∏
(N
∑
J=1
i=1
)
g(riJ )
(4.3)
を用いた計算を、幾つかの密度や量子パラメターで行った。その結果、対称化された試行
関数 ΦSNJ (R) と対称化されていない試行関数 ΦNJ (R) から計算されるエネルギーは誤差
の範囲内で一致した。そのため、本研究では計算時間がより短い ΦNJ (R) を用いて計算
をおこない、それらの結果をボース粒子のエネルギーと見なすことにする。
13
第4章
4.2
温度ゼロにおける気体-固体相転移
結果
Solid
(a)
eg,s / ε
20
η
25
−6/5
30
η=0.10
0.18
0.40
0.60
0.80
1.00
η=0.10
0.18
0.40
0.60
0.80
1.00
15
1.4
図9
Gas
η
3/10
1.6
v/σ
1.8
3
Gas
(b)
スケールされたエネルギーの有効体積依存性。クローズドシンボルは固体相の計
3
vfr / σ
3
算結果であり、オープンシンボルは気体相の計算結果である。
2
拡散モンテカルロ法を用いて、IPL ポテンシャルが働く系の気体−固体相転移を調べ
1
た。前章のスケーリング則の議論により、温度ゼロにおいては１粒子あたりのエネルギー
Solid
e と体積 v は、次のようにスケールされる事が分かっている。
0
0e
0.5
n
ẽ = η n−2
η
3
ṽ = η n−2
v1
σ3
(4.4)
そのため固体相、気体相どちらの場合でも、量子パラメター η と体積 v を様々な値に変
えてシミュレーションを行なった計算結果はスケールされるはずである。図 9 は、ポテ
ンシャルの指数が n = 12 の場合の計算結果をプロットしたものである。固体相、気体相
どちらにおいても、スケーリング則が上手く成り立っている。（固体相と気体相のエネル
ギーがかなり近いため、曲線が非常に見づらくなっている。そのため、共存圧力 Pc と転
移体積からのズレ ∆v = v − vml を使って、共通接線からの差分をプロットしたしたのが
図 10 である。）
よって、固体相、気体相のスケールされたエネルギーをそれぞれ多項式でフィッティン
グし、マクスウェルの規則によって、転移点を求めることにする。得られた接点のうち、
体積が大きい凝固体積 vfr と体積が小さい融解体積 vml とし、共通接線の傾きである共存
圧力 Pc をスケールされた物理量で表すと、
3
η 10
vml
= 1.55(1)
σ3
3
η 10
vfr
= 1.61(1)
σ3
η− 2
3
Pc σ 3
= 24.82(32)
(4.5)
となる。なお、括弧内の数字は誤差を表す。この結果から、スケールされていない元々の
転移体積は、量子パラメター η に対して、vml, fr ∝ η −3/10 のように振る舞うため、η が
14
温度ゼロにおける気体-固体相転移
Solid
(a)
Gas
η=0.10
0.18
0.40
0.60
24
23
22
1.00
η
−6/5
(eg,s − pc∆v) / ε
第4章
1.4
0.80
η
3/10
1.6
v/σ
Gas
(b)
3
1.8
3
3
図 10 固体相と気体相のエネルギーを、共通接線 pc ∆v との差分でプロット。実線は
vfr / σ
固体相の、鎖線は気体相のエネルギー曲線である。
2
大きく量子性が強い物質ほど転移体積が小さくなることが分かる。つまり、温度ゼロにお
1
いては強い零点運動によって固体相が不安定化することを示している。なおこの結果は、
Solid
0
用いたポテンシャルは異なるが、
Nosanow らが示した結果と定性的に一致している。
0
0.5
η
1
以上は、ポテンシャルの指数が n = 12 の場合についての結果であったが、指数を大き
くしていくと Kalos 達が計算した剛体球（n → ∞）の計算結果に近づくことが予想され
るため、指数を大きくした場合についても同様の計算をおこなった。その結果、図 11 に
示すように、指数が大きくなるほど量子パラメターを大きくしたときの転移体積の変化は
緩やかになり、温度ゼロにおける量子剛体球の転移体積 v/σ 3 = 4.0 に近づいて行くよう
になる。なお、指数を大きくした際に、体積変化が緩やかになることは、スケーリング則
vml, fr ∝ η −3/(n−2) が示す振る舞いと一致している。
15
第4章
温度ゼロにおける気体-固体相転移
3
reduced Volume V/(Nσ )
4
n=50
3
n=30
2
n=12
0
0.5
quantum parameter η
1
図 11 IPL ポテンシャルで冪乗の指数を変えた場合の気体−固体相転移相図。丸シン
ボルは n = 50、三角シンボルは n = 30、四角シンボルは n = 12 の場合で、クローズ
ドシンボルは融解体積、オープンシンボルは凝固体積を表す。
4.3
考察
スケーリング則の議論で、温度ゼロにおける特徴的な物理量は、有効体積だけであるこ
とを述べた。そのため、固体−気体の相転移は有効体積が、ある値になった際に起こると
3
予想され、前節で示した実際の計算でこのことが確かめられた。また有効体積は η 10
v
σ3
と表されるため、量子パラメターが大きくなると、スケールされていない元々の体積は小
さくなる事が分かった。この振る舞いは温度ゼロにおいて、エネルギーが有効体積のみの
関数である事に起因しており、もし特徴的な物理量がこの他にも現れた際には、また異
なった振る舞いをする事が考えられる。そのため次の章では、有限温度における計算を行
なうことによって、有効体積の他に有効温度という新しい特徴的な物理量が現れた系の固
体−気体相転移を調べてゆくことにする。
16
第5章
第5章
5.1
有限温度における気体-固体相転移
有限温度における気体-固体相転移
手法：経路積分モンテカルロ法
有限温度における計算では、経路積分モンテカルロ法 (PIMC) を用いているが、ここで
は本章で重要となる粒子の統計性について述べるにとどめる。手法の詳細については付録
F、または Ceperley17) のレビュー論文を参照していただきたい。
粒子が区別できるというのは、シミュレーション上において全ての粒子にナンバリング
をおこない、位置 ri にいる i 番目の粒子は虚時間 β = M τ 経過後に必ず元の位置 ri に戻
ることを意味している。そのため、区別できる粒子の経路は図 12 の赤実線のようになっ
ている。
しかし、ボース粒子の場合、波動関数が粒子の入れ替えに対して対称であるため、系の
状態 |Ri は、
|Ri =
1 ∑
|r1 , r2 , · · · , rN i
N!
(5.1)
P
を意味する。ここでは、位置のすべての置換について和をとっている。つまり、最後の座
標 |RM i と最初の |R0 i の対応する位置を同一視するだけでなく、両方の結びついた座標
で粒子位置の入れ替えも許すような虚時間に沿った周期境界条件を課すことになる。この
粒子の入れ替えのイメージが図 13 である。元々の経路が実線で表され、粒子の交換をあ
るタイムステップにおいて入れると、点線の経路にそれぞれ変更される。実際のシミュ
レーションでは、全てのタイムステップにおいて粒子の交換を試す必要はなく、また粒子
間の距離が離れている場合の交換は採択確率がかなり小さくなるため、近接粒子について
のみ入れ替えを考える。
図 12 区別出来る粒子の経路変更；元の経
図 13 ボース粒子の経路変更；i 番目の粒
路が赤実線であり、モンテカルロによって
子の元の経路を赤実線、j 番目の粒子の元
点線の経路に変更される。
の経路を青実線で示す。
17
第5章
5.2
有限温度における気体-固体相転移
Helmholz の自由エネルギーの計算方法
マクスウェルの規則に従い、気体−固体の相転移点を求めるためには、各状態の自由エ
ネルギーが分かっていなければならない。しかしながら、PIMC のシミュレーションで直
接得られる物理量は、内部エネルギー E と圧力 P だけであるため、自由エネルギー F を
求めるには工夫が必要となる。このような場合によく使われる方法として、熱力学積分と
いうものがある 18) 。カノニカルアンサンブルにおいて、圧力は自由エネルギーを体積で
１回微分した、
(
P =−
)
∂F
∂V
(5.2)
N,T
で表され、密度 ρ = N/V と、１粒子当たりの自由エネルギー f = F/N を用いて書き直
すと
P
=
ρ2
(
∂f
∂ρ
)
(5.3)
N,T
となる。そしてこの式を密度 [ρ0 , ρ] の区間で積分すると、
∫
ρ
f (ρ) = f (ρ0 ) +
dρ0
ρ0
P (ρ0 )
ρ02
(5.4)
となり、求めたい f (ρ) の値はリファレンス点 ρ0 から圧力を積分する事で計算できる。な
おシミュレーションにおいて、各密度における圧力は virial 定理を用いて計算する事が出
来る。
1
P =
3V
*
∑
dv(rij )
2K −
rij
drij
i<j
+
(5.5)
ここで K は運動エネルギー、v(rij ) は粒子間の相互作用を表し、h· · · i はアンサンブル平
均である。
以上のように熱力学積分をおこなえば、自由エネルギーは計算可能になる。しかしなが
ら、リファレンス点をどこにとるかという問題が残る。気体相の場合、理想気体となる希
薄極限をリファレンス点として扱うのが便利である。十分希薄な領域においては、理想気
体の状態方程式 P = ρkB T が成り立っており、式 (5.4) は
∫
ρ
fid (ρ) = fid (0) +
0
0
0 ρ kB T
dρ
ρ02
(5.6)
と表される。なお fid は理想気体の１粒子あたりの自由エネルギーである。
[
3
fid (ρ) = kB T ln ρ + ln
2
(
2π~2
mkB T
)
]
−1
(5.7)
そして式 (5.4) で ρ0 = 0 とした式から式 (5.6) を引くと、
∫
ρ
f (ρ) = fid (ρ) +
0
18
dρ0
P (ρ0 ) − ρ0 kB T
ρ02
(5.8)
第5章
有限温度における気体-固体相転移
となり、ある密度における気体相の自由エネルギー f (ρ) を希薄極限から圧力を積分する
ことによって計算することができる。
気体相については希薄極限から積分できたが、固体相ではさらなる工夫が必要となる。
固体相ではリファレンス点 ρ0 における自由エネルギー f (ρ0 ) を、アインシュタイン結晶
から計算するという方法をとる 18, 19) 。
アインシュタイン結晶とは格子点の周りを、調和振動子的に運動する粒子群のことであ
り、アインシュタインによって提唱された、固体のフォノン結晶のモデルである。このア
インシュタイン結晶からの計算をするために、次のような仮想ポテンシャルを導入する。
[
N
)2
k ∑(
(0)
Uλ (R) = (1 − λ)U (R) + λ U (R(0) ) +
ri − ri
2 i=1
]
(0)
(5.9)
(0)
(0)
ここで R = (r1 , r2 , · · · , rN ) は N 個の粒子の位置座標であり、R(0) = (r1 , r2 , · · · , rN )
は各格子点に粒子が存在している事を表す。k は格子点周りを振動するときのバネ定数で
あり、任意の定数である。λ はスイッチングパラメターと呼ばれ、λ = 0 のときに、仮想
ポテンシャルは Uλ (R) = U (R) となり、元々の系の相互作用になる。また λ = 1 のとき
には Uλ (R)=U (R(0) ) +
k
2
∑N
i=1 (ri
(0)
− ri )2 となり、格子点周りでの調和振動子ポテン
シャルを表す。つまり、λ = 0 で元々の斥力ポテンシャルが作用する系を表し、λ = 1 の
ときに、アインシュタイン結晶になっている系を再現するようになる。なお U (R(0) ) は、
すべての粒子がそれぞれの格子点に存在しているときの総ポテンシャルエネルギーであ
り、シミュレーション中では定数として扱われる。
では、このスイッチングパラメター λ を導入したことにより、分配関数や自由エネル
ギーがどうなるかみてみる。分配関数 Z で λ に依存している部分に注目すると、
∫
Z(N, V, T, λ) ∝
dr N exp(−βUλ )
(5.10)
と書かれる。また Helmholtz の自由ネルギーの λ 微分を考えると
(
∂F (λ)
∂λ
)
=−
N,V,T
∫
1 ∂
1
∂Z(N, V, T, λ)
ln Z(N, V, T, λ) = −
β ∂λ
βZ(N, V, T, λ)
∂λ
dr N (∂Uλ /∂λ) exp[−βUλ ]
∫
dr N exp[−βUλ ]
∂Uλ
=
∂λ λ
(5.11)
=
となる。ここで h· · · iλ は仮想ポテンシャル Uλ が働く系でのアンサンブル平均を表す。
さらに、この式を λ で積分してやると、
∫
λ=1
F (λ = 1) − F (λ = 0) =
dλ
λ=0
∂Uλ
∂λ
(5.12)
λ
となり、もし λ = 1 の自由エネルギーが分かっていれば、h∂Uλ /∂λiλ を区間 [0,1] の範囲
で積分することで、求めたい λ = 0 の自由エネルギーを計算することが出来る。
19
第5章
有限温度における気体-固体相転移
実際に使う式 (5.9) の仮想ポテンシャルで、この積分を書き直すと、
∫
*
λ=0
F = FEin +
dλ
λ=1
+
N
)2
k ∑(
(0)
ri − ri
− [U (R) − U (R(0) ]
2 i=1
(5.13)
λ
となり、アインシュタイン結晶の自由エネルギー FEin は、N 個の３次元量子調和振動子
として計算でき、~ω = (~2 k/m)1/2 を使って書くと、
(
)
3
FEin = U (R(0) ) + N ~ω + 3N kB T ln 1 − e−β~ω
2
(5.14)
と表される。ここで、仮想ポテンシャルの λ = 1 の場合において、相互作用に定数項
U (R(0) ) が含まれているため、自由エネルギー FEin にも U (R(0) ) を含めなければなら
ない。
以上のように、ある密度 ρ0 において仮想ポテンシャルを用いて、アインシュタイン結
晶から固体相の自由エネルギー f (ρ0 ) を求める事が出来るようになった。そして式 (5.4)
を用いて、任意の密度 ρ の自由エネルギー f (ρ) が計算可能になった。しかし、この固体
相の計算手法にもいくつかの注意すべき点がある。
一つ目は、粒子と格子点が１対１対応してナンバリングされているため、粒子はすべて
区別出来るもとのして扱われてしまう事である。もともと粒子が区別出来る場合の計算で
あれば問題はないが、ボース粒子の場合には不適切な計算をおこなってしまい、異なる結
果が得られる可能性がある。しかしながら、後の 5.5 節で詳しく述べるが、ボース粒子系
における固体相の内部エネルギーと圧力は、区別出来る粒子系の場合とほぼ一致するた
め、自由エネルギーの計算において問題は生じない。
二つ目は、仮想ポテンシャルにおいて、バネ定数が任意である事である。どのようなバ
ネ定数を用いても、式 (5.13) の被積分関数の値を計算で求めることは可能だが、最適化し
free energy fs*
0
< ∂Uλ / ∂λ >
-20
-40
30.7
30.6
-60
0
0.5
100
1
200
300
time slice M
λ
図 14
T ∗ = 0.5, η = 0.6, ρ∗ = 0.9 のとき
の被積分関数のスイッチングパラメター λ
T ∗ = 0.5, η = 0.6, ρ∗ = 0.9 のとき
の自由エネルギーのタイムスライス M 依
依存性
存性
図 15
20
第5章
有限温度における気体-固体相転移
ていないバネ定数を使うと、λ = 1 近傍で被積分値の変化が激しくなり、積分全体の精度
が低下してしまうのである。式 (5.13) の積分において、被積分値の変化を押さえるため
には、 k2
∑N
i=1 (ri
(0)
− ri )2 − U (R) の値を出来るだけ小さくするのがよい。これは、アイ
ンシュタイン結晶の相互作用と、元々の系の相互作用が出来るだけ近くなるようにバネ定
数を設定するということであり、実際には λ = 0 と λ = 1 における粒子の平均移動距離
が出来るだけ等しくなるようにする。
*
N
∑
(0)
(ri − ri )2
i=1
+
+
*N
∑
(0) 2
≈
(ri − ri )
i=1
λ=0
(5.15)
λ=1
このようにバネ定数を設定したときの被積分関数の λ 依存性をプロットしたものが図 14
である。この方法を用いても、λ = 1 近傍における激しい変化が見られるが、最適化して
いないバネ定数を使うとさらに変化が激しくなり、次に述べる数値積分において区分点の
数を非常に多くしなければならなくなる。
以上は、この手法自身の注意すべき点であったが、その他にも数値計算をおこなう際
に注意すべき点がある。積分にはガウス・ルジャンドルの数値積分法を用いる。この方法
では、N 個の区分点を用いた場合、2N − 1 次の高次の多項式まで近似でき、一般によ
く使われている積分の台形公式にくらべ、よい精度で計算することが出来る。図 14 では
N = 20 の区分点について計算したものであり、ガウス・ルジャンドル法では区分点は等
間隔ではなく、積分区間の両端近くでは細かく区分され、中央付近ではやや荒くなってい
る。そのため今回のように、λ = 1 近傍で変化が激しく、中央付近で変化が少ない場合に
は、有効的に計算することができる。なお、自由エネルギーの計算では N = 12, 16, 20, 24
などの区分数で計算を行い、N = 20 であれば問題がないことを確認している。
またこれらの計算には、経路積分モンテカルロ法を用いているため、虚時間のタイムス
ライス M 依存性があり、注意する必要がある。図 15 は、そのタイムスライス依存性を詳
しく調べたものである。通常の Lennard-Jones ポテンシャルなどでは M ≤ 120 程度で
十分であるが、アインシュタイン結晶の調和振動子ポテンシャルは、各格子点に局在して
おり、狭い範囲を正確にサンプルするためには、タイムスライスを多くしなければならな
くなる。特に λ = 1 近傍ではタイムスライス依存性に注意しなければならない。なお最
終的な自由エネルギー f (ρ0 ) の値には、タイムスライスを増やした際に、前後の数値が誤
差の範囲内で一致するような値を採用している。
21
第5章
5.3
有限温度における気体-固体相転移
synthetic data 法
前節において気体相では式 (5.8)、固体相では式 (5.4) と式 (5.13) を用いて圧力から自
由エネルギーを計算する熱力学積分について述べたが、この節では実際に積分を実行する
方法や、誤差評価について述べる。数値積分の台形公式や、ガウス・ルジャンドル法を用
いて積分値を求めることも可能ではあるが、積分区間が [0, ρ] などとなっているため、求
めたい密度 ρ において、その都度区分点を割り振る必要があり、あまり効率的ではない。
そのため本研究では圧力 P を、密度の関数としてフィッティングした後、その関数を用
いて積分を実行し、自由エネルギーを求めることにした。固体相の圧力は
Psolid =
m
∑
al (ρ − ρ0 )l
(5.16)
l=0
でフィットし、気体相の圧力は
(
Pgas = ρkB T
1+
m
∑
)
bl ρ l
(5.17)
l=1
でフィットした。固体相の ρ0 は係数 al を大きくせず、精度を保つために便宜上導入した
値であり、サンプルした密度のいずれかの値を使う。また気体相は希薄領域 ρ 1 におい
て、理想気体として振る舞うはずであるから、Pgas = ρkB T となるよう関数を設定した。
このようなフィッティング関数を用いて、熱力学積分をおこない自由エネルギーを求め
ていくが、積分等を実行しているため通常の四則演算で書かれている場合の誤差伝搬と
は、異なる誤差評価が必要となる。本研究では synthetic data20) という考え方を用いて
誤差評価をおこなう。図 16 は synthetic data 法の簡単なイメージフローである。まずシ
ミュレーションで得られるデータ群 D(0) は、それぞれ平均値を中心としたガウス分布を
していると期待される。つまり、各密度 ρi で得られる圧力のデータは、平均値 Pi を中
心として偏差 σi でガウス分布に従って広がっているはずである。このようなときに各密
図 16 synthetic data 法のイメージ
22
第5章
有限温度における気体-固体相転移
100
80
60
40
20
0
-0.04
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
0.04
図 17 synthetic data 法により生成された自由エネルギーの分布：T ∗ = 1, η = 0.1815
の固体相における圧力（赤）と、自由エネルギー（青）の分布。
度で、中央値が 0 で偏差が 1 のガウス分布によって生成される乱数 ζi を用いて、人工的
(α)
なデータ Pi
(α)
= Pi + σi ζi
を多数生成する。α は各密度で生成されたデータの番号で
あり、おおよそ Nd = 10, 000 個程度生成する。これで各密度における Nd 個のデータを
使って、synthetic data の組 D S (α) を Nd 個作ることが出来るので、それぞれの組に対
して χ2 フィッティングを行い、フィッティングパラメター aS (α) を求める。最後にパラ
メター aS (α) を使って、熱力学積分を実行すると、各 synthetic data の組に対して自由
エネルギーの値が得られる。フィッティングを行っているので、どのような密度の自由
エネルギーも計算可能だが、シミュレーションでサンプルされている密度 ρi と同じ点に
おいてのみ、自由エネルギーの値 f (ρi ) を計算している。密度 ρi において、Nd 個の自
由エネルギーの値が得られたため、そこから平均値と標準偏差を計算し誤差を評価する。
synthetic data によって生成された自由エネルギーの値もガウス分布に従っていると考え
られ、その分布をヒストグラムで描いたものが図 17 である。圧力（赤）から生成された
自由エネルギー（青）もガウス分布に従って分布していることがよく分かる。またデータ
の生成数 Nd を増やしても、ヒストグラムの形がよりキレイになり、標準偏差の精度がよ
くなるだけである事が確認できた。
23
第5章
5.4
区別出来る粒子の場合
5.4.1
区別出来る粒子の結果
有限温度における気体-固体相転移
16
Solid
free energy f*
12
8
Gas
4
1.2
1.6
2
reduced Volume v*
図 18
マクスウェルの規則：T ∗ = 1 で η = 0.1815 のときの、固体相と気体相の自由
エネルギー。誤差はシンボルよりも小さい。
まずは区別出来る粒子の場合について計算をおこなった。前節までで述べた方法で、
∗
T = 1, η = 0.1815 の場合における各状態の自由エネルギーを計算し、体積に対してプ
ロットしたものが図 18 である。（気体相における自由エネルギーの計算の際、希薄極限の
理想気体から熱力学積分をおこなうので、実際に計算したデータ点はさらに多い。）マク
スウェルの規則に従い、相転移点を求めるために、それぞれの自由エネルギーのデータを
体積または密度に対して多項式フィットする。
f (x) = a0 +
m
∑
ai xi
(5.18)
i=1
χ2 計算をしてみると、今回計算を行った温度 T と量子パラメター η の範囲では固体相は
m = 3、気体相は m = 5 が最適であった。この関数を用いて、共存圧力となる接線の傾
きと、融解、凝固の転移点となる接点を求める。
表 3 は温度 T ∗ と量子パラメター η を変えて、融解、凝固密度や共存圧力を求めた結
果である。そのなかで温度を一定にして、量子パラメターを変化させたときの凝固体積
vfr /σ 3 = 1/(ρfr σ 3 ) の変化をプロットしたものが図 19 になる。融解体積 vml は凝固体積
と同じ振る舞いを示すので、ここでは省略してある。温度一定で、質量を軽くし量子性を
強くすると（η が大きくなると）
、右上がりに凝固体積が大きくなっている。つまり、より
低密度で固化するようになるため、量子性によって固化が促進されていると言え、量子剛
体球の系で見られた振る舞いと同様な結果となっている。
24
第5章
表3
有限温度における気体-固体相転移
温度と量子パラメターを変えたときの凝固密度、融解密度、融解圧
凝固密度 ρ∗fr
融解密度 ρ∗ml
0.644 ± 0.008
0.663 ± 0.008
12.41 ± 0.52
T = 1, η = 0.6
0.521 ± 0.013
0.530 ± 0.014
10.66 ± 0.84
T ∗ = 1, η = 1.0
0.495 ± 0.038
0.505 ± 0.036
12.34 ± 2.82
T ∗ = 0.5, η = 0.1815
0.474 ± 0.006
0.483 ± 0.007
4.22 ± 0.20
T ∗ = 0.5, η = 0.6
0.437 ± 0.019
0.440 ± 0.022
6.15 ± 0.83
T ∗ = 0.5, η = 1.0
0.399 ± 0.045
0.401 ± 0.046
6.61 ± 2.09
T ∗ = 1, η = 0.1815
∗
融解圧 p∗c
vfr / σ
3
3
2
1
0
0
0.5
1
η
図 19 凝固体積の量子パラメター依存性：赤い四角シンボルは T ∗ = 1 のとき、黒い
丸シンボルは T ∗ = 0.5 の結果。青いダイヤモンドシンボルは有限温度における量子剛
体球のデータを T ∗ = 0.5 としてプロットした。また η = 0 のデータは先行研究で得ら
れている古典系 (剛体球、IPL ポテンシャル) の計算結果。また緑実線は前章の温度ゼ
ロにおける計算結果をプロットしたもの。
次に、第 3 章で述べたスケーリング則が成り立っているのか検証していく。スケーリン
グ則によれば、もとの温度 T 、圧力 P 、密度 ρ は
kB T̃ = η −6/5
kB T
,
P̃ = η −3/2
P σ3
,
ρ̃ = η −3/10 ρσ 3
(5.19)
にスケールされる。このスケールされた量によって、密度−温度相図と圧力−温度相図
を作成した。（図 20）どちらの相図においてもデータ点は少ないがスケーリングが上手く
いっているように見え、有効温度が減少すると凝固密度（または融解圧）が減少する振る
舞いを示している。また Hoover ら 21) がおこなった、本研究と同様の IPL ポテンシャル
を用いた、古典系におけるシミュレーションによると、
ρ̃fr = 0.813 T̃ 1/4 ,
P̃ = 15.95 T̃ 5/4
25
(5.20)
第5章
有限温度における気体-固体相転移
200
Solid
Solid
100
η
η
-3/2
-3/10
3
ρ fr σ
Pσ /ε
3
1
Gas
0.5
Gas
0
2
4
η
図 20
-6/5
6
0
0
8
2
4
η
kB T / ε
-6/5
6
8
kB T / ε
左側；スケーリングされた密度−温度状態相図。右側；スケーリングされた圧
力−温度状態相図。黒○は η = 0.1815、赤△は η = 0.6、青□は η = 1.0 のシミュレー
ションデータ。Hoover ら 21) によって導かれた古典系におけるスケーリング則を鎖線
で表し、前章の温度ゼロにおける結果を×でプロットした。
のスケーリング則が成り立っており、この関係式をプロットしたものが図 20 の鎖線であ
る。PIMC を用いた量子系の計算結果は、古典系の計算結果よりも下側に現れている。つ
まり、量子性によって固体相が安定化し、固体領域が広がっている事を示している。ま
た、この振る舞いは次のような見方も出来る。スケーリングの議論により、１粒子あたり
の自由エネルギーが f ≈ kB T f¯(T̃ , T̃ 3/n ṽ) で表され、２つ目の変数は η に依存しないた
め、η を変えると T̃ − ṽ 相図上で、T̃ 3/n ṽ = const. となるような曲線上を動く。そこで古
典極限として扱えるほど高温にあり、凝固密度よりも僅かに小さい密度の気体粒子を考え
る。温度を一定にしながら、粒子の質量を軽く（量子パラメターを大きく）すると、上図
の鎖線の共存曲線に沿って状態が変化してゆく。有効温度がある程度下がってくると量子
性が効いてきて、元々は気体であったものが固体になる振る舞いを示すことになる。ここ
からさらに有効温度を下げていった場合については、現状のサンプルデータではハッキリ
した事は言えないので、次項の考察で議論する。
5.4.2
区別出来る粒子の考察
有限温度の計算では、質量が軽く量子パラメターが大きな物質ほど、転移体積が大きく
なっていくという、量子性によって固化が促進される振る舞いが得られた。しかしながら
第 4 章の温度ゼロの計算結果では、量子パラメターが大きな物質ほど転移体積が小さくな
る振る舞いを示していた。この違いがどこから来るのか考えてみる。
第 3 章のスケーリングの議論で、IPL ポテンシャルが粒子間に働く系における特徴的な
物理量は、有効温度 T̃ = η −n/(n−2)
kB T
と有効体積 ṽ = η 3/(n−2)
V
N σ3
であることを示
し、そして量子パラメターの変化によって、系の状態に異なる２つの影響が現れる事を述
べた。ひとつ目は、量子パラメターが大きくなると、有効体積が大きくなり、固体相を不
26
第5章
有限温度における気体-固体相転移
安定化し、もうひとつは、量子パラメターが大きくなると、有効温度が小さくなり、固体
相を安定化する働きがある。温度ゼロにおいては、ひとつめの"有効体積の増加効果"のみ
が現れて、質量が軽く量子パラメターが大きな物質ほど固体相が不安定となり、転移体積
が小さくなっていた（図 19 の緑実線）。今回計算した有限温度における IPL ポテンシャ
ルが働く系では、有効体積の増加と有効温度の減少のふたつの効果が競合しているはずだ
が、今回計算したパラメター領域では"有効温度の減少効果"が勝ってしまい、質量が軽く
量子パラメターが大きな物質ほど転移体積が大きくなる振る舞いを示したのではないかと
思われる。その理由に、今回の振る舞いが剛体球系における振る舞いと似ている事が挙
げられる。剛体球の系は、今回もちいた IPL ポテンシャルで n → ∞ とした場合に対応
している。その際、有効体積は ṽ →
T̃ →
η −1 kB T
V
N σ3
となり量子パラメターに依存せず、有効温度は
と表され、剛体球の系では"有効温度の減少効果"のみが現れる。よって、
この有効温度の減少効果のみがある、剛体球系の計算 7, 8) と同様の振る舞いが、今回の計
算でも見られるため、計算したパラメター領域では有効温度の減少効果が勝っていると考
えられる。
では、どのようなパラメター領域であったならば、有限温度で"有効体積の増加効果"が
勝っている振る舞いを見ることが出来るのか議論してみる。温度ゼロのときには、エネル
ギースケールとしては零点エネルギーしかなく、有効体積の増加効果のみ現れたが、有限
温度にすると温度という新たなエネルギースケールが出て来たため、新たな有効温度の減
少効果が現れたと考えられる。そこで零点エネルギーと温度が同程度の大きさになった
~2 /(mv 2/3 ) ∼ kB T のときが、振る舞いが変わる境界だとすると、体積に対する境界線は
v
∼
σ3
(
kB T
)3/2
η 3/2
(5.21)
と表される。この境界線を T ∗ = kB T / = 1 の場合（赤鎖線）と T ∗ = 0.5（黒鎖線）の
場合をプロットしたのが図 19 の鎖線である。この線より上側が、"有効温度の減少効果"
が支配的な領域であり、下側が"有効体積の増加効果"が支配的な領域になる。すると、今
回計算した量子パラメターのほとんどの領域で、有効温度の減少効果が強く働いており、
量子性が強い物質ほど固化しやすくなっている振る舞いに説明がつく。T ∗ = 1, η = 1.0
のデータは境界線の下側に存在するが、境界に近いため有効体積の増加効果が顕著に見ら
れなかったと考えられる。温度 T ∗ をさらに下げると、この境界線の傾きは大きくなるた
め η < 1 の領域であっても、有効体積の増加効果が強く働いて、量子性が強くなると固化
しにくくなる振る舞いが見られる可能性がある。なお、温度ゼロにしてしまえば、境界線
の傾きは無限大となり、どんな量子パラメターの領域であっても、有効体積の増加効果が
強く働くことになる。
これまで有限温度と温度ゼロにおいて、量子性の強さが固化にどのような影響を及ぼす
のか議論してきた。次に状態相図上で、有限温度の計算結果と温度ゼロの計算結果がどう
つながるのか考える。図 20 では、温度ゼロの結果は×印で、有限温度の結果は○、△、□
でプロットされている。そこで相図上の PIMC を用いた有限温度の計算結果をそのまま
27
第5章
有限温度における気体-固体相転移
T → 0 に外挿すると、DMC で計算した温度ゼロの値とズレているように見える。このズ
レには粒子の統計性の違いが大きく関係していると考えられる。
ここまでの有限温度の計算では、粒子はすべて区別出来るものとして扱ってきた。しか
し、第 4 章でも述べたように、温度ゼロではボース粒子に対する計算がおこなわれてお
り、粒子の統計性が全く異なる。Feynman22) が行なった議論によると、統計性によらず
温度ゼロの基底状態の波動関数は節なしであるため、区別出来る粒子の基底エネルギーと
ボース粒子の基底エネルギーは同じになるはずである。そのため区別出来る粒子とボース
粒子の結果を温度ゼロにおいて繋げてもよいと考えられる。このときの振る舞いは、現状
得られているデータから推測すると、相図上の融解圧曲線に極小点が現れてくる可能性が
あり、融解圧曲線に極小点が現れるというのは、特異なことであり非常に興味深い。
しかしながら、低エネルギー励起について考えると奇妙な点が残る。通常のシミュレー
ションにおいては、サンプルに使う粒子数によってエネルギーなどに変化が現れる場合が
ある（有限サイズ効果と言われている）。本来求めるべき値は、粒子数 N に依存しない
N → ∞ の値である。本研究でも、PIMC または DMC のシミュレーションの際は粒子
数を N = 32 ∼ 256 などと変化させ、N = 108 程度を使えば有限サイズ効果は無視出来
ることを確認している。つまり相図上のデータはあらかじめ N → ∞ の処理がされてい
ると見なせる。この状態で温度ゼロに値を外挿するという行為は、粒子数の極限操作をし
た後に温度の極限操作をすることになる。しかしながら、温度ゼロでの計算は温度の極限
操作がなされた後に粒子数の極限操作がされて、順序が逆になっている。この極限操作の
順序は、ボース粒子については問題にならないが、区別出来る粒子にとっては重要な問題
となる。温度ゼロ近傍の場合を考えると、区別出来る粒子の低エネルギー励起は、粒子数
の極限操作を先にすると連続の励起状態が出来て、エントロピーに負の対数発散などを起
こしてしまう可能性がある。よって温度ゼロにおけるボース粒子の結果と有限温度におけ
る区別出来る粒子の結果を比較するのは不適切であるとも考えられる。なお Boninsegni
らによっても、区別出来る粒子とボース粒子の振る舞いが議論されている 23) 。
ここでまで議論してきたように、基底エネルギーは一致するとしても、どんなに低温で
あっても有限温度で N → ∞ とすると、区別出来るか否かで自由エネルギーまたはエン
トロピーに違いがでる可能性があり、低温極限において PIMC の区別出来る粒子の計算
結果と、温度ゼロにおける DMC のボース粒子の結果を、そのまま比較しても良いのか
はっきりしない。そのため次の節からは、有限温度においても粒子をボース統計に従うも
のとして扱った計算をおこない、その結果について議論してゆく。
28
第5章
5.5
ボース粒子の場合
5.5.1
ボース粒子の結果
有限温度における気体-固体相転移
経路積分モンテカルロ法においてボース粒子を扱うには、ふたつの粒子の経路を虚
時間に沿って交換する必要がある。この交換のプロセスが入るためボース粒子の計
算には、区別出来る粒子の数倍の計算時間が必要となる。今回は計算時間の都合上、
[T ∗ = 0.5, η = 0.6] の場合と [T ∗ = 0.5, η = 0.1815] の場合について詳しくシミュレー
ションをおこなった。図 21 は T ∗ = 0.5, η = 0.6 の気体状態において、PIMC シミュ
レーション中のエネルギーと圧力の推移をプロットしたものである。実際のシミュレー
ションステップ数は横軸の数値に 3000 を掛けたもので、最初の約 1 万ステップは粒子を
区別出来るものとしてサンプリングし、その後にボース粒子の交換効果を取り入れてい
る。ρ∗ = (N σ 3 )/V = 0.33 と ρ∗ = 0.25 の気体状態の場合において、内部エネルギーは
粒子の交換効果をいれると、区別出来る粒子のエネルギーよりも大きく下がる。しかし圧
力は、ρ∗ = 0.25 の場合はエネルギーと同様に交換効果をいれると下がるが、ρ∗ = 0.33
の場合は交換効果を入れても変化は見られなかった。この振る舞いは次のように理解さ
れる。圧力は P = −∂F/∂V で表されるため、ρ∗ = 0.33 近傍のボース粒子の自由エネル
ギーは、区別出来る粒子の自由エネルギーを平行移動させたものになっている。一方で、
ρ∗ = 0.25 近傍でのボース粒子の自由エネルギーの変化は、区別出来る粒子の変化よりも
緩やかになっていることを意味している。
このように、粒子の交換効果を取り入れたことによる熱力学量の変化は密度によって
∗
∗
異なっている。区別出来る粒子 Edist
とボース粒子 Ebose
の気体状態の内部エネルギー差
は、密度に対して図 22 の左図のような振る舞いを見せる。理想気体と見なせるほど希薄
だと、区別出来る粒子とボース粒子のエネルギー差はほとんどなく、そこから密度が大き
2.78
1.24
5.2
2.74
1.22
3.2
4.8
0
50
100
1.2
3
2.72
0
150
100
200
Monte Carlo steps
Monte Carlo steps
図 21 PIMC シミュレーション中の気体状態の内部エネルギー（赤実線）と圧力（青
実線）の推移。エネルギーは左縦軸、圧力は右縦軸に対応する。左側は T ∗ = 0.5, η =
0.6, ρ∗ = 0.33 の場合。右側は T ∗ = 0.5, η = 0.6, ρ∗ = 0.25 の場合。
29
pressure P*
5
energy E*
energy E*
2.76
pressure P*
3.4
第5章
2
entropy S / N kB
E*dist - E*bose
0.4
0.2
0
有限温度における気体-固体相転移
0.2
1
0
0
0.4
density ρ*
0.2
0.4
density ρ*
∗
とボース粒子
図 22 T ∗ = 0.5, η = 0.6 の気体相の場合。左側；区別出来る粒子 Edist
∗
の内部エネルギー差。右側；区別出来る粒子（黒実線）とボース粒子（赤実線）
Ebose
のエントロピー。
くなるとボース粒子のエネルギーが減少し始める。密度がある程度大きくなると、粒子の
交換効果によるエネルギーの低下は頭打ちになり、エネルギーの変化量は密度が変わって
もほぼ一定となった。
それぞれの密度において、粒子の交換効果を入れた際のエネルギーと圧力をシミュレー
ションで求めることが出来れば、5.2 節で述べた熱力学積分を行うことによって、転移点
を求めるのに必要な自由エネルギーが計算できる。その結果、T ∗ = 0.5, η = 0.6 の気体
相では、ボース粒子にしたことによって内部エネルギーが下がり、最終的には自由エネル
ギーも下がることが計算で分かった。自由エネルギーを計算した後に、F = E − T S の関
係式から１粒子あたりのエントロピーを各密度で見積もったものが図 22 の右図である。
黒実線は区別出来る粒子のエントロピーで、赤実線はボース粒子のエントロピーである。
内部エネルギーのときと同様、十分希薄であればどちらの場合でも、理想気体と見なせて
エントロピーの大きさは同じになる。そして密度が大きくなると、ボース粒子のエントロ
ピーは区別出来る粒子のエントロピーよりも低くなった。これは区別出来る粒子よりも
ボース粒子のほうが、状態数が少ないためエントロピーが低くなっていると理解できる。
なお粒子の交換効果を入れた場合のエントロピーの減少は、自由エネルギーを増加させる
方向に働くが、内部エネルギーの減少による自由エネルギーの減少の方が大きいため、最
終的に自由エネルギーは減少することになる。
Fdist − Fbose = Edist − Ebose − T (Sdist − Sbose )
(5.22)
以上は T ∗ = 0.5, η = 0.6 の気体相に、粒子の交換効果を入れてボース粒子として扱っ
た場合の計算結果であったが、固体相においてもボース粒子の計算を行なったところ、内
部エネルギーと圧力は区別出来る粒子と同じ値になった。つまり、粒子の交換効果は気体
相に対しては、自由エネルギーを低下させる働きがあるが、固体相に対しては影響を及ぼ
30
第5章
有限温度における気体-固体相転移
T ∗ = 0.5, η = 0.6 における転移密度と融解圧
表4
凝固密度 ρ∗fr
融解密度 ρ∗ml
融解圧 Pc
ボース粒子
0.472 ± 0.011
0.482 ± 0.011
7.85 ± 0.30
区別出来る粒子
0.437 ± 0.019
0.440 ± 0.022
6.15 ± 0.83
80
Solid
3
Pσ /ε
60
η
-3/2
40
20
Gas
0
0
2
η
-6/5
4
kB T / ε
図 23 ボース粒子系における圧力−温度相図；黒●は T ∗ = 0.5, η = 0.1815 の赤▲
は T ∗ = 0.5 η = 0.6 のボース粒子の計算結果。×印は DMC による温度ゼロのボース
粒子の計算結果である。比較するために区別出来る粒子の場合をオープンシンボルで
プロットしてある。
さない。また、この時の転移密度や融解圧は表 4 のようになり、量子統計性を考えてやる
事で気体相が安定化し、より高密度で固化することが分かった。
次に、温度は同じでも量子パラメターが小さな T ∗ = 0.5, η = 0.1815 の場合について計
算した。その結果、このパラメターでは固体相、気体相のどちらにおいても、粒子の交換
効果を入れてもエネルギーや圧力に変化は見られなかった。つまり、ボース粒子の転移密
度や融解圧は、区別出来る粒子の場合とまったく同じになる。これらの結果を区別出来る
粒子の計算で作成した、スケールされた圧力−温度相図上にプロットしたものが図 23 で
ある。ボース粒子の計算は、T ∗ = 0.5, η = 0.1815、T ∗ = 0.5, η = 0.6 と温度ゼロの計算
結果の３点しかないが、これらの点を繋げるとボース粒子の融解圧曲線には極小点がある
ことが分かる。極小点 P̃min がある場合、系の圧力が P̃min < P̃ < P̃ (T = 0) となるよう
に設定し温度を上げてゆくと、最初は気体であったものが固体になり、さらに温度を上げ
るとまた気体に戻るというリエントラントな振る舞いを示す。この極小点が何を意味する
のかについては次項で議論するが、今回のボース粒子系の計算で、"量子統計性を考えて
やると、区別出来る粒子に比べボース粒子の気体相は安定化し、また圧力−温度相図上の
融解圧曲線に極小点が現れる"ことが分かった。
31
第5章
5.5.2
有限温度における気体-固体相転移
ボース粒子の考察
区別出来る粒子の計算では、極小点があるのか定かではなかったが、有限温度において
もボース粒子の計算を行うことで、融解圧曲線に極小点があることが分かった。一般に融
解圧曲線上ではクラジウス・クラペイロンの式
dP
Sgas − Ssolid
=
dT
Vgas − Vsolid
(5.23)
が成り立つことが知られている。右辺は、気体から固体への凝固に伴うエントロピー変
化 ∆S と体積変化 ∆V の比で表されている。融解圧曲線上では、気体相の体積の方が
固体相よりも大きいため常に ∆V > 0 となっている。通常の物質では気体相のエントロ
ピーの方が大きいため ∆S > 0 であり、融解圧曲線は温度の増加関数である。しかし、融
解圧曲線に極小点がある場合、極小点より高温側では dP/dT > 0 となり、低温側では
dP/dT < 0 となっているはずである。この低温領域では、∆V > 0 であるため ∆S < 0
となって、固体相のエントロピーの方が気体相のエントロピーよりも大きくなる振る舞い
を示す。なお、温度ゼロ近傍では熱力学第三法則に従うように、∆S ' 0 となり融解圧曲
線の傾きはゼロになっているはずである。
現実の物質でも融解圧曲線に極小点が現れ、固体相のエントロピーの方が気体相のエン
トロピーよりも大きくなるものがある。Straty ら 24) によって、ボース粒子であるヘリウ
ム４の融解圧曲線が詳細に調べられ、T = 0.775 K 付近に極小点があることが実験的に分
かっている。また統計性は違うが、フェルミ粒子であるヘリウム３においてより顕著な極
小点があることが分かっている。（低温の液体相ではフェルミ粒子であることが効いてエ
ントロピーが小さくなっており、また固体相は局在スピン S = 1/2 の集まりと見なせる
のだが、スピン間に働く相互作用が小さいため、低温ではスピンの自由度が残ってしまい
N ln 2 という大きなエントロピーを持つ。そのためヘリウム３では固体相のエントロピー
の方が大きくなり極小点が現れている。）
今回計算した IPL ポテンシャルが働く系でも、融解圧曲線に極小点が現れたが、この系
でなぜ固体相のエントロピーの方が大きくなるのか考えてみる。低温における固体のエネ
ルギー励起は、フォノンのモードがほとんどである。そのため、このフォノンに関連する
物理量を計算することによって、温度ゼロ近傍の低温においてエントロピーが増加する原
因を考察する。
今回は、温度ゼロにおけるバルクモジュラス K とシアモジュラス µ を計算し、検証す
ることにした。単一原子の結晶では、フォノンは１つの縦波と２つの横波のモードから成
り立っており、縦波の速度 vl はバルクモジュラス、横波の速度 vt はシアモジュラスを用
いて
√
vl =
√
K
,
ρ
vt =
µ
ρ
(5.24)
と表される。バルクモジュラス K は系全体に一様な力を加えたときの弾性率であり、体
32
第5章
有限温度における気体-固体相転移
積を変化させたときの圧力変化
K = −V
dP
dP
=ρ
dV
dρ
(5.25)
で計算出来る。またシアモジュラス µ は、ある方向にひずみを加えたときのの弾性率であ
り、とくに z 0 = z + x の変位を与えたときのエネルギー変化 E 0 = E + 12 µ2 から求め
る事が出来る。バルクモジュラスとシアモジュラスも、スケーリングの議論によって
K̃ = η −3/2
K
,
σ3
µ̃ = η −3/2
µ
σ3
(5.26)
となるはずで、図 24 は異なる量子パラメター η でこのバルクモジュラスとシアモジュラ
スを温度ゼロにおいて計算した結果である。バルクモジュラスとシアモジュラスどちらと
も、スケールされた体積によって同一曲線を描く振る舞いをしている。このなかで、固体
相のバルク・シアどちらのモジュララスも転移点に向かって急激に減少している。特に固
体相のシアモジュラスは、気体相のバルクモジュラスよりもかなり小さくなっているのが
特徴的である。またシアモジュラスは転移点近傍で、ゼロに近い値をとっているため、こ
の相転移は２次転移に近いものであると考えられる。式 (5.24) により、小さくなったシ
アモジュラスは、フォノンの横波速度が遅くなったことを意味している。一般にフォノン
の速度が遅くなると、状態密度が大きくなり、エントロピーの増加を引き起こすため、こ
の横波のソフト化が固体相のエントロピーを大きくしているのではないかと考えられる。
今回の検証では、簡単のために温度ゼロにおける計算をおこなったが、温度ゼロ近傍の低
温でも同様な横波のソフト化が起こり、固体相のエントロピーが気体相よりも大きくなる
可能性がある。今後は、極小点より低温の領域で同様の計算をするなどして、固体相のエ
ントロピーの方が大きくなる原因を探って行きたい。
Gas
3
µ /(ε/σ )
Solid
η
200
−3/2
0.60
η
−3/2
3
0.30
K/(ε/σ ) , η
0.1815
100
0
1
1.5
η
3/10
v/σ
3
2
図 24 温度ゼロにおける、バルクモジュラス K とシアモジュラス µ。オープンシンボ
ルは固体相のバルクモジュラス、クローズドシンボルは気体相のバルクモジュラスを表
し、十字シンボルなどは固体相のシアモジュラスを表す。なお、凝固・融解体積を点線
で示してある。
33
第6章
まとめ
本研究では、まずスケーリングの議論によって零点運動には、相反する２つの効果があ
ることが明らかとなった。ひとつ目は、温度および密度を一定にして質量を軽くすると、
有効体積が大きくなることで固化を妨げる効果である。もうひとつは、質量が軽くなると
有効温度が減少する効果で、これは固化を促進する。一般には、両者の競合により、固化
に対する零点運動の効果が決まる。
また数値計算の結果、温度ゼロでは零点運動による有効体積の増加効果のみが見られ、
質量が軽い物質ほど高密度で固化し、量子効果によって固化は抑制されることが確かめら
れた。一方で、今回計算した有限温度のパラメター領域は、熱エネルギーが零点エネル
ギーより大きい有限温度領域であったため、零点エネルギーを大きくすると、有効温度の
減少効果が勝り、固化が促進されることが明らかとなった。もし、同様な計算をさらに低
温で行なった場合に、零点エネルギーを大きくすると、有効温度の減少効果が支配的な領
域から、有効体積の増加が支配的な領域にクロスオーバーする事が期待される。
有限温度においてボース粒子の統計性を取り入れると、気体相の自由エネルギーは区別
出来る粒子の場合よりも低下し、気体相が安定化することが明らかとなった。またこの際
に、圧力−温度相図において、融解圧曲線に弱い極小が現れることを見出した。これは固
体相が気体相に比べて高エントロピーの状態にある事を示唆する。この振る舞いを理解す
る為に、固体相におけるずり弾性を調べたところ、相境界近傍においてこれがソフト化し
ていることが分かった。これにより、固体相の横波フォノンの速度が減少し、エントロ
ピーの相対的な増加の原因となっている可能性が考えられる。
謝辞
本研究に際し、テーマの示唆から論文執筆に至るまで指導教員である平島大教授に非常
にお世話になりました。私の研究スピードが遅いため、いろいろご迷惑をおかけしてしま
いましたが、最後まで指導していただけた事に深く感謝いたします。私の研究は理論的な
ものであるため、実験の知識が乏しくまた先行研究についても多くをしらないため、名古
屋大学、極低温量子物性研究室の和田教授、松下助教、檜枝講師には様々なことを教えて
いただきました。ありがとうございました。経路積分モンテカルロ法のプログラムコード
は、Konkuk University（韓国）の Yongkyung Kwon 教授が研究で使われているコード
を、使用させていただきました。本研究でも使用し、改良することを快諾してくださった
ことを、心より感謝いたします。最後に、絶えず様々な分野の知識と刺激を与えてくださ
り、大学院生活を楽しく有意義なものにしてくださった物性理論研究室の先生方及び院生
の皆様には大変お世話になりました。
34
付録 A
付録 A
Monte Carlo integration
Monte Carlo integration
実軸上の区間 [a, b] 上の滑らかな関数 f の積分
∫
b
I=
(A.1)
dx f (x)
a
を計算したいとする。この問題に関する標準的な数値解法は、等間隔の xi の集合上で（た
だし点が等間隔ではないガウス積分を除く）関数の値を求め、和
b−a∑
I=
wi f (xi )
N i=1
N
(A.2)
を計算するということである。重み wi は f に依存しない。この重みにより、この方法の
精度が決まる。このような方法は、通常被積分関数の多項式近似に基づいており、精度
σ は積分点の間隔 h のベキで表せる（σ ∼ hk ∼ n−k で k は正の整数）。モンテカルロ
(MC) 積分でも式 (A.2) を使うが、wi はすべて 1 に等しくここでは xi をランダムに選ぶ。
もしランダムな座標 xi が [a, b] に一様に分布していて、N が十分大きいなら、式 (A.2)
の和は、厳密な積分値に近くなる。そのときの分散は
*(
σ2 =
b−a∑
fi
N i=1
N
)2 +
(*
−
b−a∑
fi
N i=1
N
+)2
(A.3)
になる。h i は、ランダム座標 xi のすべての可能な列の実現について平均を表す。右辺の
最後の項で平均をとると、区間 [0, 1] 上での関数 f の平均 f¯ の 2 乗になる。右辺第 1 項の
和を、i = j の場合と i 6= j の場合に分け、少し計算すると、
σ2 =
(b − a)2 ¯2 ¯2
(f − f )
N
(A.4)
が得られる。ここで、上線は区間 [0, 1] 上での関数の平均を表す。誤差は、[0, 1] 上での
√
f の分散に比例することがわかる。σ ∼ 1/ N という事実は、中心極限定理からも予想
できる。このスケーリングは、標準的な求積法に比べ、明らかに悪い。この事を標準的
な求積法の場合と比較してみる。標準的な求積法では xi に関して等間隔な値を用いて、
N −k (k ≥ 1) の誤差になる。しかし、d 次積分に対しては、もし積分領域が 1 辺 L の超立
方体だとすると、この中には N = (L/h)d 個の点があり、したがって結果の誤差は N −k/d
でスケールされる。また MC 積分の誤差は次元に依存しない（中心極限定理が次元に依ら
ない）ため、O(N −1/2 ) のままである。よってこの誤差を標準的な方法の誤差と比較する
と、d > 2k のときにはオーダー k のアルゴリズムよりも、MC 積分の方が効率的である。
このように高次元で MC 積分が優れている理由は、ある意味で規則的な格子点よりも
ランダム分布の方が、一様性が高いためである。例えば、積分領域の中の直方体区域を考
えてみよう。積分点が一様に分布するということは、この直方体区域中の点の数が、体積
35
付録 A
Monte Carlo integration
に近似的に比例することを意味する。直方体の辺が標準的な積分法で使うような格子点の
軸に平行に向くように置くとする。このとき区域の大きさを増やすに従って、区域の面が
積分点の列を通るたびに、内部の点の数が階段状に増加することは明らかである。このよ
うな格子点数のステップは極端に起こりにくいで、この観点でランダム分布はより一様で
ある。
さらにモンテカルロ積分の誤差を減らすため、いくつかの方法が考案されている。その
一つに、重み付きサンプリングモンテカルロと呼ばれるものがある。実際の多くの場合で
は、積分区域内の異なる領域の積分への寄与は大きく異なる。単純な MC 積分は関数を
一様にサンプリングするので、もし被積分関数への主要な寄与が、積分領域内の小さな領
域から主に来るなら、その領域での関数のサンプリングには少しの MC 点しか使われて
おらず、大きな統計誤差が出る。この効果は式 (A.4) で、区域上での f の大きな分散とし
て現れる。誤差へのこのような寄与を減らすために、f の絶対値が大きくなる領域で、点
の密度を高くする。つまり、領域 [a, b] 上で関数 ρ(x) が f の形に近い、すなわち f /ρ が
ほぼ定数であるものを用いる。ここで、ρ は次式のように規格化されている。
∫
b
(A.5)
dxρ(x) = 1
a
さらに f の積分を
∫
∫
b
b
dxf (x) =
a
a
[
f (x)
dxρ(x)
ρ(x)
]
(A.6)
と書くと、括弧の中の関数は十分変動がな少くなる（ρ として、f に形が似たものを選ん
だため）。この関数の前に付いた重み ρ(x) は、分布 ρ(x) に従うランダムな点 xi を選べ
ば、積分に含めることができ、そうするとモンテカルロによる積分結果は、式 (A.2) で与
えられ、結果の誤差がかなり減るのは、この場合の式 (A.4) の分散を求めれば分かる。

(∫ [
)2 
]2
]
∫ b[


b
1
f (x)
f (x)
2
σ =
ρ(x)dx −
ρ(x)dx

N  a ρ(x)
ρ(x)
a
(A.7)
実際に f /ρ が近似的に定数になるように ρ を選ぶことができるなら、式 (A.5) の結果、
中括弧の中の式は「単純なサンプリング」の場合よりも、はるかに小さくなる。これが重
み付きサンプリングモンテカルロである。ここで注目すべきは、この方法を用いる場合は
関数 f についてあらかじめ知識が必要であるという事である。このようなあらかじめの
知識を必要としない適応型モンテカルロ法も考案されている。また一般的に MC 積分で
は、乱数発生の相関に敏感ではなく、相関により、点 xi が生成される順番は影響される
が、その和は影響されない。実際、（疑似）ランダム性を満たそうとしなくても、高次元
の区域を一様に満たして、積分に適しているような人工的な数列を発生させることができ
る。そのような方法を「準モンテカルロ法」と呼んでいる。
36
付録 B
付録 B
B.0.3
Monte Carlo simulations
Monte Carlo simulations
Importance sampling
カノニカルアンサンブルでの、ある物理量 Q の平均を求めることを考えよう。平均は
∑
hQi =
µ
∑
Qµ e−βEµ
µ
e−βEµ
(B.1)
で与えられ、ボルツマン因子 exp[−βEµ ] に従って状態に重みを付けて、全ての状態 µ
について和（積分）を取らなくてはいけない。式 (B.1) は系のサイズがとても小さい場
合には扱いやすいが、系のサイズが大きい場合には、ある部分領域に対して平均をと
るのが、実際の計算時においては現実的である。この様な系のサイズが大きいとき、モ
ンテカルロ法はある確率分布 pµ に従ってこの部分領域を選択する。選択した部分領域
M {µ1 , · · · , µM } において hQi の推定量 QM は
∑M
QM =
−1 −βEµi
i=1 Qµi pµi e
∑M −1 −βEµ
j
j=1 pµj e
(B.2)
で与えられ、M の数が増加するとより正確になり、hQi の良い推定量になる（M → ∞
のとき QM = hQi となる）。しかし、部分領域 M について平均を取るようにしても、そ
の領域全てを等確率でサンプリングすると、重要でない部分と重要な部分が同じ頻度で計
算され、効率が悪い。そのため、重要な部分は重点的にサンプリングされ、重要でない部
分は多くサンプリングされないようにすれば効率が良くなるはずである。このような考え
方は Importance sampling（重み付きサンプリング）と呼ばれ、多くのシミュレーション
アルゴリズムに組み込まれている。
Importance sampling は次のように行われる。M の全ての状態が等確率で選ばれるの
ではなく、ボルツマン因子に比例した pµ = Z −1 e−βEµ の確率で、ある状態 µ が選択され
るようにする。このとき hQi の推定量は式 (B.2) より
QM
M
1 ∑
=
Qµi
M i=1
(B.3)
で与えられる。ここで注目すべきは、推定量 QM においてボルツマン因子がキャンセル
して消えていることである。こうすることにより、状態の数が少ないところで系がシミュ
レーションの時間の大半を過ごしているときに役に立つ。（例えば、低温において状態が
基底状態に落ち着いているときなど。）ここで残る問題は、どのようにして正確にボルツ
マン確率によって状態を生成するかである。この問題を解決するには、統計的に独立な状
態を生成させることはやめて、"Markov process" によって状態を生成するという方法が
ある。
37
付録 B
B.0.4
Monte Carlo simulations
Markov processes
多くの MC 計算において、マルコフ過程 (markov process) を用いて状態を生成してい
る。マルコフ過程とは、統計的に独立な状態を生成するのではなく、前の状態に依存した
確率分布で新しい状態を生成するのである。例えば、系がある状態 µ であるときに、次の
新しい状態 ν を生成する（系が新しい状態に遷移する）ことを考えよう。ただし、「初期
状態が µ であるときはいつでも、ある特定の状態を必ず生成する」ようなことは起こらな
い。（ボルツマン確率に従って、状態を生成したいのだから。
）状態 µ から状態 ν を生成す
る確率は、遷移確率T (µ → ν) と呼ばれ、マルコフ過程に従って状態を生成する際は、以
下の二つの条件を満たしていなければならない。
(i) 遷移確率 T (µ → ν) はシミュレーション時間中に変化しない。
(ii) 遷移確率 T (µ → ν) は、系がどのような状態を経由して遷移したかに依らず、注目
した二つの状態 µ と ν にのみ依存する。
この二つの条件は、状態 µ が与えられたときは、マルコフ過程によって状態 ν が生成され
る確率はいつでも同じであるということを意味している。またマルコフ過程において、遷
移確率 T (µ → ν) は必ず次の状態を生成しなければならないので以下のように規格化され
ている。
∑
T (µ → ν) = 1
(B.4)
ν
しかしここで注意しなければいけないのは、遷移確率 T (µ → µ) は必ずしもゼロである必
要はない事である。つまり、マルコフ過程において系の状態がある状態 µ に留まることは
あり得るのである。
MC シミュレーションにおいて、マルコフ過程は繰り返し用いられ、状態のマルコフ連
鎖をつくる。状態 µ から始まり、新しい状態 ν をつくり、その状態をもとにまた新しい状
態 λ をつくることを繰り返す。さらにシミュレーションでは正確にボルツマン確率で状
態を生成できるように、ボルツマン因子に比例した分布を持つ系の状態のマルコフ連鎖を
つくり、この分布が連鎖内の位置や初期状態に無関係なようにしなければらない。（つま
り、選ばれた初期条件を忘れるだけの時間を必要とし、平衡状態になっている。）ある条
件下では、マルコフ連鎖はこのような不変な分布を、少なくとも長期的には持つ事ができ
る。それには上に挙げたマルコフ過程であるための条件の他に、エルゴート性や詳細釣り
合いの条件を満たしていなければならない。
B.0.5
Ergodicity
エルゴート性 (ergodicity) の条件とは、「アンサンブルに含めたい全ての状態は、マル
コフ過程によって他の状態から、有限のステップで到達することができる」ことである。
この条件は、正確にボルツマン確率に従って状態を生成できる平衡状態に到達するために
必要なのである。例えば、ある状態 ν がボルツマン分布にしたがって、ゼロではない確率
38
付録 B
Monte Carlo simulations
pν で現れるという状況（平衡状態）を考えたいとする。このとき、その他の状態 µ から
どれだけ時間をかけても状態 ν にたどり着けないとすると、その系では ν が現れず平衡
状態にいたらない。（マルコフ連鎖の中で状態 ν になる確率はゼロになり、有限の出現確
率であって欲しい pν がゼロになる。）
つまりエルゴート性の条件を満たすマルコフ過程では、ある状態からある状態への遷移
確率がゼロであっても良いが、選んだ二つの状態の間に少なくとも一つはゼロではない遷
移確率を持った経路がなくてはならない。多くの MC アルゴリズムでは、このエルゴー
ト性の条件は常に満たされており、あまり気にする必要はない。
B.0.6
Detailed balance
次にもう一つの詳細釣り合い (detailed balance) の条件について考える。この条件は系
が平衡状態に達した後、ボルツマン確率分布で状態を生成できることを保証するために必
要である。では、ボルツマン分布 exp(−βH) に基づき、状態のマルコフ連鎖を生成する
ときの平衡状態について考えよう。与えられた定常分布 pµ （いまの場合はボルツマン分
布）を与えるような遷移確率 T (µ → ν) を見つける必要がある。状態を記憶することに
なるコンピュータのメモリは有限なので、可能な状態の数 N は有限であると仮定する。
よって状態の出現する確率を表すベクトル p
~ には N 個の成分を持つが、状態間の遷移確
率を表す行列 T̂ は N 2 個の要素を持つ。このことは、この問題において多くの異なる解
があり、解を見つけるときには自由度があるという事を意味する。
ある時間、つまりマルコフステップ t において、状態 µ が出現する確率を与える関数を
pµ (t) とする。あるステップから次へ行くときこの関数の変化は、次の２つの過程に支配
される。① 時間 t での状態 µ から t + 1 での別の状態 ν に行く過程。これにより pµ は減
少する。② 時間 t でのある状態 ν から時間 t + 1 で µ に行く過程。これにより pµ は増加
する。この仕組みは、以下の方程式
pµ (t + 1) − pµ (t) = −
∑
T (µ → ν)pµ (t) +
∑
ν
T (ν → µ)pν (t)
(B.5)
ν
と書ける。このとき定常状態であるためには、pµ (t + 1) = pµ (t) が要求される。した
がって
∑
T (µ → ν)pµ (t) =
∑
ν
T (ν → µ)pν (t) (B.6)
ν
が得られる。となる。また式 (B.4) を用いて上式 (B.6) を簡単にすると、
pµ (t) =
∑
T (ν → µ)pν (t)
(B.7)
ν
となり、この式を満たす遷移確率によって、確率分布 pµ はマルコフ過程のなかで平衡状
態になっている。（以降は定常状態を考えるので、p の t 依存性を省略する。）ここで注意
しなければならないのは、解に自由度があるため式 (B.6) の一般解を求めるのが難しいこ
とである。ところが、全ての状態の組 µ, ν に対して、
T (µ → ν)pµ = T (ν → µ)pν
39
(B.8)
付録 B
Monte Carlo simulations
の特殊解があり、この解を詳細釣り合いの解と呼んでいる。
以上より、式 (B.8) を満たす遷移確率の組を選ぶことによってマルコフ過程における、
ある状態の出現する確率 pµ を作ることがができるようになった。（平衡状態においてボ
ルツマン分布によって状態を出現させられる。）また式 (B.8) を変形させて、ボルツマン
分布に従う平衡状態の場合、
T (µ → ν)
pν
=
= e−β(Eν −Eµ )
T (ν → µ)
pµ
(B.9)
と書くこともでき、式 (B.9) と式 (B.4) は遷移確率を選ぶ際の制約となる。もしこれらの
条件とエルゴート性の条件を満たしていれば、マルコフ過程の平衡分布はボルツマン分布
になる。条件に合った遷移確率を選ぶ方法は、それぞれのアルゴリズムの肝心な部分であ
り、シミュレーションの精度が良いものや、所要時間が短いアルゴリズムを用いるとよ
い。しかし、すべての条件を満たしているアルゴリズムを用いてシミュレーションを行っ
たとしても、平衡状態になるまでの時間がどれだけ必要かは、それぞれのシミュレーショ
ン上で検証しなければならないことに注意しよう。
B.0.7
Acceptance ratios
以上のことより、式 (B.9) と式 (B.4) を満たす遷移確率を選ぶことにより、マルコフ連
鎖を用いてボルツマン分布に従って新しい状態を生成できるようになったが、このままの
形の遷移確率ではあまり実用的ではない。そこで採択確率 (acceptance ratios) と呼ばれ
るものを導入する。マルコフ過程について述べた際、遷移確率 T (µ → µ) はゼロではなく
ても良く、状態が変わらない場合も許されていたことを思い出そう。この状態が変わらな
い場合は、式 (B.9) によりどのような遷移確率の形を選んだとしても、詳細釣り合いの条
件を満たしていることが分かる。0 < T (µ → µ) < 1 になるよう注意しながら、式 (B.4)
を満たすように T (µ → ν) を調整する自由度がある。さらに T (µ → ν) を調整すること
ができたら、式 (B.9) を使って同時に T (ν → µ) を調整することができる。このことから
遷移確率次のように書き換えると便利である。
T (µ → ν) = g(µ → ν)A(µ → ν)
(B.10)
g(µ → ν) は試行ステップ確率 (selection probability) と呼ばれ、状態 µ が与えられたと
き、新しい状態 ν を作り出す確率であり、A(µ → ν) は採択確率 (acceptance ratio) であ
る。新しい状態が作られたとき、その状態を採択し状態を変化させるかどうかを、採択確
率で決める。このとき、採択確率は 0 ≤ A(µ → ν) ≤ 1 でなければならない。もし採択確
率をすべての遷移についてゼロにすると、すべての新しい状態は棄却され T (µ → µ) = 1
になり状態は常に変化しない。（つまり実際にシミュレーションを行うときは、採択確率
をゼロにしてはいけない。）また式 (B.9) を書き換えると
g(µ → ν)A(µ → ν)
T (µ → ν)
=
= e−β(Eν −Eµ )
T (ν → µ)
g(ν → µ)A(ν → µ)
40
(B.11)
付録 B
Monte Carlo simulations
となり、A(µ → ν)/A(ν → µ) はゼロから無限大のどれかの値を取ることができ、
g(µ → ν) に関しては我々はどのような値でも選ぶことができる。また式 (B.4) による制
約は、系がマルコフ過程によってできた状態に最終的に落ち着いていれば満たされてい
る。このようにデザインされたアルゴリズムを用いることによって、系は平衡状態に達
し、正しいボルツマン確率に従って状態を効率的に生成することができる。ただし注意す
べき点として、アルゴリズムをデザインするときに、採択確率が小さくなるようにしてし
まうと状態があまり遷移しなくなり、シミュレーションの時間の無駄使いになってしま
う。好ましいのは、状態空間を素早く動き、異なる状態を広くサンプリングでき、効率的
なアルゴリズムを用いることである。この問題の簡単な解決方法は、採択確率をできるだ
け 1 に近づけることであり、最も良い方法は遷移確率 T (µ → ν) の状態 µ と ν に対する
依存性を、できるだけ g(µ → ν) に組み込み、採択確率 A(µ → ν) に入れないようにする
ことである。
41
付録 C
付録 C
C.0.8
Metropolis algorithm
Metropolis algorithm
Metropolis algorithm
Metropolis algorithm は 1953 年に Nicolas Metropolis らによって、剛体球気体のシ
ミュレーションで導入された。Metropolis algorithm はモンテカルロシミュレーション
のなかでもっとも多く用いられている。ここでは、Metropolis algorithm を Ising model
で説明する際、single-spin-flip dynamics と呼ばれる、スピンを一つずつフリップさせる
手法をとることにする。（注意：Metropolis algorithm は single-spin-flip dynamics であ
ることを意味しない。Metropolis algorithm を特徴づけるのは、採択確率の選び方であ
り、このアルゴリズムでも複数のスピンをフリップさせることがある。）
Metropolis algorithm では、試行ステップ確率 g(µ → ν) をすべての可能な新しい状態
ν に対して等しく設定し、その他はゼロとする。例えば、シミュレーションを行う系に N
個のスピンがあるとすると、状態 µ から到達可能な新しい状態 ν は N 個あり、それらの
状態はすべて等しい試行ステップ確率を持つということである。その場合の試行ステップ
確率は、
g(µ → ν) =
1
N
(C.1)
となる。さらにこのときの詳細釣り合いの条件式 (B.9) は、
T (µ → ν)
g(µ → ν)A(µ → ν)
A(µ → ν)
=
=
= e−β(Eν −Eµ )
T (ν → µ)
g(ν → µ)A(ν → µ)
A(ν → µ)
(C.2)
と表され、この式を満たすように採択確率 A(µ → ν) を選ばなくてはいけない。もっと
も簡単な選び方としては、A(µ → ν) = A0 e− 2 β(Eν −Eµ ) があるが、この選び方はエネル
1
ギーの差分 Eν − Eµ に対して採択確率が指数関数的に小さくなるので、シミュレーショ
ンのなかで多くの時間を状態の棄却に費やしてしまい非効率である。このように採択確
率をある関数形にしてしまいがちだが、式 (C.2) は採択確率のペアに対する条件である
ため、条件さえ満たしていればどのような形でもよいのである。さらにモンテカルロシ
ミュレーションでは、採択確率ができるだけ 1 に近づくようにアルゴリズムをデザイン
するのがよい。では状態 µ のほうが状態 ν のエネルギーよりも低い場合を考えて見よう
(Eµ < Eν )。このとき採択確率のペアのうち大きい方（つまり確率 1）を A(ν → µ) にす
ると、式 (C.2) を満たすために A(µ → ν) は e−β(Eν −Eµ ) になる。
{
A(µ → ν) =
e−β(Eν −Eµ )
1
if Eν − Eµ > 0
otherwise
(C.3)
つまり、新しい状態が現在の状態よりも低いまたは等しいエネルギーを持つ場合は、常
に新しい状態への遷移を受け入れ、反対にエネルギーが高くなる場合には、上で与えら
れた確率で状態の遷移を受け入れる。これが single-spin-flip dynamics をもちいた Ising
model の Metropolis algorithm である。
42
付録 C
Metropolis algorithm
実際に Metropolis algorithm を計算に用いるとき、状態をどうやって確率 A(µ → ν) ≤
1 で採択して、確率 1 − A(µ → ν) で棄却するのかという問題がある。これは、0 と 1 の間
の一様な乱数 r を生成させることによって解決できる。それは r < A(µ → ν) なら状態を
採択し、そうでないなら棄却する。同じ確率 A(µ → ν) を使って、何度もこの手順を繰り
返せば、全試行回数に対して A(µ → ν) の割合で状態が採択されることは明らかである。
C.0.9
Problems in Metropolis algorithm
ここでは、Metropolis algorithm の問題点を議論する。Ising model の平衡状態では、
スピンは大きさ ξ の cluster をつくる。ξ は相関長 (correlation length) とよばれ、転移温
度 Tc で発散する。このとき、reduced temperature とよばれる無次元のパラメーターを
定義する。
t=
T − Tc
Tc
(C.4)
この t を使って相転移点近傍の ξ を表すと、
ξ ∼ |t|−ν
(C.5)
となる。（ν は臨界指数 (critical exponent) と呼ばれている正の量。）臨界指数は空間次
元や秩序変数の対称性だけから決まる普遍性を持ち、個々の体系（Ising model では、交
換相互作用の強さや格子の形状）によらない。また ν は MC シミュレーションに用いる
algorithm にもよらない。さらに、相関時間 (correlation time) τ と呼ばれる格子あたり
の MC ステップ数によって測られる量を定義する。
τ ∼ |t|−zν
(C.6)
ここでもちいた z は dynamic exponent と呼ばれ、臨界指数 ν などのような普遍性を
持った量ではなく、用いたシミュレーション algorithm に依存する。z の値が大きけれ
ば、相転移点に近づくにつれて相関時間は急激に大きくなり、シミュレーションに時間が
かかるようになる。逆に z の値が小さければシミュレーションは早く終わるようになる。
式 (C.5) と式 (C.6) をまとめると、
τ ∼ ξz
(C.7)
となり、転移点近傍で相関長が大きく発散していくにつれ、どの程度相関時間が長くなる
か教えてくれる。しかし実際には、有限の大きさを持つ系では相関長は転移点近傍では発
散しない。系の大きさが Ld （d は系の次元）のとき、相関長はカットされて式 (C.7) は
τ ∼ Lz
(C.8)
となる。これで実際に z をシミュレーションで見積もることのできる表式になった。z を
Metropolis algorithm を使って計算するとおおよそ z = 2.17 になり、これは２次元 Ising
model に対する MC アルゴリズムの中でも大きな値である。つまり他のアルゴリズムに
比べ時間がかかるということである。まず Metropolis algorithm の問題点ををスピンの
43
付録 C
Metropolis algorithm
観点から見てみる。転移点近傍では相関長が大きくなり、スピンは同じ向きに揃いだして
domains をつくる。（注意：domain はスピンの向きが揃った領域のことであり、cluster
はスピン間の相関がある領域のことである。）強磁性的なモデルの場合、domain の中心
近くのスピンは周りを同じ向きのスピンに囲まれているので、フリップする可能性はかな
り低く、それによって domain 全体をフリップさせることが難しくなる。domain の端の
スピンならばフリップする可能性があり、これが domain 全体をフリップさせようとした
場合の主要過程となる。しかし、転移点近傍では domain はかなり大きくなっており、端
のスピンをフリップさせているだけではかなりの時間を必要とする。もう一つの大きな理
由は、Metropolis algorithm の採択確率の選び方である。このアルゴリズムでは採択確
率がボルツマン因子に比例しており、温度が低くなると採択確率が極端に小さくなって
しまう。そのような場合、新しい状態が生成されにくくなり、平衡状態に達するのに非
常に時間がかかってしまい、シミュレーション効率が悪いのである。以上が Metropolis
algorithm の問題点である。
44
付録 D
付録 D
Variational Monte Carlo
Variational Monte Carlo
変分モンテカルロ法 (VMC) は変分原理を用いて、適切なパラメターを持った変分関数
を設定し、変分エネルギーの最小値を数値的に求める手法である。系の変分エネルギー
は、S 個の変分パラメターに依存する試行多体波動関数 ψT を用いて、
∫
hEi =
dR ψT∗ (R)HψT (R)
∫
dR |ψT |2 (R)
(D.1)
と表される。ここでは、N 個の粒子のすべての位置座標を R = r1 , · · · , rN と書く。この
変分エネルギーが最小となる波動関数を求めたいのだが、現実の系では多体波動関数は、
配置空間の広い領域で小さな値をとると考えられ、一様に分布したランダムな点をサンプ
リングする単純な方法では効率が悪い。そのため、波動関数が大きな値をとると考えられ
る配置空間の領域を、重点的にサンプリングするためメトロポリスアルゴリズムを使うと
よい。試行関数を用いて、局所エネルギー
EL (R) =
HψT (R)
ψT (R)
(D.2)
を定義する。EL (R) は粒子の位置に依存し、ψT が厳密な波動関数に近づくにつれて、
EL は R に応じた変化が弱くなり、ハミルトニアンの厳密な固有関数であれば定数とな
る。そこで局所エネルギーを用いて、変分エネルギーを書き直すと、
∫
hEi =
dR |ψT (R)|2 EL (R)
∫
dR |ψT (R)|2
(D.3)
となるため、メトロポリスアルゴリズムによるサンプリングには、次の定常分布 ρ(R) を
使う。
ρ(R) = ∫
|ψT (R)|2
dR |ψT (R)|2
(D.4)
以上をまとめると、変分モンテカルロ法では、変分パラメターを変化させながら、定常
分布 ρ(R) を用いたメトロポリス・モンテカルロ計算により局所エネルギー EL が最小と
なる場合を探ることになる。この方法は非常にシンプルであり、変分パラメターの数が少
なければ計算時間が短くなるため、様々な研究に使われている。しかし、得られるエネル
ギーが試行関数 ψT に大きく依存するため、波動関数の選び方に注意する必要がある。
次に紹介する拡散モンテカルロ法は、試行関数の選び方に原理的に依存しないため、よ
り正確な計算が可能となる。
45
付録 E
付録 E
Diffusion Monte Carlo
Diffusion Monte Carlo
シュレーディンガー方程式と拡散方程式の類似性を用いて、粒子の古典的拡散過程のシ
ミュレーションを行なうことにより、相互作用する量子力学的粒子の特性を計算する方法
を、拡散モンテカルロ法 (DMC) と呼ぶ。虚時間のシュレーディンガー方程式は、虚時間
を τ = it/~ として
−
]
∂Φ(R, τ ) [
= −D∇2 + V (R) Φ(R, τ )
∂τ
(E.1)
と表される。ここで D = ~2 /2m であり、この方程式は 3N 次元空間における拡散方程
式の類似形になっている。通常の拡散方程式はグリーン関数を用いて解く事ができるが、
虚時間のシュレーディンガー方程式ではポテンシャル項が存在するため、そのままではグ
リーン関数全体の規格化ができない。そこで規格化因子として exp(τ ET ) を導入する。こ
の場合に、方程式は
−
]
∂Φ(R, τ ) [
= −D∇2 + V (R) − ET Φ(R, τ )
∂τ
(E.2)
と修正され、グリーン関数をもちいて解くことが出来る。つまり、ランダムウォーカーの
拡散過程として数値解を求めることができ、この場合 Φ(R, τ ) は拡散するウォーカーの密
度分布と見なされる。そしてこの拡散過程は、拡散定数 D の拡散プロセスと、通常の拡
散方程式にはない項である、[V (R) − ET ]Φ(R, τ ) の効果を反映する、分岐プロセスから
構成される。
波動関数の時間発展は、グリーン関数 G を用いて
∫
0
Φ(R , τ + ∆τ ) =
dR G(R, R0 , ∆τ )Φ(R, τ )
(E.3)
と書かれ、グリーン関数は詳細釣り合いの条件を満たす遷移確率となる。さらに詳細釣り
合いの条件を満たしているため定常状態が存在し、その場合の拡散方程式は、
[
]
−D∇2 + V (R) Φ(R, τ ) = ET Φ(R, τ )
(E.4)
となり、時間に依らないシュレーディンガー方程式を表す。またグリーン関数は、時間間
隔 ∆τ が十分に小さい場合、鈴木・トロッター分割を用いて、
0
G(R, R , ∆τ ) ' (4πD∆τ )
−3N/2
{
}
{
}
[R0 − R]2
V (R) + V (R0 )
exp −
×exp −∆τ [
− ET ]
4D∆τ
2
(E.5)
と表される。
シミュレーションでは、グリーン関数の拡散部分、つまり運動エネルギーによる部分で
決まる遷移確率で、ウォーカーを新しい位置に動かす。そしてポテンシャル項の効果は、
位置 R0 に来たウォーカーに対して重み因子として働く。しかし、この手続きはあまり効
率的ではない。なぜなら、重みは小さいのに、もっと有利なところに移動した場合と、同
46
付録 E
Diffusion Monte Carlo
様な計算量が必要となるためである。そのため、各ウォーカーが動く領域の重要性に比例
して、不利な場合はウォーカーを消滅させ、有利な場合はウォーカーを新しく生成すると
いう分岐過程を入れることで、計算をより効率的に出来るようになる。
虚時間のシュレーディンガー方程式を解く際、ランダムな拡散過程を用いるとシミュ
レーション効率が悪くなる場合がある。例えば、ポテンシャルがある条件下で発散してし
まう場合などは、分岐プロセスにおいて過度な生成・消滅が行なわれることになり、定常
分布への収束が遅くなる。この問題を改善する方法が、Importance sampling である。こ
の方法では、適当な試行関数 ΨT (R) を選び、
(E.6)
f (R, τ ) = Φ(R, τ )ΨT (R)
と密度分布を置換する。これを用いると、式 (E.2) は、
−
∂f (R, τ )
= −D∇2 f (R, τ ) + (EL (R) − ET )f (R, τ ) + D∇ · [f (R, τ ) F (R)] (E.7)
∂τ
と書き直される。ここで、
EL (R) =
HΨT (R)
,
ΨT (R)
F (R) = 2
∇ΨT (R)
ΨT (R)
(E.8)
である。式 (E.7) では第１、第３項が拡散プロセスに相当し、第２項は分岐プロセスに使
われる。
この形であれば、真の波動関数に近い試行関数を選ぶ事が出来れば、EL (R) が大幅に
平坦化されるため、効率的にサンプリングを行なうことができる。
なお、Importance sampling を用いた場合、密度分布の時間発展は
∫
0
f (R , τ + ∆τ ) =
dR G(R, R0 , ∆τ )f (R, τ )
(E.9)
となり、グリーン関数は
0
−3N/2
G(R, R , ∆τ ) ' (4πD∆τ )
{
}
[R0 − R − D∆τ F (R)]2
exp −
4D∆τ
{
}
EL (R) + EL (R0 )
× exp −∆τ [
− ET ]
(E.10)
2
と表される。
以上が拡散モンテカルロ法の大まかな説明であるが、この数値計算では鈴木・トロッ
ター分割を用いているため、虚時間の時間間隔 ∆τ に依存した誤差が残る。正しい結果を
得るには、複数の ∆τ の値について計算結果を得た後に、外挿によって ∆τ → 0 のときの
値を求めなければならない事に注意してほしい。
より詳しい説明は Reynolds15) らの論文を参照していただきたい。
47
付録 F
付録 F
Path Integral Monte Carlo
Path Integral Monte Carlo
経路積分モンテカルロ法 (PIMC) は Feynman の経路積分の概念に基づいており、ここ
では経路積分についての簡単な定式化を与え、モンテカルロ法によって何を計算するのか
を述べる。
有限温度での量子系の分配関数 Z は、
Z(β) = Tr(e
−βH
∫
)=
dR R|e−βH |R
(F.1)
と書かれる。ここで、R は N 個の粒子の位置座標を表し、β = 1/(kB T ) である。式 (F.1)
はハミルトニアンの指数関数を含んでおり、時間発展演算子のトレースを直接計算するの
は難しい。そこで、指数が十分小さければ短時間近似を行い計算をする事が可能であるた
め、β を M 個に分割することを考える。このとき演算子は、τ = β/M として
)M
(
eβH = eτ H
(F.2)
と表される。この τ はタイムステップと呼ばれ、PIMC において非常に重要なパラメター
である。分割された時間発展演算子を用いて、分配関数は
∫
dR0 R0 |e−βH |R0 =
∫
dR0 dR1 dR2 · · · dRM −1
× R0 |e−τ H |R1 R1 |e−τ H |R2 · · · RM −1 |e−τ H |R0
と書ける。時間発展演算子の間に、M − 1 個の単位演算子
∫
(F.3)
dRi |Ri i hRi | を挿入した。
また行列要素の積で、最初と最後の状態が等しいことから、τ 方向に周期境界条件がある
ことが分かる。短時間の時間発展演算子の行列要素は、
R|e−τ H |R0 =
0 2
1
e−τ V (R) e−(R−R ) /(2τ )
3N/2
(2πτ )
(F.4)
であるから、式 (F.3) に代入すると、
∫
dR0 R0 |e−βH |R0 =
1
(2πτ )3N M/2
{
∫
dR0 dR1 dR2 · · · dRM −1
× exp −
M
−1 [
∑
m=0
(Rm+1 − Rm )2
+ τ V (Rm )
2τ
]}
(F.5)
となる。この式では RM = R0 である。
上式の指数関数の引数に含まれる、和の個々の項は座標 Ri の古典的多粒子系の離
散化時間でのラグランジアンになっている。すると、和は作用を意味し、積分は座標
R0 , · · · , RM のすべての可能な組についての和となる。このような組は配置空間内の経路
を表し、世界線あるいはビーズなどと呼ばれている。
48
付録 F
Path Integral Monte Carlo
N 個の３次元粒子の量子分配関数は、虚時間方向に M 個のバネでつながれた粒子の、
３次元 N 粒子古典的分配関数と解釈することができる。よって、量子分配関数をもとめ
るには、単に古典系で標準的なモンテカルロシミュレーションをすればよい。シミュレー
ションの簡単な流れは次のようになる。「ランダムに選ばれた、ある時間スライス m̃ にお
いて、N 個の粒子の中から１個をランダムに選ぶ。それにランダムな変位を加え、採択
確率が p = exp[−τ (Hold − Hnew )] のメトロポリス判定をおこない、新しい状態を生成
する。そして、これを平衡状態に達するまで繰り返す。」以上が、PIMC の基本概念とシ
ミュレーション方法である。
最後に、経路積分においてエネルギーがどのように書けるのか示す。有限温度での量子
系の物理量 A の期待値は、
hAi =
Tr(Ae−βH )
Tr(e−βH )
(F.6)
で表されるので、経路積分におけるエネルギーの期待値は、
Tr(He−βH )
∂
=−
ln Z(β)
Z
∂β
= hU1 i + hU2 i
hEi =
U1 =
M
−1
∑
3
1
−
(Rm+1 − Rm )2
2τ
2M τ 2 m=0
M −1
1 ∑
U2 =
V (Rm )
M m=0
(F.7)
(F.8)
(F.9)
(F.10)
となる。
今回使った PIMC のコードには、計算精度を上げつつ計算時間を短くする様々な工夫
がなされている。詳しくは Ceperley の論文を見て頂きたい。
49
引用、参考文献
引用、参考文献
1)
2)
3)
4)
5)
6)
7)
8)
9)
10)
11)
12)
13)
14)
15)
16)
17)
18)
19)
20)
21)
22)
23)
24)
L. H. Nosanow, L. J. Parish, and F. J. Pinski, Phys. Rev. B 11 (1975) 191.
小池 幸生、修士論文（平成 15 年度、名古屋大学理学研究科）
B. J. Alder and T. E. Wainwright, J. Chem. Phys. 33 (1960) 1439.
W. G. Hoover and F. H. Ree, J. Chem. Phys. 49 (1968) 3609.
J-P. Hansen, D. Levesque, and D. Schiff, Phys. Rev. A 3 (1971) 776.
M. H. Kalos, D. Levesque, and L. Verlet, Phys. Rev. A 9 (1974) 2178.
K. J. Runge and G. V. Chester, Phys. Rev. B 38 (1988) 135.
L. M. Sesé, J. Chem. Phys. 126 (2007) 164508; L. M. Sesé and L. B. Bailey, ibid, 126
(2007) 164509.
H. Freimuth, H. Wiechert, H. P. Schildberg, and H. J. Lauter, Phys. Rev. B 42 (1990)
587.
T. Momoi, Phys. Rev. B 67 (2003) 014529.
D. S. Hirashima, T. Momoi, and T. Takagi, J. Phys. Soc. Jpn. 72 (2003) 1446.
H. Kitamura, S. Tsuneyuki, T. Ogitsu, and T. Miyake, Nature 404 (2000) 259.
T. Cui, Y. Takada, Q. Cui, Y. Ma, and G. Zou, Phys. Rev. B 64 (2001) 024108.
T. E. Markland, J. A. Morrone, B. J. Berne, K. Miyazaki, E. Rabani, and D. R.
Reichman, Nature Phyics 7 (2011) 134.
P. J. Reynolds, D. M. Ceperley, B. J. Alder, and W. A. Lester, Jr., J. Chem. Phys. 77
(1982) 5593.
C Cazorla, G E Astrakharchik, J Casulleras, and J Boronat, New J. Phys. 11 (2009)
013047.
D. M. Ceperley, Rev. Mod. Phys. 67 (1995) 279.
D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic, San Diego,
2002)
D. Frenkel, Phys. Rev. Lett. 56 (1986) 858.
W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes
in FORTRAN. The Art of Scientific Computing (Cambridge Universiity Press, Cambridge, 1992).
W. G. Hoover, M. Ross, K. W. Johnson, D. Henderson, J. A. Barker, and B. C. Brown,
J. Chem. Phys. 52 (1970) 4931; W. G. Hoover, S. G. Gray, and K. W. Johnson, ibid.,
55 (1971) 1128.
R. P. Feynman, Statistical Mechanics: A set of Lecture (Westview Press, Boulder, 1972)
M. Boninsegni, L. Pollet, N. Prokof’ev, and B. Svistunov, Phys. Rev. Lett. 109 (2012)
025302.
G. C. Straty and E. D. Adams, Phys. Rev. Lett. 17 (1966) 290.
50