Perturbation Calculation of Helium Atom Using Monte Carlo Integration

This is the 12th day of the Julia Advent Calendar 2023 🎄

Introduction

In this article, we will learn about the perturbation method using the electronic state calculation of a helium atom as a case study. To keep the calculations light and avoid analytical derivations, we will use Antique.jl for wavefunctions and MonteCarloIntegration.jl for integration. For a formal explanation of perturbation theory, please refer to Professor Muto's lecture notes from Tokyo Institute of Technology.

Package Usage

Before dealing with the helium atom $\mathrm{He}$, let's first learn how to use the packages using the helium atom ion $\mathrm{He}^+$ as a subject. $\mathrm{He}^+$ is what is known as a hydrogen-like atom, composed of a nucleus with a $+2$ charge and a single electron with a $-1$ charge. The eigenvalues and eigenfunctions for hydrogen-like atoms are available in Antique.jl. First, run the following to install it.

Input

using Pkg; Pkg.add("Antique")

Since $\mathrm{He}^+$ has an atomic number (nuclear charge) of $Z=2$, we specify Z=2 following this example.

Input

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

In a hydrogen-like atom where the nuclear mass is sufficiently large, the eigenvalues for the principal quantum number $n \geq 1$ are given by

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

For the ground state of $\mathrm{He}^+$, i.e., $Z=2, n=1$, it becomes

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

This can be calculated by specifying n=1 in the function E() of the previously defined module He⁺ as follows:

Input

He⁺.E(n=1)

Output

-2.0

Regarding energy units, using the conversion factors from the 2018 CODATA recommended values, it can be converted to the following value:

Input

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")

Output

E = -2.0 Eₕ
  = -54.422772491976 eV

Next, let's check the normalization condition of the wavefunction. We will use MonteCarloIntegration.jl for integration, so install it using the following command.

Input

using Pkg; Pkg.add("MonteCarloIntegration")

Here is an example of calculating the following integral with this package:

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

Input

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)

Output

<ψ|ψ> = 0.9983502045775035 ± 0.0007877199112562297

Since we obtained approximately 1, we can consider the calculation to be correct.

Explanation of the program

Here is some supplementary information about the program. First, r[1] corresponds to $r$, r[2] to $\theta$, and r[3] to $\varphi$. With this explained, the rest should be self-evident, but while we should ideally integrate over $0\leq r \lt \infty, 0\leq \theta \lt \pi, 0\leq \varphi \lt 2\pi$, we cannot sample $r$ to infinity, so we cut it off at approximately $r=100$. Also, He⁺.ψ(r...) is a shorthand for He⁺.ψ(r[1],r[2],r[3]) using the splat syntax.

Hamiltonian

Derivation of the Hamiltonian

We treat the helium atom as a three-body system consisting of two electrons and one nucleus (alpha particle = two neutrons and two protons). Considering both electrons and the nucleus as point masses and assuming only Coulomb forces for interactions, the non-relativistic Hamiltonian is expressed as follows:

$$\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}}$$

Converting this to dimensionless form using atomic units (for details, see 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}$$

Here, the Hamiltonian is made dimensionless using the Hartree energy $E_\mathrm{h}$, with the relationship $\hat{H}_\mathrm{tot} / E_\mathrm{h} = \hat{H}_\mathrm{tot}^\ast$. Furthermore, while $m_\mathrm{\alpha}/m_\mathrm{e}=7294.299~541~42(24)$ is available from the 2018 CODATA recommended values, we will ignore nuclear quantum effects ($m_\mathrm{\alpha}\rightarrow\infty$) this time. Furthermore, omitting these asterisks hereafter, we obtain the following Hamiltonian.

$$\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}}$$

In this note, we will consider this by dividing it into the unperturbed Hamiltonian

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

and the perturbation term

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

Unperturbed Energy and Wavefunctions

Perturbation theory starts with the premise that the unperturbed problem has already been solved, so we first prepare the eigenvalues and eigenfunctions of the unperturbed Hamiltonian $\hat{H}_0$ (ignoring electron-electron repulsion). Looking closely, the unperturbed Hamiltonian $\hat{H}_0$ is the sum of the Hamiltonians for two helium ions $\mathrm{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}$$

If we let these eigenvalues be $E_1$ and $E_2$, and their respective eigenfunctions be $\psi_1(\pmb{r}_1)$ and $\psi_2(\pmb{r}_2)$, we can immediately expect from analogies such as examples of 2D potential wells or 2D harmonic oscillators that the eigenvalues of the unperturbed Hamiltonian $\hat{H}_0$ will be $E_1 + E_2$ and its eigenfunctions will be $\psi_1(\pmb{r}_1) \psi_2(\pmb{r}_2)$. In fact,

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

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

Using these two equations, we can show that:

$$\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}$$

Thus, it can be shown that $\psi_1 \psi_2$ is an eigenfunction of the unperturbed Hamiltonian $\hat{H}_0$ with the eigenvalue $E_1 + E_2$.

Note on wavefunctions

However, simply satisfying the Schrödinger equation is not enough. In a system containing two or more electrons (fermions), the wavefunction must be antisymmetric with respect to the exchange of their spatial and spin coordinates. This is an axiom independent of the Schrödinger equation and is not automatically satisfied just by solving the Schrödinger equation; therefore, calculations must be performed using basis functions that satisfy this symmetry in advance. Discussions on this can be found in Section 2.1.3 of ["Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory" by A. Szabo and N. S. Ostlund] (Japanese translation by Kimio Ohno, Takeo Sakai, and Yuji Mochizuki, University of Tokyo Press, 1987). That said, in this case, we are limited to $n_1=n_2=1$ and considering the case where the spatial part is symmetric and the spin part is antisymmetric (parahelium), so $\psi_1(\pmb{r}_1) \psi_2(\pmb{r}_2)$ is acceptable. While it feels like the normalization factor might be off, a similar explanation can be found on Wikipedia.

Now, let's determine the specific energy. Since it is a simple sum, it can be calculated as follows:

Input

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")

Output

E = -4.0 Eₕ
  = -108.845544983952 eV

Next, let's check the normalization condition of the wavefunction. The wavefunction is defined as a function of six variables as follows:

Input

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

Let's check the normalization condition using Monte Carlo integration as before.

Input

res = vegas(x -2
        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)

Output

<ψ|ψ> = 1.0004144350516748 ± 0.0017535335574600287

It has been confirmed that the result is approximately 1.

First-order Energy

The unperturbed wavefunction $\psi_1 \psi_2$ prepared earlier is an exact eigenfunction of the unperturbed Hamiltonian, but it is not an eigenfunction of the Hamiltonian once the perturbation term is added. However, assuming it is a reasonable approximation, we calculate the expectation value of the perturbation term (electron-electron repulsion) using $\psi_1 \psi_2$. This corresponds to the calculation of the first-order energy approximation described in 14.3.2 of these lecture notes.

Since the wavefunction is defined in spherical coordinates, we convert it to Cartesian coordinates before calculating the distance between the electrons (there might be a way to calculate it directly from spherical coordinates).

Input

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

By inserting this into the integration code used earlier to verify the normalization condition, we can calculate the expectation value of the perturbation term (electron-electron repulsion).

Input

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")

Output

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

Adding this to the previous $-4 E_{\mathrm{h}}$, the total energy including the perturbation term is calculated as follows:

Input

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

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

Output

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

While this is higher than the currently most accurate known energy $-2.903~724~377 034~119~598~311~159~245~194~4~E_{\mathrm{h}}$, we can confirm it is not an entirely unreasonable value.

Considerations from the variational principle

From the perspective of the variational principle, while the unperturbed wavefunction is not completely off the mark for the Hamiltonian including the perturbation term, it is likely not a very accurate wavefunction. Note that the variational principle implies that "the lower the energy, the more accurate the energy is." Since other expectation values are not guaranteed to be more accurate, be aware that energy is not an absolute criterion for judging the overall quality of a wavefunction. For more details, please refer to this explanation.

Conclusion

Perturbation theory is not an approach that ends after considering only the first-order approximation; it is an approach that involves adding higher-order corrections. In this article, we shall be satisfied with having confirmed the validity of calculating the expectation value of the perturbation term using unperturbed wavefunctions (the lowest-order perturbation calculation).

References

• Professor Muto's lecture notes from Tokyo Institute of Technology
• Antique.jl
• MonteCarloIntegration.jl
• IUPAC Greenbook 3.2.9
• 2018 CODATA recommended values
• A. Szabo and N. S. Ostlund, "Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory" (Translated by Kimio Ohno, Takeo Sakai, and Yuji Mochizuki, University of Tokyo Press, 1987)
• H. Nakashima, H. Nakatsuji, Phys. Rev. Lett. 101, 240406 (2008)
• J. F. Pérez-Torres, J. Chem. Educ. 96, 4, 704–707 (2019)

https://gist.github.com/ohno/d63d79ca816a3614b201bc5728399fb7