⚛️

モンテカルロ積分を用いたヘリウム原子の摂動計算

2023/12/12に公開

Julia

Julia Advent Calendar 2023 12日目です🎄

はじめに

ここではヘリウム原子の電子状態計算を題材として摂動法を学びます. 解析的な計算を避けて気軽に計算できるよう波動関数にはAntique.jlを使用し, 積分にはMonteCarloIntegration.jlを使用します. なお, 摂動論のちゃんとした説明は東京工業大学の武藤先生の講義ノートを参照してください.

パッケージの使い方

ヘリウム原子 $\mathrm{He}$ を扱う前に, まずはヘリウム原子イオン $\mathrm{He}^+$ を題材としてパッケージの使い方を学びましょう. $\mathrm{He}^+$ は $+2$ の電荷をもった核と $-1$ の電荷をもった電子が1つずつ集まってできている, いわゆる水素様原子（すいそようげんし）です. 水素様原子の固有値・固有関数はAntique.jlで利用できます. まず下記を実行してインストールしましょう.

入力

using Pkg; Pkg.add("Antique")

$\mathrm{He}^+$ は原子番号（核の電荷）が $Z=2$ なので, こちらを例にしてZ=2を指定します.

入力

using Antique
He⁺ = antique(:HydrogenAtom, Z=2, Eₕ=1.0, a₀=1.0, mₑ=1.0, ℏ=1.0)

核の質量が十分に大きいときの水素様原子では, 主量子数 $n \geq 1$ に対する固有値は

E_n = -\frac{Z^2}{2 n^2} E_{\mathrm{h}}

で与えられるので, $\mathrm{He}^+$ の基底状態, つまり $Z=2, n=1$ では

E_1 = -\frac{2^2}{2 1^2} E_{\mathrm{h}} = -2 E_{\mathrm{h}}

となります. これは先ほど定義したモジュールHe⁺の関数E()にn=1を指定して次のように計算できます.

入力

He⁺.E(n=1)

出力

-2.0

エネルギーの単位について, 2018年CODATA推奨値の換算係数を用いると次のような値に換算できます.

入力

Eₕ2eV = 27.211386245988 # eV # https://physics.nist.gov/cgi-bin/cuu/Value?hrev
println("E = ", He⁺.E(n=1), " Eₕ")
println("  = ", He⁺.E(n=1) * Eₕ2eV, " eV")

出力

E = -2.0 Eₕ
  = -54.422772491976 eV

次に, 波動関数の規格化条件を確認しましょう. 積分にはMonteCarloIntegration.jlを使用するので次のコマンドを実行してインストールします.

入力

using Pkg; Pkg.add("MonteCarloIntegration")

このパッケージで次の積分を計算する例です.

\int |\psi_n(r,\theta,\varphi)|^2 r^2 \sin\theta~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi = 1

入力

using MonteCarloIntegration
res = vegas(r -> abs(He⁺.ψ(r...))^2 * r[1]^2 * sin(r[2]), [0,0,0], [100,π,2π], atol=1e-5)
println("<ψ|ψ> = ", res.integral_estimate, " ± ", res.standard_deviation)

出力

<ψ|ψ> = 0.9983502045775035 ± 0.0007877199112562297

おおよそ1が得られましたので, 正しく計算できると考えられます.

プログラムの説明

プログラムについて補足します. まずr[1]は $r$ , r[2]は $\theta$ , r[3]は $\varphi$ に対応します. これだけ説明すれば後は察せられると思いますが, 本来ならば $0\leq r \lt \infty, 0\leq \theta \lt \pi, 0\leq \varphi \lt 2\pi$ を計算するところを, $r$ に関しては無限遠までサンプリングすることはできないので適当に $r=100$ 程度までで打ち切っています. また, He⁺.ψ(r...)はsplat記法を用いたHe⁺.ψ(r[1],r[2],r[3])の略記です.

ハミルトニアン

ハミルトニアンの導出

ヘリウム原子を2つの電子と1つの原子核（アルファ粒子＝中性子2つと陽子2つ）からできた3体系として扱います. 電子と核はいずれも質点と考え, 相互作用はクーロン力のみとすれば, 非相対論的なハミルトニアンは次式で書き表されます.

\hat{H}_\mathrm{tot}
=
- \frac{\hbar^2}{2 m_\mathrm{\alpha}} \nabla_{\alpha}^2
- \frac{\hbar^2}{2 m_\mathrm{e}} \nabla_{1}^2
- \frac{\hbar^2}{2 m_\mathrm{e}} \nabla_{2}^2
- \frac{2 e^2}{4 \pi \varepsilon r_1}
- \frac{2 e^2}{4 \pi \varepsilon r_2}
+ \frac{  e^2}{4 \pi \varepsilon r_{12}}
これを原子単位（詳細はIUPAC Greenbook 3.2.9を参照）で無次元化すると

\hat{H}_\mathrm{tot}^\ast =
- \frac{1}{2 m_\mathrm{\alpha}/m_\mathrm{e}} \nabla_{\alpha}^{2\ast}
- \frac{1}{2} \nabla_{1}^{2\ast}
- \frac{1}{2} \nabla_{2}^{2\ast}
- \frac{2}{r_1^\ast}
- \frac{2}{r_2^\ast}
+ \frac{1}{r_{12}^\ast}
となります. ここで, ハミルトニアンはハートリーエネルギーE_\mathrm{h}を用いて無次元化されており, \hat{H}_\mathrm{tot} / E_\mathrm{h} = \hat{H}_\mathrm{tot}^\astという関係で結ばれます. また, 2018年CODATA推奨値よりm_\mathrm{\alpha}/m_\mathrm{e}=7294.299~541~42(24)が利用できますが, 今回は核の量子効果を無視（m_\mathrm{\alpha}\rightarrow\infty）します. さらに, 以降はこれらのアスタリスクを省略するものとすると, 次のハミルトニアンが得られます.

\hat{H} = \hat{H}_0 + \hat{H}' = - \frac{1}{2} \nabla_{1}^{2} - \frac{1}{2} \nabla_{2}^{2} - \frac{2}{r_1} - \frac{2}{r_2} + \frac{1}{r_{12}}

このノートでは, これを非摂動ハミルトニアン

\hat{H}_0 = - \frac{1}{2} \nabla_{1}^{2} - \frac{1}{2} \nabla_{2}^{2} - \frac{2}{r_1} - \frac{2}{r_2}

と摂動項

\hat{H}' = \frac{1}{r_{12}}

に分けて考えていきます.

非摂動のエネルギーと波動関数

摂動論では非摂動の問題は既に解けているという前提から始まりますので, まずは非摂動（電子間反発を無視した）ハミルトニアン $\hat{H}_0$ の固有値と固有関数を用意します. よく見ると非摂動ハミルトニアン $\hat{H}_0$ は2つのヘリウムイオンHe⁺のハミルトニアン

\hat{h}_1 = - \frac{1}{2} \nabla_{1}^2 - \frac{2}{r_1}

\hat{h}_2 = - \frac{1}{2} \nabla_{2}^2 - \frac{2}{r_2}

の和になっています. これらの固有値をそれぞれ $E_1$ と $E_2$ , 固有関数をそれぞれ $\psi_1(\pmb{r}_1)$ と $\psi_2(\pmb{r}_2)$ とすると, 2次元の井戸型ポテンシャルや2次元の調和振動子の例などの類推から直ちに非摂動ハミルトニアン $\hat{H}_0$ の固有値は $E_1 + E_2$ , 固有関数は $\psi_1(\pmb{r}_1) \psi_2(\pmb{r}_2)$ であると予想されます. 実際,

\hat{h}_1 \psi_1 = E_1 \psi_1

\hat{h}_2 \psi_2 = E_2 \psi_2

の2式を用いれば

\begin{aligned} \hat{H}_0 \psi_1 \psi_2 &= \left( \hat{h}_1 + \hat{h}_2 \right) \psi_1 \psi_2 \\ &= \hat{h}_1 \psi_1 \psi_2 + \hat{h}_2 \psi_1 \psi_2 \\ &= \psi_2 \hat{h}_1 \psi_1 + \psi_1 \hat{h}_2 \psi_2 \\ &= \psi_2 E_1 \psi_1 + \psi_1 E_2 \psi_2 \\ &= E_1 \psi_1 \psi_2 + E_2 \psi_1 \psi_2 \\ &= \left( E_1 + E_2 \right) \psi_1 \psi_2 \end{aligned}

となり, 非摂動ハミルトニアン $\hat{H}_0$ について $\psi_1 \psi_2$ が固有関数となり, $E_1 + E_2$ が固有値となることが示せます.

波動関数についての注意

しかし, Schrödinger方程式を満足するだけではダメです. 電子（フェルミ粒子）が2つ以上含まれる系では, それらの空間座標およびスピン座標の置換に対して波動関数は反対称にならなければいけません. これはSchrödinger方程式とは独立した公理であり, Schrödinger方程式を満足するだけは自動的に満足されないため, あらかじめ対称性を満足するような基底関数のみを用いるなどして計算を行わなければいけません. この辺りの話はA.ザボ, N.S.オストランド著, 大野公男, 阪井健男, 望月祐志訳『新しい量子化学上　電子構造の理論入門』(東京大学出版会, 1987)の2.1.3に書いてあります. とはいえ, 今回は $n_1=n_2=1$ に限定しており, 空間部分が対称でスピン部分が反対称（パラヘリウム）の場合を考えているので $\psi_1(\pmb{r}_1) \psi_2(\pmb{r}_2)$ で問題ありません. 規格化因子が間違っている気がしますが, Wikipediaに同様の説明が書いてあります.

さて, 具体的なエネルギーを求めましょう. 単純な和ですので, 下記のように計算できます.

入力

Eₕ2eV = 27.211386245988 # eV # https://physics.nist.gov/cgi-bin/cuu/Value?hrev
println("E = ", (He⁺.E(n=1) + He⁺.E(n=1)), " Eₕ")
println("  = ", (He⁺.E(n=1) + He⁺.E(n=1)) * Eₕ2eV, " eV")

出力

E = -4.0 Eₕ
  = -108.845544983952 eV

次に波動関数の規格化条件を確認しましょう. 波動関数は6変数関数として次のように定義します.

入力

Ψ(r₁,θ₁,φ₁,r₂,θ₂,φ₂) =  He⁺.ψ(r₁,θ₁,φ₁) * He⁺.ψ(r₂,θ₂,φ₂)

先ほどと同様にモンテカルロ積分を用いて規格化条件を確認しましょう.

入力

res = vegas(x ->
        x[1]^2 * sin(x[2]) *
        x[4]^2 * sin(x[5]) *
        abs(Ψ(x...))^2 ,
        [0,0,0,0,0,0],
        [50,π,2π,50,π,2π],
        atol=1e-5
    )

println("<ψ|ψ> = ", res.integral_estimate, " ± ", res.standard_deviation)

出力

<ψ|ψ> = 1.0004144350516748 ± 0.0017535335574600287

おおよそ1となることが確かめられました.

第一近似のエネルギー

先ほど用意した非摂動波動関数 $\psi_1 \psi_2$ は非摂動ハミルトニアンの厳密な固有関数になっていますが, 摂動項を加えたハミルトニアンの固有関数にはなっていません. しかし, まあ近い形になっているだろうということで $\psi_1 \psi_2$ を使って摂動項（電子間反発）の期待値を計算します. こちらの講義ノートの14.3.2に書いてある第一近似のエネルギーの計算に相当します.

波動関数が球面座標で定義されているのでデカルト座標に変換してから電子間距離を計算します.（球面座標から直接計算する方法もあるかもしれません.）

入力

x(r,θ,φ) = r*sin(θ)*cos(φ)
y(r,θ,φ) = r*sin(θ)*sin(φ)
z(r,θ,φ) = r*cos(θ)

function r₁₂(r₁,θ₁,φ₁,r₂,θ₂,φ₂)
    return sqrt(
        (x(r₁,θ₁,φ₁) - x(r₂,θ₂,φ₂)) ^ 2 + 
        (y(r₁,θ₁,φ₁) - y(r₂,θ₂,φ₂)) ^ 2 + 
        (z(r₁,θ₁,φ₁) - z(r₂,θ₂,φ₂)) ^ 2
    )
end

これを先ほどの規格化条件の確認をした積分のコードに挟み込むと, 摂動項（電子間反発）の期待値が計算できます.

入力

res = vegas(x ->
        x[1]^2 * sin(x[2]) *
        x[4]^2 * sin(x[5]) *
        conj(Ψ(x...)) *
        1 / r₁₂(x...) *
        Ψ(x...) ,
        [0,0,0,0,0,0],
        [50,π,2π,50,π,2π],
        rtol = 1e-6,
        atol = 1e-6
    )

∫ψψV₁₂ψψdr₁dr₂ = res.integral_estimate

println("<V₁₂> = ", ∫ψψV₁₂ψψdr₁dr₂, " Eₕ")
println("      = ", ∫ψψV₁₂ψψdr₁dr₂ * Eₕ2eV, " eV")

出力

<V₁₂> = 1.2495678829756485 Eₕ
      = 34.0024743042319 eV

これを先ほどの $-4 E_\mathrm{h}$ に足すと, 摂動項も含めた全エネルギーが次のように計算できます.

入力

E = He⁺.E(n=1) + He⁺.E(n=1) + ∫ψψV₁₂ψψdr₁dr₂

println("<H> = ", E, " Eₕ")
println("    = ", E * Eₕ2eV, " eV")

出力

<H> = -2.7504321170243515 Eₕ
    = -74.8430706797201 eV

これは現在知られている最も正確なエネルギー $-2.903~724~377 034~119~598~311~159~245~194~4~E_\mathrm{h}$ と比べると高いですが, 全くおかしな値ではないことが確認できます.

変分原理からの考察

変分原理の観点から言えば, 摂動項も含めたハミルトニアンにとって非摂動波動関数は的外れではないものの, あまり正確な波動関数にはなっていないであろうということです. なお, 変分原理の意味するところは「エネルギーが低くければ低くいほどエネルギーが正確である」という意味です. 他の期待値も正確になるとは限りませんので, 波動関数の良し悪しを判断する絶対的な基準ではないということに注意してください. 詳しくはこちらの解説をご覧ください.

おわりに

摂動論は第一近似だけ考えて終わりではなく, より高次の補正を加えていくアプローチです. ここでは, 非摂動波動関数を用いて摂動項の期待値を計算すること（最低次の摂動計算）の正当性が確認できたということで満足するものとします.