水素分子のエネルギーに現れる積分の計算

なぜ気体の水素は$\mathrm{H}{\phantom{X}}_2$なのか

今年は量子力学誕生100周年の年です。温故知新ということでこのような記事を書いてみようかと。

水素は気体として存在する際には1原子としては存在せず$\mathrm{H}{\phantom{X}}_2$になります。なぜでしょうか。それは原子2つの距離の関数として水素原子2体系（陽子2つ+電子2つ）のエネルギーを描くと、2原子の距離が原子の大きさ程度となるときにエネルギーが最小になるからです。

陽子は電子より2000倍も重いです。よって電子が運動する時間スケールと比較し、陽子の運動は非常に遅いです。そこで陽子の運動を無視します。これをBorn-Oppenheimer近似と呼びます。この近似では陽子は静止したプラス電荷のsourceとして扱われるため、系のエネルギーを陽子2体の距離の関数として表すことができます。電子が陽子間に存在すると陽子の正電荷により電子同士の斥力が遮蔽されてエネルギーが低くなります。よって電子は陽子同士の間に高い確率で存在するほうがエネルギー的に得であり、これが糊のような（？）役目をすることで水素原子2体は近い場所に存在することになります。

この記事では、Born-Oppenheimer近似および後に説明するHeitler-Londonの方法において、水素原子2体系のエネルギーに現れるすべての積分を解析的に実行します（指数積分だけ残りますが）。この計算は初等的な量子化学の教科書やネット上の記事にもある程度載っているのですが、「交換積分」と呼ばれる積分の解析的な計算はあまり載っていません。これをネット上に記しておくことはそれなりに意味があるのではないかと思います。交換積分の計算は丁寧に書いたので冗長に感じるかもしれませんがご容赦ください。

主な参考文献は[1][2][3]です。特にこの記事のメインである$V_{EX}^{(12)}$の計算は[3]を参考にしており、計算結果が当該文献と一致していることを確認しています。[4]はポテンシャルの計算と球関数に関して議論している数学書であり[3]で引用しています。[5][6]はSlaterの量子化学に関する教科書であり、関連事項に関する詳しい議論がなされていると思います（[6]に関しては読んでいません...）。[7]はスタンダードな量子力学の教科書です。[8]は楕円体座標(spheroidal coordinate)に関するWikipediaの記事です。

Formalism

まずは水素原子2体系のエネルギーの表式を導出します。

Born-Oppenheimer近似 (BO近似)

改めてBorn-Oppenheimer近似（以下BO近似と記す）に関して説明します。

まずnotationに関してまとめておきます：

水素原子2体系、すなわち陽子2つ+電子2つの系のHamiltonianは以下です：
$$\begin{array}{rl}& \hat{H}:={\hat{H}}_{\text{p}}+{\hat{H}}_{\text{e}}+{\hat{H}}_I,\\ \\ & \{\begin{array}{ll}{\displaystyle {\hat{H}}_{\text{p}}:=-\frac{{\hslash }^2}{2m_p}{\mathrm{\Delta }}_A-\frac{{\hslash }^2}{2m_p}{\mathrm{\Delta }}_B}& \text{...陽子の運動項}\\ {\displaystyle {\hat{H}}_{\text{e}}:=-\frac{{\hslash }^2}{2m_e}{\mathrm{\Delta }}_1-\frac{e^2}{r_{A1}}-\frac{{\hslash }^2}{2m_e}{\mathrm{\Delta }}_2-\frac{e^2}{r_{B2}}}& \text{...電子1,2の運動項+陽子Aと電子1のクーロン力+陽子Bと電子2のクーロン力}\\ {\displaystyle {\hat{H}}_I:=\frac{e^2}{r_{12}}-\frac{e^2}{r_{A2}}-\frac{e^2}{r_{B1}}}& \text{...電子1,2、陽子Aと電子2、陽子Bと電子1の間のクーロン力}\end{array}\end{array}$$
陽子の運動は電子の運動と比較して非常に遅いため、BO近似では${\hat{H}}_{\text{p}}$を無視します。改めてBO近似の下でのHamiltonian ${\hat{H}}_{\text{BO}}$を以下のように定義します：
$$\begin{array}{rl}& {\hat{H}}_{\text{BO}}:={\hat{H}}_{\text{e}}+{\hat{H}}_I,\\ \\ & \{\begin{array}{l}{\displaystyle {\hat{H}}_{\text{e}}:=-\frac{{\hslash }^2}{2m_e}{\mathrm{\Delta }}_1-\frac{e^2}{r_{A1}}-\frac{{\hslash }^2}{2m_e}{\mathrm{\Delta }}_2-\frac{e^2}{r_{B2}}}\\ {\displaystyle {\hat{H}}_I:=\frac{e^2}{r_{12}}-\frac{e^2}{r_{A2}}-\frac{e^2}{r_{B1}}}\end{array}\end{array}$$
この近似では、系の波動関数は電子$1,2$に関する波動関数で表されます。

Heitler-Londonの方法

さて、あとは${\hat{H}}_{\text{BO}}$のエネルギー固有値と固有関数を求めればよいのですが、2つの離れた電荷のsource（=2つの陽子）の下で電子2体の波動関数を解析的に求めるのは難しいです。そこで電子2体の基底状態（=エネルギー最小の状態）の波動関数に関して以下に示すような単純なものを仮定します。

陽子2つが非常に遠くに離れていれば、系のエネルギーが最小になるのは2つの水素原子が両方とも基底状態である場合です。水素原子の基底状態（1s状態）の波動関数は以下で与えられます：

ただし電子2体の場合、両者は区別がつかないという条件から、波動関数はラベル1,2の入れ替えに対し対称または反対称となります。ゆえに陽子A,Bが非常に遠くに離れている場合の系の基底状態の波動関数は

となります。そしてHeitler-Londonの方法では、実際の水素分子の波動関数もEq.(1)で表されることを仮定します（※この方法に関して、「まとめ」の章もご参照ください）。このような取り扱いは定性的にも定量的にもそれなりに良いことが知られています。そこで本記事では電子2体の波動関数をEq.(1)だと仮定します。そして水素分子のエネルギーをこの波動関数を使って計算します。

電子2体のスピン状態

ここでEq.(1)の波動関数の対称性と電子2体のスピンの関係に関して言及しておきます。電子はフェルミオンなので、電子$1,2$の波動関数は$1,2$のラベルの入れ替えで反対称になります。ただしここで言う波動関数は空間の波動関数とスピンの波動関数の直積です：
$$\begin{array}{r}{\psi }_{\mathrm{t}\mathrm{o}\mathrm{t}}(1,2)={\psi }_{\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{e}}({\mathit{r}}_1,{\mathit{r}}_2)\otimes {\psi }_{\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}(1,2)\end{array}$$
ここで空間の波動関数${\psi }_{\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{e}}$とはすなわちEq.(1)の波動関数です。一方でスピン波動関数${\psi }_{\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}$は

のように与えられます。そして${\psi }_{\mathrm{t}\mathrm{o}\mathrm{t}}(1,2)$が$1,2$の入れ替えで反対称なので

• ${\psi }_{\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{e}}$が対称な場合${\psi }_{\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}$は反対称$\to$スピン1重項状態
• ${\psi }_{\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{e}}$が反対称な場合${\psi }_{\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}$は対称$\to$スピン3重項状態

になります。そこでEq.(1)の2つの波動関数を改めて以下のように呼び直します：

ここでsubscriptの$s$はスピン1重項（$\text{spin }\underset{―}{s}\text{inglet}$）を、$t$はスピン3重項（$\text{spin }\underset{―}{t}\text{riplet}$）を表します。このように、粒子の統計性（=フェルミオンかボソンか）は、関係のないように見える自由度の間にも相関をもたらします。

エネルギーの具体的な表式

系のエネルギーはHamiltonianの期待値をEq.(2)(3)の波動関数により計算すればよいです。スピン1重項状態に対するエネルギーを$E_s$、3重項に対するエネルギーを$E_t$とすると
$$\begin{array}{}\text{(4)}& E_s(R)& =\int d^3{\mathit{r}}_1\int d^3{\mathit{r}}_2{u}_s^{\ast }({\mathit{r}}_1,{\mathit{r}}_2){\hat{H}}_{\text{BO}}{u}_s({\mathit{r}}_1,{\mathit{r}}_2),\text{(5)}& E_t(R)& =\int d^3{\mathit{r}}_1\int d^3{\mathit{r}}_2{u}_t^{\ast }({\mathit{r}}_1,{\mathit{r}}_2){\hat{H}}_{\text{BO}}{u}_t({\mathit{r}}_1,{\mathit{r}}_2)\end{array}$$
です。$R$は陽子間の距離です。

ここで${\hat{H}}_e$を
$$\begin{array}{rl}& {\hat{H}}_e={\hat{H}}_{A1}+{\hat{H}}_{B2},\\ \\ & \{\begin{array}{l}{\displaystyle {\hat{H}}_{A1}:=-\frac{{\hslash }^2}{2m_e}{\mathrm{\Delta }}_1-\frac{e^2}{r_{A1}},}\\ {\displaystyle {\hat{H}}_{B2}:=-\frac{{\hslash }^2}{2m_e}{\mathrm{\Delta }}_2-\frac{e^2}{r_{B2}}}\end{array}\end{array}$$
のようにわけておきます。$u_A({\mathit{r}}_1),u_B({\mathit{r}}_2)$をそれぞれ陽子$A$を核とする電子$1$の水素原子の基底状態の波動関数と電子$B$を核とする電子$2$の水素原子の基底状態の波動関数とすると
$$\begin{array}{r}{\hat{H}}_{A1}{u}_A({\mathit{r}}_1)=E_{1s}{u}_A({\mathit{r}}_1),{\hat{H}}_{B2}{u}_B({\mathit{r}}_2)=E_{1s}{u}_B({\mathit{r}}_2)\end{array}$$
が成立します。ここで$E_{1s}$は水素の基底状態である1s状態のエネルギーです。このことからEq.(4)(5)において${\hat{H}}_e$の寄与は$E_{1s}$に置き換わります。また$u_A,u_B$は規格化条件
$$\begin{array}{r}\int d\mathit{r}{u}_A^{\ast }(\mathit{r})u_A(\mathit{r})=\int d\mathit{r}{u}_B^{\ast }(\mathit{r})u_B(\mathit{r})=1\end{array}$$
を満たします。さらにEq.(2)(3)に現れる規格化因子$N_s,N_t$は$u_s,u_t$の規格化条件
$$\begin{array}{r}\int d\mathit{r}{u}_s^{\ast }(\mathit{r})u_s(\mathit{r})=\int d\mathit{r}{u}_t^{\ast }(\mathit{r})u_t(\mathit{r})=1\end{array}$$
より
$$\begin{array}{rl}& N_s(R)=\frac{1}{\sqrt{2(1+S(R{)}^2)}},N_t(R)=\frac{1}{\sqrt{2(1-S(R{)}^2)}},\\ & S(R):=\int d^3\mathit{r}{u}_A^{\ast }(\mathit{r})u_B(\mathit{r})\end{array}$$
となります。これらの量は陽子間距離$R$に依存することに注意してください。

以上を用いてEq.(4)(5)を計算すると

となります[7]。Eq.(6)はCoulomb積分、Eq.(7)は交換積分と呼ばれます。

楕円体座標（spheroidal coordinate）

さて$S,V_C,V_{EX}$を計算するのですが、それには楕円体座標(spheroidal coordinate)を使うのがよいです。この座標ではまず2次元における座標（ここではこれを$y\text{-}z$平面とする）を、図1の$A,B$からの距離$r_A,r_B$により表します。ここでは$A:(-R/2,0),B:(R/2,0)$とします（そしてこれらは静止した2つの陽子の位置とします）。そして

$$\begin{array}{r}\xi :=\frac{r_A+r_B}{R},\eta :=\frac{r_A-r_B}{R}\end{array}$$
を定義し、$z,y$を$\xi ,\eta$で表すと
$$\begin{array}{r}z=\frac{R}{2}\xi \eta ,y=\pm \frac{R}{2}\sqrt{({\xi }^2-1)(1-{\eta }^2)}\end{array}$$
となります。

楕円体座標の設定（2次元）

3次元の位置は、上の2次元平面を$z$軸を中心にしてぐりんと$\phi$だけ回して表します（図2）。改めて$xyz$軸を図2のように定義します。図1の平面は図2の点線で表された平面であるとし、また図1の$y$の値を図2では$r$とします。$\xi ,\eta ,\phi$を用いて$x,y,z$を表すと
$$\begin{array}{r}\{\begin{array}{l}{\displaystyle x=r\cos \phi ,}\\ {\displaystyle y=r\sin \phi ,}\\ {\displaystyle z=\frac{R}{2}\xi \eta ,}\end{array}r=\frac{R}{2}\sqrt{({\xi }^2-1)(1-{\eta }^2)}\end{array}$$
となります。

楕円体座標の設定（3次元）。$A,B$および青点が乗る平面は点線で示されている。この平面は$z$軸を回転軸として$x$軸から$\phi$だけ傾いている。

以下ではこの座標を用いて積分を行います。そのためにJacobianを計算しておきましょう。計算は単純なので、結果だけ記しておきます：
$$\begin{array}{r}dxdydz={\left(\frac{R}{2}\right)}^3({\xi }^2-{\eta }^2)d\xi d\eta d\phi \end{array}$$
また各変数がとりうる値の範囲は
$$\begin{array}{r}1\le \xi <+\mathrm{\infty },-1\le \eta \le 1,0\le \phi <2\pi \end{array}$$
です。

改めて、楕円体座標は以下のような座標系です[2]：

ちなみに公式3は[2]での定義と比較し$x$軸と$y$軸が入れ替わっています。世の定義では公式3のように定義しているものが多いみたいなのでそうしておきました[8]。この記事ではその違いが問題になることはないです。

積分を実行する

それでは$S,V_C,V_{EX}$に現れる積分を実行していきます。$V_C,V_{EX}$を
$$\begin{array}{rl}& V_C(R)=V_C^{(12)}(R)+V_C^{(B1)}(R)+V_C^{(A2)}(R),\\ & V_{EX}(R)=V_{EX}^{(12)}(R)+V_{EX}^{(B1)}(R)+V_{EX}^{(A2)}(R),\\ \\ & \{\begin{array}{ll}{\displaystyle V_C^{(12)}(R)}& {\displaystyle =\int d^3{\mathit{r}}_1\int d^3{\mathit{r}}_2{u}_A^{\ast }({\mathit{r}}_1)u_B^{\ast }({\mathit{r}}_2)\left(\frac{e^2}{r_{12}}\right)u_A({\mathit{r}}_1)u_B({\mathit{r}}_2),}\\ {\displaystyle V_C^{(B1)}(R)}& {\displaystyle =\int d^3{\mathit{r}}_1\int d^3{\mathit{r}}_2{u}_A^{\ast }({\mathit{r}}_1)u_B^{\ast }({\mathit{r}}_2)(-\frac{e^2}{r_{B1}})u_A({\mathit{r}}_1)u_B({\mathit{r}}_2),}\\ {\displaystyle V_C^{(A2)}(R)}& {\displaystyle =\int d^3{\mathit{r}}_1\int d^3{\mathit{r}}_2{u}_A^{\ast }({\mathit{r}}_1)u_B^{\ast }({\mathit{r}}_2)(-\frac{e^2}{r_{A2}})u_A({\mathit{r}}_1)u_B({\mathit{r}}_2),}\\ {\displaystyle V_{EX}^{(12)}(R)}& ={\displaystyle \int d^3{\mathit{r}}_1\int d^3{\mathit{r}}_2{u}_B^{\ast }({\mathit{r}}_1)u_A^{\ast }({\mathit{r}}_2)\left(\frac{e^2}{r_{12}}\right)u_A({\mathit{r}}_1)u_B({\mathit{r}}_2),}\\ {\displaystyle V_{EX}^{(B1)}(R)}& ={\displaystyle \int d^3{\mathit{r}}_1\int d^3{\mathit{r}}_2{u}_B^{\ast }({\mathit{r}}_1)u_A^{\ast }({\mathit{r}}_2)(-\frac{e^2}{r_{B1}})u_A({\mathit{r}}_1)u_B({\mathit{r}}_2),}\\ {\displaystyle V_{EX}^{(A2)}(R)}& ={\displaystyle \int d^3{\mathit{r}}_1\int d^3{\mathit{r}}_2{u}_B^{\ast }({\mathit{r}}_1)u_A^{\ast }({\mathit{r}}_2)(-\frac{e^2}{r_{A2}})u_A({\mathit{r}}_1)u_B({\mathit{r}}_2)}\end{array}\end{array}$$
のように分解しておきます。ただし対称性から
$$\begin{array}{r}{V}_C^{(B1)}=V_C^{(A2)},V_{EX}^{(B1)}=V_{EX}^{(A2)}\end{array}$$
です。

$S(R)$の計算

最も簡単なのは
$$\begin{array}{r}S(R)=\int d^3\mathit{r}{u}_A^{\ast }(\mathit{r})u_B(\mathit{r})\end{array}$$
の計算です。
$$\begin{array}{rl}S(R)& =\int d^3\mathit{r}{u}_A(\mathit{r})u_B(\mathit{r})\\ & =\left(\frac{1}{\pi a_0^3}\right)\int d^3\mathit{r}\exp (-\frac{r_A+r_B}{a_0})\end{array}$$
であり、
$$\begin{array}{r}\xi =\frac{r_A+r_B}{R},\eta =\frac{r_A-r_B}{R}\end{array}$$
とすると（図1の$A,B$の位置に陽子$A,B$が存在する。青い点の座標を$\mathit{r}$とする）
$$\begin{array}{r}S(R)=\left(\frac{1}{\pi a_0^3}\right)\times {\left(\frac{R}{2}\right)}^3\times {\int }_1^{\mathrm{\infty }}d\xi {\int }_{-1}^{1}d\eta {\int }_0^{2\pi }d\phi ({\xi }^2-{\eta }^2)\exp (-\frac{R}{a_0}\xi )\end{array}$$
となります。これは簡単に積分できて

です。

$V_C^{(A2)}(=V_C^{(B1)})$の計算

次に簡単なのは$V_C^{(A2)}$の計算です。この積分は${\mathit{r}}_1$に関しては$u_A$のnormalizationから1となります。残りの${\mathit{r}}_2$に関しては上と同じように
$$\begin{array}{r}\xi =(r_{A2}+r_{B2})/R,\eta =(r_{A2}-r_{B2})/R\end{array}$$
と座標を設定することで
$$\begin{array}{rl}{V}_C^{(A2)}& =\int d^3{\mathit{r}}_1{u}_A^{\ast }({\mathit{r}}_1)u_A({\mathit{r}}_1)\int d^3{\mathit{r}}_2{u}_B^{\ast }({\mathit{r}}_2)(-\frac{e^2}{r_{A2}})u_B({\mathit{r}}_2)\\ & =\int d^3{\mathit{r}}_2{u}_B^{\ast }({\mathit{r}}_2)(-\frac{e^2}{r_{A2}})u_B({\mathit{r}}_2)\\ & =-e^2\frac{2}{R}\left(\frac{1}{\pi a_0^3}\right)\times {\left(\frac{R}{2}\right)}^3{\int }_1^{\mathrm{\infty }}d\xi {\int }_{-1}^{1}d\eta {\int }_0^{2\pi }d\phi ({\xi }^2-{\eta }^2)\frac{1}{\xi +\eta }\exp (-\frac{R}{a_0}(\xi -\eta ))\\ & =-\frac{e^2}{2R}{D}^3{\int }_1^{\mathrm{\infty }}d\xi {\int }_{-1}^{1}d\eta (\xi -\eta )\exp (-D(\xi -\eta ))\end{array}$$
を得ます。これも簡単に計算できて

となります。

$V_{EX}^{(A2)}(=V_{EX}^{(B1)})$の計算

$V_{EX}^{(A2)}(=V_{EX}^{(B1)})$の計算も簡単です。
$$\begin{array}{r}{V}_{EX}^{(A2)}=\int d{\mathit{r}}_1\int d{\mathit{r}}_2{u}_B^{\ast }({\mathit{r}}_1)u_A^{\ast }({\mathit{r}}_2)(-\frac{e^2}{r_{A2}})u_A({\mathit{r}}_1)u_B(\mathit{r})\end{array}$$
ですが、$r_{A2}$は${\mathit{r}}_1$の座標を含まないので、${\mathit{r}}_1$の積分を実行すると$S$がfactor outされます：
$$\begin{array}{r}=-e^{2}S(R)\times \int d{\mathit{r}}_2{u}_A({\mathit{r}}_2)u_B({\mathit{r}}_2)\frac{1}{r_{A2}}\end{array}$$
あとは上の計算と同じように${\mathit{r}}_1$の楕円体座標を設定します：
$$\begin{array}{r}\xi =\frac{r_{A2}+r_{B2}}{R},\eta =\frac{r_{A2}-r_{B2}}{R}\end{array}$$
この座標のもとで$V_{EX}^{(A2)}$は
$$\begin{array}{r}{V}_{EX}^{(A2)}=\frac{2}{a_0^3}{\left(\frac{R}{2}\right)}^2\times {\int }_1^{\mathrm{\infty }}d\xi {\int }_{-1}^{1}d\eta (\xi -\eta )\exp (-\frac{R}{a_0}\xi )\end{array}$$
となります。この積分は今までと同様簡単に積分できて

となります。

$V_C^{(12)}$の計算

次に$V_C^{(12)}$の計算を行います。この積分の被積分関数の分母には$r_{12}$が存在し、これは${\mathit{r}}_1,{\mathit{r}}_2$の両方に依存するため今までのように${\mathit{r}}_1,{\mathit{r}}_2$のどちらかを単純に積分することはできません。しかしちょっとした考察により積分ができます。

まず図3のように座標を設定します。この図では$1,A,B$がなす平面と$2,A,B$がなす平面が一致しているように見えますが、ホントは（一般には）ずれています。

電子$1,2$に対する楕円体座標。$A,B$にそれぞれ陽子が存在する。

$V_C^{(12)}(R)$は
$$\begin{array}{r}{V}_C^{(12)}(R)=\int d^3{\mathit{r}}_1\int d^3{\mathit{r}}_2{u}_A^{\ast }({\mathit{r}}_1)u_B^{\ast }({\mathit{r}}_2)\left(\frac{e^2}{r_{12}}\right)u_A({\mathit{r}}_1)u_B({\mathit{r}}_2)\end{array}$$
ですが、先に${\mathit{r}}_2$の積分を行います。被積分関数で${\mathit{r}}_2$に依存するのは$r_{12},r_{B2}$（図3のオレンジ色の実線）です。${\mathit{r}}_1$は後で積分するため固定されています。図3の$r_{B1},r_{12}$を
$$\begin{array}{r}{r}_{B1}\to R(\text{fixed}),r_{12}\to r_{A2}\end{array}$$
と対応させると、$V_C^{(12)}(R)$の${\mathit{r}}_2$に関する積分は
$$\begin{array}{r}\int d^3{\mathit{r}}_2\frac{1}{r_{A2}}{e}^{-2r_{B2}/a_0}\end{array}$$
となり、これは$V_C^{(A2)}$の計算に現れた積分と一致します。この積分の結果より
$$\begin{array}{rl}\frac{e^2}{\pi a_0^3}\int \frac{1}{r_{12}}{e}^{-2r_{B2}/a_0}& =-(V_C^{(A2)}(R)\text{において}R\to r_{B1}\text{としたもの})\\ & =\frac{e^2}{r_{B1}}(1-e^{-2r_{B1}/a_0}-\frac{r_{B1}}{a_0}{e}^{-2r_{B1}/a_0})\end{array}$$
となります。よって
$$\begin{array}{r}{V}_C^{(12)}(R)=\frac{e^2}{\pi a_0^3}\int d^3{\mathit{r}}_1{e}^{-2r_{A1}/a_0}\times \frac{e^2}{r_{B1}}(1-e^{-2r_{B1}/a_0}-\frac{r_{B1}}{a_0}{e}^{-2r_{B1}/a_0})\end{array}$$
あとは今までと同様に楕円体座標を
$$\begin{array}{r}\xi =(r_{A1}+r_{B1})/R,\eta =(r_{A1}-r_{B1})/R\end{array}$$
と設定して積分すれば
$$\begin{array}{r}{V}_C^{(12)}(R)=\frac{e^4}{2a_0}{D}^2{\int }_1^{\mathrm{\infty }}d\xi {\int }_{-1}^{1}d\eta (\xi +\eta )\{1-e^{-D(\xi -\eta )}-\frac{D}{2}(\xi -\eta )e^{-D(\xi -\eta )}\}{e}^{-D(\xi +\eta )}\end{array}$$
この積分は（面倒ですが）困難なく積分できて

を得ます。

$V_{EX}^{(12)}$の計算

この計算が本記事のメインです。
$$\begin{array}{r}{V}_{EX}^{(12)}(R)=\int d^3{\mathit{r}}_1\int d^3{\mathit{r}}_2{u}_B^{\ast }({\mathit{r}}_1)u_A^{\ast }({\mathit{r}}_2)\left(\frac{e^2}{r_{12}}\right)u_A({\mathit{r}}_1)u_B({\mathit{r}}_2)\end{array}$$
の積分を実行します。$V_C^{(12)}$同様、座標を図3のように設定します。そして
$$\begin{array}{rl}& {\xi }_1:=(r_{A1}+r_{B1})/R,{\eta }_1:=(r_{A1}-r_{B1})/R,r_1:=\frac{R}{2}\sqrt{({\xi }_1^2-1)(1-{\eta }_1^2)},{\phi }_1:1,A,B\text{がなす平面の}x\text{軸からの角度}\\ & {\xi }_2:=(r_{A2}+r_{B2})/R,{\eta }_2:=(r_{A2}-r_{B2})/R,r_2:=\frac{R}{2}\sqrt{({\xi }_2^2-1)(1-{\eta }_2^2)},{\phi }_2:2,A,B\text{がなす平面の}x\text{軸からの角度}\end{array}$$
とします。$r_{12}$をこの座標で具体的に書き下すと
$$\begin{array}{r}{r}_{12}=\sqrt{r_1^2+r_2^2-2r_1{r}_2\cos ({\phi }_1-{\phi }_2)+\frac{R^2}{4}({\xi }_1{\eta }_1-{\xi }_2{\eta }_2{)}^2}\end{array}$$
となります。

$V_{EX}^{(12)}$の場合、被積分関数に${\mathit{r}}_2$に依存する変数が$r_{12},r_{B2}$だけでなく$u_A^{\ast }({\mathit{r}}_2)$からもたらされる$r_{A2}$も存在するので、$V_C^{(12)}$のようには積分することができません。

電磁気学においてこういう計算ではよく多重極展開を行いますが、実は楕円体座標でも多重極展開に類似した(?)以下の公式が成立します[3][4]：

この展開を用いると、$V_{EX}^{(12)}(R)$の被積分関数の中で$\phi$依存性をもつのは展開に現れる$\cos$の部分だけであり、${\int }_0^{2\pi }d{\phi }_1{\int }_0^{2\pi }d{\phi }_2$を行うと$\nu =0$の寄与のみが残ります。そして
$$\begin{array}{r} P_{\tau }^0=P_{\tau },Q_{\tau }^0=Q_{\tau },D_{\tau 0}=\frac{2\tau +1}{2}(P_{\tau },Q_{\tau }\text{はそれぞれ第1種および第2種のルジャンドル多項式})\end{array}$$
なので（Appendix参照のこと）、件の積分は
$$\begin{array}{rl}{V}_{EX}^{(12)}(R)& =\frac{e^2{D}^6}{16}{\int }_1^{\mathrm{\infty }}{e}^{-D{\xi }_1}d{\xi }_1{\int }_{-1}^1({\xi }_1^2-{\eta }_1^2)d{\eta }_1{\int }_1^{\mathrm{\infty }}{e}^{-D{\xi }_2}d{\xi }_2{\int }_{-1}^1({\xi }_2^2-{\eta }_2^2)\frac{1}{{\tilde{r}}_{12}}d{\eta }_2,\\ \\ \frac{1}{{\tilde{r}}_{12}}& =\frac{2}{a_{0}D}\sum _{\tau =0}^{\mathrm{\infty }}\frac{2\tau +1}{2}(P_{\tau }({\xi }_1)Q_{\tau }({\xi }_2)P_{\tau }({\eta }_1)P_{\tau }({\eta }_2)\theta ({\xi }_2-{\xi }_1)\\ & +P_{\tau }({\xi }_2)Q_{\tau }({\xi }_1)P_{\tau }({\eta }_1)P_{\tau }({\eta }_2)\theta ({\xi }_1-{\xi }_2))\end{array}$$
となります。以上から$V_{EX}^{(12)}(R)$は
$$\begin{array}{rl}{V}_{EX}^{(12)}(R)=& \frac{e^2{D}^6}{16}\times \frac{2}{a_{0}D}\sum _{\tau =0}^{\mathrm{\infty }}\frac{2\tau +1}{2}\times {\int }_1^{\mathrm{\infty }}{e}^{-D{\xi }_1}d{\xi }_1{\int }_{-1}^1({\xi }_1^2-{\eta }_1^2)d{\eta }_1\\ & \times [{\int }_1^{{\xi }_1}{e}^{-D{\xi }_2}d{\xi }_2{P}_{\tau }({\xi }_2)Q_{\tau }({\xi }_1)P_{\tau }({\eta }_1)\times {\int }_{-1}^1({\xi }_2^2-{\eta }_2^2)P_{\tau }({\eta }_2)d{\eta }_2\\ & +{\int }_{{\xi }_1}^{\mathrm{\infty }}{e}^{-D{\xi }_2}d{\xi }_2{P}_{\tau }({\xi }_1)Q_{\tau }({\xi }_2)P_{\tau }({\eta }_1)\times {\int }_{-1}^1({\xi }_2^2-{\eta }_2^2)P_{\tau }({\eta }_2)d{\eta }_2]\end{array}$$
になります。ここで${\eta }_2$の積分を実行すると
$$\begin{array}{rl}{\int }_{-1}^1({\xi }_2^2-{\eta }_2^2)P_{\tau }({\eta }_2)d{\eta }_2& =\{\begin{array}{ll}{\displaystyle 2({\xi }_2^2-\frac{1}{3})}& \text{for }\tau =0\\ {\displaystyle -\frac{4}{15}}& \text{for }\tau =2\\ 0& \text{other }{\tau }^{\prime }s\end{array}\\ \\ & =2({\xi }_2^2-\frac{1}{3}){\delta }_{\tau 0}-\frac{4}{15}{\delta }_{\tau 2}\end{array}$$
となるので、のこるのは$\tau =0,2$のみです。$Q_0,Q_2,P_0,P_2$は
$$\begin{array}{rl}{Q}_0(x)& =\log \frac{x+1}{x-1},Q_2(x)=-3x+P_2(x)\log \frac{x+1}{x-1},\\ P_0(x)& =1,P_2(x)=\frac{3x^2-1}{2}\end{array}$$
です（Appendix参照）。これらを$V_{EX}^{(12)}(R)$に代入して
$$\begin{array}{rl}{V}_{EX}^{(12)}(R)=\frac{e^2{D}^5}{8a_0}\times \{& \frac{1}{2}{\int }_1^{\mathrm{\infty }}{e}^{-D{\xi }_1}d{\xi }_1{\int }_{-1}^1({\xi }_1^2-{\eta }_1^2)d{\eta }_1{\int }_1^{{\xi }_1}{e}^{-D{\xi }_2}d{\xi }_2\times 2({\xi }_2^2-\frac{1}{3})\times P_0({\xi }_2)Q_0({\xi }_1)P_0({\eta }_1)\\ & +\frac{5}{2}{\int }_1^{\mathrm{\infty }}{e}^{-D{\xi }_1}d{\xi }_1{\int }_{-1}^1({\xi }_1^2-{\eta }_1^2)d{\eta }_1{\int }_1^{{\xi }_1}{e}^{-D{\xi }_2}d{\xi }_2\times (-\frac{4}{15})\times P_2({\xi }_2)Q_2({\xi }_1)P_2({\eta }_1)\\ & +\frac{1}{2}{\int }_1^{\mathrm{\infty }}{e}^{-D{\xi }_1}d{\xi }_1{\int }_{-1}^1({\xi }_1^2-{\eta }_1^2)d{\eta }_1{\int }_{{\xi }_1}^{\mathrm{\infty }}{e}^{-D{\xi }_2}d{\xi }_2\times 2({\xi }_2^2-\frac{1}{3})\times P_0({\xi }_1)Q_0({\xi }_2)P_0({\eta }_1)\\ & +\frac{5}{2}{\int }_1^{\mathrm{\infty }}{e}^{-D{\xi }_1}d{\xi }_1{\int }_{-1}^1({\xi }_1^2-{\eta }_1^2)d{\eta }_1{\int }_{{\xi }_1}^{\mathrm{\infty }}{e}^{-D{\xi }_2}d{\xi }_2\times (-\frac{4}{15})\times P_2({\xi }_1)Q_2({\xi }_2)P_2({\eta }_1)\}\\ =\frac{e^2{D}^5}{8a_0}\times [& \frac{1}{2}{\int }_1^{\mathrm{\infty }}{e}^{-D{\xi }_1}d{\xi }_1{\int }_{-1}^1({\xi }_1^2-{\eta }_1^2)d{\eta }_1{\int }_1^{{\xi }_1}{e}^{-D{\xi }_2}d{\xi }_2\times 2({\xi }_2^2-\frac{1}{3})\times \log \frac{{\xi }_1+1}{{\xi }_1-1}\\ & +\frac{5}{2}{\int }_1^{\mathrm{\infty }}{e}^{-D{\xi }_1}d{\xi }_1{\int }_{-1}^1({\xi }_1^2-{\eta }_1^2)d{\eta }_1{\int }_1^{{\xi }_1}{e}^{-D{\xi }_2}d{\xi }_2\times (-\frac{4}{15})\times \frac{3{\xi }_2^2-1}{2}\{-3{\xi }_1+P_2({\xi }_1)\log \frac{{\xi }_1+1}{{\xi }_1-1}\}\times \frac{3{\eta }_1^2-1}{2}\\ & +\frac{1}{2}{\int }_1^{\mathrm{\infty }}{e}^{-D{\xi }_1}d{\xi }_1{\int }_{-1}^1({\xi }_1^2-{\eta }_1^2)d{\eta }_1{\int }_{{\xi }_1}^{\mathrm{\infty }}{e}^{-D{\xi }_2}d{\xi }_2\times 2({\xi }_2^2-\frac{1}{3})\times \log \frac{{\xi }_2+1}{{\xi }_2-1}\\ & +\frac{5}{2}{\int }_1^{\mathrm{\infty }}{e}^{-D{\xi }_1}d{\xi }_1{\int }_{-1}^1({\xi }_1^2-{\eta }_1^2)d{\eta }_1{\int }_{{\xi }_1}^{\mathrm{\infty }}{e}^{-D{\xi }_2}d{\xi }_2\times (-\frac{4}{15})\times \frac{3{\xi }_1^2-1}{2}\{-3{\xi }_2+P_2({\xi }_2)\log \frac{{\xi }_2+1}{{\xi }_2-1}\}\times \frac{3{\eta }_1^2-1}{2}]\end{array}$$
を得ます。