Variational computations of Helium

- General procedure

Quantum Mechanics
 K. Konishi, G. Paffuti

Study analytically  the variational problem for s-states of Helium with a low dimensional basis.  Use both Hylleraas and Pekerris coordinates.

ClearAll["Global`*"];
Off[General::spell, General::spell1,InterpolatingFunction::"dmval"];

Formulation of the problem

General features

Variational principle states that the computation of lowest eigenvalue of the Schrödinger equation can be formulated as the search of the minimum of a quadratic form on a unit sphere of Hilbert space, i.e. on vectors of unit norm. The principle is of great numerical relevance when combined with Rietz method: a trial wave function is defined as linear combination of a finite set of base functions, possibly depending on some parameter λ. For each λ the minimum problem is equivalent to a finite matrix formulation of Schrödinger equation and an eigenvalue E(λ) is found. The subsequent minimization with respect to λ gives the minimum in this approximate Hilbert space. We make two observations which will be relevant in the numerical work

If we choose N functions for the basis as a subset of a complete basis (when N → ∞),  we are guaranteed to obtain the exact solution in the limit, and the differnce in norm between our approximation and the exact solution wil tend to zero in this limit.

For large enough bases the choice of parameters λ should be irrelevant, in principle, but in practice it can make a lot of difference on the convergence rate.

Whitin the stated approximation, the variational approach will provide an estimate also for the excited states: the higher eigenvalues of the finite matrix. In principles the optimization of λ parameters must be done on each state separately.

In Quantum Mechanics the extra parameter λ is usually a scale, i.e. has dimension of a length, or 1/length depending on the formalism. Quite intuitively, each problem has a best length scale of its own.

Helium

In this case the variational problem lies in finding the minimum of

ℰ =(( M - L)/N) ;

where

M =(1/2) ∫(((((Subscript[∇, 1])Ψ))^2)+ ((((Subscript[∇, 2])Ψ))^2)) ;             L = ∫((Z/(Subscript[r, 1]))+(Z/(Subscript[r, 2])) -(1/(Subscript[r, 12])) )(Ψ^2) ;                    N = ∫(Ψ^2).

Let us suppose we choose a wave function depending on a scale λ, ψ(λ x). By the change of variables x → x/λ this dependence cancels in the wave function and reappears in the kinetc term and in the potential term, which are homogeneous under scale transformation.  The quadratic form ℰ takes the form

ℰ =(((λ^2)  M - λ L)/N) ;

while quadratic forms M, L, N are always given by previuous formulae, computed with λ=1. At the final stage ℰ will be minimized with respect to λ. We add that clearly ℰ is the mean energy on the state ψ.

Let us consider a basis of functions (Subscript[f, i]). Our trial function will have the form

Ψ = (∑_(i=1)^N)(Subscript[c, i])(Subscript[f, i])

Quadratic forms become quadratic forms of coefficients  (we adopt Einstein's convention on sums from now on, all sums are from 1 to N; all (Subscript[f, i])are taken real for simplicity):

M = (Subscript[M, ij])(Subscript[c, i])(Subscript[c, j]);         L = (Subscript[L, ij])(Subscript[c, i])(Subscript[c, j]);          N = (Subscript[N, ij])(Subscript[c, i])(Subscript[c, j]);

(Subscript[M, ij])=∫(Subscript[f, i]) M (Subscript[f, j]) =(1/2) ∫{((Subscript[∇, 1])(Subscript[f, i]))((Subscript[∇, 1])(Subscript[f, j]))+((Subscript[∇, 2])(Subscript[f, i]))((Subscript[∇, 2])(Subscript[f, j]))};
(Subscript[L, ij])= ∫(Subscript[f, i])L (Subscript[f, j])= ∫(Subscript[f, i])((Z/(Subscript[r, 1]))+(Z/(Subscript[r, 2])) -(1/(Subscript[r, 12])) )(Subscript[f, j]) ;   (Subscript[N, ij])= ∫(Subscript[f, i])N (Subscript[f, j])= ∫(Subscript[f, i])(Subscript[f, j]) ;

A trivial computation shows that

((δℰ)/(Subscript[δc, i]))= 0  ⇒  ((λ^2)(Subscript[M, ij])-λ (Subscript[L, ij]))(Subscript[c, j])= ℰ[λ] (Subscript[N, ij])(Subscript[c, j])

i.e. the extremum appears as an eigenvalue for the truncated Schrödinger equation (as expected). The only unusual thing is the matrix N in the right hand side, due to the possible non orthogonality of base functions. We have only to minimize the eigenvalue found in solving (). Once we find a solution if we need the wave function we used (), after performing the inverse of the scale transformation used before, i.e. we have to write:

Ψ = (∑_(i=1)^N)(Subscript[c, i])(Subscript[f, i])[λ x]

Hylleraas coordinates

The first variational computation of Helium spectrum has been made by Hylleraas [1]. We limit to S wave functions. As discussed in notebook [15] the three particles nucleus + 2 electrons determine an invariant triangle, fixed by its edges (Subscript[r, 1]), (Subscript[r, 2]), (Subscript[r, 12]),  respectively, the distances between the two electrons and the nucleus and their relative distance.

The wave function for s - states depends only on these variables. Hylleraas  variables are a linear combination of these variables

s = (Subscript[r, 1]) + (Subscript[r, 2]) ; t = (Subscript[r, 1])-(Subscript[r, 2]) ; u = (Subscript[r, 12]);               0≤ t ≤ u ≤ s ≤ ∞

The bounds come from the  known inequalities on lengths of the edges of a triangle. In these variables a straightforward computation gives, apart for an overall factor which cancels in the ratio ℰ :

M =(∫_0^∞)ds (∫_0^s)du(∫_0^u)dt  (   s ((u^2)-(t^2)) (((∂Ψ)/(∂s)) ((∂Ψ)/(∂u))+((∂Ψ)/(∂s)) ((∂Ψ)/(∂u)))+
t ((s^2)-(u^2)) (((∂Ψ)/(∂t)) ((∂Ψ)/(∂u))+((∂Ψ)/(∂t)) ((∂Ψ)/(∂u))  )+ ((s^2)-(t^2)) u (((∂Ψ)/(∂s)) ((∂Ψ)/(∂s))+((∂Ψ)/(∂t)) ((∂Ψ)/(∂t))+((∂Ψ)/(∂u)) ((∂Ψ)/(∂u)))  );
L = (∫_0^∞)ds (∫_0^s)du(∫_0^u)dt ( 4 Z u s - ((s^2)-(t^2))) (Ψ^2) ; 
N = (∫_0^∞)ds (∫_0^s)du(∫_0^u)dt u ((s^2)-(t^2)) (Ψ^2) ;

The reader can check at least some of these terms, or follow the notebook [15] for their derivation within Mathematica. The base functions will be ψ[i][s,t,u] (we suppose to have reabsorbed the scale factor λ)

Ψ = ∑ C[i] ψ[i][s,t,u]

and the matrices associated at the quadratic forms are

M[i,j]=(∫_0^∞)ds (∫_0^s)du(∫_0^u)dt (s ((u^2)-(t^2)) (((∂ψ(i))/(∂u)) ((∂ψ(j))/(∂s))+((∂ψ(i))/(∂s)) ((∂ψ(j))/(∂u)))+t ((s^2)-(u^2)) (((∂ψ(i))/(∂u)) ((∂ψ(j))/(∂t))+((∂ψ(i))/(∂t)) ((∂ψ(j))/(∂u)))+((s^2)-(t^2)) u (((∂ψ(i))/(∂s)) ((∂ψ(j))/(∂s))+((∂ψ(i))/(∂t)) ((∂ψ(j))/(∂t))+((∂ψ(i))/(∂u)) ((∂ψ(j))/(∂u)))) ;
L[i,j] =(∫_0^∞)ds (∫_0^s)du(∫_0^u)dt ( 4 Z u s - ((s^2)-(t^2))) ψ(i) ψ(j) ;
N[i,j] =(∫_0^∞)ds (∫_0^s)du(∫_0^u)dt u ((s^2)-(t^2))  ψ(i) ψ(j) ;

For completeness, and for comparison with other notebooks, we add that the first correction to the nonrelativistic Hamiltonian used so far is due to the finite mass of the nucleus. The corrected Hamiltonian at first order m/M takes the form (in atomic units, with m = ℏ = e = 1)

H = (1/(2μ))(((Subscript[p, 1])^2)+ ((Subscript[p, 2])^2)) + V + (1/M)(Subscript[p, 1]).(Subscript[p, 2])

where μ is the reduced mass of the electron,  μ = m M/(m+M), and the last term, absent in the two body problem, is named mass polarization correction. These two terms give rise to a correction in eigenvalues:

(Subscript[δE, mass]) = ((μ/m)-1) E;           (Subscript[δ, Pol]) = (m/M) ((⟨ψ|K |ψ⟩)/(⟨ψ|ψ⟩)).

The first term is trivial and comes from the fact that eigenvalues are proportional to m and we have to make the substitution m → μ. The second term is the matrix element of (Subscript[p, 1]).(Subscript[p, 2]). Also this operator have to be translated in Hylleraas variables, the result is

K=(∫_0^∞)ds (∫_0^s)du(∫_0^u)dt 
 (-2s ((u^2)-(t^2)) (((∂Ψ)/(∂u)) ((∂Ψ)/(∂s)))+u ((s^2)+(t^2)-2 (u^2)) ((((((∂Ψ)/(∂s))))^2) -(((((∂Ψ)/(∂t))))^2) )-2t ((s^2)-(u^2)) (((∂Ψ)/(∂u)) ((∂Ψ)/(∂t)))-u ((s^2)-(t^2)) (((((∂Ψ)/(∂u))))^2)   ) ;

The relative matrix is (see notebook [15] ):

K[i,j]=(∫_0^∞)ds (∫_0^s)du(∫_0^u)dt  (-s ((u^2)-(t^2)) (((∂ψ(i))/(∂u)) ((∂ψ(j))/(∂s))+((∂ψ(i))/(∂s)) ((∂ψ(j))/(∂u)))+
u ((s^2)+(t^2)-2 (u^2)) (((∂ψ(i))/(∂s)) ((∂ψ(j))/(∂s))-((∂ψ(i))/(∂t)) ((∂ψ(j))/(∂t)))-t ((s^2)-(u^2)) (((∂ψ(i))/(∂u)) ((∂ψ(j))/(∂t))+((∂ψ(i))/(∂t)) ((∂ψ(j))/(∂u)))-u ((s^2)-(t^2)) ((∂ψ(i))/(∂u)) ((∂ψ(j))/(∂u))  );

Hylleraas basis

Hylleraas choice for the basis is

(Subscript[φ, ijk])(s,t,u)= (s^i)(t^(2j+p))(u^k) Exp[-s/2]

The choice p = 0, 1 gives rise to symmetric or antisymetric functions, alowing the description of paraHelium and orthoHelium. The great advantage of this basis is that all matrix elements are trivially computed, as the most general integral appearing in the computation has the form

(∫_0^∞)ds (∫_0^s)du(∫_0^u)dt  (s^α) (t^β) (u^γ) Exp[-s] = (((2+α+β+γ)!)/((1+β)(2+β+γ)));

Computing integrals

We wil hav to do a lot of integrals, each one elementary. The standard routine Integrate can be used, both for generic arguments and for specific arguments:

risExact[a_,b_,c_]:= (((2+a+b+c)!)/((1+b)(2+b+c)))

Integrate[s^a t^b u^c Exp[-s],{s,0,∞},{u,0,s},{t,0,u},Assumptions-> {a≥ 0,b≥ 0,c≥ 0}]//Timing

Integrate[s^2 t^2 u^1 Exp[-s],{s,0,∞},{u,0,s},{t,0,u}]//Timing

It is faster to use indefinite integrals for t and u and perform the last integral in s, which gives just a factorial, by substitution:

tInt[x_]:= ((#1/.{t-> u})- (#1//.{t-> 0}))& [Together[Integrate[x,t]]];
uInt[x_]:= ((#1/.{u-> s})- (#1/.{u-> 0}))& [Integrate[x,u]];

s integral is the definition of factorials:

(∫_0^∞)ds   (s^k)  Exp[-s] = k! ;

sInt[x_]:= Together[Expand[x]//.{s^a_. Exp[-s]-> a!, Exp[-s]-> 1}];

Then we can define:

myInt[x_]:= sInt[uInt[tInt[x]]];

The computation is faster of at least a factor of 30:

myInt[s^2 t^2 u^1 Exp[-s]]//Timing

This function will be improved below.

First computation

As a warm up we take a simple example with only three basis functions

ψ[1] = Exp[-s/2];
ψ[2] = u ψ[1];
ψ[3] = (t^2)   ψ[1];

dim = 3;
Z=2;

Let us define our matrices as written above and compute them using the function myInt[] (we compute more than we need, we are not using symmetry of matrices);

Mij[i_,j_]:= 
u (s^2-t^2) (D[ψ[i],s] D[ψ[j],s]+D[ψ[i],t] D[ψ[j],t]+D[ψ[i],u] D[ψ[j],u])+s (u^2-t^2) (D[ψ[i],s] D[ψ[j],u]+D[ψ[i],u] D[ψ[j],s])+t (s^2-u^2) (D[ψ[i],t] D[ψ[j],u]+D[ψ[i],u] D[ψ[j],t])

matM=Table[myInt[Mij[i,j]],{i,dim},{j,dim}];

Lij[i_,j_]:= ( 4 Z u s - (s^2-t^2))ψ[i] ψ[j];
matL=Table[myInt[Lij[i,j]],{i,dim},{j,dim}];

Nij[i_,j_]:=u (s^2-t^2)ψ[i] ψ[j];
matN=Table[myInt[Nij[i,j]],{i,dim},{j,dim}];

Kij[i_,j_]:= 
u(s^2+t^2-2u^2)(D[ψ[i],s] D[ψ[j],s]-D[ψ[i],t] D[ψ[j],t]) -
u(s^2-t^2)D[ψ[i],u] D[ψ[j],u] - 
s(u^2-t^2)(D[ψ[i],s] D[ψ[j],u]+D[ψ[i],u] D[ψ[j],s]) -
t(s^2-u^2)(D[ψ[i],t] D[ψ[j],u]+D[ψ[i],u] D[ψ[j],t]);
matK=Table[myInt[Kij[i,j]],{i,dim},{j,dim}];

Here is the result :

MatrixForm/@{matM,matL,matN}

Eigenvalues and eigenvectors 1

We will analyse several possibilities to find the eigenvalues. The first and most simple is just to decide which state we want to study and use the Mathematica function Eigenvalues

calcEne[λ_Real,state_Integer]:=Sort[N[Eigenvalues[{λ^2 matM - λ matL,matN}]]][[state]]

The folowing graph shows the dependence on λ of the eigenvalue

Plot[calcEne[x,1],{x,2,5}]

To find the minimum the simplest option is to use Mathematica function FindMinimum

minEigenvalue = FindMinimum[calcEne[x,1],{x,3,5}]

{ene,λ} = {minEigenvalue[[1]],x/.minEigenvalue[[2]]}

We will find below a quite accurate value for the ground state energy

eneAccurate= -2.903724377034119859;

We see that even with such a small basis we have a correct eigenvalue,  up to four digits.

Having found the optimal λ value we can now compute eigenvectors

{eigH,eigV} = Eigensystem[{λ^2 matM - λ matL,matN}];
ord = Ordering[eigH];
{eigH,eigV} = {eigH[[ord]],eigV[[ord]]};

eigV[[1]]

Our set of basis functions is

base = Table[ψ[i],{i,1,dim}];

Then, after performing the rescaling :

Simplify[(eigV[[1]].base/.{s-> k s,t-> k t,u-> k u})/eigV[[1,1]]]

We have to put k = λ, we leave in symbolic form to clarify its role.

Eigenvalues and eigenvectors 2

A second possibility to compute the eigenvalues is to compute zeros of the matrix H-E, where H is the effective Hamiltonian and E here stands for ℰ N. it would be prefereable to compute numerically the determinant, here we use NSolve to extract the roots:

calcEneDet[λ_Real,state_Integer]:=Block[{x},Sort[x/.NSolve[Det[λ^2 matM - λ matL - x matN],x]][[state]]];

minEig = FindMinimum[calcEneDet[x,1],{x,3,5}]

which is identical to previous result.

Eigenvalues and eigenvectors 3: arbitrary precision

Even if in this case there is no physical motivation, we want to show how one can use Mathematica ability to obtain results with very high precision. We have computed exact matrices, the elements are rational numbers not just machine numbers. If we want to work with high precision we have to use input parameters with the wanted precision. As an example our matrix matM is exact, but let us compare the following results:

Sort[N[Eigenvalues[matM]]][[1]]//InputForm

Sort[N[Eigenvalues[matM],40]][[1]]//InputForm

Sort[N[Eigenvalues[1.0matM],20]][[1]]//InputForm

Even if we require a 20 digit precision the third result has the same precision as the first because the parameter 1.0 has the precision of Machine Numbers, in Mathematica:

$MachinePrecision

If we want a precise output we have to fix the precision of the parameter, with the function SetPrecision

SetPrecision[1.0,40]

Unfortunately we cannot just put the parameter λ with a high precision in our functions because the generalized form of Eigenvalues accepts only machine  numbers:

λ1 = SetPrecision[2.0,20];

temp=N[Eigenvalues[{λ1^2 matM - λ1 matL,matN}],20];

We can in this simple case invert the matrix N and transform the problem in a usual eigenvalue problem:

matNinv = N[Inverse[matN],20];

matMHP =  matNinv.matM; matLHP = matNinv.matL;

We define an high precision version of our matrices

calcEneHP[λ_Real,state_Integer,p_]:= (λHP = SetPrecision[λ,p];
Sort[N[Eigenvalues[λHP^2 matMHP - λHP matLHP],p]][[state]])

and comput the eigenvalues

minEig = FindMinimum[calcEneHP[x,1,20],{x,3,5},WorkingPrecision->20,PrecisionGoal->12]//Timing

This difficulty is avoided if we use the determinant to find the eigenvalues

calcEneDetHP[λ_Real,state_Integer,p_]:= calcEneDet[SetPrecision[λ,p],state];

minEig = FindMinimum[calcEneDetHP[x,1,20],{x,3,5},WorkingPrecision->20,PrecisionGoal->12]//Timing

Second computation

In this section we will extend the basis of functions and collect some results to be used for check computations performed in other notebooks.

We will classify the basis by the maximum power in each variables, s, t ,u. For t we use the index j in (t^(2j)). We introduce also a number, parity, with possible values 0,1 to allow generation of basis functions both for even and odd states under the electron exchange.

Z = 2; parity = 0;
maxs = 1; maxt =2; maxu = 2;

It is convenient to collect these definitions in a function wich returns the list of basis functions and the dimension of the base

getBase[maxs_,maxt_,maxu_]:=(
ind=0;
Do[Do[Do[
ind++;
  ψ[ind]= s^i1* t^(2 j1+parity)* u^k1 * Exp[-s/2]/(i1+2*j1+k1)!;
,
{i1,0,maxs}],{j1,0,maxt}],{k1,0,maxu}];
dim = ind;
base = Table[ψ[k],{k,1,dim}];
{base,dim})

We have introduced a numerically convenient factor in the definition of base vectors, just to avoid numbers which are too big.

{base,dim} = getBase[maxs,maxt,maxu];

We repeat here for the reader the definition of the matrices and write a function which compute them. As the execution time grows we use the

symmetry properties to compute only about half of the matrices, other elements are obtained by (Subscript[M, ji])= (Subscript[M, ij]).

Mij[i_,j_]:= u (s^2-t^2) (D[ψ[i],s] D[ψ[j],s]+D[ψ[i],t] D[ψ[j],t]+D[ψ[i],u] D[ψ[j],u])+s (u^2-t^2) (D[ψ[i],s] D[ψ[j],u]+D[ψ[i],u] D[ψ[j],s])+t (s^2-u^2) (D[ψ[i],t] D[ψ[j],u]+D[ψ[i],u] D[ψ[j],t]);

Lij[i_,j_]:= ( 4 Z u s - (s^2-t^2))ψ[i] ψ[j];

Nij[i_,j_]:=u (s^2-t^2)ψ[i] ψ[j];

Kij[i_,j_]:= u(s^2+t^2-2u^2)(D[ψ[i],s] D[ψ[j],s]-D[ψ[i],t] D[ψ[j],t]) -
u(s^2-t^2)D[ψ[i],u] D[ψ[j],u] - 
s(u^2-t^2)(D[ψ[i],s] D[ψ[j],u]+D[ψ[i],u] D[ψ[j],s]) -
t(s^2-u^2)(D[ψ[i],t] D[ψ[j],u]+D[ψ[i],u] D[ψ[j],t]);

getMatrix[X_]:= Module[{mat},
mat = 0*IdentityMatrix[dim];
Table[mat[[i,j]]=myInt[X[i,j]];,{i,1,dim},{j,i,dim}];
mat + Transpose[mat] - DiagonalMatrix[Tr[mat,List]]];

matM1 = getMatrix[Mij];//Timing

matL1 = getMatrix[Lij];
matN1 = getMatrix[Nij];

Improvement of the function myInt

In the Hylleraas basis we can read with the function Exponent the power of various variables in ψ and apply directly the formula

(∫_0^∞)ds (∫_0^s)du(∫_0^u)dt  (s^α) (t^β) (u^γ) Exp[-s] = (((2+α+β+γ)!)/((1+β)(2+β+γ)));

to compute integrals. We need exponents for each monomial in an expression, the simple thing is to transform the expression in a list, by applyng List. For example

Apply[List,(s^2)Exp[-s] + s t Exp[-s]]

The exception is the case of single negative monomial. In this case the Head of the expression is Times and the produced list contains a sputious -1:

FullForm[-(s^2) t Exp[-s]]

Apply[List,-(s^2) t Exp[-s]]

To avoid this problem we just add 1, without exponential, to the expression:

Apply[List,-(s^2) t Exp[-s]+1]

When we will use Exponent we just drop the first result, {0, 0, 0} corresponding to the added number 1. As an example

espr= - s^2 t Exp[-s] + 3 s t Exp[-s];

Drop[Exponent[#,{s,t,u}]&/@Apply[List,Expand[espr+1]],1]

and also

espr2= - s^2 t Exp[-s] ;
Drop[Exponent[#,{s,t,u}]&/@Apply[List,Expand[espr2+1]],1]

The integral now is easily defined as

myInt[x_]:=Module[{cf,pow},pow=((((2+#[[1]]+#[[2]]+#[[3]])!)/((1+#[[2]])(2+#[[2]]+#[[3]]))))&/@ Drop[Exponent[#,{s,t,u}]&/@Apply[List,Expand[x+1]],1];
cf = Drop[Apply[List,Expand[x+1]]/.{E^c_.-> 1,s->1,t-> 1,u-> 1},1];
cf.pow]

This procedure is much faster than the previous one :

matM = getMatrix[Mij];//Timing

matL = getMatrix[Lij];
matN = getMatrix[Nij];
matK = getMatrix[Kij];

and give the same results:

{matM== matM1, matL == matL1, matN== matN1}

Eigenvalues and eigenvectors

minEig = FindMinimum[calcEne[var,1],{var,3,5}]//Timing

or with high precision

matNinv =Inverse[ matN];//Timing

matMHP =  matNinv.matM; matLHP = matNinv.matL;

prec=20;

(minEig =N[ FindMinimum[calcEneHP[x,1,prec],{x,3,5},WorkingPrecision->4*prec,PrecisionGoal->prec],prec])//Timing

The previous routine with determinants is too slow :

TimeConstrained[FindMinimum[calcEneDetHP[x,1,20],{x,3,5},WorkingPrecision->20,PrecisionGoal->12],10]

The reason is that in the instruction

calcEneDet[λ_Real, state_Integer] := Block[{x}, Sort[x /. NSolve[Det[λ^2 matM - λ matL - x matN], x]][[state]]];

the determinant was computed symbolically. But let us suppose that we are only interested in one eigenvalue, as is usually the case, and we have a rough estimate of it. We can compute numerically the determinant with

calcDet[x_?NumberQ,ham_,norma_,prec_]:=
Module[{det},
 det=Det[ham - SetPrecision[x,prec]* norma];
precdet = precdetnew;
precdetnew =Precision[det];
det]

then compute the eigenvalue with FindRoot and apply a minimization procedure.

calcEneDet[λ_?NumberQ,x0_,prec_]:=Block[{x },
x/.N[FindRoot[calcDet[x,λ^2 matM - λ matL,matN,prec],{x,x0,x0+1/10},
WorkingPrecision-> 2.5prec, AccuracyGoal-> prec],prec]]

Consider an example. A rough estimate for energy is -3.0. We use 30 digit precision

ene0 = -3.0;
prec=30;

(A comparison of the=N[FindMinimum[calcEneDet[var,SetPrecision[ene0,prec],prec],{var,3,5},WorkingPrecision->2.5*prec,PrecisionGoal->prec],prec])//Timing

This method is faster than the previous one, and we have worked with an even higher precision.

The global variable precdet tell us with wich precision has been effectively computed the determinant, than give us also a control on the results:

precdet

A comparison of the two results is :

{minEigDet[[1]]-minEig[[1]], (var/.minEigDet[[2]]) - (x/.minEig[[2]]) }

We can be quite sure that the first 19 digits for the eigenvalue are correct, within this approximation (the correctnes with respect the true value is another problem).

Third computation: some checks and some results

Some checks

Our first task is to check our results against other notebooks. For efficiency in intensive computational routines the previous matrices have been introduced manually, i.e. writing explicitely the matrix elements. It is better to check for typos. Here we report the rows and columns from 10 to 15 for the previous matrices computed with maxs = 1; maxt =2; maxu = 2; in notebook [17]

checkmatM ={{23200,39168,85664,(5859/2),9828,21747},{39168,84480,197120,(17712/5),(64152/5),(189266/5)},{85664,197120,495936,(32868/5),(132132/5),(1237379/15)},{(5859/2),(17712/5),(32868/5),1248,3328,2768},{9828,(64152/5),(132132/5),3328,(31552/3),(27680/3)},{21747,(189266/5),(1237379/15),2768,(27680/3),(62480/3)}};

checkmatL ={{138930,213772,(1496404/3),23280,85360,125376},{213772,397208,993020,(125376/5),(501504/5),196320},{(1496404/3),993020,(7944160/3),(250752/5),(1086592/5),458080},{23280,(125376/5),(250752/5),10260,30780,20762},{85360,(501504/5),(1086592/5),30780,102600,(228382/3)},{125376,196320,458080,20762,(228382/3),114228}};

checkmatN ={{99840,134400,336000,18018,72072,87516},{134400,225792,602112,(80784/5),(350064/5),120120},{336000,602112,1705984,(175032/5),(816816/5),300300},{18018,(80784/5),(175032/5),9600,32000,15680},{72072,(350064/5),(816816/5),32000,(352000/3),62720},{87516,120120,300300,15680,62720,77760}};

checkmatK ={{(2080/3),576,1168,-(189/2),-168,891},{576,1056,2464,-216,-(3168/5),(2794/5)},{1168,2464,(17696/3),-(2376/5),-(8316/5),(11011/15)},{-(189/2),-216,-(2376/5),192,512,32},{-168,-(3168/5),-(8316/5),512,(3968/3),(320/3)},{891,(2794/5),(11011/15),32,(320/3),(3040/3)}};

Let us compare with the semi - analytical results previously produced :

{matM[[Range[10,15],Range[10,15]]]== checkmatM,
matL[[Range[10,15],Range[10,15]]]== checkmatL,
matN[[Range[10,15],Range[10,15]]]== checkmatN,
matK[[Range[10,15],Range[10,15]]]== checkmatK}

This gives some confidence on the matrices computed in other notebooks.

Large basis

Now let us see why the method exposed up to this point is not the best choice if one has to do computations with larger basis and to repeat some of them many times. In the following example dim = 120, i.e. our matrices are 120 × 120.

Z = 2; parity = 0;
maxs = 3; maxt =4; maxu = 5;

{base,dim} = getBase[maxs,maxt,maxu];

matM1 = getMatrix[Mij];//Timing

The time needed to perform the previous operation clearly depends on the computer, but in any case it is quite long for repeated computations.

It is better to generate an analytical expression for the matrix elements. A variant of this expression is the one used in notebook [17]. The variable shT will denote parity.

Remove[i,j,i1,j1,k1,i2,j2,k2,shT];

φ[i]=  s^i1 t^(2j1+shT) u^k1 Exp[-s/2]/(i1+2j1+k1)!;
φ[j]=  s^i2 t^(2j2+shT) u^k2 Exp[-s/2]/(i2+2j2+k2)!;

esp=Simplify[(u (s^2-t^2) (D[φ[i],s] D[φ[j],s]+D[φ[i],t] D[φ[j],t]+D[φ[i],u] D[φ[j],u])+s (u^2-t^2) (D[φ[i],s] D[φ[j],u]+D[φ[i],u] D[φ[j],s])+t (s^2-u^2) (D[φ[i],t] D[φ[j],u]+D[φ[i],u] D[φ[j],t]))];

If the reader wants to see the expressions below he can take away the semi - colons

expressionM=Simplify[myInt[esp]];

esp=Simplify[( 4 valoreZ* u s - (s^2-t^2))φ[i] φ[j]];
expressionL=Simplify[myInt[esp]];

esp=Simplify[u (s^2-t^2)φ[i] φ[j]];
expressionN=Simplify[myInt[esp]];

esp=Simplify[u(s^2+t^2-2u^2)(D[φ[i],s] D[φ[j],s]-D[φ[i],t] D[φ[j],t]) -
u(s^2-t^2)D[φ[i],u] D[φ[j],u] - 
s(u^2-t^2)(D[φ[i],s] D[φ[j],u]+D[φ[i],u] D[φ[j],s]) -
t(s^2-u^2)(D[φ[i],t] D[φ[j],u]+D[φ[i],u] D[φ[j],t])];
expressionK=Simplify[myInt[esp]];

Finally we write a function which computes the matrices using these analytical expressions :

computeMatrices[Zval_,parity_,maxs_,maxt_,maxu_]:=
Module[{matM,matL,matN,matK,espM,espN,espL,espK},maxinds = (maxs+1); maxindt = (maxt+1);maxindu = (maxu+1);
dim = maxinds*maxindt*maxindu;
matM = matL = matN = matK =  0*IdentityMatrix[dim];
espM = expressionM/.{valoreZ->Zval,shT-> parity};
espL = expressionL/.{valoreZ->Zval,shT-> parity};
espN = expressionN/.{valoreZ->Zval,shT-> parity};
espK = expressionK/.{valoreZ->Zval,shT-> parity};
	

Do[  Do[ Do[Do[Do[Do[
i=i1+maxinds*j1+maxinds*maxindt*k1+1;
j=i2+maxinds*j2+maxinds*maxindt*k2+1;
 If[i ≥ j,
matM[[i,j]]= espM;
matL[[i,j]]= espL;
matN[[i,j]]= espN;
matK[[i,j]]= espK;
] ;
,
{k2,0,maxu}],{j2,0,maxt}],{i2,0,maxs}],
{k1,0,maxu}],{j1,0,maxt}],{i1,0,maxs}];
matM =matM+Transpose[matM] - DiagonalMatrix[Tr[matM,List]];
matL =matL+Transpose[matL] - DiagonalMatrix[Tr[matL,List]];
matN =matN+Transpose[matN] - DiagonalMatrix[Tr[matN,List]];
matK =matK+Transpose[matK] - DiagonalMatrix[Tr[matK,List]];
{matM,matL,matN,matK}];

The computation is almost an order of magnitude faster

{matM,matL,matN,matK}=computeMatrices[2,0,maxs,maxt,maxu];//Timing

while the result is identical :

matM1== matM

This method is the one effectively used in the notebook [17]. The ground state energy is

(minEig = FindMinimum[calcEne[x,1],{x,3,5}])//InputForm//Timing

minEig[[1]]-eneAccurate

We give also an example for the computation of ground state of orthoHelium, we have just to change the parameter parity :

{matM,matL,matN,matK}=computeMatrices[2,1,maxs,maxt,maxu];//Timing

(minEigOrtho = FindMinimum[calcEne[x,1],{x,3,5}])//InputForm

t - dependent exponents

For higher excited states line 1s-ns one expect a quite different scale between the two electrons. In the Hylleraas basis this feature is not well reproduced as the exponential depends in a symmetric way on the two coordinates. A possible alternative is to have exponential which depends on t. In order to have even and odd states these exponent must appear in the form Cosh[t] and Sinh[t].

We give here an example with a small basis and we will compute the ground state of orthoHelium.

parity = 1;
signParity = If[parity== 0,1,-1];
dim = 6;  Z=2;

As a basis we choose

ψ[1]=Exp[-s/2] (Exp[c t/2]+signParity * Exp[-c t/2])/2;
ψ[2]=s ψ[1];
ψ[3]=t Exp[-s/2] (Exp[c t/2]- signParity * Exp[-c t/2])/2;
ψ[4]=u ψ[1];
ψ[5]=u ψ[2];
ψ[6]=u ψ[3];

We can define our integration procedure along the same lines folowed in the previous sections

tInt[x_]:= ((#1/.{t-> u})- (#1//.{t-> 0}))& [Together[Integrate[x,t]]];
uInt[x_]:= ((#1/.{u-> s})- (#1/.{u-> 0}))& [Integrate[x,u]];

After t integration the generic integral in s is of the form (s^n)Exp[α s] and we can use

Integrate[Exp[α s] s^n, {s, 0, Infinity},
          Assumptions -> Re[n] > -1 && Re[α] < 0]

to define

sInt[x_]:=Together[(Expand[x]//.{E^y_:>E^Factor[y]})//.{s^e_. Exp[s f_.]->e! (-f)^(-e-1),Exp[s f_.]->-1/f}]

myInt[x_]:= sInt[uInt[tInt[x]]];

We generate the matrices with the usual function getMatrix:

matM=getMatrix[Mij];//Timing

matL = getMatrix[Lij];

matN = getMatrix[Nij];

Up to this point the matrices depend on the unknown parameter c. To compute the eigenvalues we use a numerical value for c. An alternative procedure, which can

be testd by the reader is to compute also matrices for numerical values of  c,  inside the following function

calcEneExp[λ_Real,valc_Real,stato_Integer]:=Sort[N[Eigenvalues[{λ^2 (matM/.{c-> valc}) - λ (matL/.{c-> valc}),(matN/.{c-> valc})}]]][[stato]]

The search for the minimum is the standard one, but this time we have to minimize with respect two parameters, λ and c

minEigExp = FindMinimum[calcEneExp[x,y,1],{{x,2.5},{y,0.4}},PrecisionGoal-> 8]//InputForm

This result is compatible with the one found in the previous section in Hylleraas basis.

Triangular variables (Pekeris)

The Hylleraas variables are intuitive and have the remarkable feature of allowing an analytical computation of the matrix elements. They have a drawback: the domain of the three variables s, u, t is not independent and  it is impossible, or at least not simple, to construct simple orthogonal basis.  As a consequence,  the matrices in this basis are full matrices, not sparse, and this is very bad from the numerical point of view: both memory requirements and numerical rounding errors cause severe limitations. In fact  also in our elementary considerations we have seen that it is difficult to have a basis larger than 100 - 150 elements. While these dimensions suffice for a reasonable precise determination of the ground state energy, the method appears not well suited for excited states. A physical reason related to these numerical considerations lies in the too symmetric form of the exponential decrease of functions: Exp[-((Subscript[r, 1])+ (Subscript[r, 2]))]. It is expected that for exited states the characteristic lenghts for the two electrons are different and a better approximation would be a t-dependent exponential: it is possible to mimic an exponential with a polynomial, as implicitly done in the Hylleraas variables, but this requires a high degree and again we go into troubles with memory eand precision.

A simple way out for these problems has been pointed out by Pekeris, [2]. The first observation is that if one uses more symmetrical variables, the so called triangular variables

x = r1 + r2 -r12 (= s - u);      y = r1-r2 +r12 (= u - t);     z = -r1+r2+r12 (= u + t);

by inequalities on the lenghts of the edges of a triangle all these variables have the same, uncorrelated, range

0 ≤ x, y, z ≤ ∞ .

Being independent it is now possible to choose a really factorized basis, which, taking into account the possibilty of even and odd terms under exchage of the electrons has the general form

(Subscript[φ, ijk])(x,y,z)= (Subscript[F, i])[α x]((Subscript[F, j])[β y](Subscript[F, k])[γ z]± (Subscript[F, j])[β z](Subscript[F, k])[γ y]).

Under the exchange  of electrons (orbital) variables ( 1 ↔ 2 )

x → x ; y → z ; z → y

then ± signs in () refer to symmetric (paraHelium) and antisymmetric (orthoHelium) states. α, β, γ are three parameters which are at our disposal. One of them is just a scale, as λ parameter in the Hylleraas variables, let us consider β as this scale parameter, then for α = a β, γ = c β, β is reabsorbed by a scale transformation and effective dependence is on (a x) and (c z). A priori these two parameters can be also variational parameters but we make a simpler choice. Assume that the F functions have an exponential behavior, then in the base function we have

Exp[- a x - y - c z ] = Exp[ - {( a+1-c)(Subscript[r, 1]) + (a-1-c)(Subscript[r, 2]) + (-a + 1 +c) (Subscript[r, 12])]

If we want to avoid the θ dependent exponentials we have to choose  a = c+1. We make the simplest choice and put c=1, a=2. The only parameter left  is the overall scale, as in the Hylleraas case.

Now we come to the essential point: being x,y,z three independent variables we can now choose as base functions some component of a complete set in the interval [0,∞]. The natural choice is to consider the Lagerre polynomials, i.e. basis of the form

(Subscript[F, i])(x)= Exp[-x/2] LaguerreL[i,x].

Now we have just to repeat the steps performed with the Hylleraas variables

δ (((λ^2)  M - λ L)/N) = 0

Quadratic forms

The computation is straightforward as these coordinates are a linear combination of the Hylleraas variables. We leave it as an exercise. In notebook [15] we took the hard way, i.e. we directly obtained the quadratic forms by a change of variables from Cartesian coordinates. The result was:

M=(∫_0^∞)d𝓍 (∫_0^∞)d𝓎(∫_0^∞)d𝓏 (𝓍 ((𝓎^2)+(𝓏^2)+𝓍 (𝓎+𝓏)) (((((∂Ψ)/(∂𝓍))))^2) -2𝓍 𝓎 (𝓍+𝓎) (((∂Ψ)/(∂𝓎)) ((∂Ψ)/(∂𝓍)))- 2𝓍 𝓏 (𝓍+𝓏) (((∂Ψ)/(∂𝓏)) ((∂Ψ)/(∂𝓍)))+
𝓎 (2 (𝓍^2)+2 (𝓎+𝓏) 𝓍+𝓏 (𝓎+𝓏)) (((((∂Ψ)/(∂𝓎))))^2) +𝓏 (2 (𝓍^2)+2 (𝓎+𝓏) 𝓍+𝓎 (𝓎+𝓏)) (((((∂Ψ)/(∂𝓏))))^2) );
L=(∫_0^∞)d𝓍 (∫_0^∞)d𝓎(∫_0^∞)d𝓏  ((1/4) ((((𝓎-𝓏))^2)-(((2 𝓍+𝓎+𝓏))^2)+4 Z (𝓎+𝓏) (2 𝓍+𝓎+𝓏)) ((Ψ )^2));
N=(∫_0^∞)d𝓍 (∫_0^∞)d𝓎(∫_0^∞)d𝓏 ((1/2) ((𝓍+𝓎) (𝓍+𝓏) (𝓎+𝓏))(Ψ^2)) ;
K=(∫_0^∞)d𝓍 (∫_0^∞)d𝓎(∫_0^∞)d𝓏 ( 2 𝓍 𝓎 𝓏 (((((∂Ψ)/(∂𝓍))))^2) +2𝓍 𝓎 (𝓍+𝓎) (((∂Ψ)/(∂𝓎)) ((∂Ψ)/(∂𝓍)))+
2𝓍 𝓏 (𝓍+𝓏) (((∂Ψ)/(∂𝓏)) ((∂Ψ)/(∂𝓍)))-2 𝓍 𝓎 (𝓍+𝓎+𝓏) (((((∂Ψ)/(∂𝓎))))^2) -2𝓎 𝓏 (𝓎+𝓏) (((∂Ψ)/(∂𝓏)) ((∂Ψ)/(∂𝓎)))-2 𝓍 𝓏 (𝓍+𝓎+𝓏) (((((∂Ψ)/(∂𝓏))))^2) );

The corrsponding matrix elements, on basis functions ψ (i) are evidently

M[i,j]=(∫_0^∞)d𝓍 (∫_0^∞)d𝓎(∫_0^∞)d𝓏 (𝓍 ((𝓎^2)+(𝓏^2)+𝓍 (𝓎+𝓏)) ((∂ψ(i))/(∂𝓍)) ((∂ψ(j))/(∂𝓍))-𝓍 𝓎 (𝓍+𝓎) (((∂ψ(i))/(∂𝓎)) ((∂ψ(j))/(∂𝓍))+((∂ψ(i))/(∂𝓍)) ((∂ψ(j))/(∂𝓎)))-
𝓍 𝓏 (𝓍+𝓏) (((∂ψ(i))/(∂𝓏)) ((∂ψ(j))/(∂𝓍))+((∂ψ(i))/(∂𝓍)) ((∂ψ(j))/(∂𝓏)))+
𝓎 (2 (𝓍^2)+2 (𝓎+𝓏) 𝓍+𝓏 (𝓎+𝓏)) ((∂ψ(i))/(∂𝓎)) ((∂ψ(j))/(∂𝓎))+𝓏 (2 (𝓍^2)+2 (𝓎+𝓏) 𝓍+𝓎 (𝓎+𝓏)) ((∂ψ(i))/(∂𝓏)) ((∂ψ(j))/(∂𝓏)));
L[i,j] =(∫_0^∞)d𝓍 (∫_0^∞)d𝓎(∫_0^∞)d𝓏  ((1/4) ((((𝓎-𝓏))^2)-(((2 𝓍+𝓎+𝓏))^2)+4 Z (𝓎+𝓏) (2 𝓍+𝓎+𝓏)) ψ(i) ψ(j) );
N[i,j] =(∫_0^∞)d𝓍 (∫_0^∞)d𝓎(∫_0^∞)d𝓏 ((1/2) ((𝓍+𝓎) (𝓍+𝓏) (𝓎+𝓏)) ψ(i) ψ(j)) ;
K[i,j] =(∫_0^∞)d𝓍 (∫_0^∞)d𝓎(∫_0^∞)d𝓏 ( 2 𝓍 𝓎 𝓏 ((∂ψ(i))/(∂𝓍)) ((∂ψ(j))/(∂𝓍))+𝓍 𝓎 (𝓍+𝓎) (((∂ψ(i))/(∂𝓎)) ((∂ψ(j))/(∂𝓍))+((∂ψ(i))/(∂𝓍)) ((∂ψ(j))/(∂𝓎)))+
𝓍 𝓏 (𝓍+𝓏) (((∂ψ(i))/(∂𝓏)) ((∂ψ(j))/(∂𝓍))+((∂ψ(i))/(∂𝓍)) ((∂ψ(j))/(∂𝓏)))-2 𝓍 𝓎 (𝓍+𝓎+𝓏) ((∂ψ(i))/(∂𝓎)) ((∂ψ(j))/(∂𝓎))-𝓎 𝓏 (𝓎+𝓏) (((∂ψ(i))/(∂𝓏)) ((∂ψ(j))/(∂𝓎))+((∂ψ(i))/(∂𝓎)) ((∂ψ(j))/(∂𝓏)))-2 𝓍 𝓏 (𝓍+𝓎+𝓏) ((∂ψ(i))/(∂𝓏)) ((∂ψ(j))/(∂𝓏))) .

We write them in Mathematica form:

Clear["Mij*","Lij*","Nij*","Kij*"]

Mij[i_,j_]:=
𝓏(2 (𝓍^2)+ 2(𝓎+𝓏)𝓍 + 𝓎(𝓎+𝓏))D[ψ[i],𝓏] D[ψ[j],𝓏]+
𝓎(2 (𝓍^2)+ 2(𝓎+𝓏)𝓍 + 𝓏(𝓎+𝓏))D[ψ[i],𝓎] D[ψ[j],𝓎]+
𝓍 (𝓎^2 + 𝓏^2 + 𝓍(𝓎+𝓏))D[ψ[i],𝓍] D[ψ[j],𝓍] -
𝓍 𝓎(𝓍+𝓎)( D[ψ[i],𝓍] D[ψ[j],𝓎]+ D[ψ[i],𝓎] D[ψ[j],𝓍])-
𝓍 𝓏 (𝓍 + 𝓏)(D[ψ[i],𝓍] D[ψ[j],𝓏]+ D[ψ[i],𝓏] D[ψ[j],𝓍]);

Lij[i_,j_]:= (1/4)( (𝓎-𝓏)^2 -(2𝓍+𝓎+𝓏)^2 + 4(𝓎+𝓏)(2𝓍+𝓎+𝓏)Z)ψ[i]ψ[j];
Nij[i_,j_]:= (1/2)( (𝓍+𝓎)(𝓍+𝓏)(𝓎+𝓏))ψ[i]ψ[j];

Kij[i_,j_]:= 
2𝓍 𝓎 𝓏  D[ψ[i],𝓍] D[ψ[j],𝓍] -2 𝓍 𝓎(𝓍+𝓎+𝓏) D[ψ[i],𝓎] D[ψ[j],𝓎]-
2 𝓍 𝓏(𝓍+𝓎+𝓏)D[ψ[i],𝓏] D[ψ[j],𝓏] +
𝓍 𝓎(𝓍+𝓎)( D[ψ[i],𝓍] D[ψ[j],𝓎]+ D[ψ[i],𝓎] D[ψ[j],𝓍]) -
𝓎 𝓏 (𝓎+𝓏)(D[ψ[i],𝓎] D[ψ[j],𝓏]+ D[ψ[i],𝓏] D[ψ[j],𝓎])+
𝓍 𝓏(𝓍+𝓏)(D[ψ[i],𝓍] D[ψ[j],𝓏]+ D[ψ[i],𝓏] D[ψ[j],𝓍]);

Matrices

Now we have to give an order to our base vectors, which are (with the choice α = 2; β = 1 = γ ):

Exp[- 𝓍 - 𝓎/2 - 𝓏/2] LaguerreL[i1,2𝓍]LaguerreL[i2,𝓎]LaguerreL[i3,𝓏].

We choose to give an ordering based on the degree of the polynomials. Our free parameters will be ω, the maximum degree, and all polynomials of total degree less than ω will be taken. This allows simpler control on the dimension of the basis. The following function just writes the base functions looping on the indices of the Laguerre polynomials, and taking into account the exchange symmetry:

getBase[ω_,parity_]:=
(signparity = If[parity== 0,1,-1];
ind=0;
Do[
Do[
Do[
If[i1+i2+i3 ≤ ω ,
ind=ind+1;
ψ[ind]= (1/2) Exp[-𝓍-(𝓎/2)-(𝓏/2)](LaguerreL[i1,2𝓍]LaguerreL[i2,𝓎]LaguerreL[i3,𝓏]+
signparity*LaguerreL[i1,2𝓍]LaguerreL[i2,𝓏]LaguerreL[i3,𝓎]);
];,
{i1,0,ω-i2-i3}],{i2,i3+parity,ω-i3}],{i3,0,ω}];
dim = ind;
base = Table[ψ[k],{k,1,dim}];
{base,dim});

Note : For odd states the term with i2=i3 will be identically 0, this is why in the loop there is the instruction {i2, i3 + parity, ω - i3}, for parity = 1 the zero term is avoided.

Now there is again a trick to perform quickly integrals. Once expanded products of polynomials and exponents are a sum of terms of the form

(x^a) (y^b)(z^c)Exp[-2x - y - z]

and their integral in x, y, z is trivial

(∫_0^∞)𝕕x(∫_0^∞)𝕕y(∫_0^∞)𝕕z  (x^a) (y^b)(z^c)Exp[-2x - y - z] = ((a! b! c!)/(2^(a+1))).

We can again perform all integrals just by substitution with the function

myInt[x_]:=(pow=((1/(2^((#[[1]])+1)))(#[[1]])!*(#[[2]])!*(#[[3]])!)&/@ Drop[Exponent[#,{𝓍,𝓎,𝓏}]&/@Apply[List,Expand[x+aa]],1];
cf = Drop[Apply[List,Expand[x+aa]]/.{E^c_.-> 1,𝓍->1,𝓎-> 1,𝓏-> 1,aa-> 0},1];
cf.pow)

The function getMatrix will work again :

getMatrix[X_]:= Module[{mat},
mat = 0*IdentityMatrix[dim];
Table[mat[[i,j]]=myInt[X[i,j]];,{i,1,dim},{j,i,dim}];
mat + Transpose[mat] - DiagonalMatrix[Tr[mat,List]]];

An example

Let us choose ω = 3 and even states for Helium:

ω = 3; parity = 0; Z=2;

We get the base functions :

{base,dim}=getBase[ω,parity];

This is the dimension of our basis :

dim

We can now produce matrices (unfortunately the computation is not very fast)

(matM = getMatrix[Mij];
matL = getMatrix[Lij];
matN = getMatrix[Nij];
matK = getMatrix[Kij];
)//Timing

A numerical check

As in the case of Hylleraas variables we use this notebook to check some other computations. In this case the question is quite important: in notebook [17] matrices are generated in a completely different way: one can use recursion relations of Laguerre polynomials and their orthogonality properties to write explicitly the matrix elements. A good feature of the method is that it is written in analytic form, i.e. it can be used for arbitrary ω;  a drawback is that it contains very long sequences of terms. matM for example has 747 terms and it is clear that some check to its structure is welcome. Here we report the matrices computed with this alternative procedure for ω = 3.

matM1 = {{4,(1/2),-(1/2),0,-(9/4),-1,(1/4),(1/2),(1/2),0,-(5/2),(1/2),1},{(1/2),14,-1,-(3/2),-4,-(23/4),-2,(5/2),(3/2),0,(1/2),-(7/2),0},{-(1/2),-1,29,-(9/2),(5/4),-10,-(35/4),0,5,0,0,1,0},{0,-(3/2),-(9/2),49,0,(15/4),-18,0,0,0,0,0,0},{-(9/4),-4,(5/4),0,(15/2),2,-(1/2),-(19/4),-2,(3/4),-(9/4),-1,-(13/4)},{-1,-(23/4),-10,(15/4),2,(45/2),2,-7,-(49/4),(15/4),-4,-(23/4),(7/4)},{(1/4),-2,-(35/4),-18,-(1/2),2,(85/2),(5/4),-16,0,(5/4),-10,0},{(1/2),(5/2),0,0,-(19/4),-7,(5/4),21,(23/4),-(99/8),-7,1,(7/8)},{(1/2),(3/2),5,0,-2,-(49/4),-16,(23/4),54,-15,1,-9,-2},{0,0,0,0,(3/4),(15/4),0,-(99/8),-15,(155/4),3,0,-(27/4)},{-(5/2),(1/2),0,0,-(9/4),-4,(5/4),-7,1,3,49,(9/2),-(29/2)},{(1/2),-(7/2),1,0,-1,-(23/4),-10,1,-9,0,(9/2),105,-14},{1,0,0,0,-(13/4),(7/4),0,(7/8),-2,-(27/4),-(29/2),-14,(163/4)}};
matL1 ={{(27/2),-2,-(1/2),0,-(25/2),(3/2),0,4,0,0,3,0,0},{-2,(33/2),-2,-(3/2),(3/2),-(31/2),3,0,4,0,0,3,0},{-(1/2),-2,(33/2),0,0,3,-(37/2),0,0,0,0,0,0},{0,-(3/2),0,(27/2),0,0,(9/2),0,0,0,0,0,0},{-(25/2),(3/2),0,0,(99/4),-(5/2),-(1/4),-(41/2),(3/2),6,-(37/2),(3/2),5},{(3/2),-(31/2),3,0,-(5/2),(117/4),-4,(3/2),-(47/2),0,(3/2),-(43/2),0},{0,3,-(37/2),(9/2),-(1/4),-4,(129/4),0,3,0,0,3,0},{4,0,0,0,-(41/2),(3/2),0,(207/4),-4,-(171/4),6,0,-(49/4)},{0,4,0,0,(3/2),-(47/2),3,-4,(237/4),(9/4),0,6,(3/4)},{0,0,0,0,6,0,0,-(171/4),(9/4),(369/4),0,0,(9/2)},{3,0,0,0,-(37/2),(3/2),0,6,0,0,(183/2),-8,-53},{0,3,0,0,(3/2),-(43/2),3,0,6,0,-8,(213/2),3},{0,0,0,0,5,0,0,-(49/4),(3/4),(9/2),-53,3,(345/4)}};
matN1 ={{4,-(5/2),(1/2),0,-(19/4),2,-(1/4),(3/2),-(1/2),0,(9/2),-(1/2),-1},{-(5/2),10,-7,(3/2),2,-(37/4),5,-(1/2),(5/2),0,-(1/2),(11/2),0},{(1/2),-7,19,-(27/2),-(1/4),5,-(61/4),0,-1,0,0,-1,0},{0,(3/2),-(27/2),31,0,-(3/4),9,0,0,0,0,0,0},{-(19/4),2,-(1/4),0,(21/2),-4,(1/2),-(33/4),3,(9/4),-(63/4),3,(33/4)},{2,-(37/4),5,-(3/4),-4,(39/2),-10,3,-(59/4),-(3/4),3,-(89/4),-(3/4)},{-(1/4),5,-(61/4),9,(1/2),-10,(63/2),-(1/4),7,0,-(1/4),7,0},{(3/2),-(1/2),0,0,-(33/4),3,-(1/4),19,-(31/4),-(129/8),13,-1,-(139/8)},{-(1/2),(5/2),-1,0,3,-(59/4),7,-(31/4),36,6,-1,15,2},{0,0,0,0,(9/4),-(3/4),0,-(129/8),6,(137/4),-3,0,(51/4)},{(9/2),-(1/2),0,0,-(63/4),3,-(1/4),13,-1,-3,55,-(29/2),-(91/2)},{-(1/2),(11/2),-1,0,3,-(89/4),7,-1,15,0,-(29/2),87,8},{-1,0,0,0,(33/4),-(3/4),0,-(139/8),2,(51/4),-(91/2),8,(297/4)}};
matK1 ={{0,(5/2),-(1/2),0,-(5/4),-2,(1/4),(1/2),(1/2),0,(7/2),(1/2),-1},{(5/2),-2,7,-(3/2),1,-(19/4),-4,-(3/2),(3/2),0,(1/2),(13/2),0},{-(1/2),7,-7,(27/2),-(3/4),4,-(39/4),0,-3,0,0,1,0},{0,-(3/2),(27/2),-15,0,-(9/4),9,0,0,0,0,0,0},{-(5/4),1,-(3/4),0,-(1/2),4,-(1/2),-(3/4),-3,(3/4),-(13/4),-4,(19/4)},{-2,-(19/4),4,-(9/4),4,-(3/2),10,2,-(17/4),-(9/4),-1,-(43/4),-(1/4)},{(1/4),-4,-(39/4),9,-(1/2),10,-(11/2),-(3/4),6,0,-(3/4),0,0},{(1/2),-(3/2),0,0,-(3/4),2,-(3/4),-1,(27/4),-(15/8),-3,1,-(13/8)},{(1/2),(3/2),-3,0,-3,-(17/4),6,(27/4),-6,(9/2),1,3,-(3/2)},{0,0,0,0,(3/4),-(9/4),0,-(15/8),(9/2),-(9/4),3,0,-(15/4)},{(7/2),(1/2),0,0,-(13/4),-1,-(3/4),-3,1,3,1,(37/2),-(13/2)},{(1/2),(13/2),1,0,-4,-(43/4),0,1,3,0,(37/2),-3,0},{-1,0,0,0,(19/4),-(1/4),0,-(13/8),-(3/2),-(15/4),-(13/2),0,-(17/4)}};

For historical reasons they differ by an inessential overall factor 2 from the one generated here (this factor cancels in the variational calculation). Let us check the consistency:

{matM1== 2 matM, matL1 == 2matL, matN1 == 2matN, matK1 == 2 matK}

Eigenvalues

Computation of eigenvalues proceeds along the standard ways.

calcEne[λ_Real,state_Integer]:=Sort[N[Eigenvalues[{λ^2 matM - λ matL,matN}]]][[state]]

We can see the λ dependence

Plot[calcEne[x,1],{x,1,3}]

and compute the minimum

NumberForm[minimum = FindMinimum[calcEne[x,1],{x,1.5,2.5}],10]

References

1) E.A. Hylleraas, Z. Phys. 48, 469 (1928).

2) C.L. Pekeris, Phys. Rev. 112, 1649 (1958).

3) H.A. Bethe and E.E. Salpeter, Quantum Mechanics of One- and Two-electron atoms, Springer (1977).