少数電子系の高精度計算(1)

最近、"精密計算化学"みたいな分野が存在してほしい気持ちが、私の中で高まってきた(精密化学という名前は、fine chemistryの訳語に既に使われてる)。電子数が1〜4の原子や分子では、Hylleraasが始めて、James & Coolidgeらに拡張された変分計算が(Born-Oppenheimer近似を仮定した限りでは)最高精度なので、最初に、この方法を見ていく。

歴史

17世紀末〜18世紀前半に、Keplerの法則と(二体問題における)万有引力の法則の"等価性"が示された。また、概算では、地球と月に及ぼしている力と、地球が地表付近の物体に及ぼす力を比較した結果、同じ種類の力が働いてると推測された。万有引力の法則は、適用範囲が広くて新しい予測が沢山出るので、どれくらい強固な理論か試験されることになった。

簡単に出る予測として、地表付近の重力が高度変化する時の一次係数の値がある。変化率は、3メートルにつき100万分の一程度という割合(0.308mGal/m)なので、重力の精密測定が必要になる。実際は、地球の密度が一様でないので、高度補正が、どれくあい綺麗に出るかは分からないし、17〜18世紀に、測定した人がいたかも不明。一応、17世紀の技術でも、振り子の周期測定などから、重力の高度依存性を検出できただろう。1670年代に、Jean Richerは、振り子の周期が場所によって異なるのを発見してるから発想はありえたと思うけど、ボルダの振り子やケーターの振り子が開発されたのは18世紀末〜19世紀初頭。まあ当時はそんな測定しても、実用的な意味もなかったろう。

別の予測として、Keplerの法則で説明しきれない、月の軌道の精密計算が目標となった。これは実用的価値があったけど、三体問題なので計算は難しく、ニュートンオイラーも解決できなかった。計算精度が低かったために、一時期は、万有引力の法則が厳密には正しくない可能性や、惑星や月の間に、複数種類の力が働いてる可能性も議論された。結局、1740年代後半に計算法の問題だと判明して、1750年頃、当時のエライ人たちは、(当時の測定水準では)月の軌道からは万有引力の法則の逸脱を検出できないと合意した。

この時点でも、月の軌道計算の精度は、実用的に十分な水準には達してなかったようで、1750年以後も、オイラーは計算法の改良に取り組んだ。オイラーも、これといった成果も残せず、月の軌道の精密計算は、19世紀にも課題として持ち越された。計算法の改良は、ハレー彗星の回帰時期予測にも適用されて、ハレーの予測を改善したらしい。ハレー彗星の回帰は、1758年末〜1759年に観察されたけど、予測と観測には、まだ大きな差があった。



量子力学の場合、Keplerの法則と似た立ち位置に、Rydbergの公式(やBohr模型)があった。水素原子に対する量子力学の予測は、それ単体で見れば、簡単なものを難しいものに置き換えただけなので、新しい予測として、ヘリウム原子や水素分子で、基底状態エネルギーの計算を試みた人がいた。

1927年1月出版の論文で、オランダの物理学者Øyvind Burrauは、水素分子イオンH2+に対して、シュレディンガー方程式を解いて、基底状態エネルギーを計算した。水素分子イオンは、一電子系なので、電子相関はない。
Berechnung des Energiewertes des Wasserstoffmolekel-Ions (H{2/+}) im Normalzustand
http://gymarkiv.sdu.dk/MFM/kdvs/mfm%201-9/mfm-7-14.pdf

論文には、1922年に、PauliとK.F. Niessenが(独立に)、古典力学に"量子化条件"を考慮することで、エネルギーの計算を行ったことが書いてある(計算結果がどうだったのかは書いてない)。Born-Oppenheimer近似の下では、エネルギーは核間距離に依存し、更に、このポテンシャルエネルギーを使って原子核に対するシュレディンガー方程式を解くことになる。普通、ポテンシャルエネルギー曲線は解析的な形では書けないので、二原子分子ではMorseポテンシャルのようなモデルを使ったり、もっと単純にエネルギー極小値点近傍で調和振動子近似を適用することが、しばしば行われる(Morseの論文は、1929年に出てる)。

Burrauの論文によると、エネルギーの極小値は、-16.29(eV)で、これに極小点近傍のゼロ点振動分0.07eVを加味して、水素分子イオンのイオン化エネルギーを、16.22(eV)と算出した。1926年に、James FranckとPascual Jordanが測定した結果によれば、水素分子イオンのイオン化エネルギーは、16.1(eV)だったそうだ。0.1(eV)程度の差があるけど、測定誤差の範囲である。

正確な実測値が、どのくらいだったか覚えてないので、数値を確認しとくと、2017年の論文

Fundamental Transitions and Ionization Energies of the Hydrogen Molecular Ions with Few ppt Uncertainty
https://doi.org/10.1103/PhysRevLett.118.233001

には、QEDで、H2+のイオン化エネルギーを、12桁精度まで計算した値(自分で確認はしてないので合ってるかは知らない)が書いてあって、131058.121993(cm-1)=16.24913(eV)程度らしい。1920〜1930年代頃の論文を見てると、1Ry=13.53(eV)という数値が使われてて、今と単位の基準が違ったか、単位電荷の測定精度が低かったかしてたっぽい。単純比較はできないものの、Burrauの結果は、かなりよく測定値を再現してたと思われる。

Burrauの論文には、1926年に水素分子の解離エネルギーが測定されて、4.4(eV)だったことも書いてある。これも、一応、数値を確認しとくと、
(H2の解離エネルギー)=(H2のイオン化エネルギー)+(H2+のイオン化エネルギー)-2*(Hの基底状態エネルギー)
で、H2+のイオン化エネルギーは上に書いた通りの値を採用して、H2のイオン化エネルギーは、

Hydrogen@NIST Webbook
https://webbook.nist.gov/cgi/cbook.cgi?ID=C1333740&Mask=20#Ion-Energetics

の数値だと、15.4259(eV)程度。従って、H2の解離エネルギーは、16.24913+15.4259-13.598433*2=4.478164(eV)≒432(kJ/mol)程度と計算できて、1926年の測定値4.4(eV)は、そこそこ正しそう。

1927年4月付けの論文で、Edward Condonという人は、水素分子イオンの計算と、ヘリウムの第一イオン化ポテンシャルの測定値を使って、水素分子の解離エネルギーを、4.4(eV)と見積もれると書いている。Condonの議論は言葉で書かれてる部分が多くて、妥当なのかは分からない。この議論を信用するなら、ヘリウムの基底状態エネルギーが正確に計算できれば、水素分子の解離エネルギーも、当時の測定誤差と同程度の精度で計算できることになる。Edward Condonは、Franck-Condon原理で知られる人と同一人物だと思う。



1927年に、Albrecht Unsöld、John Slater、Georg W. Kellnerという人らが、独立に、シュレディンガー方程式に基づいて、ヘリウムの基底状態のエネルギーを計算した。ニ電子系では、電子相関があって、電子相関を考慮すれば、高い精度で計算できるかは当時まだ未知だった(まだ電子相関という言葉はなかったけども)。

ヘリウムの第2イオン化エネルギーは、水素原子と同様に計算できるので、基底状態のエネルギーが分かれば、第1イオン化エネルギーも分かる。3人の計算法は全員異なっていて、一番精度の高い計算をしたのはSlater。論文の日付は1927年4月26日になってる。
The Structure of the Helium Atom: I
https://doi.org/10.1073/pnas.13.6.423
論文には、"摂動論"だと書いてるけど、現代の人が摂動論と聞いて想像する方法とは違うように見える(実は等価だったりするのかもしれないけど、よく考えてない)。方法はともかく、結果を見ると、第1イオン化エネルギーの計算値と測定値の差は0.13eVとなっている。論文には、計算値と測定値の一致は、very remarkableだと書いてある。当時は、Theodore Lymanの1924年の論文に記載されてた198298(cm-1)=24.586(eV)が実験値として採用されてた(Slaterの論文には、おそらく、1Ry=13.53eVとされてたことと関連して、24.586ではなく、24.48という数値が書いてある)

Lymanは、水素原子の離散スペクトルの内、遠紫外領域にある(Schumann plateというものを使って測定してたらしく、当時の論文には、Schumann regionとか書いてあったりする)ライマン系列の測定に成功した後、1915年には、極紫外領域にあるヘリウム原子の離散スペクトルのために、測定方法の改良を始めている。



1927年、HeitlerとLondonは、6月30日付けの論文で、水素分子の基底状態に対する簡単な計算を行った。Heitler-London理論は、0次の摂動波動関数(つまり無摂動波動関数)を使って、1次までの摂動エネルギーを計算した近似である。定量的には精度が低く、水素分子の解離エネルギーの計算値は測定値より1.5eVくらい小さい。論文では、ポテンシャルエネルギー曲線のおおまかな形状を再現することが強調されている。

ヘリウムで、Heitler-London理論と同じように、基底状態のエネルギーを、一次の摂動項まで考慮して計算すれば、-2.75(au)になって、測定値−2.903(au)との差は、4(eV)強。多分、この計算は、Albrecht Unsöldがやったのが最初だと思う。1929年に、Slaterは、一般の多電子原子に対する同様の計算を行い、1次の摂動項まで考慮した基底状態エネルギーと、Hundの規則が整合的であることを示した。

摂動論が期待通りに機能するなら、高次の摂動まで考慮することで、水素分子でも、多電子原子でも、エネルギーと波動関数を精密に計算できるはず。こういう方針で証明や計算が行われたことがあるのかは分からないけど、多分、収束は遅い。

一般論として、二次の摂動エネルギーを計算するには、一次の摂動波動関数が必要になり、ヘリウムの場合は、厳密に解析的な関数で書ける(この事実は、多分1950年代末に発見された)。水素分子の場合には、多分、一次の摂動波動関数すら解析的には計算できないので、無摂動ハミルトニアンの固有関数で展開することになる(通常、無摂動ハミルトニアンは、連続スペクトルを持ち、それを無視するのは正しくない気がするけど、扱い方が分からないので、固有値と固有関数のみが計算に使われることになるのだと思う)。

一次の摂動波動関数が厳密に分かるヘリウムの場合、二次までの摂動エネルギーを計算に入れても、まだ測定値との差が1(eV)以上残る。Slaterが、1927年の論文で、ヘリウムの計算を、どうやったのか詳細を見てないけど、計算値と測定値の差が0.13eVだったのは、割と健闘してたと言える。

Hetiler-London流の摂動計算は、精度が低かったけど、比較的すぐに多原子分子に拡張された。1910年代には共有結合の定性的モデルは存在していたけど、Heitler-Londonの理論の延長上に、共有結合量子力学的解釈が行われた。実際のところ、無摂動ハミルトニアンを、どう取るかは自由なので、原子軌道を特別な基底と思う先験的な理由は別にないはずだけど、現在でも化学結合の定性的説明が、このような方法で与えられる。普通、数式とかは使わずに説明されてるので、それが摂動論的描像なのは分かりにくい。



イギリスのDouglas Hartreeは、1928年に出版された論文(受領は1927年)
The Wave Mechanics of an Atom with a Non-Coulomb Central Field. Part II. Some Results and Discussion
https://doi.org/10.1017/S0305004100011920
の中で、Self-Consistent Field法とHartreeが命名した方法によって、原子やイオンの基底状態エネルギーの計算を実行した。特に、ヘリウムの第一イオン化エネルギーを、1.835(Ry)=24.97(eV)と計算してる。Hartreeの論文を受けて出版されたJohn Gauntの論文では、1.75(Ry)=23.81(eV)と計算するべきだろうと書いてあって、測定値との差は、0.8(eV)弱。Slaterの計算には及んでないが、上記の摂動計算よりはいい。

Hartreeの方法は、色々と修正や改良が加えられたものの、根幹部分は、Hartree-Fock近似として、現在も、計算化学で一番基本的な方法として生き残っている。けど、以下の1938年の論文を読むと、Hartreeの方法は、しばらくの間、多電子原子にのみ使われてて、多原子分子で使える方法にはなってなかったらしい。
Self-consistent field for molecular hydrogen
https://doi.org/10.1017/S0305004100020089
現在では、Hartree-Fock法は、分子軌道法の一部のように扱われるけど、1930年代の時点では、2つの方法は区別されていた。分子が扱えない方法に、分子軌道も何もない。分子軌道法は、1929年のLennard-Jonesの論文に始まると大体書いてある。定量性は低いけど、今でも教科書に書かれてるヒュッケル法などが、当時の分子軌道法の代表例。現在、計算化学で、分子軌道の係数の初期値を決めるのにも、拡張ヒュッケル法というものを使う。


1928年に出版された論文で、Shou-Chin Wang(王守竞)という人は、水素分子に対して、単純な非線形試行関数を用いた変分計算を示した。Wangの使った試行関数は、ただ一つのパラメータしか含んでなかった。が、Heitler-Londonの理論よりも、いい結果を与えた。

The Problem of the Normal Hydrogen Molecule in the New Quantum Mechanics
https://doi.org/10.1103/PhysRev.31.579

ヘリウムについて、同様の変分計算は、1927年のKellnerの論文に書かれている。
Die Ionisierungsspannung des Heliums nach der Schrödingerschen Theorie
https://doi.org/10.1007/BF01391720
この変分計算は、今でも教科書で見かけることがあるが、測定値と計算値との差は、1.5(eV)程度になる。Kellnerは、パラメータの多い試行関数でも計算を行ったが、複数パラメータを扱う時には、線形パラメータについてのみ変分計算を行い、非線形パラメータは決め打ちした値を使っているっぽい。最終的に得られた計算値は、Slaterが同年に与えたものより測定値との差が大きい。

1928年に、Egil Hylleraasは、線形試行関数を用いて、ヘリウムの基底状態の変分計算を行った。ここで使われた試行関数は、Kellnerの試行関数で、非線形パラメータを1に固定したものと見れる。論文に書いてある結果は、Kellnerのものより大分良くて、Slaterの計算と同程度の精度を達成している。KellnerやHylleraasの論文には、"Ritzの方法"と書いてあるけど、Ritzの1909年頃の論文でも、線形試行関数が使われてて、それに沿ったものと見ることが出来る。



1928年頃、ニ電子系でさえ、測定値と測定誤差の範囲で一致する計算値を出せてない。この惨状を、当時の物理学者や化学者が、どう思ってたのか分からないけど、Hylleraasは、ヘリウムの基底状態の計算精度に満足してなかったらしい。

Hylleraasは1929年の論文で、扱いやすい形の非線形試行関数を使って変分計算を行い、ヘリウムの基底状態のエネルギーを誤差0.01(eV)未満の精度で求めた。当時の測定精度を考えれば、これは十分と言える。

観測値とよく合う計算値が出たけど、Hylleraasの計算は、変分関数の極小値を出してるだけなので、真の極小値が、これよりずっと小さい可能性を否定しておかないと、理論の検証として十分とは言えない。基底状態エネルギーの下限を求める簡単な方法は、例えば、D.H.Weinsteinという人の1932年の以下の論文に書かれている。
A Lower Limit for the Ground State of the Helium Atom
https://doi.org/10.1103/PhysRev.40.797
論文の末尾に、この問題を示唆したのは、ポーリングだと書いてある。但し、ここでは、Hylleraasの論文で最良の結果を与えた試行関数ではなく、単純な波動関数を使ってるので、ヘリウム原子の基底状態エネルギーが、-3.1(hartree)よりは、大きいとしか言ってない。同じ方法で、Hylleraas試行関数を使って、下限の評価を行った論文は、1956年に出たものが最初と思われる。

1930年代には、Hylleraasの変分計算が、正しい値に収束することを疑う人もいたらしい。1937年には、以下のような論文が書かれている。
On the Convergence of the Hylleraas Variational Method
https://doi.org/10.1103/PhysRev.51.855

大分後の1977年に、Hylleraasの方法で、ヘリウムの基底状態エネルギーが厳密解に収束する証明をしたと主張する論文が出ている(以下のIIのTheorem10)。

The convergence of the Rayleigh-Ritz Method in quantum chemistry I
https://doi.org/10.1007/BF00548026

The convergence of the Rayleigh-Ritz Method in quantum chemistry II
https://doi.org/10.1007/BF00548027

尤も、収束するというだけなら、1927年のKellnerの基底や、1928年のHylleraasの基底でも収束はすると思う(IIのTheorem9は、Hylleraasの1928年の基底でも収束することを述べている)ので、現実的には、収束の早さも大切。この論文は、他にも色々書いてあって、GTO(Gauss-Type Orbital)基底が、正規直交基底をなすことの証明を、私は始めて見た。


1933年になって、Hubert JamesとAlbert Coolidgeという人が、Hylleraasの方法を拡張して、水素分子に適用した。James & Coolidgeは、計算値と測定値の差を、0.03(eV)と報告している。当時、水素分子の解離エネルギーの測定は、複数のグループで行われたが、総合的に見て、4.4〜4.5(eV)の範囲という程度の信頼性しかなかったようなので、測定誤差と区別が付かない程度の精度で計算できてたことになる。

水素分子に関して、James & Coolidgeの試行関数と、Shou-Chin Wangの試行関数を見比べると、1パラメータしかない時でも形が違っていて、Wangの試行関数の方が、いい結果を与える。

水素分子の励起状態の計算も、複数の方法で行われてたが、James & Coolidgeの方法の延長上に計算を行ったのは、1935年出版の
On the Two‐Quantum Σ‐States of the Hydrogen Molecule
https://doi.org/10.1063/1.1749607
だと思われる。James & Coolidgeも、水素分子の後、1930年代の内に、リチウム原子の基底状態と水素分子の励起状態を、同様の変分法で計算した。


結局、1930年代の電子状態計算法には、大雑把に言って、いくつかの半定量的方法、Hartree-Fockの方法(現在のHartree-Fock法と同一ではない)、Hylleraas型の変分計算の3種類があった。半定量的方法は、それ以前の化学の"理論"に比べれば、見かけ上の定量性は有してるけど、測定値を再現したり予測することは期待できない。1930年代には、一般の多原子分子を扱える方法は、半定量的なものしかなかった。現在では使うこともないと思うけど、1920~30年代に行われた半定量的計算は、今でも、量子化学の教科書には書いてある。現在流布している化学結合の定性的理解は、こうした方法に基づいて作られたものだし、新しい経験則を発見するのにも貢献してきた。

Hartree-Fockの方法は、現在のHartree-Fock法の礎になったけど、1930年代には、一般の多原子分子に使うことはできなかった。現在のHartree-Fock法に相当するアルゴリズムを、いつ誰が完成させたのかは、知らない。GaussianとGamessは1970年台にリリースされ、1970年代末には、post-Hartree-Fock法という言葉が使われるようになってる。タンパク質やDNAのような高分子の計算は現在でも難しいけど、化学者が必要とする多くの化合物の化学物性を、そこそこ定量的に計算できるようになっている。

Hylleraas流の試行関数を使った変分計算は、最大で4電子系までしか行われてないようである。Hylleraasの試行関数は、電子間距離を変数に含むので、電子数の二乗オーダーで、変数の数が増えることになるし、互いに独立でもなくなる。なので、電子数が多い時には、それほど効率的でなくなるのかもしれない。

もう少し電子数が多い場合には、CI計算が行われる。例えば、1974年には、以下の論文が出てる。

Configuration-Interaction Study of Atoms. I Correlation Energies of B,C,N,O,F, and Ne
https://doi.org/10.1103/PhysRevA.9.17

Hylleraas型計算は、化学者が知りたい化合物の化学物性を知るのには使えないけど、それでも電子数の少ない系では、今でも一番精度の高い計算法で、21世紀になっても、時々論文が出ている。


エネルギーを何桁くらいまで計算すればいいか考える。愚かな人類が、2〜3桁も合ってれば、よさそうと考えがちなのは、何を測定するにしても、2〜3桁より高い精度で測定するのは大変なことが多いからだろうけど、基本的な物理定数の測定誤差に由来する不定性を何桁も超えて計算しても、あまり意味がない。測定値と比較することを考えれば、相対論的補正より細かい精度を要求するかは、一つの基準になりそうに思える。他にも、原子核が点電荷でないことから来る補正とかもあるけど、評価しづらいので保留。

てことで、相対論的補正が、どれくらいなのか概算する。原子核電荷Zとした時、Bohr速度は、v=\alpha Z cになる。\alphaは微細構造定数。\beta = v/c = Z \alphaとおくと、相対論的補正は
m_{e} c^2 \left( \dfrac{1}{\sqrt{1 - \beta^2}} - 1 - \dfrac{1}{2}\beta^2 \right)
くらいのオーダーだろうと思われる。例としてヘリウムの場合を考えると、Z=2だから、\betaは、2/137くらいで、相対論的補正は、0.01eVオーダーになる。リチウム原子の場合は、0.05eVくらい。

1eVは約96500(J/mol)で、100(kJ/mol)弱なので、0.01eVは、1(kJ/mol)のオーダー。リチウムの0.05eVは、1(kcal/mol)に近い。これは化学的精度と言われる数kcal/molくらいの水準に近いけど、最内殻電子の話で、通常の化学物性とかで影響が大きいのは外殻電子で、外殻電子では相対論的補正の影響は、もっと小さいはず。

現在の計算では、原子単位系を使うのが一般的で、原子単位系だと、0.05eVは、0.05/(13.6*2) ≒ 0.0018(hartree)なので、1.8ミリハートリーくらい。相対論的補正を無視する水準なら、ミリハートリーかサブミリハートーリー精度の計算が出来てれば、かなりいいだろうと考えられる。

常温300(K)が、エネルギー換算で2.5(kJ/mol)になって、ほぼ1ミリハートリー相当。これより弱いかもしれないけど知りたい場合がある相互作用として、分子間力がある。分子間力の大きさを概算する方法は思いつかないけど、文献値だと、He-Heで、30マイクロハートリーとかのオーダーになる。次に弱いNe-Neや、H2-H2で、0.1ミリハートリー程度。

超微細構造を見たい場合は、MHzオーダーまで計算の対象になり、参考値として、水素原子の2P状態の微細構造が10GHzオーダー、ラムシフトが1GHzオーダー。周波数では、1(kJ/mol)は、2.5(THz)くらい。1MHzは、0.1ナノハートリー程度。

というわけで、大雑把な目安として、
・単一分子の化学物性を見たいなら、ミリハートリー(これは、いわゆる化学的精度と同一オーダー)
・分子間相互作用を見たいなら、マイクロハートリー
・微細構造、超微細構造を見たいなら、ナノハートリー
くらいの計算精度が欲しいと考えられる。

現在のところ、超微細構造を見たい理由が、理論の検証以外にあるのか分からない。原子数の多い多原子分子とかで超微細構造まで見える精度の計算をしたい人は、あんまりいないだろう(簡単に計算できるなら知りたいかもしれないけど)けど、ヘリウム原子やリチウム原子だと、QED補正まで計算してる人もいる。そういう計算をする時は、非相対論的近似として、Hylleraas型試行関数を使った変分計算を出発点にしてることがある。



あと、無視されてることもあるけど、重心運動の自由度を分離した時に、核質量が有限であることから来る補正として、ヘリウムの場合は、mass polarization termというのが出る(何で、こんな名前なのかは分からない)。水素原子の場合は、mass polarizationはないけど、核質量が有限であることから、電子質量ではなく、換算質量を使うようにしないといけない。

水素原子の基底状態の電子エネルギーは、量子力学の教科書とか読めば、-13.6eVとよく書いてある。もうちょっと細かく見ると、0.5(Eh)=13.60569…(eV)で、Dirac方程式だと、13.60587…(eV)くらい。水素原子のイオン化エネルギーを、直接高精度で実験的に決定するのは難しいのか、あんまりデータがない。
Hydrogen atom@NIST Webbook
https://webbook.nist.gov/cgi/cbook.cgi?ID=C12385136&Units=CAL&Mask=20
には、イオン化エネルギーとして、13.59844(eV)という値が載ってるけど、methodがevaluationとなってるので、多分、計算値(?)。0.5(Eh)と比べると、0.26ミリハートーリー程度の差がある。電子質量でなく、換算質量を使ってシュレディンガー方程式の計算結果に放り込むと、13.59828…(eV)になって、記載値との差は6マイクロハートリー程度まで縮まる。

換算質量を使った方が文献値に近付くのが偶然でないことは、正電荷粒子が陽子より軽いポジトロニウムや、核質量が陽子より重い重水素などの水素様原子の基底状態エネルギーを調べると、はっきりと確認できるはず。逆に、NIST Webbookで見ると、重水素三重水素の実験値は、13.603eVとなってて、核質量が重い場合、電子質量をそのまま使った計算値に近づいている。外場近似で、原子核の質量補正を扱う方法は、arXiv:hep-ph/0002158("Theory of Light Hydrogen Atoms")によると、1967年に提示されたそうで、何やかんやで、α^4までの項をQEDで計算すると、基底状態のエネルギーは、-13.598434…(eV)とかになって、記載値と近付く。

ヘリウムの場合も、電子質量ではなく、換算質量を使うようにしないといけないけど、それに加えて、mass polarizationという項が出現する。Hylleraasの1929年の論文では、Massenkorrektion bei Heliumという項で、換算質量による影響を評価してるっぽいけど、mass polarizationは考慮されてない。

以下の論文によると、1933年、Betheが、ヘリウムに対して、mass polarizationの評価を行ったことが書いてある(論文を読む限り、Betheの計算には、ミスがあったのかもしれない)。

Lower Bound to the Ground-State Energy and Mass Polarization in Helium-Like Atoms
https://doi.org/10.1103/PhysRev.103.112

この論文の数値を信じるなら、mass polarizationの影響は、0.02(ミリハートリー)程度なのだと思われる。この論文は基底状態エネルギー下限の評価も行ってて、1932年のD.H. Weinsteinの論文に書いてあるのと同じ方法だけど、Hylleraas試行関数に適用したため、下限値を大きく改善している。

核質量補正を忘れると、数学的練習問題を解いてるだけになってしまう。


ヘリウム原子(Hylleraasの計算)

既に書いた通り、Hylleraasも、最初は、Kellnerと同様、線形試行関数を使った。その後、非線形試行関数を使って、ヘリウム原子の変分計算を実行した。後者の論文は、以下。
Neue Berechnung der Energie des Heliums im Grundzustande, sowie des tiefsten Terms von Ortho-Helium
https://doi.org/10.1007/BF01375457

1956年に、木下東一郎が、試行関数を少し拡張して、また電子計算機を使用して、精度を上げた。この論文も、よく引用されている。この人は、電子の異常磁気モーメントの高精度計算がライフワークみたいな人だったっぽい。
Ground State of the Helium Atom
https://doi.org/10.1103/PhysRev.105.1490

Ground State of the Helium Atom. II
https://doi.org/10.1103/PhysRev.115.366

木下の論文では、mass polarizationによる補正と相対論的補正の計算もされている。


ヘリウムの電子数は2個なので、6次元空間上のシュレディンガー方程式を考えることになるけど、回転対称性を利用して、変数を3つまで減らせる。1つ目、2つ目の電子の原点からの距離をr_1,r_2として、電子間距離を、r_{12}とする。Hylleraasは、s=r_{1}+r_{2},t=-r_{1}+r_{2},u=r_{12}という変数を使った。この変数で、シュレディンガー方程式を書き直すとかは面倒だけど、機械的な計算でできる。

Solving the electron-nuclear Schrödinger equation of helium atom and its isoelectronic ions with the free iterative-complement-interaction method
https://doi.org/10.1063/1.2904562

とかに、換算質量を使用し、mass polarizationも含む形で書いてある。ここでは、とりあえず、Hylleraasの結果を再現するために、換算質量も使わず、mass polarizationも考慮しない。

Hylleraasの使った試行関数は
\psi(s,t,u) = e^{-\alpha s} \displaystyle \sum_{p,q,r=0}^{\infty} c_{p,q,r} s^{p} t^{q} u^{r}
という形のものと本質的には同じ(但し、qが奇数の時、係数は0)。uの奇数次の項は、原理的にはなくてもいい(つまり収束自体はする)けど、ないと収束は遅くなる。

非線形パラメータは一個だけで、このまま計算してもいいけど、Hylleraasは、多分、当時まだ電子計算機がなかったこともあって、もう少し工夫して
\phi(s,t,u) = e^{-\dfrac{s}{2}} \displaystyle \sum_{p,q,r=0}^{\infty} c_{p,q,r} s^{p} t^{q} u^{r}
\psi(s,t,u) = \phi(ks,kt,ku)
という形の関数にしている。こうすると、変分関数は、kの二乗に比例する項と一乗に比例する項に分離する。仮に、
k^2 A - k B =\mathrm{min.}
という形になってたら、k = \dfrac{B}{2 A}となるので、非線形パラメータを決定しなくていい。そして、変分関数は
k^2 A - k B = \dfrac{-B^2}{4A}
に帰着する。A,Bは、線形パラメータc_{p,q,r}に依存する関数。実際に計算すると、結局、4次の同次式/4次の同次式という形になって、それほど簡単でもないけど。

Hylleraasが当時、どういう方法で、非線形最適化を計算したのか、論文にも詳細が書かれてなくて、ハッキリとはしないが、何らかの勾配法を使ったと思われる。線形パラメータは、定数項以外は、0に近いので、比較的見つけやすかったのだと思う。

変分関数の計算には
\phi_{i} = e^{-\dfrac{s}{2}} s^{p_i} t^{q_i} u^{r_i}
\phi_{j} = e^{-\dfrac{s}{2}} s^{p_j} t^{q_j} u^{r_j}
として、以下の積分を計算する必要がある。
M_{ij} = \displaystyle \int_{0}^{\infty} ds \int_{0}^{s} du \int_{0}^{u} dt \left( u(s^2-t^2)(\dfrac{\partial \phi_i}{\partial s} \dfrac{\partial \phi_j}{\partial s} + \dfrac{\partial \phi_i}{\partial t} \dfrac{\partial \phi_j}{\partial t} + \dfrac{\partial \phi_i}{\partial u} \dfrac{\partial \phi_j}{\partial u}) + 2s(u^2-t^2)\dfrac{\partial \phi_i}{\partial s} \dfrac{\partial \phi_j}{\partial u} + 2t(s^2-u^2)\dfrac{\partial \phi_i}{\partial t} \dfrac{\partial \phi_j}{\partial u} \right)
L_{ij} = \displaystyle \int_{0}^{\infty} ds \int_{0}^{s} du \int_{0}^{u} dt (2 s u - \dfrac{s^2 - t^2}{4})\phi_i \phi_j
N_{ij} = \displaystyle \int_{0}^{\infty} ds \int_{0}^{s} du \int_{0}^{u} dt u(s^2-t^2) \phi_i \phi_j
これらの積分値は全部解析的に出るけど、書くと長いので省略。


以下に、Hylleraasの計算を追試するための実装を作った。非線形最適化は、scipy.optimizeに丸投げ。Hylleraasの論文の式(25)には、-1.4516という数値が書いてあって、Hylleraasは、4(Ry)が1になるように、エネルギーを無次元化している(式(2)の直後)ので、-2.9032(hartree)程度に相当する。一方、同じ条件(2次以下の6項を使用)で、私の計算結果は、-2.903 329(hartree)となり、合ってるっぽい。

4次までの項を全部入れると、-2.903 713(hartree)になる(scipyが収束してるか怪しいけど)。厳密値は、-2.903 724...(hartree)になるらしいけど、Hylleraasの試行関数では、項数を増やしても、なかなか改善しない。木下論文には、39項使用した最良の結果として、2.903 722 5という数値が書いてある。

上の計算値は、核質量補正が入ってないけど、
Ionization Potential of the Helium Atom
https://arxiv.org/abs/physics/0105108
には、非相対論的エネルギーの値として、2.903 304 557 727 940 23(1)というのが記載されている。核質量補正ありだと、こんなものになるのだろう。

Atomic Spectra Database
https://dx.doi.org/10.18434/T4W30F
で実験値を調べると、2.903 385 87...(hartree)となるので、核補正なしでも、誤差1ミリハートーリー未満は達成できている。

import math
import numpy as np
import scipy.optimize

"""
\int_{0}^_{\infty} ds \int_{0}^{s} du \int_{0}^{u} dt e^{-s} a^{p} t^{q} u^{r}
"""
def T(p,q,r):
   if p<0 or q<0 or r<0:return 0
   return math.factorial(p+q+r+2)/((q+1)*(q+r+2))


def N(I,J):
   ret = T(I[0]+J[0]+2 , I[1]+J[1] , I[2]+J[2]+1)
   ret -= T(I[0]+J[0] , I[1]+J[1]+2 , I[2]+J[2]+1)
   return ret


def L(I,J):
   ret = 2*T(I[0]+J[0]+1 , I[1]+J[1] , I[2]+J[2]+1)
   ret -= T(I[0]+J[0]+2 , I[1]+J[1] , I[2]+J[2])/4
   ret += T(I[0]+J[0] , I[1]+J[1]+2 , I[2]+J[2])/4
   return ret


def M(I,J):
   pi,qi,ri = I
   pj,qj,rj = J
   alpha = 0.5
   ret = alpha*alpha*( T(pi+pj+2 , qi+qj , ri+rj+1) - T(pi+pj , qi+qj+2 , ri+rj+1))
   ret -= alpha*(pi+pj)*(T(pi+pj+1 , qi+qj , ri+rj+1) - T(pi+pj-1 , qi+qj+2 , ri+rj+1))
   ret += pi*pj*(T(pi+pj , qi+qj , ri+rj+1) - T(pi+pj-2 , qi+qj+2 , ri+rj+1))
   ret += qi*qj*(T(pi+pj+2 , qi+qj-2 , ri+rj+1) - T(pi+pj , qi+qj , ri+rj+1))
   ret += ri*rj*(T(pi+pj+2 , qi+qj , ri+rj-1) - T(pi+pj , qi+qj+2 , ri+rj-1))
   ret -= 2*alpha*rj * (T(pi+pj+1,qi+qj,ri+rj+1) - T(pi+pj+1,qi+qj+2,ri+rj-1))
   ret += 2*pi*rj*(T(pi+pj,qi+qj,ri+rj+1) - T(pi+pj,qi+qj+2,ri+rj-1))
   ret += 2*qi*rj*(T(pi+pj+2,qi+qj,ri+rj-1) -T(pi+pj,qi+qj,ri+rj+1))
   return ret


## -- calculate helium ground state energy from Hylleraas-type trial function
class Hyl:
   def __init__(self , indices):
      n = len(indices)
      self.Lm = np.zeros( (n,n) )
      self.Mm = np.zeros( (n,n) )
      self.Nm = np.zeros( (n,n) )
      for i,I in enumerate(indices):
         for j,J in enumerate(indices):
              self.Lm[i,j] = L(I,J)
              self.Mm[i,j] = M(I,J)
              self.Nm[i,j] = N(I,J)
   def __call__(self , cs): 
       Lval = np.dot(np.dot(self.Lm , cs) , cs)
       Mval = np.dot(np.dot(self.Mm , cs) , cs)
       Nval = np.dot(np.dot(self.Nm , cs) , cs)
       return -4*Lval*Lval/(Mval*Nval)  #-- returns in hartree


def generate_indices(max_sum):
   ret = []
   for loc_sum in range(max_sum+1):
       for p in range(loc_sum+1):
           rest_sum = loc_sum - p
           for q in range(rest_sum+1):
               r = rest_sum - q
               if q%2==0:ret.append( (p,q,r) )
   return ret


## see Hyleraas1929 (22)-(25)
if __name__=="__main__":
   indices =[ (0,0,0) , (0,0,1),(0,2,0),(1,0,0),(2,0,0),(0,0,2)]
   #indices = generate_indices(4)
   h = Hyl(indices)
   x0 = np.zeros(len(indices))
   x0[0]=1
   res = scipy.optimize.minimize(h , x0,  options={"maxiter":10000} )
   print(h(res.x))