# Solving the Helium Eigenvalue Problem—where to start?

**URL:** <https://discourse.julialang.org/t/solving-the-helium-eigenvalue-problem-where-to-start/137783>\
**Category:** Specific Domains\
**Tags:** question, quantum, physics, eigenvalues\
**Created:** [June 24, 2026, 5:37pm UTC](https://discourse.julialang.org/t/solving-the-helium-eigenvalue-problem-where-to-start/137783 "2026-06-24T17:37:36Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![ducksoverip](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ducksoverip/32/31967_2.png) [@ducksoverip](https://discourse.julialang.org/u/ducksoverip)\
**Post date:** [June 24, 2026, 5:37pm UTC](https://discourse.julialang.org/t/solving-the-helium-eigenvalue-problem-where-to-start/137783/1 "2026-06-24T17:37:36Z")

</div>

I have inherited Fortran code (the same one in [this paper](http://dx.doi.org/10.1103/PhysRevE.81.056705)) that solves the time-independent Schrodinger equation for the ground state of the helium atom. This is generally a very difficult problem because the Schrodinger equation for helium is a non-separable 6-D PDE:

\hat{H}\psi(\vec{r\_1}, \vec{r\_2}) = E \psi(\vec{r\_1}, \vec{r\_2}),

where E is a constant corresponding to the total energy of the system, and \hat{H} is an operator defined as follows:

\nabla^2\_1 + \nabla^2\_2 - \frac{2}{r\_1} - \frac{2}{r\_2} + \frac{1}{|\vec{r\_1} - \vec{r\_2}|}

The Fortran code that I have reduces this to a 2-D problem by expanding \psi in terms of coupled spherical harmonics, then summing over the different parameters of those harmonics—this is referred to as the method of partial waves. The two radial coordinates r\_1 and r\_2 are then mapped to a variable finite-element grid. Finally, the program solves the problem by solving the time-_dependent_ Schrodinger equation (TDSE) instead (\hat{H}\psi(\vec{r\_1}, \vec{r\_2}) = i \frac{\partial}{\partial t} \psi(\vec{r\_1}, \vec{r\_2})), but with an imaginary time parameter. This is because if you naively solve the TDSE by separating the derivative, you obtain the following:

\psi(\vec{r\_1}, \vec{r\_2}, t + dt) = e^{i \hat{H}dt} \psi(\vec{r\_1}, \vec{r\_2}, t)

If you change t to \tau = it, then the matrix exponential above becomes real and decaying. If \psi is a convenient test function (say, a Gaussian) that contains a mix of all possible eigenfunctions, then since e^{i \hat{H}d\tau} returns e^{i E\_n d\tau} when applied to any given eigenfunction \psi\_n, those states with higher E\_n decay to zero more quickly as the wavefunction is propagated, leaving only the lowest-E\_n state, which is the ground state. In-code, this propagation is done via the split-operator method, Fourier-transforming the system such that each operator matrix is diagonal when applied, except for the coupling term \frac{1}{|\vec{r\_1} - \vec{r\_2}|} . The code I have currently takes a couple minutes to run on 100 cores on an HPC cluster and returns both the right ground state energy (-2.903 au) and the ground-state wavefunction.

So, if I have code that works, runs quickly, and I know how to use it, why am I trying to reinvent the wheel? I’m glad you asked, imaginary curious interlocutor. There’s several reasons:

1. The code is nigh-unmaintainable. There’s minimal comments, obscure-at-best variable names, and this is on its 4th generation of being handed down, so I can’t easily ask questions about how it works or what effects making changes would have.
2. The code is 20 years old. This is no bad thing in itself, but in the intervening decades, as computers have gotten better and people have put more and more time into efficiently solving difficult PDEs, I’d like to explore what new methods are available to me.
3. I want to understand how my own code works. This is more about my own learning than practicality, but I’d like to move beyond plugging-and-chugging someone else’s code and towards understanding how to actually solve these kinds of problems, so that I can better extend it and explain it.

I know that Julia has many options for solving PDEs, but it’s hard for me to figure out which methods are right for solving this kind of problem, especially since the choices made for the code I have were determined by constraints 20 years ago. Finite differences were too inaccurate and spectral methods took too long, but I don’t know if that’s changed. For example, MethodOfLines.jl says it can handle spherical Laplacians, but I can’t find further documentation on that. None of the packages I’ve looked at (MethodOfLines, Ferrite, Gridap) talk about coupled terms like the one I’m dealing with. I would greatly appreciate some direction from one of the fine linear algebra experts on this forum. Thank you!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 24, 2026, 6:01pm UTC](https://discourse.julialang.org/t/solving-the-helium-eigenvalue-problem-where-to-start/137783/2 "2026-06-24T18:01:30Z")

</div>

> [@ducksoverip](#):
>
> The Fortran code that I have reduces this to a 2-D problem by expanding \psi 𝜓 in terms of coupled spherical harmonics, then summing over the different parameters of those harmonics—this is referred to as the method of partial waves.

In numerics this would typically be called a “spectral” method. (It’s not really “2d” in the sense that the number of spherical harmonics that you need is still determined by the 3d resolution, but you can get exponentially fast convergence with the number of harmonics.)

> [@ducksoverip](#):
>
> If you change t to \tau = it, then the matrix exponential above becomes real and decaying. If \psi is a convenient test function (say, a Gaussian) that contains a mix of all possible eigenfunctions, then since e^{i \hat{H}d\tau} returns e^{i E\_n d\tau} when applied to any given eigenfunction \psi\_n, those states with higher E\_n decay to zero more quickly as the wavefunction is propagated, leaving only the lowest-E\_n state, which is the ground state.

My understanding is that “imaginary time evolution” like this has mostly been superseded by implicitly restarted Lanczos, LOBPCG, and similar Krylov-type methods.

> [@ducksoverip](#):
>
> In-code, this propagation is done via the split-operator method, Fourier-transforming the system such that each operator matrix is diagonal when applied, except for the coupling term \frac{1}{|\vec{r\_1} - \vec{r\_2}|}

You can still use this in a Krylov method in order to do your matrix–vector multiplications. There are packages like FastSphericalHarmonics.jl which can help.

> [@ducksoverip](#):
>
> Finite differences were too inaccurate and spectral methods took too long, but I don’t know if that’s changed.

The method of “partial waves” _is_ a spectral method.

---

<div class="post-metadata">

**Author:** ![ducksoverip](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ducksoverip/32/31967_2.png) [@ducksoverip](https://discourse.julialang.org/u/ducksoverip)\
**Post date:** [June 24, 2026, 6:11pm UTC](https://discourse.julialang.org/t/solving-the-helium-eigenvalue-problem-where-to-start/137783/3 "2026-06-24T18:11:06Z")

</div>

Thank you for the suggestions! With the current code, 6 partial waves is enough for convergence to the ground state, so that tracks with what you’re saying.

Also, my fault for not proofreading/checking my sources thoroughly enough—I meant to say _global_ spectral methods were considered too computationally demanding relative to the level of resolution needed, hence the use of finite elements (in this case, spanned by Legendre polynomials, so the grid points correspond to the roots of the polynomials for easy quadrature). I do think part of the difficulty is that I never had a proper class on linear algebra or PDE solving, so I’m very much picking things up as I go.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 24, 2026, 6:18pm UTC](https://discourse.julialang.org/t/solving-the-helium-eigenvalue-problem-where-to-start/137783/4 "2026-06-24T18:18:17Z")

</div>

> [@ducksoverip](#):
>
> _global_ spectral methods were considered too computationally demanding relative to the level of resolution needed, hence the use of finite elements (in this case, spanned by Legendre polynomials, so the grid points correspond to the roots of the polynomials for easy quadrature).

That sounds like a pseudo-spectral collocation method? i.e. you use a spectral basis for the unknowns, and you enforce the equations at a grid of sample points, carefully chosen (usually related roots of the spectral basis or its derivatives).

---

<div class="post-metadata">

**Author:** ![ducksoverip](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ducksoverip/32/31967_2.png) [@ducksoverip](https://discourse.julialang.org/u/ducksoverip)\
**Post date:** [June 24, 2026, 6:39pm UTC](https://discourse.julialang.org/t/solving-the-helium-eigenvalue-problem-where-to-start/137783/5 "2026-06-24T18:39:30Z")

</div>

That sounds right, though I haven’t encountered that exact phrasing before.

---

<div class="post-metadata">

**Author:** ![lazarusA](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lazarusa/32/6571_2.png) [@lazarusA](https://discourse.julialang.org/u/lazarusA)\
**Post date:** [June 24, 2026, 6:48pm UTC](https://discourse.julialang.org/t/solving-the-helium-eigenvalue-problem-where-to-start/137783/6 "2026-06-24T18:48:14Z")

</div>

> [@ducksoverip](#):
>
> [this paper](http://dx.doi.org/10.1103/PhysRevE.81.056705)

Some years ago I did something here:

> **[GitHub - lazarusA/LightMatterInteraction: Time evolution in quantum mechanics, a lot of...](https://github.com/lazarusA/LightMatterInteraction)**
>
> Time evolution in quantum mechanics, a lot of helpful function to perform light-matter interaction.

which is the Julia version of the fortran code used in

> **[Perspectives for analyzing non-linear photo-ionization spectra with deep...](https://pubs.rsc.org/en/content/articlehtml/2020/fd/d0fd00117a)**

maybe there is something useful there, if not, apologies for the noise.

Edit: There is not a a lot of documentation or even a proper README, but it works, maybe I should bring this back to life :D.

---

<div class="post-metadata">

**Author:** ![ducksoverip](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ducksoverip/32/31967_2.png) [@ducksoverip](https://discourse.julialang.org/u/ducksoverip)\
**Post date:** [June 27, 2026, 9:44pm UTC](https://discourse.julialang.org/t/solving-the-helium-eigenvalue-problem-where-to-start/137783/7 "2026-06-27T21:44:22Z")

</div>

So, thinking about the problem a bit more and doing some more research, I’ve determined the following:

- Yes, the method used in the [paper](https://www.academia.edu/download/95215177/fulltext.pdf) is a pseudospectral collocation method (per [Wikipedia](https://en.wikipedia.org/wiki/Pseudo-spectral_method), apparently pseudospectral is a synonym for discrete variable representation, which is indeed what the Fortran code uses. Seems I’m one of today’s [lucky ten thousand](https://xkcd.com/1053/).)
- I inadvertently left a detail out of my description; after expanding \psi and 1/|\vec{r\_1} - \vec{r\_2}| in terms of coupled spherical harmonics, the entire equation is integrated over both pairs of angular coordinates, which (because of the orthogonality properties of the coupled harmonics) yields the following, where V\_{F0}(i,j) is the expanded and integrated form of 1/|\vec{r\_1} - \vec{r\_2}|, and i and j index bundles of the integer partial wave parameters L, l\_1, and l\_2:

-\frac{1}{2}\frac{\partial^2}{\partial r\_1^2} -\frac{1}{2}\frac{\partial^2}{\partial r\_2^2} + \frac{l\_1(l\_1 + 1)^2}{2r\_1^2} + \frac{l\_2(l\_2 + 1)^2}{2r\_2^2} - \frac{2}{r\_1} - \frac{2}{r\_2} \\+ \sum\_j V\_{F0}(i,j) = i \frac{\partial}{\partial t} \psi\_i(r\_1, r\_2)

- The above coupled PDEs are what are properly discretized on a 2D grid of Legendre polynomial roots and solved via reverse imaginary time propagation.
- If I’m purely interested in the ground-state eigenvalues and eigenvectors of this problem, then if I can construct the operators as a matrix, I can use KrylovKit.jl or ArnoldiMethod.jl to solve the eigenvalue problem for a given set of partial waves, making sure that I use enough partial waves to converge on the correct ground-state energy.

So, with this in mind, I think I can ask some more specific questions, namely:

1. Is pseudospectral collocation (what the paper calls finite element-discrete variable representation) the only viable approach, or is it worth trying a different method?
2. Based on the answer to the above, what Julia packages will help me specify the grid and construct the operators for this problem? I’m not really sure what the operators will look like if I’m not doing the split-operator + reverse imaginary time method.

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [June 28, 2026, 4:31pm UTC](https://discourse.julialang.org/t/solving-the-helium-eigenvalue-problem-where-to-start/137783/8 "2026-06-28T16:31:42Z")

</div>

> [@ducksoverip](#):
>
> 1. Based on the answer to the above, what Julia packages will help me specify the grid and construct the operators for this problem? I’m not really sure what the operators will look like if I’m not doing the split-operator + reverse imaginary time method.

Took me a minute to think through but I think the split-operator method just drops out fully.

### split-operator method

Split-operator method approximates the time-evolution operator if the Hamiltonian can be separated into 2 non-commuting parts H=H\_A + H\_B that are simple on their own:

U(t) = exp(-i H t) = exp(-iH \Delta t)^N\approx [exp(-iH\_A \Delta t) exp(-iH\_B \Delta t)]^N

In code this looks something like:

```julia-auto
function U(t, psi, N)
    diagonal_H_A = ...
    diagonal_H_B = ...

    diagonal_U_A = @. exp(-im*diagonal_H_A*t/N)
    diagonal_U_B = @. exp(-im*diagonal_H_B*t/N)

    for _ in 1:N
        psi .*= diagonal_U_A
        psi = fft(psi)
        psi .*= diagonal_U_B
        psi = ifft(psi)
    end
return psi
end

```

### your case

However for KrylovKit you just require the action of H itself and you can actually compute this in a similar way but exact:

H |\psi\rangle = H\_A |\psi\rangle + H\_B|\psi\rangle

You compute each part of the sum the same way that you would compute it in the split-operator method. So say |\psi\rangle is given in real space where H\_A is diagonal and H\_B is diagonal in momentum space, then you’d do roughly:

```julia-auto
function H(psi)
    diagonal_H_A = ...
    diagonal_H_B = ...

    psi_A = psi .* diagonal_H_A
    psi_momentum_space = fft(psi)
    psi_momentum_space .*= diagonal_H_B
    psi_B = ifft(psi_momentum_space)
    return psi_A + psi_B
end

```

Then simply use `KrylovKit.eigsolve(H, psi_0)` done. All you need to do really is construct the diagonals of H\_A and H\_B. I am not sure whether there are packages that can help you out-of-the-box. Perhaps WignerSymbols.jl is useful if you need Clebsch-Gordon coefficients.

---

<div class="post-metadata">

**Author:** ![ducksoverip](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ducksoverip/32/31967_2.png) [@ducksoverip](https://discourse.julialang.org/u/ducksoverip)\
**Post date:** [June 29, 2026, 7:18pm UTC](https://discourse.julialang.org/t/solving-the-helium-eigenvalue-problem-where-to-start/137783/9 "2026-06-29T19:18:53Z")

</div>

Thanks for the insight. The only tricky bit would be diagonalizing the electron-electron repulsion (1/|\vec{r\_1} - \vec{r\_2}|) but in principle there are routines to do that so long as I can construct the matrix. I’ll do some more looking into the matrix-construction end.