続・君はあの相転移を見たか

これが、見たかった―――。

始まりは、的確な指摘から―――。

　前回こんな記事を書きました。
君はあの相転移を見たか - 幸福の物理
　いろいろと反響があって嬉しかったのですが、なかでも関心したのがとある物理学科の友人の指摘。

　「それってただ温度に対応する粒子の分布を示しただけである。それだけでは相転移が見えたとは言えぬ。」
　確かにとある友人の言う通りである。というわけで今回はちゃんと物理量を測定して検証しました。

レッツ観測！

プログラムを改良

　従来のプログラムでは無駄な計算が多かったために計算に非常に時間がかかっていた。したがって、フラグをたてるといった枝刈りを行った。
　また、前回では周期的境界条件で求めていたが、それは考えている系は壁のようなもので覆われているという仮定に不適であると判断。ゆえに境界条件を壁に変更した。
　さらに、前回導入したGLFWGLFW、始めました。 - 幸福の物理をGLUTの代わりに用いてみた。

実験内容

　求める物理量はズバリ「比熱」！もし、ある(逆)温度で相転移が起こるならば、その前後における(逆)温度に対応する比熱の成長の仕方に急激な違いが見られることが知られています。また、相転移点で比熱が不連続になる(系が無限に大きいならば発散する)ことも知られています。
　そこで、今回は以下の条件で物理量を測定しました。

• 測定するものは「逆温度」と「比熱」である。
• 比熱は系全体のエネルギーのゆらぎに逆温度をかけたものである。

\begin{eqnarray}
C = k_{B} \beta^2 (\langle E^2 \rangle - \langle E \rangle^2)
\end{eqnarray}

• 系の大きさは62×62(縦横60のち両端2マスは壁、のこり60マスは粒子が存在しうる領域)
• 粒子数は1080
• 壁との相互作用はゼロ
• 相互作用ポテンシャルは1、重力加速度は0.01
• モンテカルロ法を100万回実行し平衡状態に達した後、全体のエネルギーを観測する。これを100回行い、エネルギーのゆらぎを求めて比熱を算出する。これを各逆温度において実行する。

実験結果

　これは逆温度を0.01刻みでシミュレートした結果です。

　そして、これは逆温度を0.001刻みでシミュレートした結果です。

つまり・・・どういうことだってばよ！？

　0.001刻みはさすがにやり過ぎた感じがあります。これ、計算が終わるまで15時間くらい待たされました。こんなに時間がかかるならばもう少し計算の工夫を考えればよかったかもしれません。
　それはさておき、どうやら\( \beta = 1.5 \)で相転移が生じているみたいです。

　フィッティング(とはいってもラフなものですが)するとこんな感じでしょうか。

青線が今回のシミュレート結果のフィッティング、緑が理論値(厳密には予想)です。ちゃんと相転移が表れているようです。

グラフだけじゃわかんねぇよ・・・。

　各温度では、粒子はいったいどのような様子でしょうか。実際に見てみましょう。
　以下、図の青四角は壁を、赤四角は粒子をそれぞれ表します。

　この図は\( \beta = 0.0 \)のときの平衡状態の1例です。

青い壁に囲まれた空間全体に一様に分布していますね。完全に気体の状態といっていいでしょう。

　この図は\( \beta = 0.7 \)のときの平衡状態の1例です。

こころなしか粒子同士がくっつき始めているように見えます。この時点ではまだ気体の状態ですね。

　この図は\( \beta = 1.5 \)のときの平衡状態の1例です。

一気に様子が変わりました。底の方に粒子が溜まり始めていますね。ただ、空間上部にもまだ粒子が存在しています。これが気液平衡状態ですね。

　この図は\( \beta = 3.0 \)のときの平衡状態の1例です。

殆どの粒子は底に溜まっています。ほぼ液体の状態であるといっていいでしょう。\( beta = 0.7 \)とほぼ同じ比熱であるのに、それぞれの状態は全く異なります。おもしろいですね。

　この図は\( \beta = 5.0 \)のときの平衡状態の1例です。

完全に液体の状態になりました。

まだまだ深い！？

　ながながと見てきましたが、まだまだ深い話が有るようです。粒子数が変化する場合はどうなるのか。平均場近似じゃない近似だったらどうなるのか。話題は絶えません。これからもっと調べていきたいです。

(((՞ةڼ◔)))＜幸福のゆらぎィ～！