scipyでMetropolisアルゴリズムを書いた

matplotlibのmoviewriterを使って動画を生成しようとしている。
どうもFFmpegというのをダウンロードしないといけないらしい。

色々弄った結果、三角格子状の反強磁性イジングモデルメトロポリスアルゴリズムのムービーを作成できるようになったから満足。

長距離力が働く系の熱力学、統計力学[執筆中]

前回の記事、「水素型ポテンシャルの分配関数」(http://d.hatena.ne.jp/tomosoeoc/20150404/1428122715)では、
水素型ポテンシャルの分配関数が体積にどのように依存するかについて分析しました。
そこで、体積が小さい領域では分配関数がVの1/2乗に比例することを根拠に、この系は熱力学的にやや不自然な振る舞いをすると結論しました。
この点について友人から幾つか疑問点が提示されたので、それについて考えてみたいと思います。

ポイントは
・「体積が小さい領域」とは、Vの1/2乗の項がVに線形な項に比べて無視できない領域のことを言う。
この領域の大きさは、「ミクロな大きさ」とは限らない。(なぜなら1/rは長距離力だから)
実際、電荷などのパラメーターが変わればVの1/2乗の項の係数は変わるし、「体積が小さい領域」の体積はいくらでも大きくなりうる。
・長距離力の存在下では、系は相加性を満たさない。(清水熱力学より)
・V^1/2乗という寄与は、Vが無限大に飛ぶ極限を取れば消えてしまうので、系の挙動に影響はないかのように見えるが、
実際には系の様々なパラメーターのゆらぎの大きさと同じなので、そこに影響しうる、気がする。(まだ真面目に考えていない)

水素原子型ポテンシャルの分配関数

Twitter上で友人に教えて貰った問題が面白かったので、ここに書いてみたいと思います。

水素原子型のポテンシャル -1/rを考えると、
E<0の領域で固有値は-1/n^2に比例するので、
ナイーブに分配関数を評価すると発散してしまいます。
何故発散するのか?どのようにすれば発散は防げるのか?というのが問題です。

結論から言えば箱に入れて、ポテンシャルの範囲を有限の領域に制限すれば良いです。
まず結果から説明しましょう。

半径Rで体積Vの箱に入れた場合を想定します。
このとき、ポテンシャルの最大値は-1/Rとなり、無限大の箱の中に入っていた時(E=0)よりも小さいです。
そのため、束縛状態の数も減少します。
大雑把には、ポテンシャルの最大値よりも低いエネルギーを持つ束縛状態は生き残り、
そうでない状態は自由粒子のスペクトルの中に取り込まれます。

生き残った束縛状態の数を数えるにはまず、許されるnの数を数えましょう。

  • 1/n^2 < (定数) *1/R

を満たすnの数を数えればよいので、結局
nはRの2分の1乗に比例することがわかります。
角運動量の自由度も考えると、各nに対してn^2個の状態があります。
よって、許される束縛状態の数はおよそn^3に比例することがわかります。
nが十分に大きい範囲からの寄与を考えることにすると、結局分配関数は
Rの2分の3乗に比例することがわかります。


自由粒子の領域も考慮した上でこれらの結果をまとめてみましょう。
初歩の統計力学から、自由粒子の領域の分配関数はおよそ体積に比例することがわかる。
一方で、上記の計算から、水素原子の束縛状態の分配関数は体積の2分の1乗に比例することがわかります。
合わせると、体積が小さな領域では分配関数は束縛状態に支配され、
体積が大きい領域では分配関数は自由粒子状態に支配されるようです。
普通、分配関数は体積に比例するので、体積が小さな領域での振る舞いが非常に非熱力学的であることもわかります。


最後に、箱のなかに入れる操作が妥当な理由について考えましょう。
そもそもこのポテンシャルはE>0に連続領域を持っているので、この部分の分配関数をまともに評価したければ、最初から箱に入れておくしかありません。
それにともなってE<0への領域も変更を受け、結果的にVの2分の1乗の依存性が出るのは上で見たとおりです。
よって、今回の発散は、本来は箱に入れるべきポテンシャルを箱に入れなかったために起こった人為的なものだと考えられます。

実用上の問題について何点か。
水素原子の束縛状態部分の分配関数を議論したい場合は、わざわざ箱に入れて計算せずとも、
分配係数の計算を適当なnで打ち切ってしまえば十分でしょう。よほど大きな変化がない限り、熱的にアクセス出来るエネルギーの範囲は変わらないので、これでまっとうな答えが出るはずです。
また、以上の議論から「温度Tにある水素原子のアンサンブルのエネルギー期待値を求めよ」という問題は、体積に依存する非常に厄介な問題だということもわかります。

グランドカノニカル分布

ふと友人から言われて勉強しなおしたので一言。

グランドカノニカル分布では、確率は粒子数に対してexp(N\beta\mu)の形で線形に依存している。この時、\muの値には全くN依存性がない。
状況設定をよく把握せずにこの式を解釈しようとしたらここで引っかかってしまった。つまりこういうミスだ。

ある系の粒子数を変動させることを考える。この時、最初の粒子数から次の粒子数に移るときは、化学ポテンシャルの値はある初期値\mu_0をとっているだろう。
しかし、粒子数を変動させるに連れて、系の様子もどんどん変化していくのだから、\muの値も徐々に変わっていくはずだ。
よって、確率が粒子数に比例するのはおかしくて、本当は積分系で書かれるべきだ。

このミスが起こったのは、グランドカノニカル分布が「粒子浴」と隣接しているという状況にあることを失念していたからだと思う。
\muの値は系の特質ではなくて、非常に巨大な粒子浴の性質なのだ。あまりに巨大なので、系に出入りする程度の粒子数では化学ポテンシャルは変化せず、
一定として扱ったとしても問題はなくなっている。

グランドカノニカル分布は、導出の際の設定と、実際の応用の仕方がカノニカル分布の場合以上に直感に反しているので、
こうした勘違いが起こりやすいのではないかと思う。

日記

これからはもう少し頻繁に自分の思考を記録していこうと思う。
この二年間、自分と向き合うことを放置していたつけが回ってきているような気もするから。

まあ何より、あとで読み返した時に面白いだろう。

研究室

ここ一ヶ月ほどまともに研究室に行って実験をしていない。
週4日課題と戦うだけで心身ともに疲れ果てて、金曜〜土曜は実験を行う気力が残っていない事が多いのだ。
春休みを使って実験をしようとも計画していたのだが、結局風邪を引いて寝込んでしまった。
面白い現象が見える兆候は色々とあるのに、自分の責任でプロジェクトが止まっているのは申し訳ない。

好みの問題

量子力学の数学的構造がよくわかったとして、頭のなかがすっきりするだけで、
計算ができるようにもならないし、実験の役にも立たない気がするけれど、
「頭のなかがすっきりするだろう」という希望だけで勉強してしまう。