調和振動子の分布

個人的な疑問が２つ解消されたのでメモ。

コンテンツ

• 単体の単振動
• カノニカル集団での単振動

それぞれ位置とエネルギーの分布関数を調べます。

単体の単振動

調和振動子があります。ポテンシャル

U(x)= \frac{1}{2}kx^2 

に従って単振動しています。普通の古典力学の範囲でこのポテンシャルでの運動方程式を解くと、座標は

x = A\sin(\omega t + \alpha)

で、速度は座標を使って

v = \omega A\cos(\omega t + \alpha) = \pm\omega \sqrt{A^2 - x^2}

とかけます。符号が適当ですが、今回の議論では使わないので適当にさせてください。

位置の分布関数

さて、時刻 $t$ が一様分布する場合、位置はどういう分布をしているでしょうか。これが第１の疑問。単純に知らなかっただけですが、調べてみました。

答えを言うと逆正弦分布です。

P(x \mid E) = \frac{1}{Z\sqrt{A^2 - x^2}}

導出

座標の微小区間 $x\sim x+dx$にいる時間 $dt$ は速度で割ればいいので、

dt = \frac{dx}{\omega \sqrt{A^2 - x^2}}

とかけます。時刻は一様分布としたので、この座標にいる時間はその座標にいる確率と考えます。すると逆正弦分布が得られます。

または、時間に関して何らかの期待値を計算する積分を書いたとき、この微小時間の関係を使って変数変換したら、位置に関する期待値の式になるので、その内の確率分布部分とも考えられるのではないでしょうか。

また、調和振動子が動ける範囲は限られています。全エネルギー

E = \frac{1}{2}kx^2 + \frac{1}{2}mv^2

は一定で、 $x = A$ で位置エネルギーが最大の時 $v=0$ となりますので、$A = \sqrt{2E/k}$ ということになります。また$x$ が存在する範囲は $-A \le x \le A$ になります。これを踏まえると、確率分布関数をしっかり書くと

P(x \mid E) = 
\left\{\begin{array}{lr}
\frac{1}{Z\sqrt{A^2 - x^2}} &\left(-A \le x \le A\right)\\
0&(\mathrm{otherwise})
\end{array} \right.

となります。 ここで

Z = \int_{-A}^{A}\frac{dx}{\sqrt{A^2 - x^2}} = \pi

は規格化のための係数です。

位置エネルギー

また、位置エネルギーは座標の２乗に比例するので、倍角公式 $\cos 2\theta = 1 - 2\sin^2 \theta$ を使えば

U(t)= \frac{1}{2}kx(t)^2 = \frac{1}{2}kA^2\sin^2(\omega t + \alpha) = \frac{1}{4}kA^2\left(1 - \cos(2\omega t + 2\alpha)\right)

とかけますので、これも(中心がずれた)逆正弦分布になります。

\frac{dU}{dt} = \frac{1}{2}\omega kA^2\sin(2\omega t + 2\alpha) = \omega\sqrt{\left(\frac{1}{4}kA^2\right)^2 - \left(U - \frac{1}{4}kA^2\right)^2}

というポテンシャルエネルギーの時間変化（位置に対する速度を考えました）から、

P(U) = \left\{\begin{array}{l}
\frac{1}{\pi\sqrt{\left(E/2\right)^2 - \left(U - E/2\right)^2}} \ \ \ \  \left(0 \le U \le E\right)\\
0
\end{array} \right.

という感じになります。普通の逆正弦分布から座標を $kA^2/4 = E/2$ ずらした感じですね。

実験

実際にPythonで確かめてみます。$A = 1$、$E = 1$、$k = 2$ とします。

import numpy as np
import matplotlib.pyplot as plt 
t = np.random.uniform(0, 2*np.pi, size=2000)
xt = np.sin(t)
x = np.linspace(-1, 1, 202)[1:-1]
fx = 1./np.pi/np.sqrt(1. - x**2)
count, bins, hist = plt.hist(xt, bins=40)
dx = bins[1] - bins[0]
plt.plot(x, fx*xt.size*dx)
plt.show()

Ut = xt**2
U = np.linspace(0, 1, 202)[1:-1]
fU = 1.0/np.pi/np.sqrt((0.5)**2 - (U - 0.5)**2)
count, bins, hist = plt.hist(Ut, bins=40)
dU = bins[1] - bins[0]
plt.plot(U, fU*Ut.size*dU)
plt.show()

結果は以下のようになります。

• 座標

• ポテンシャルエネルギー

再現できているかなと思います。

カノニカル集団での分布関数

一方、分子動力学シミュレーションを行い、調和振動子ポテンシャルに従う化学結合の結合長の分布を調べてみると、正規分布のような形をしています。そしてポテンシャルエネルギーはラプラス分布のように見えます。
熱浴下での運動とはいえ、単体だと一番分布が薄かったところが今度は一番存在することになります。これはどういう仕組みなんでしょうか。僕には直感的でなかったので、これが第二の疑問です。
これを確かめます。

位置の分布関数

位置の確率分布関数を全エネルギーで周辺化してみます。

P(x) = \int_{E_0}^\infty dEP(x\mid E)P(E)  

第１項 $P(x\mid E)$ は全エネルギーを固定したときの分布で、逆正弦分布です。
第２項はこの振動子がもつエネルギーの分布です。さて、分子動力学は統計集団を取り扱いますが、その部分系はカノニカル集団とできます。この振動子がカノニカルだとすると、エネルギーはボルツマン分布 $P(E) = \exp(-\beta E)$ になります。

積分範囲の下限 $E_0$ は $x$ がいるところの位置エネルギー $U(x)$ です。もしこれ未満の$E$ という状況を考えると、運動エネルギーが負ということになり物理的に考えない範囲になります。また逆に考えて、上で定義した $E$に対する $x$ の範囲を思い出してみると、そこからはみ出した $x$ があり、これは $P(x \mid E) = 0$ となって積分範囲から除外されます。

ちゃんと代入すると、

P(x) = \int_{\frac{1}{2}kx^2}^\infty dE\frac{\exp(-\beta E)}{\pi\sqrt{\frac{2E}{k} - x^2}}  

ここまで書くと多分普通に積分できますが、今回はMaximaに任せてみます。wxMaximaを立ち上げ、

integrate(exp(- β*E)/(π*sqrt(2*E/k - x^2)), E, k*x^2/2, inf);

と入力、Shift+Enterを押すと、$\beta$ と $k$ の正負を聞かれるので両方 positive ( + Shift+Enter) と入力すると、

\frac{\sqrt{k}\, {{\% e}^{-\frac{k\, {{x}^{2}} \beta }{2}}}}{\sqrt{2}\, \sqrt{\pi }\, \sqrt{\beta }}

という出力が得られます。これを整理すると、

\sqrt{\frac{k}{2\pi\beta}}\exp \left(-\frac{k \beta x^2}{2} \right)

となり、確かに正規分布することがわかりました。一応規格化しておくと

P(x) = \sqrt{\frac{k\beta}{2\pi}}\exp\left(-\frac{k\beta x^2}{2}\right)

となります。$P(x)dx$ が確率としての意味を持つので、$P(x)$ は位置の逆数の次元を持ちます。

• 追記

投稿直後に気づいたのですが、ボルツマン因子を $P(E) = \beta\exp(-\beta E)$ とちゃんと規格化しておけば $P(x)$ も規格化された状態で計算されましたね。

位置エネルギーの分布関数

またポテンシャルエネルギーの分布ですが、同様に考えると

P(U) =  \int_{E_0}^\infty dEP(U\mid E)P(E) 
\sim \int_{U}^\infty dE\frac{\exp(-\beta E)}{\sqrt{\left(E/2\right)^2 - \left(U - E/2\right)^2}}

となります。これもMaximaに投げることにします。

integrate(exp(-β*E)/√((E/2)²-(U - E/2)²), E, U, ∞)

を入力、$\beta$ について聞かれますので positive と入力すると

\frac{\sqrt{\pi}\, {{\% e}^{-U \beta }}}{\sqrt{U}\, \sqrt{\beta }}

が得られます。整理すると

\sqrt{\pi} \frac{\exp(-\beta U )}{\sqrt{\beta U}}

となりました。微妙にラプラス分布とはちがいますね。なお、規格化のためには係数を $\sqrt{1/\pi}$ とするとよいです。

これはまあ $k=1$ のときの$\chi$2乗分布

f(x) = \sqrt{\frac{1}{2\pi}}\frac{\exp(-x/2)}{\sqrt{x}}

の形をしています。 $x \rightarrow 2\beta U$ とすると係数が4倍違うのですが、まあ今回は無視します。正規分布する座標 $x$ の２乗が位置エネルギーですので、Wikipediaにある$\chi^2$ 分布の定義からも 位置エネルギーの分布は $\chi^2$ になるのが妥当そうです。

確認

別件で SPC/E flexibleモデルの水のMDを行っていて、そのトラジェクトリを使ってみます。
このモデルでは $k = 345000\ \mathrm{kJ/mol/nm^2}$です。MDは $300\ \mathrm{K}$で行いました。

一方で、熱浴下の振動子は熱浴との相互作用により周波数が変わることが知られています。何となく、 $k = 312000\ \mathrm{kJ/mol/nm^2}$ という数字を受信したので、これを使ってみます。

このパラメータを分布関数の式に代入しますと、以下のようなプロットを得ます。

青い線が $k = 345000$、緑の線が $k = 312000$ としたときの線です。これから、変調した周波数のほうがよくフィットすることが分かります（もっとちゃんと言うなら尤度なりp値なりを確認すればいいかなと思います）。これは周波数の変調は、カノニカル分布下でエネルギーが調和振動子に分配されていることからは示せないことが言えるのではと思います。実際、調和振動子熱浴モデルにより説明されます。
https://ocw.kyoto-u.ac.jp/ja/graduate-school-of-science-jp/course-chemical-statics/pdf/lect05.pdf

ポテンシャルエネルギーのほうの確認は面倒になったので省略します。

その他

カノニカル集団に現れる$\chi^2$分布として、今回はポテンシャルエネルギーを示しました。他には分子の並進運動における速度の分布がマクスウェル分布になることが知られていますが、これは自由度３の$\chi^2$分布を適宜変数変換したものになります。