永久双極子の相互作用エネルギー

はじめに

永久双極子の相互作用が距離の6乗に反比例するという話があります。Lennard-Jonesエネルギーの$1/r^6$ の項を説明する物理の一つになるらしい（他には誘起双極子とか）。
これを確認したいと思います。

Keesom相互作用という名前もあるらしいです。

注意

たぶん計算間違いしています。ただ結果は時々見かけるものと同じで、 $1/r^6$ に比例することは導けましたので、いったん公開します。

ということで、以下の2点を募集しますので、もしお時間ある方がいましたら考えていただけたらありがたいです。

間違い探し

熱平均を取るところで、テイラー展開の一般項に間違いがあると思います。

もっとスマートな数学

$U(r, \Omega)$ を計算するところで、今回は高校数学に毛が生えたようなベクトル演算を使いましたが、もっとスマートな数学がある気がしています。

話の流れ

まず、2つの双極子の相互作用エネルギーが位置・角度の関数としてどう書けるかを確かめます。
続いて、角度に関してボルツマン分布でエネルギーの熱平均を調べてみます。

準備

ベクトル $\vec{x}$ を引数にとり実数を返す関数 $U(\vec{x}) = f(\vec{x})$ のテイラー展開を思い出します。

まず、位置 $\vec{x}$ での勾配はベクトルで

$$\vec{G}(\vec{x}) = \nabla f(\vec{x})$$

とかけます。またヘシアン $H(\vec{x})$ は行列で、行 $i$ 列 $j$ の成分は

$$H_{ij}(\vec{x}) = \frac{\partial^2}{\partial x_i\partial x_j}f(\vec{x})$$

とかきます。すると微小ベクトル $\vec{a}$ を与えると

$$f(\vec{x} + \vec{a}) \approx f(\vec{x}) + \vec{G}(\vec{x})^\mathrm{T}\vec{a} + \frac{1}{2}\vec{a}^\mathrm{T}H(\vec{x})\vec{a}$$

具体例

$$f(\vec{x}) = \frac{1}{\sqrt{\left(\vec{x} + \vec{a}\right)^2}}$$

$$\vec{G}(\vec{x}) = \frac{\vec{x} + \vec{a}}{\sqrt{\left(\vec{x} + \vec{a}\right)^2}^3}$$

ヘシアンについては、対角成分

$$H_{ii} = -\frac{1}{\sqrt{\left(\vec{a} + \vec{x}\right)^2}^3} + \frac{3}{2} \frac{(a_i + x_i)(2a_i + 2x_i)}{\sqrt{\left(\vec{a} + \vec{x}\right)^2}^5}$$

と非対角成分

$$H_{ij} = \frac{3}{2} \frac{(a_i + x_i)(2a_j + 2x_j)}{\sqrt{\left(\vec{a} + \vec{x}\right)^2}^5}$$

を合わせてまとめると、ヘシアンは

$$H_{ij}= -\frac{\delta_{ij}}{\sqrt{\left(\vec{a} + \vec{x}\right)^2}^3} + \frac{3(a_i + x_i)(a_j + x_j)}{\sqrt{\left(\vec{a} + \vec{x}\right)^2}^5}$$

となります。なので、 $f(\vec{x} + \vec{a})$ の $\vec{x}$ 周りのテーラー展開は

$$\frac{1}{\sqrt{\left(\vec{x} + \vec{a}\right)^2}} \approx \frac{1}{|\vec{x}|} -\frac{\vec{x}\cdot\vec{a}}{|\vec{x}|^3} -\frac{1}{2}\frac{\vec{a}^2}{|\vec{x}|^3} +\frac{3}{2}\frac{\left(\vec{x}\cdot\vec{a}\right)^2}{|\vec{x}|^5}$$

と書けます。

双極子相互作用（絶対零度）

要は熱揺動による分布を考えない、位置を固定した場合。

状況設定

双極子が2つあって、一つは $\vec{\mu}_1$ 、もう一つは $\vec{\mu}_2$ 。大きさを評価するときは $|\vec{\mu}_i| = \mu_i$ とベクトル記号をとったものとします。記号の枯渇、というよりも自分の発想の貧困。。。そして、双極子1から2を見た時のベクトルを $\vec{r}$ とします。大きさ $r = |\vec{r}|$ は双極子間の距離になります。

さて、双極子のままではエネルギーを具体的に計算するすべを持ちません、というかそれをこれからやります。なので、双極子を点電荷の形で近似します。具体的には、１つの双極子に対して2つの点電荷を少し距離をごくわずか離して双極子の方向に並べます。2つの点電荷は正負が逆で大きさが同じです。これで遠くから見たら電荷は0だが電場を発生する双極子として扱えます。双極子 $i$ に対して2つの電荷を $\pm q_i$ とします。$-q_i$ からみた $+q_i$ の位置ベクトルを$\vec{h}_i$ と書きます。符号用に $\sigma_i = \pm 1$ と書きます。

ベクトルでは角度に関する平均は取りにくいので、角度表記もしてみます。$\vec{r}$ の方向に単位ベクトル $\vec{u} = \vec{r}/r$ をとります。これと $\vec{\mu}_i (i = 1, 2)$ との角度を $\theta_i$ とします。また、$\vec{u}$ の方向を軸にしたときの回転で、 $\vec{\mu}_1$ と $\vec{\mu}_2$ の回転角度の差を $\phi$ とします。これらの角度をまとめて $\Omega = (\theta_1, \theta_2, \phi)$ と表記します。

角度の自由度を2つもつ双極子が2つで自由度は4つかと思いきや、 $\vec{u}$ 周りに2つが同じように回転した系は同じものとみなし、1自由度引きまして、角度については自由度3になります。($\phi = \phi_2 - \phi_1$ のイメージです。)

ほしいもの

$$U(r, \Omega)$$

計算

点電荷の相互作用

図で分かるように、この状況では点電荷が4つあります。点電荷同士の相互作用を考えますが、ここで同じ双極子同士の相互作用、例えば $+q_1$ と $-q_1$ の相互作用はないものとします。
すると $(+q_1, +q_2)$ $(+q_1, -q_2)$ $(-q_1, +q_2)$ $(-q_1, -q_2)$ の4つの相互作用を考えることになります。双極子相互作用のエネルギーはこれの合計です。
$\sigma_1q_1$ から $\sigma_2q_2$ への位置ベクトルを

$$\vec{r}_{\sigma_1\sigma_2} = \vec{r} - \sigma_1\frac{1}{2}\vec{h}_1 + \sigma_2\frac{1}{2}\vec{h}_2$$

とかくと、相互作用エネルギーは

$$U = \sum_{\sigma_1 = \pm}\sum_{\sigma_2 = \pm} \frac{(\sigma_1q_1)(\sigma_2q_2)}{4\pi\epsilon}\frac{1}{\sqrt{\left(\vec{r}_{\sigma_1\sigma_2}\right)^2}}$$

となります。

展開

さて、 $|\vec{r}| >> |\vec{h}_i|$ でした。 $- \sigma_1\frac{1}{2}\vec{h}_1 + \sigma_2\frac{1}{2}\vec{h}_2$ を微小量 $\vec{a}$ と思って展開すると、

$$U \approx \sum_{\sigma_1 = \pm}\sum_{\sigma_2 = \pm} \frac{(\sigma_1q_1)(\sigma_2q_2)}{4\pi\epsilon} \left( \frac{1}{r} -\frac{\vec{r}\cdot\vec{a}}{r^3} -\frac{1}{2}\frac{\vec{a}^2}{r^3} +\frac{3}{2}\frac{(\vec{r}\cdot\vec{a})^2}{r^5} \right)$$

となります。さて、各項ごとに見ていきます。

$$\sum_{\sigma_1 = \pm}\sum_{\sigma_2 = \pm} \frac{(\sigma_1q_1)(\sigma_2q_2)}{4\pi\epsilon r} = \left(+1 -1 -1 +1\right)\frac{q_1q_2}{4\pi\epsilon r}$$

より、 $1/r$ の項は消えます。同様に、

$$\vec{r}\cdot\vec{a} = -\frac{1}{2}\sigma_1\vec{r}\cdot\vec{h}_1 +\frac{1}{2}\sigma_2\vec{r}\cdot\vec{h}_2$$

$$\vec{a}^2 = \frac{1}{4}\left(\vec{h}_1^2 + 2\sigma_1\sigma_2\vec{h}_1\vec{h}_2 + \vec{h}_2^2\right)$$

$$\left(\vec{r}\cdot\vec{a}\right)^2 = \frac{1}{4}\sigma_1^2\left(\vec{r}\cdot\vec{h}_1\right)^2 -\frac{1}{2}\sigma_1\sigma_2\left(\vec{r}\cdot\vec{h}_1\right)\left(\vec{r}\cdot\vec{h}_2\right) +\frac{1}{4}\sigma_2^2\left(\vec{r}\cdot\vec{h}_2\right)$$

より、$\sum_{\sigma_1 = \pm}\sum_{\sigma_2 = \pm}$ で0にならない項、つまり $\sigma_1$ と $\sigma_2$ がともに偶数乗となっているのは $\vec{a}^2$ と $\left(\vec{r}\cdot\vec{a}\right)^2$ の2番目の項のみとなります。
これらを持ってくると$\sum_{\sigma_1 = \pm}\sum_{\sigma_2 = \pm}1 = 4$ を考え、

$$U = \frac{q_1q_2}{4\pi\epsilon} \left( -\frac{\vec{h}_1\cdot\vec{h}_2}{r^3} +\frac{3\left(\vec{r}\cdot\vec{h}_1\right)\left(\vec{r}\cdot\vec{h}_2\right)}{r^5} \right)$$

とまとめられ、さらに $q_i \vec{h}_i = \vec{\mu}_i$ を思い出すと、

$$U = \frac{1}{4\pi\epsilon r^3} \left( -\vec{\mu}_1\cdot\vec{\mu}_2 +\frac{3\left(\vec{r}\cdot\vec{\mu}_1\right)\left(\vec{r}\cdot\vec{\mu}_2\right)}{r^2} \right)$$

とかけます。これを続く節で角度 $\Omega$ を使って表してみます。

角度

角度の定義をベクトルの言葉で書き直します。ここで双極子1から2への方向の単位ベクトルを $\vec{u} = \vec{r}/r$ としたことを思い出します。

双極子を結んだ方向 $\vec{u}$ と双極子のなす角 $\theta_i$

$$\vec{\mu}_i\cdot\vec{u} = \mu_i \cos\theta_i$$

双極子を結んだ方向 $\vec{u}$ 周りの回転角の差 $\phi$

双極子方向を向き $\vec{u}$ と垂直なベクトルをクロス積で得ます。

$$\left(\vec{\mu}_1\times\vec{u}\right)\cdot \left(\vec{\mu}_2\times\vec{u}\right) = \mu_1\mu_2\sin\theta_1\sin\theta_2\cos\phi$$

スカラー4重積

さらにここでベクトルの積に関する定理を持ってきます。

$$(\vec{A}\times\vec{B})\cdot(\vec{C}\times\vec{D})=(\vec{A}\cdot\vec{C})(\vec{B}\cdot\vec{D}) - (\vec{B}\cdot\vec{C})(\vec{A}\cdot\vec{D})$$

まとめ

スカラー4重積の左辺と $\phi$ 定義の左辺が同じ形をしていますので放り込んでやると、$\vec{u}\cdot\vec{u} = 1$ と $\theta_i$ の定義より

$$\left(\vec{\mu}_1\times\vec{u}\right)\cdot \left(\vec{\mu}_2\times\vec{u}\right) = \vec{\mu}_1\cdot\vec{\mu}_2 - (\vec{\mu}_1\cdot\vec{u} )(\vec{\mu}_1\cdot\vec{u}) =\vec{\mu}_1\cdot\vec{\mu}_2 - \mu_1\mu_2\cos\theta_1\cos\theta_2$$

となり、$\phi$ の定義の右辺を使いますと

$$\vec{\mu}_1\cdot\vec{\mu}_2 = \mu_1\mu_2\cos\theta_1\cos\theta_2 + \mu_1\mu_2\sin\theta_1\sin\theta_2\cos\phi$$

$$\frac{3\left(\vec{r}\cdot\vec{\mu}_1\right)\left(\vec{r}\cdot\vec{\mu}_2\right)}{r^2} = 3\left(\vec{u}\cdot\vec{\mu}_1\right)\left(\vec{u}\cdot\vec{\mu}_2\right) = 3\mu_1\mu_2\cos\theta_1\cos\theta_2$$

とベクトル内積を角度で書けましたので。まとめると $U(r, \Omega)$ をベクトルなしで書けます。

結論

$$U(r, \Omega) = -\frac{\mu_1\mu_2}{4\pi\epsilon r^3}\left(2\cos\theta_1\cos\theta_2 - \sin\theta_1\sin\theta_2\cos\phi\right)$$

熱平均

状況

位置は固定しているが、双極子はぐるぐる回転しているとします。で、その角度はボルツマン分布しているとします。

表記を楽にします。

$$f(\Omega) = 2\cos\theta_1\cos\theta_2 - \sin\theta_1\sin\theta_2\cos\phi$$

$$g(r) = - \frac{\mu_1\mu_2}{4\pi \epsilon r^3}$$

とすると

$$U(r, \Omega) = g(r)f(\Omega)$$

です。

ほしいもの

$$\left<U\right>(r) = \frac{\int_0^\pi\sin\theta_1d\theta_1\int_0^\pi \sin\theta_2d\theta_2\int_0^{2\pi} d\phi U(r, \Omega)\exp(-\beta U(r, \Omega))} {\int_0^\pi\sin\theta_1d\theta_1\int_0^\pi \sin\theta_2d\theta_2\int_0^{2\pi} d\phi \exp(-\beta U(r, \Omega))}$$

ここで、 $\theta_i$ については緯度に関する積分と同じように、微小体積 $d\Omega$ は $\theta_i$ に依存します。 なので $\sin\theta_i d\theta_i$ が積分にかかってきています。 $\phi$ は経度と同じようなものなので$d\phi$ だけになります。

方針

ボルツマン因子 $\exp\left(-\beta U\right)$ をテイラー展開して $\int$ と $\sum$ を入れ替え、2次以上を落とします。$\exp$ を展開してるから絶対収束してると思います。

計算

ここから間違いがあるかもしれない。

展開

$$\exp\left(-\beta f g\right) = 1 - \beta f g + \frac{1}{2!}(\beta f g)^2 - \frac{1}{3!}(\beta f g)^3 + \cdots$$

これを使うと分母は $\int$ と $\sum$ を入れ替え

$$分母 = \sum_{n=0}^\infty\frac{1}{n!}(-\beta g)^n\int_0^\pi \sin\theta_1d\theta_1\int_0^\pi\sin\theta_2d\theta_2\int_0^{2\pi}d\phi(f(\Omega))^n = \sum_{n=0}^\infty\frac{1}{n!}(-\beta g)^n I_n$$

と書けます。ここで角度の関数についての第$n$項の積分を記号$I_n$にしました。

$$I_n = \int_0^\pi \sin\theta_1d\theta_1\int_0^\pi\sin\theta_2d\theta_2\int_0^{2\pi}d\phi(f(\Omega))^n$$

一方で分子は分母を $\beta$ で微分するものとしてもいいですが、単に $U = fg$ をかけているだけなので、同じように展開してしまいます。

$$分子 = -\frac{1}{\beta}\sum_{n = 1}^\infty \frac{1}{(n -1)!}(-\beta g)^n I_n$$

一般項

ここから頑張って計算します。 $(f(\Omega))^n$ の係数を考えます。

$$I_n = \int_0^\pi \sin\theta_1d\theta_1\int_0^\pi\sin\theta_2d\theta_2\int_0^{2\pi}d\phi\sum_{i=0}^n\ _nC_i(2\cos\theta_1\cos\theta_2)^{n-i}(-\sin\theta_1\sin\theta_2\cos\phi)^i$$

積分と和を入れ替えつつ、まとめます。

$$I_n = \sum_{i=0}^n\ _nC_i(-1)^i2^{n-i} \left(\int_0^\pi d\theta_1\sin^{i+1}\theta_1\cos^{n-i}\theta_1\right) \left(\int_0^\pi d\theta_2\sin^{i+1}\theta_2\cos^{n-i}\theta_2\right) \left(\int_0^{2\pi} d\phi \cos^i\phi\right)$$

三角関数の累乗の積分はベータ関数になります。

$$I_n = \sum_{i=0}^n\ _nC_i (-1)^i2^{n-i}\left[\frac{(-1)^{n-i} + 1}{2}B\left(\frac{i+2}{2}, \frac{n - i + 1}{2}\right)\right]^2 \left[\frac{(-1)^i + 1}{2}B\left(\frac{1}{2}, \frac{i + 1}{2}\right)\right]$$

なんか $(-1)^i + 1$ は $i$ の偶奇で 0 になるかどうかが決まります。まず $i$ は偶数のみ、 $n - i$ も偶数のみという条件から $n$ も偶数のみに限られます。これから、仮変数を書き換えますと

$$I_{2n} = \sum_{i=0}^n\ _{2n}C_{2i}2^{2n - 2i}\left[B\left(i + 1, n - i + \frac{1}{2}\right)\right]^2\left[B\left(\frac{1}{2}, i + \frac{1}{2}\right)\right]$$

で、分解するべくガンマ関数で書いてやると

$$I_{2n} = \sum_{i=0}^n\ _{2n}C_{2i} 2^{2n-2i} \left[\frac{\Gamma\left(i + 1\right)\Gamma\left(n - i + \frac{1}{2}\right)}{\Gamma\left(n + 1 + \frac{1}{2}\right)}\right]^2 \left[\frac{\Gamma\left(\frac{1}{2}\right)\Gamma\left(i + \frac{1}{2}\right)}{\Gamma\left(i + 1\right)}\right]$$

となりますが、引数が半整数か整数だから飛び階乗と $\pi$ で書けて、さらに軽く整理すると

$$I_{2n} = \sum_{i=0}^n\ _{2n}C_{2i} 2^{2n-2i} (i!) \left[2^{i+1}\frac{(2n - 2i - 1)!!}{(2n + 1)!!}\right]^2 \left[\pi\frac{(2i - 1)!!}{2^i}\right]$$

で、思いつく限りできるだけ整理すると

$$I_{2n} = \frac{2^{2n + 4}\pi}{(2n + 1)(2n + 2)}\frac{\left[(n + 1)!\right]^2}{(2n + 2)!} \sum_{i=0}^n\frac{(2n - 2i)!}{\left[(n-i)!\right]^2}$$

打ち切り

とりあえず$300 \mathrm{K} = 2.4942 \mathrm{kJ/mol} = 208.5 \mathrm{cm}^{-1}$ ということで、無理があると思うのですが、
$\beta U(r, \Omega) << 1$ として 1次の項までで打ち切ると

$$I_0 = \frac{2^4\pi}{1*2}\frac{1}{2!}\frac{0!}{0!} = 4\pi$$

$$I_2 = \frac{2^6\pi}{3*4}\frac{(2!)^2}{(4)!}\left[\frac{2!}{1^2} + \frac{0!}{0!}\right]=\frac{8}{3}\pi$$

$$分子 \approx gI_1 - \beta g^2 I_2 = -\beta\left(-\frac{\mu_1\mu_2}{4\pi\epsilon r^3}\right)^2\frac{8}{3}\pi$$

$$分母 \approx I_0 - \beta gI_1= 4\pi$$

結論

以上をまとめると

$$\left<U\right> = -\frac{2}{3k_BT}\left(\frac{\mu_1\mu_2}{4\pi\epsilon}\right)^2\frac{1}{r^6}$$

となり、 $1/r^6$ に比例することがわかりました。

間違い可能性箇所

$I_0$ が一般項から計算すると $4\pi$ ですが、これは一般項ではなくそのまま計算すると

$$\left(\int_0^\pi d\theta \sin\theta \right)^2\int_0^{2\pi}d\phi = \left(\left[-\cos\theta\right]_0^\pi\right)^2\left[\phi\right]_0^{2\pi}=2^2\times 2\pi = 8\pi$$

です。おそらく一般項の $2^{2n + 4}$ の指数が1足りないのかなと思いますが、自分では見つけられませんでした。

あとがき

これ計算してまとめるのにどんだけ時間かかったか……しかも間違えてるし……もうやりたくない。