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

なぜ気体の水素は$\ce{H2}$なのか

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

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

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

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

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

Formalism

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

Born-Oppenheimer近似 (BO近似)

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

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

水素原子2体系、すなわち陽子2つ+電子2つの系のHamiltonianは以下です：
\begin{align} &\hat H:=\hat H_{\text{p}}+\hat H_{\text{e}}+\hat H_I, \\ \\ & \begin{cases} \displaystyle \hat H_{\text{p}}:=-\frac{\hbar^2}{2m_p}\Delta_A -\frac{\hbar^2}{2m_p}\Delta_B &\text{...陽子の運動項} \\ \displaystyle \hat H_{\text{e}}:=-\frac{\hbar^2}{2m_e}\Delta_1-\frac{e^2}{r_{A1}} -\frac{\hbar^2}{2m_e}\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{cases} \end{align}
陽子の運動は電子の運動と比較して非常に遅いため、BO近似では$\hat H_{\text{p}}$を無視します。改めてBO近似の下でのHamiltonian $\hat H_{\text{BO}}$を以下のように定義します：
\begin{align} &\hat H_{\text{BO}}:=\hat H_{\text{e}}+\hat H_I, \\ \\ & \begin{cases} \displaystyle \hat H_{\text{e}}:=-\frac{\hbar^2}{2m_e}\Delta_1-\frac{e^2}{r_{A1}} -\frac{\hbar^2}{2m_e}\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{cases} \end{align}
この近似では、系の波動関数は電子$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{align} \psi_{\rm tot}(1,2)=\psi_{\rm space}(\boldsymbol{r}_1,\boldsymbol{r}_2)\otimes\psi_{\rm spin}(1,2) \end{align}
ここで空間の波動関数$\psi_{\rm space}$とはすなわちEq.(1)の波動関数です。一方でスピン波動関数$\psi_{\rm spin}$は

のように与えられます。そして$\psi_{\rm tot}(1,2)$が$1,2$の入れ替えで反対称なので

• $\psi_{\rm space}$が対称な場合$\psi_{\rm spin}$は反対称$\to$スピン1重項状態
• $\psi_{\rm space}$が反対称な場合$\psi_{\rm spin}$は対称$\to$スピン3重項状態

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

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

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

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

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

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

となりますIgi。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{align} \xi:=\frac{r_A+r_B}{R}, \quad \eta:=\frac{r_A-r_B}{R} \end{align}
を定義し、$z,y$を$\xi,\eta$で表すと
\begin{align} z=\frac{R}{2}\xi\eta, \quad y=\pm\frac{R}{2}\sqrt{(\xi^2-1)(1-\eta^2)} \end{align}
となります。

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

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

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

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

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

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

積分を実行する

それでは$S, V_C, V_{EX}$に現れる積分を実行していきます。$V_C, V_{EX}$を
\begin{align} &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{cases} \displaystyle V_C^{(12)}(R)& \displaystyle=\int d^3\boldsymbol{r}_1\int d^3 \boldsymbol{r}_2 \ u^*_A(\boldsymbol{r}_1)u^*_B(\boldsymbol{r}_2) \left(\frac{e^2}{r_{12}}\right) u_A(\boldsymbol{r}_1)u_B(\boldsymbol{r}_2), \\ \displaystyle V_C^{(B1)}(R)& \displaystyle=\int d^3\boldsymbol{r}_1\int d^3 \boldsymbol{r}_2 \ u^*_A(\boldsymbol{r}_1)u^*_B(\boldsymbol{r}_2) \left(-\frac{e^2}{r_{B1}}\right) u_A(\boldsymbol{r}_1)u_B(\boldsymbol{r}_2), \\ \displaystyle V_C^{(A2)}(R)& \displaystyle=\int d^3\boldsymbol{r}_1\int d^3 \boldsymbol{r}_2 \ u^*_A(\boldsymbol{r}_1)u^*_B(\boldsymbol{r}_2) \left(-\frac{e^2}{r_{A2}}\right) u_A(\boldsymbol{r}_1)u_B(\boldsymbol{r}_2), \\ \displaystyle V_{EX}^{(12)}(R) &=\displaystyle \int d^3\boldsymbol{r}_1\int d^3 \boldsymbol{r}_2 \ u^*_B(\boldsymbol{r}_1)u^*_A(\boldsymbol{r}_2) \left(\frac{e^2}{r_{12}}\right) u_A(\boldsymbol{r}_1)u_B(\boldsymbol{r}_2), \\ \displaystyle V_{EX}^{(B1)}(R) &=\displaystyle \int d^3\boldsymbol{r}_1\int d^3 \boldsymbol{r}_2 \ u^*_B(\boldsymbol{r}_1)u^*_A(\boldsymbol{r}_2) \left(-\frac{e^2}{r_{B1}}\right) u_A(\boldsymbol{r}_1)u_B(\boldsymbol{r}_2), \\ \displaystyle V_{EX}^{(A2)}(R) &=\displaystyle \int d^3\boldsymbol{r}_1\int d^3 \boldsymbol{r}_2 \ u^*_B(\boldsymbol{r}_1)u^*_A(\boldsymbol{r}_2) \left(-\frac{e^2}{r_{A2}}\right) u_A(\boldsymbol{r}_1)u_B(\boldsymbol{r}_2) \end{cases} \end{align}
のように分解しておきます。ただし対称性から
\begin{align} V_C^{(B1)}=V_C^{(A2)}, \quad V_{EX}^{(B1)}=V_{EX}^{(A2)} \end{align}
です。

$S(R)$の計算

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

です。

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

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

となります。

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

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

となります。

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

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

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

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

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

を得ます。

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

この計算が本記事のメインです。
\begin{align} V_{EX}^{(12)}(R) = \int d^3\boldsymbol{r}_1\int d^3 \boldsymbol{r}_2 \ u^*_B(\boldsymbol{r}_1)u^*_A(\boldsymbol{r}_2) \left(\frac{e^2}{r_{12}}\right) u_A(\boldsymbol{r}_1)u_B(\boldsymbol{r}_2) \end{align}
の積分を実行します。$V_C^{(12)}$同様、座標を図3のように設定します。そして
\begin{align} &\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)}, \ \varphi_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)}, \ \varphi_2:2,A,B\text{がなす平面の}x\text{軸からの角度} \end{align}
とします。$r_{12}$をこの座標で具体的に書き下すと
\begin{align} r_{12}=\sqrt{r_1^2+r_2^2-2r_1r_2\cos(\varphi_1-\varphi_2)+\frac{R^2}{4}(\xi_1\eta_1-\xi_2\eta_2)^2} \end{align}
となります。

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

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

この展開を用いると、$V_{EX}^{(12)}(R)$の被積分関数の中で$\varphi$依存性をもつのは展開に現れる$\cos$の部分だけであり、$\int_0^{2\pi}d\varphi_1\int_0^{2\pi}d\varphi_2$を行うと$\nu=0$の寄与のみが残ります。そして
\begin{align}　 P_\tau^0=P_\tau, \ Q_\tau^0=Q_\tau, \ D_{\tau 0}=\frac{2\tau+1}{2} \quad (P_\tau, Q_\tau\text{はそれぞれ第1種および第2種のルジャンドル多項式}) \end{align}
なので（Appendix参照のこと）、件の積分は
\begin{align} V_{EX}^{(12)}(R) &=\frac{e^2D^6}{16} \int_1^\infty e^{-D\xi_1}d\xi_1 \int_{-1}^1(\xi_1^2-\eta_1^2)d\eta_1 \int_1^\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}^\infty\frac{2\tau+1}{2} \big( P_\tau(\xi_1)Q_\tau(\xi_2)P_\tau(\eta_1)P_\tau(\eta_2)\theta(\xi_2-\xi_1) \\ &\hspace{4.5cm}+P_\tau(\xi_2)Q_\tau(\xi_1)P_\tau(\eta_1)P_\tau(\eta_2)\theta(\xi_1-\xi_2) \big) \end{align}
となります。以上から$V_{EX}^{(12)}(R)$は
\begin{align} V_{EX}^{(12)}(R) =& \frac{e^2D^6}{16} \times \frac{2}{a_0D} \sum_{\tau=0}^\infty \frac{2\tau+1}{2} \times \int_1^\infty e^{-D\xi_1}d\xi_1 \int_{-1}^1 (\xi_1^2-\eta_1^2)d\eta_1 \\ &\times\Bigg[ \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 \\ &\quad+\int_{\xi_1}^{\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 \Bigg] \end{align}
になります。ここで$\eta_2$の積分を実行すると
\begin{align} \int_{-1}^1(\xi_2^2-\eta_2^2)P_\tau(\eta_2)d\eta_2 &= \begin{cases} \displaystyle 2\left(\xi_2^2-\frac{1}{3}\right) & \text{for } \tau=0 \\ \displaystyle -\frac{4}{15} & \text{for }\tau=2 \\ 0 & \text{other } \tau's \end{cases} \\ \\ &= 2\left(\xi_2^2-\frac{1}{3}\right)\delta_{\tau 0} -\frac{4}{15}\delta_{\tau 2} \end{align}
となるので、のこるのは$\tau=0,2$のみです。$Q_0,Q_2,P_0,P_2$は
\begin{aligned} Q_0(x)&=\log\frac{x+1}{x-1},\quad Q_2(x)= -3x+P_2(x)\log\frac{x+1}{x-1}, \\ P_0(x)&=1, \quad P_2(x)=\frac{3x^2-1}{2} \end{aligned}
です（Appendix参照）。これらを$V_{EX}^{(12)}(R)$に代入して
\begin{aligned} V_{EX}^{(12)}(R) = \frac{e^2D^5}{8a_0} \times \Bigg\{ &\frac{1}{2} \int_1^\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\left(\xi_2^2-\frac{1}{3}\right) \times P_0(\xi_2)Q_0(\xi_1)P_0(\eta_1) \\ &+\frac{5}{2} \int_1^\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\left(-\frac{4}{15}\right)\times P_2(\xi_2)Q_2(\xi_1)P_2(\eta_1) \\ &+\frac{1}{2} \int_1^\infty e^{-D\xi_1}d\xi_1 \int_{-1}^1(\xi_1^2-\eta_1^2)d\eta_1 \int_{\xi_1}^\infty e^{-D\xi_2}d\xi_2 \times 2\left(\xi_2^2-\frac{1}{3}\right)\times P_0(\xi_1) Q_0(\xi_2)P_0(\eta_1) \\ &+\frac{5}{2} \int_1^\infty e^{-D\xi_1}d\xi_1 \int_{-1}^1(\xi_1^2-\eta_1^2)d\eta_1 \int_{\xi_1}^\infty e^{-D\xi_2}d\xi_2 \times\left(-\frac{4}{15}\right)\times P_2(\xi_1)Q_2(\xi_2)P_2(\eta_1) \Bigg\} \\ = \frac{e^2D^5}{8a_0} \times \Bigg[ &\frac{1}{2} \int_1^\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\left(\xi_2^2-\frac{1}{3}\right) \times \log\frac{\xi_1+1}{\xi_1-1} \\ &+\frac{5}{2} \int_1^\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\left(-\frac{4}{15}\right)\times \frac{3\xi_2^2-1}{2} \left\{ -3\xi_1+P_2(\xi_1)\log\frac{\xi_1+1}{\xi_1-1}\right\} \times\frac{3\eta_1^2-1}{2} \\ &+\frac{1}{2} \int_1^\infty e^{-D\xi_1}d\xi_1 \int_{-1}^1(\xi_1^2-\eta_1^2)d\eta_1 \int_{\xi_1}^\infty e^{-D\xi_2}d\xi_2 \times 2\left(\xi_2^2-\frac{1}{3}\right) \times \log\frac{\xi_2+1}{\xi_2-1} \\ &+\frac{5}{2} \int_1^\infty e^{-D\xi_1}d\xi_1 \int_{-1}^1(\xi_1^2-\eta_1^2)d\eta_1 \int_{\xi_1}^\infty e^{-D\xi_2}d\xi_2 \times\left(-\frac{4}{15}\right)\times \frac{3\xi_1^2-1}{2} \left\{ -3\xi_2+P_2(\xi_2)\log\frac{\xi_2+1}{\xi_2-1} \right\} \times\frac{3\eta_1^2-1}{2} \Bigg] \end{aligned}
を得ます。

さて$\eta_1$はすぐ積分できます。以下の公式
\begin{align} \int_{-1}^1(\xi_1^2-\eta_1^2)d\eta_1&=2\xi_1^2-\frac{2}{3}, \\ \int_{-1}^1(\xi_1^2-\eta_1^2) \frac{3\eta_1^2-1}{2}d\eta_1 &=-\frac{4}{15} \end{align}
を用いれば
\begin{align} V_{EX}^{(12)}(R) = \frac{e^2D^5}{8a_0} \Bigg[&\int_1^\infty e^{-D\xi_2} \left(\xi_1^2-\frac{1}{3}\right) \log\frac{\xi_1+1}{\xi_1-1}d\xi_1 \times \int_1^{\xi_1}e^{-D\xi_2} \times 2 \left( \xi_2^2-\frac{1}{3} \right) \\ &+\frac{5}{2}\times\left(-\frac{4}{15}\right)^2 \times \int_1^\infty e^{-D\xi_1} \left\{ -3\xi_1+P_2(\xi_1) \log\frac{\xi_1+1}{\xi_1-1} \right\} d\xi_1 \times \int_1^{\xi_1} e^{-D\xi_2} \frac{3\xi_2^2-1}{2} d\xi_2 \\ &+2\int_1^\infty e^{-D\xi_1} \left(\xi_1^2-\frac{1}{3}\right)d\xi_1 \times \int_{\xi_1}^\infty e^{-D\xi_2} \left( \xi_2^2-\frac{1}{3} \right) \log\frac{\xi_2+1}{\xi_2-1} d\xi_2 \\ &+\frac{5}{2}\times\left(-\frac{4}{15}\right)^2 \int_1^\infty e^{-D\xi_2} \frac{3\xi_1^2-1}{2}d\xi_1 \times \int_{\xi_1}^\infty e^{-D\xi_2} \left\{ -3\xi_2+P_2(\xi_2) \log\frac{\xi_2+1}{\xi_2-1} \right\} d\xi_2 \Bigg] \\ \\ = \frac{e^2D^5}{8a_0} \Bigg[& \frac{5}{2} \times \left(\frac{4}{15}\right)^2 \left\{ \int_1^\infty (-3)\xi_1e^{-D\xi_1} d\xi_1 \int_1^{\xi_1} e^{-D\xi_2} \frac{3\xi_2^2-1}{2} d\xi_2 +\int_1^\infty e^{-D\xi_2} \frac{3\xi_1^2-1}{2}d\xi_1 \int_{\xi_1}^\infty (-3)\xi_2e^{-D\xi_2}d\xi_2 \right\} \\ &+\int_1^\infty e^{-D\xi_1} \left( \xi_1^2-\frac{1}{3} \right)d\xi_1 \times \Bigg\{ \log\frac{\xi_1+1}{\xi_1-1} \times 2 \int_1^{\xi_1} e^{-D\xi_2} \left(\xi_2^2-\frac{1}{3}\right) d\xi_2 \\ &\hspace{4.6cm} +2\int_{\xi_1}^\infty e^{-D\xi_2} \left(\xi_2^2-\frac{1}{3}\right) \log\frac{\xi_2+1}{\xi_2-1}d\xi_2 \\ &\hspace{4.6cm} +\frac{3}{2}\times\frac{5}{2} \times\left(-\frac{4}{15}\right)^2 \log\frac{\xi_1+1}{\xi_1-1} \int_1^{\xi_1} e^{-D\xi_2}\frac{3\xi_2^2-1}{2}d\xi_2 \\ &\hspace{4.6cm} +\frac{3}{2}\times\frac{5}{2} \times\left(-\frac{4}{15}\right)^2 \int_{\xi_1}^\infty e^{-D\xi_2} \frac{3\xi_2^2-1}{2} \log \frac{\xi_2+1}{\xi_2-1} d\xi_2 \Bigg\} \Bigg] \\ =\frac{e^2D^5}{10a_0} \Bigg[& -\Bigg\{ \int_1^\infty \xi_1 e^{-D\xi_1} d\xi_1 \int_1^{\xi_1} e^{-D\xi_2} \left( \xi_2^2-\frac{1}{3} \right) d\xi_2 + \int_1^\infty e^{-D\xi_1} \left(\xi_1^2-\frac{1}{3}\right) \int_{\xi_1}^\infty \xi_2e^{-D\xi_2}d\xi_2 \Bigg\} \\ &+3\int_1^\infty e^{-D\xi_1} \left(\xi_1^2-\frac{1}{3}\right) d\xi_1 \times \Bigg\{ \log\frac{\xi_1+1}{\xi_1-1} \int_1^{\xi_1}e^{-D\xi_2} \left(\xi_2^2-\frac{1}{3}\right)d\xi_2 +\int_{\xi_1}^\infty e^{-D\xi_2} \left(\xi_2^2-\frac{1}{3}\right) \log\frac{\xi_2+1}{\xi_2-1} d\xi_2 \Bigg\} \Bigg] \end{align}
となります。

上式の最終行の各積分を実行していきましょう。

最初の中括弧の積分

これは簡単に実行できます：
\begin{align} \int_1^\infty \xi_1e^{-D\xi_1}d\xi_1 \int_1^{\xi_1} e^{-D\xi_2} \left( \xi_2^2-\frac{1}{3} \right) d\xi_2 &= \int_1^\infty e^{-D\xi_1} \left(\xi_1^2-\frac{1}{3}\right) \int_{\xi_1}^\infty \xi_2e^{-D\xi_2}d\xi_2 \\ &=\frac{e^{-2D}}{24D^5} (8D^3+24D^2+30D+15) \end{align}

後半の中括弧の積分

\begin{align} \int_1^\infty e^{-D\xi_1} \left(\xi_1^2-\frac{1}{3}\right) d\xi_1 \times \Bigg\{ \log\frac{\xi_1+1}{\xi_1-1} \int_1^{\xi_1}e^{-D\xi_2} \left(\xi_2^2-\frac{1}{3}\right)d\xi_2 +\int_{\xi_1}^\infty e^{-D\xi_2} \left(\xi_2^2-\frac{1}{3}\right) \log\frac{\xi_2+1}{\xi_2-1} d\xi_2 \Bigg\} \tag{8} \end{align}
を計算します。
\begin{align} \int_a^\infty d\xi_2e^{-D\xi_2} \left(\xi_2^2-\frac{1}{3}\right)\log(\xi_2+1) \end{align}
を計算するには部分積分を用います。$e^{-D\xi_2}(\xi_2^2-1/3)$の不定積分は
\begin{aligned} \int d\xi_2 e^{-D\xi_2}\left( \xi_2^2-\frac{1}{3}\right) &=\int e^{-D\xi_2}(\xi_2^2-\frac{1}{3})d\xi_2 \\ &= -e^{-D\xi_2} \left( \frac{\xi_2^2}{D} +\frac{2}{D^2}\xi_2 +\frac{2}{D^3} \right) +\frac{1}{3}e^{-D\xi_2} \frac{1}{D} \end{aligned}
であり、また
\begin{align} \int_a^\infty e^{-D\xi_2}\frac{1}{\xi_2+1}d\xi_2&=-e^D Ei(-D(a+1)), \\ Ei(x)&:=-\text{p.v.}\int_{-x}^\infty \frac{e^{-t}}{t}dt \end{align}
です（$Ei$は指数積分と呼ばれます。$\text{p.v.}$はCauchyの主値積分）。これらを用いて計算すれば
\begin{aligned} \int_a^\infty e^{-D\xi_2} \left(\xi_2^2-\frac{1}{3}\right) \log(\xi_2+1)d\xi_2 &= e^{-aD} \left\{ \left( \frac{a^2}{D} +\frac{2a}{D^2} +\frac{2}{D^3} \right) -\frac{1}{3D} \right\} \log(1+a) \\ &+\frac{e^{-aD}}{D^3} \left\{ D(a-1)+3 \right\} \\ &+\left( -\frac{2}{3}\frac{1}{D}+\frac{2}{D^2}-\frac{2}{D^3} \right) e^D Ei(-D(a+1)) \end{aligned}
となります。同様に$a>1$のとき
\begin{aligned} \int_a^\infty e^{-D\xi_2} \left(\xi_2^2-\frac{1}{3}\right) \log(\xi_2-1)d\xi_2 &= e^{-aD} \left\{ \left( \frac{a^2}{D} +\frac{2a}{D^2} +\frac{2}{D^3} \right) -\frac{1}{3D} \right\} \log(a-1) \\ &+\frac{e^{-aD}}{D^3} \left\{ D(a+1)+3 \right\} \\ &-\left( \frac{2}{3}\frac{1}{D}+\frac{2}{D^2}+\frac{2}{D^3} \right) e^D Ei(-D(a-1)) \end{aligned}
となります。以上より
\begin{aligned} \text{Eq.}(8)= & \int_1^\infty e^{-D\xi_1} \left(\xi_1^2-\frac{1}{3}\right) \log(\xi_1+1) \Big\{ -e^{-D\xi_1} A(\xi_1,D) +\frac{1}{3}e^{-D\xi_1}\frac{1}{D} +B(D) \Big\} \\ - & \int_1^\infty e^{-D\xi_1} \left(\xi_1^2-\frac{1}{3}\right) \log(\xi_1-1) \Big\{ -e^{-D\xi_1} A(\xi_1,D) +\frac{1}{3}e^{-D\xi_1}\frac{1}{D} +B(D) \Big\} \\ &+\int_1^\infty e^{-D\xi_1} \left(\xi_1^2-\frac{1}{3}\right) \Bigg[ e^{-D\xi_1}\left(A(\xi_1,D)-\frac{1}{3D}\right)\log(1+\xi_1) \\ & \hspace{3.5cm} +\frac{e^{-D\xi_1}}{D^3} (D\xi_1-D+3) \\ & \hspace{3.5cm} +C(D)Ei(-D(\xi_1+1)) \\ & \hspace{3.5cm} -\Big\{ \left( e^{-D\xi_1}A(\xi_1,D) -\frac{1}{3}e^{-D\xi_1}\frac{1}{D} \right)\log(\xi_1-1) \\ &\hspace{3.5cm} +\frac{e^{-D\xi_1}}{D^3}D\xi_1 +\frac{e^{-D\xi_1}}{D^2} +\frac{3}{D^3}e^{-D\xi_1} \\ &\hspace{3.5cm} -B(D)Ei(-D(\xi_1-1)) \Big\} \Bigg], \\ \\ &A(\xi_1,D):= \frac{1}{D}\xi_1^2+\frac{2}{D^2}\xi_1+\frac{2}{D^3},\quad B(D):=e^{-D} \left( \frac{2}{3}\frac{1}{D} +\frac{2}{D^2} +\frac{2}{D^3} \right), \quad C(D):=e^D \left( -\frac{2}{3}\frac{1}{D} +\frac{2}{D^2} -\frac{2}{D^3} \right) \end{aligned}
を得ます。いくつかの項が打ち消し合いますが、整理すると
\begin{aligned} =&B(D) \int_1^\infty e^{-D\xi_1} \left(\xi_1-\frac{1}{3}\right) \log(\xi_1+1)d\xi_1 \\ -&B(D) \int_1^\infty e^{-D\xi_1} \left(\xi_1-\frac{1}{3}\right) \log(\xi_1-1)d\xi_1 \\ +&B(D) \int_1^\infty e^{-D\xi_1} \left(\xi_1-\frac{1}{3}\right) Ei(-D(\xi_1-1))d\xi_1 \\ +&C(D) \int_1^\infty e^{-D\xi_1} \left(\xi_1-\frac{1}{3}\right) Ei(-D(\xi_1+1))d\xi_1 \\ -&\frac{2}{D^2} \int_1^\infty e^{-2D\xi_1} \left(\xi_1^2-\frac{1} {3}\right)d\xi_1 \end{aligned}
となります。以上から$V_{EX}^{(12)}(R)$は
\begin{aligned} V_{EX}^{(12)}(R)= \frac{e^2D^5}{10a_0} \Bigg[ &-\frac{e^{-2D}}{12D^5} (8D^3+24D^2+30D+15) -\frac{6}{D^2} \int_1^\infty e^{-2D\xi_1} \left(\xi_1^2-\frac{1}{3}\right)d\xi_1 \\ &+3 \Bigg\{ B(D)\int_1^\infty e^{-D\xi_1} \left(\xi_1^2-\frac{1}{3}\right)\log(\xi_1+1)d\xi_1 \\ &\quad\quad+B(D)\int_1^\infty e^{-D\xi_1}\left(\xi_1^2-\frac{1}{3}\right) (Ei(-D(\xi_1-1)-\log(\xi_1-1))d\xi_1 \\ &\quad\quad+C(D)\int_1^\infty e^{-D\xi_1}\left(\xi_1^2-\frac{1}{3}\right) Ei(-D(\xi_1+1))d\xi_1 \Bigg\} \Bigg] \\ \\ = \frac{e^2D^5}{10a_0} \Bigg[ & -e^{-2D} \left( \frac{2}{3}\frac{1}{D^2} +4\frac{1}{D^3} +\frac{11}{2}\frac{1}{D^4} +\frac{11}{4}\frac{1}{D^5} \right) \\ &+3 \Bigg\{ B(D)\int_1^\infty e^{-D\xi_1} \left(\xi_1^2-\frac{1}{3}\right)\log(\xi_1+1)d\xi_1 \\ &\quad\quad+B(D)\int_1^\infty e^{-D\xi_1}\left(\xi_1^2-\frac{1}{3}\right) (Ei(-D(\xi_1-1)-\log(\xi_1-1))d\xi_1 \\ &\quad\quad+C(D)\int_1^\infty e^{-D\xi_1}\left(\xi_1^2-\frac{1}{3}\right) Ei(-D(\xi_1+1))d\xi_1 \Bigg\} \Bigg] \end{aligned}
となります。あとは
\begin{align} &\int_1^\infty e^{-D\xi_1} \left(\xi_1^2-\frac{1}{3}\right) \log(\xi_1+1)d\xi_1 \tag{9},\\ &\int_1^\infty e^{-D\xi_1}\left(\xi_1^2-\frac{1}{3}\right)(Ei(-D(\xi_1-1))-\log(\xi_1-1))d\xi_1 \tag{10},\\ &\int_1^\infty e^{-D\xi_1}\left(\xi_1^2-\frac{1}{3}\right) Ei(-D(\xi_1+1))d\xi_1 \tag{11} \end{align}
が計算できればよいです。

Eq.(9)の不定積分はすでに計算しました。これを用いてEq.(9)は
\begin{aligned} \int_1^\infty e^{-D\xi_1}\left(\xi_1^2-\frac{1}{3}\right)\log(\xi_1+1) d\xi_1 &= \frac{2}{3}\int_1^\infty d\xi_1 e^{-D\xi_1}P_2(\xi_1)\log(\xi_1+1) \\ &= e^{-D} \left( \frac{2}{3}\frac{1}{D} +\frac{2}{D^2} +\frac{2}{D^3} \right)\log2 \\ &+\frac{3}{D^3}e^{-D} \\ &+\left( \frac{2}{D^2} -\frac{2}{D^3} -\frac{2}{3}\frac{1}{D} \right) e^D E_i(-2D) \\ &=B(D)\log2+\frac{3}{D^3}e^{-D}+C(D)Ei(-2D) \end{aligned}
となります。

Eq.(11)は
\begin{align} \frac{d}{d\xi_1}Ei(-D(\xi_1+1))=\frac{e^{-D(\xi_1+1)}}{\xi_1+1} \end{align}
を用いれば、やはり部分積分で計算できます。結果は
\begin{aligned} &\int_1^\infty e^{-D\xi_1} \left(\xi_1^2-\frac{1}{3}\right) Ei(-D(\xi_1+1)) \\ &= \frac{1}{D^3} \left\{ 2e^{-D} \left(1+D+\frac{1}{3}D^2\right) Ei(-2D) +\frac{5}{4}e^{-3D} -2\left(1-D+\frac{1}{3}D^2\right) e^D Ei(-4D) \right\} \\ &=B(D)Ei(-2D)+\frac{5}{4D^3}e^{-3D}+C(D)Ei(-4D) \end{aligned}
となります。

最後のEq.(10)ですが、丸括弧の第1項と第2項それぞれの積分は発散します。しかしその差であるEq.(10)は有限です。実際指数積分は以下の展開を持ちますSugiura：

ここで$\gamma$はオイラー定数です。これを用いれば積分の下限で
\begin{align} \lim_{\xi_1\to 1+0}(Ei(-D(\xi_1-1))-\log(\xi_1-1)) =\gamma+\log D \end{align}
であり有限です。改めて、Eq.(10)は
\begin{align} \frac{d}{d\xi_1} (Ei(-D(\xi_1-1))-\log(\xi_1-1)) =\frac{e^{-D(\xi_1-1)}}{\xi_1-1}-\frac{ 1}{\xi_1-1} \end{align}
を用いると、これまでと同様に部分積分で計算できて
\begin{align} &\int_1^\infty e^{-D\xi_1}\left(\xi_1^2-\frac{1}{3}\right)(Ei(-D(\xi_1-1))-\log(\xi_1-1))d\xi_1 \\ &=\frac{2}{D^3} e^{-D} (D^2/3+D+1)(\gamma+\log D-\log2) -\frac{e^{-D}}{D^3} (D+7/4) \\ &=B(D)(\gamma+\log D-\log2)-\frac{e^{-D}}{D^3}\left(D+\frac{7}{4}\right) \end{align}
となります。

最終的な$V_{EX}^{(12)}$の表式

以上から$V_{EX}^{(12)}(R)$の積分は実行できて、以下のようになります：

これで$E_s, E_t$に現れるすべての積分が実行できました。

$E_s, E_t$のグラフ

以上の計算を用いて$E_s,E_t$のグラフを描いたのが図4です。

$E_s, E_t$のグラフ（縦軸の単位はeV）。$E_s$は青線、$E_t$はオレンジの線。横軸は陽子間距離$R[\text{Å}]$。エネルギーは$2E_{1s}$から測っている。

横軸は陽子間距離$R$（単位はオングストローム）、縦軸は$E_s$（青線、単位はeV）および$E_t$（オレンジの線、単位はeV）です。ただしエネルギーは$2E_{1s}$から測っていることに注意してください。$E_t$は常に$E_s$よりも大きい & 最小値を持たない一方、$E_s$は$R\simeq 0.8\text{Å}$付近で最小となります。すなわち水素分子はスピン1重項かつ水素原子が$0.8\text{Å}$程度の距離のときに安定となります。このときエネルギーは$-3\text{eV}$程度です。これらの値は実験値:$0.741\text{Å}, -4.747\text{eV}$とだいたい合っていますIgi。Heitler-Londonの方法は比較的よい近似であることがわかります。

スピン1重項のほうがエネルギーが低くなる理由は、空間の波動関数が対称であり、陽子間における電子の存在確率が高くなるからです。このとき陽子により電子間のクーロン斥力が遮蔽され、そのぶんエネルギーが低くなります。このようにして電子により原子同士が結合することを「共有結合」と呼びますIgi。

まとめ

本記事ではBorn-Oppenheimer近似およびHeitler-Londonの方法における水素分子のエネルギーを解析的に計算しました。特に$V_{EX}^{(12)}$の計算を書いた文献は少なく思えたので、これを詳細に記しました。

最後にHeitler-Londonの方法と変分法について述べておきます。本記事では電子2体系の波動関数をEq.(2)(3)の形だと仮定しましたが、もとのHeitler-Londonの方法では変分法で波動関数を求めていますOhiwa。電子2体系の波動関数が$u_A(\boldsymbol{r}_1)u_B(\boldsymbol{r}_2)$および$u_A(\boldsymbol{r}_2)u_B(\boldsymbol{r}_1)$の線形結合
\begin{align} C_1u_A(\boldsymbol{r}_1)u_B(\boldsymbol{r}_2)+C_2u_A(\boldsymbol{r}_2)u_B(\boldsymbol{r}_1) \end{align}
で書けるとして、エネルギーを最小にする条件でエネルギーおよび$C_1,C_2$を決定します。するとエネルギーは永年方程式の解として求まり、$E_s$または$E_t$となります。係数$C_1,C_2$は$E_s$に対して$C_1=C_2$、$E_t$に対して$C_1=-C_2$となって、結局本記事の結果と一致します。

おしまい。${}_\blacksquare$
${}$

Appendix: ルジャンドル多項式

第1種ルジャンドル多項式$P_n(x)$、第2種ルジャンドル多項式$Q_n(x)$は以下で定義されます：

これらは以下の微分方程式の特殊解です：
\begin{align} (1-x^2)y''-2xy'+n(n+1)y=0 \end{align}
$Q_n$は以下のように書けます：

ルジャンドルの陪多項式は以下で定義されます：

以上。${}_\blacksquare$