Helium - Pekeris variables

Quantum Mechanics
 K. Konishi, G. Paffuti

Compute the s-state levels of Helium using the triangular (Pekeris) variables.

Unprotect[dimIndR,indR,segnoParita,αKP,βKP,γKP,formamatM,formamatL,formamatN,formamatK];
Unprotect[x,y,z,u,s,t,i1,i2,i3,j1,j2,j3];
ClearAll["Global`*"];
Off[General::spell, General::spell1,InterpolatingFunction::"dmval"];
Protect[x,y,z,u,s,t,i1,i2,i3,j1,j2,j3];

N.B. In this notebook there are some global variables which must preserve their values, to avoid unwanted changes we protect them. αKP, βKP, γKP are scaling parameters for variables, segnoParita[] is a function which is ± 1 depending on exchange symmetry of the wave function, indR and dimIndR are a list of indices and its dimensions. Variables formamatM, formamatL, formamatN, formamatK contain the algebraic expressions of the quadratic forms. We choose to use them as global variables in order to simplify the life to the reader. Names x, y, z, s, t, u are reserved for Hylleraas and Pekeris variables respectively, so they are not at disposal of the reader.

As a general rule if the reader adds some instructions to the notebook it is better not to use short names and start the computation from the beginning, i.e. with ClearAll etc.

segnoParita[0]=1; segnoParita[1]=-1;
Protect[dimIndR,indR,segnoParita,αKP,βKP,γKP,formamatM,formamatL,formamatN,formamatK];

Formulation of the problem

In this notebook we present a formulation of the variational problem for Helium which allows easily the use of large basis, and consequently a more precise estimate of the energy levels, both for the ground and excited states. This approach has been mainly developed by Pekeris []. The main difference with respect usual Hylleraas approach lie in the choice of variables, here triangular variables, defined as

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

For notations and details see notebook [NB-10.15-Coord-Helium.nb], the discussion below is just a copy of part of that notebook.  These variables have an independent range

0 ≤ x, y, z ≤ ∞ ;

and this allows the definition of a basis function in terms of factorized functions, based coveniently on some orthogonal polynomials. Basis vectors are taken of the 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, than in basis 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 θ dependent exponentials we have to choise 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 now we can choose as basis functions some component of a complete set in the interval [0,∞]. The natural choice is to consider 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 Hylleraas variables

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

In notebooks [15] and [16] we have written the relevant quadratic forms and we report here the general structure of matrix elements:

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))/(∂𝓏))) .

K is the matrix needed to compute mass corrections.

In notebook  [16] matrix elements have been computed analytically, for each value of the dimension of the basis. Here we follow another procedure which allows a computation in symbolic form, valid for every dimension. The key point is that both derivative and multiplication by x for Laguerre polynomials can be written symbolically using the relations:

x L'[n,x] - n L[n,x] - n L[n-1,x];  
x L[n,x] = - (n+1)L[n+1,x] + (2n+1)L[n,x] - n L[n-1,x];

This implies that an arbitrary quadratic form in derivatives of base functions can be reduced to products of Laguerre polynomials and exponents. To this expression we can apply the orthogonality property

∫Exp[-x] L[n,x] L[k,x]dx = (Subscript[δ, nk]);

which allows a trivial computation of all integrals involved. This quite boring task is performed by the function funzExpand[term] which accepts an expression and applies in sequence, for each of the three variables x, y, z  the previous relations.

For the reader not interested in the details of the computation we summarize the only steps needed to perform calculations, only two function calls are needed:

The next section can be skipped by the non interested reader, for interested ones a brief description of the computation startegy and of the functions is presented

Basis and matrices

funzIndice

The function funzIndice has already been used in notebook []. It simply creates a numeration for vectors with an incremental loop. The values of the index are stored in indR[i1,i2,i3], where i1,i2,i3 are the degrees of Laguerre polynomials, te value is the integer associate dwith this combinationof values. The only tricky point is that the loop on i2 starts from i3 for even states and form i3+1 for odd states. This is because for odd states expression () vanishes when second and third indexes are equal.

funzIndice[parita_,ω_]:=Module[{ind},
ind=0;
Unprotect[dimIndR,indR];
Do[
Do[
Do[
If[i1+i2+i3 ≤ ω ,
ind=ind+1;
indR[i1,i2,i3]= ind;
];,
{i1,0,ω-i2-i3}],{i2,i3+parita,ω-i3}],{i3,0,ω}];
dimIndR = ind;
Protect[dimIndR,indR];
];

scambia

This is another utility function : just exchange second and third indexes in a matrix and add the results with ± depending on the symmetry. At the same time exchange left and right indexes to ensure hermiticity of the matrices. This trick only works for β = γ, which (see below) is our choice.

scambia[segno_,mat_]:=Module[{provv,provv1,provv2,provv3,provv4,provv5,
provv6,provvj,provvi,provvij,xx,yy,matOut},
provv1 = mat/.{j2-> xx};
provv2 =  provv1/.{j3-> j2};
provvj= provv2/.{xx-> j3};

provv3 = mat/.{i2-> yy};
provv4 = provv3/.{i3-> i2};
provvi = provv4/.{yy-> i3};

provv5 = provv1/.{i2-> yy};
provv6= provv5/.{j3-> j2,i3-> i2};
provvij = provv6/.{xx-> j3,yy-> i3};
matOut =  Expand[(1/4)(segno*(provvj + provvi)+mat+provvij)];
matOut];

computeAlgebricMatrices

This function constructs a formal expression for quadratic forms. We use rescaled variables, i.e. x → x/α, y → y/β , z → z/γ in such a way that Lagerre polynomial have a simple dependence. As explained in notebook [] we will have an overall scale, which will be called λ, while the relative values of the parameters has been chosed as α=2, β=1, γ = 1.

In the function F and G are a symbolic representation of left and right wave functions in the quadratic forms (i.e. bra and ket). With each of them we construct the linear differential corresponding at derivatives with respect Hylleraas variables s, t, u, then we report the quadratic form in Hylleraas form, i.e. with s, t, u variables. This detour with Hylleraas variables is mainly do to hystorical reasons, this was the first option we used, before producing explicit expression for quadratic forms in term of x, y, z directly. We decide to leave this procedure because the concidence with result of notebook [] is a good crossed check.

The function end with a call to funzExpand, which perform the hard work of transforming these formal expression in formal matrix elements.

The result is stored in global variables, formamatM, formamatL, formamatN, formamatK. In usual programming these variables could be passed to getMatrices, which compute matrices numerically. To avoid this additional work to the reader we have defined these variables as global, in this way the reader canalso ignore the details of the computation.

computeAlgebricMatrices[Z_Integer]:= Module[{ds,ds1,dt,dt1,du,du1},
F[x_,y_,z_]:= Exp[-x/2]L[i1,x]Exp[-y/2]L[i2,y]Exp[-z/2]L[i3,z];G[x_,y_,z_]:= Exp[-x/2]L[j1,x]Exp[-y/2]L[j2,y]Exp[-z/2]L[j3,z];
ds = αKP(∂_x) F[x,y,z];dt = γKP(∂_z) F[x,y,z]-βKP(∂_y) F[x,y,z];du = -αKP(∂_x) F[x,y,z]+βKP(∂_y) F[x,y,z]+γKP(∂_z) F[x,y,z];
ds1 = αKP(∂_x) G[x,y,z];dt1 = γKP(∂_z) G[x,y,z]-βKP(∂_y) G[x,y,z];du1 = -αKP(∂_x) G[x,y,z]+βKP(∂_y) G[x,y,z]+γKP(∂_z) G[x,y,z];
espM=u(s^2-t^2)(ds*ds1 + dt*dt1 + du*du1)+du1 *(s(u^2-t^2)ds+t(s^2-u^2)dt)+
du *(s(u^2-t^2)ds1+t(s^2-u^2)dt1);
espN=u(s^2-t^2)*F[x,y,z]*G[x,y,z];
espL = (4*Z*u*s - (s^2-t^2))F[x,y,z]*G[x,y,z];
espK = (s^2+t^2 - 2 u^2)*u*(ds*ds1 - dt*dt1) - (s^2-t^2)*u*du*du1 -
du1*(s*(u^2-t^2)*ds + t*(s^2-u^2)dt) - du*(s*(u^2-t^2)*ds1 + t*(s^2-u^2)dt1);
Unprotect[formamatM,formamatL,formamatN,formamatK];
formamatM = funzExpand[espM];formamatN = funzExpand[espN];
formamatL = funzExpand[espL];formamatK = funzExpand[espK];
Protect[formamatM,formamatL,formamatN,formamatK];]

funzExpand

This is the core of the computation. The function acts on a quadratic form.

funzExpand[termine_]:= Module[{},
reg1= {s-> (x/(αKP)) +(1/2) ((z/(γKP))+(y/(βKP))), t-> (1/2)((z/(γKP))-(y/(βKP))),u-> (1/2)((z/(γKP))+(y/(βKP)))};
esp2a=Expand[(termine/.reg1)/(E^(-x-y-z))];
esp2a = esp2a/.{αKP-> 2,βKP-> 1,γKP-> 1};
(* sostituzione termini non derivativi  *)
regLL2 = a_.*w^2* b_.*L[j_,w]*c_.* L[i_,w]-> a*b*c*(
i(i-1)KroneckerDelta[i-2,j]-4i^2KroneckerDelta[i-1,j]+
2(1+3i+3i^2)KroneckerDelta[i,j]-4(1+i)^2KroneckerDelta[i+1,j]+(i+1)(2+i)KroneckerDelta[i+2,j]);
regLL1 = a_.*w* b_.*L[j_,w]*c_.* L[i_,w]-> a*b*c*
(- (i+1)KroneckerDelta[i+1,j] + (2i+1)KroneckerDelta[i,j] -
i KroneckerDelta[i-1,j] );
regLL0 =  a_.*L[j_,w]*b_.* L[i_,w]-> a*b*KroneckerDelta[i,j];
test1a = Expand[(esp2a/.x-> w)/.regLL2];
test1a = test1a/.w-> x;
test1a = Expand[(test1a/.y-> w)/.regLL2];
test1a = test1a/.w-> y;
test1a = Expand[(test1a/.z-> w)/.regLL2];
test1a = test1a/.w-> z;
test1a = Expand[(test1a/.x-> w)/.regLL1];
test1a = test1a/.w-> x;
test1a = Expand[(test1a/.y-> w)/.regLL1];
test1a = test1a/.w-> y;
test1a = Expand[(test1a/.z-> w)/.regLL1];
test1a = test1a/.w-> z;
test1a = Expand[(test1a/.x-> w)/.regLL0];
test1a = test1a/.w-> x;
test1a = Expand[(test1a/.y-> w)/.regLL0];
test1a = test1a/.w-> y;
test1a = Expand[(test1a/.z-> w)/.regLL0];
test1a = test1a/.w-> z;
(* sostituzione termini  derivativi  *)
regdL2 = a_.*w^2* b_.*Derivative[0,1][L][j_,w]*c_.* Derivative[0,1][L][i_,w]-> a*b*c*(
i*j(KroneckerDelta[i,j] + 
KroneckerDelta[i-1,j-1]-KroneckerDelta[i-1,j]-
KroneckerDelta[i,j-1]));
regdL1 = a_.*w* b_.*Derivative[0,1][L][j_,w]*c_.* Derivative[0,1][L][i_,w]-> a*b*c*
i*KroneckerDelta[i,j] ;
regdL01 = a_.*w* b_.*L[j_,w]*c_.* Derivative[0,1][L][i_,w]-> a*b*c*i*(
KroneckerDelta[i,j] - KroneckerDelta[i-1,j]);
regdL02 = a_.*w^2* b_.*L[j_,w]*c_.* Derivative[0,1][L][i_,w]-> a*b*c*i(
(j+1)KroneckerDelta[i,j+1] + (2j+1)KroneckerDelta[i,j] - j*
KroneckerDelta[i,j-1] - (j+1)KroneckerDelta[i-1,j] -(2j+1)*
KroneckerDelta[i-1,j] +j* KroneckerDelta[i-1,j-1] );
regdL02= a_.*w^2* b_.*L[j_,w]*c_.* Derivative[0,1][L][i_,w]-> a*b*c*i(
(-1+i)KroneckerDelta[i-2,j] + (1-3*i)KroneckerDelta[i-1,j]+
(3*i+1)*KroneckerDelta[i,j]-KroneckerDelta[i+1,j]-i*KroneckerDelta[i+1,j]);
test1b = Expand[(test1a/.x-> w)/.regdL2];
test1b = test1b/.w-> x;
test1b = Expand[(test1b/.y-> w)/.regdL2];
test1b = test1b/.w-> y;
test1b = Expand[(test1b/.z-> w)/.regdL2];
test1b = test1b/.w-> z;
test1b = Expand[(test1b/.x-> w)/.regdL1];
test1b = test1b/.w-> x;
test1b = Expand[(test1b/.y-> w)/.regdL1];
test1b = test1b/.w-> y;
test1b = Expand[(test1b/.z-> w)/.regdL1];
test1b = test1b/.w-> z;
test1b = Expand[(test1b/.x-> w)/.regdL02];
test1b = test1b/.w-> x;
test1b = Expand[(test1b/.y-> w)/.regdL02];
test1b = test1b/.w-> y;
test1b = Expand[(test1b/.z-> w)/.regdL02];
test1b = test1b/.w-> z;
test1b = Expand[(test1b/.x-> w)/.regdL01];
test1b = test1b/.w-> x;
test1b = Expand[(test1b/.y-> w)/.regdL01];
test1b = test1b/.w-> y;
test1b = Expand[(test1b/.z-> w)/.regdL01];
test1b = test1b/.w-> z;
test1b];

getMatrices

After computeAlgebricMatrices has stored in global variables the algebraic forms a call to getMatrices produces their numerical values. With this distinction the slow algebraic computation must be performed only after the choice of the atom (this is essentially due to the fact that potential term contain Z, we could avoid any repetition splitting the L form, but for simplicity we leave things in this way). getMatrices must of course be called each time we change the dimension or the parity (exchange symmetry) of the state.

getMatrices[parita_,ω_]:= Module[{provv,xx,segno},
Nomega = (1/6) (1+ω) (2+ω) (3+ω);
(* compute indexes *)
funzIndice[parita,ω];
matM = Table[0,{i,1,dimIndR},{j,1,dimIndR}];
matL = matM; matN=matM; matK=matM;
segno = segnoParita[parita];
fM = scambia[segno,formamatM];
fL = scambia[segno,formamatL];
fN = scambia[segno,formamatN];
fK = scambia[segno,formamatK];

Do[
Do[
If[ Abs[i3-j3]≤ 2,
fM1=fM;
fL1=fL;
fN1=fN;
fK1=fK;
Do[
Do[
If[ Abs[i2-j2]≤ 2,
fM2=fM1;
fL2=fL1;
fN2=fN1;
fK2=fK1;
Do[
Do[
If[i1+i2+i3 ≤ ω && j1+j2+j3≤ ω,
If[Abs[i1-j1]≤ 2 ,
ind1 = indR[i1,i2,i3];
ind2= indR[j1,j2,j3];
matM[[ind1,ind2]] = fM2;
matL[[ind1,ind2]] = fL2;
matN[[ind1,ind2]] = fN2;
matK[[ind1,ind2]] = fK2;
];];,
{i1,0,ω-i2-i3}],{j1,0,ω-j2-j3}]];,
{i2,i3+parita,ω-i3}],{j2,j3+parita,ω-j3}]];,
{i3,0,ω}],{j3,0,ω}];
matM = N[SparseArray[matM]];
matL = N[SparseArray[matL]];
matN = N[SparseArray[matN]];
matK = N[SparseArray[matK]];
{matM,matL,matN,matK}];

Computation

The atom

Let first show with an example a computation for ground state of the Helium atom

Z=2;

this compute algebric forms for matrices:

computeAlgebricMatrices[Z];

parity and dimensions

Now we choose the singlet part of the spectrum (parity = 0) and fix a maximum degree. This operation starts to be time consuming, if the reader has a slow computer he can take a lower value for ω.

parity = 0; ω = 15;

{matM,matL,matN,matK}=getMatrices[parity,ω];//Timing

The dimension of matrices is

Dimensions[matM]

The reader can appreciate how sparse are matrices form the following graph :

ArrayPlot[matM]

eigenvalues

With sparse matrices we can use specially designed algorithms to compute eigenvalues. The one used in based on a recursive research around a value fixed by the option Shift, it is the Arnoldi method (see Mathematica Help). For the remaining part the function energy is the same that has been used in other notebooks.

numMaxEig = 5;

energy[state_?IntegerQ,λ_?NumberQ]:=Sort[Re[Eigenvalues[{λ^2 matM-λ matL,matN},
numMaxEig+3,Method->{Arnoldi,Shift->-3.0000,BasisSize-> 30}]]][[state]];

You can appeciate the speed (remember that we are working with a basis of 440 elements).

energy[1,2.]//Timing

It is not too time consuming to produce a plot of various eigenvalues as a function of λ :

λplot=Range[0.1,4,0.1];

Show[
Table[ListPlot[Transpose[{λplot,Evaluate[(energy[i,#1]&/@λplot)]}],
AxesOrigin-> {0,-3},
PlotRange-> {All,{-3,-1}},Joined-> True,
AxesLabel-> {"λ","eig"}],{i,1,5}]]

As usual to have a definite value we have to minimize eigenvalues with respect the parameter λ. We use again the function findParameter presented in notebook [].

findParameter[state_]:= Module[
{dE=1,maxTry=50,dEprec=(10^(-10)),dλprec =(10^(-3)),ene,dλ,λ,
λstart,λup,λdown,ene1,ene2},
λstart =2.5 Z/state;
λ = λstart;
dλ = 0.2*λ; 
λTry =1;
ene= energy[state,λstart];
While[(λTry ≤ maxTry && dλ > dλprec && dE > dEprec),
λup = λ + dλ; ene1 = energy[state,λup];
If[ene1 ≤ ene,
λ=λup;dE = Abs[ene-ene1];ene=ene1;,
λdown = λ - dλ; ene2 = energy[state,λdown]; 
If[ene2 < ene,
λ=λdown;dE = Abs[ene-ene2];ene=ene2;,
dλ=dλ/2];
]; λTry = λTry +1;];
λ]

Here we report conversion factors, atomic masses etc. to get numbers in atomic units or in (cm^(-1)).

unitRy = 109737.3156; unitcm = 2*unitRy;
mAtom = {1,4.002602,6.941,9.012182};
mElectronAtomicUnits = 0.510998902/931.494013;
massRatio = mAtom[[Z]]/(mElectronAtomicUnits + mAtom[[Z]]);
(* for helium we use a more accurate value, taken from NIST *)
ratioαPartElectron = 7294.2995363; 
If[Z==2,massRatio = ratioαPartElectron/(ratioαPartElectron+1);];

This function, clacEneNR is almost identical to the analogous function in notebook  [16]. It needs as input a state number and a value for the scale parameter λ., produce a list with non relativistic energy, and the two mass corrections terms. A third optional argument in input, a list of states, gives as output the same values for the list, all values being computed with the optimal value for the λ specified by the variable state. This third variable canbe used to check if and how one obtains a bad result taking a single λ value for all sates.

calcEneNR[state_,λ_,listIn_:{}]:= Module[{de1,de2,ene,eigH,eigV,list},
If[Length[listIn]== 0,list = {state};,list=listIn;];
{eigH,eigV}=Re[Eigensystem[{λ^2 matM-λ matL,matN},state+8,Method->{Arnoldi,Shift->-3.0000,BasisSize-> 30}]];
order = Ordering[eigH];
eigH = eigH[[order]];
eigV = eigV[[order]];
ene = eigH[[list]];
(* these are mass corrections *)
de1=ene*(massRatio - 1);
de2 = Table[ λ^2*((eigV[[j]].(matK.eigV[[j]]))/((eigV[[j]].(matN.eigV[[j]]))))*(mElectronAtomicUnits/(mAtom[[Z]])),{j,list}];
If[Length[list]== 1,Flatten[{ene,de1,de2}],{ene,de1,de2}]]

Example

Here we compute the energy and corrections for ground state :

λGround = findParameter[1];
NumberForm[energyGroundState=calcEneNR[1,λGround],20]

Ground and excited states

We continue our computation for excited states, this would be quite difficult with Hylleraas variables. We print the execution times in such a way that the reader can judge if he have to lower the parameter ω.

(eneStates=Table[λ=findParameter[i]; calcEneNR[i,λ],{i,1,5}];)//Timing

NumberForm[TableForm[eneStates],10]

thisComputation[ω]=Flatten[eneStates[[All,1]]]

In the literature few years ago some very accurate results on the variational estimate of energy levels was presented, see for example [1]. Here we report their results for low levels, for both singlet (ene1S) and triplet (ene3S) configurations

ene1SLiterature = 
{-2.903724377034119859,-2.145974046054412,-2.06127198974087,-2.03358671688,-2.02117685115};

ene3SLiterature = -{2.175229378236791300,2.068689067472454,2.036512083098,2.02261887230,2.01537745299};

In the next two subsections we perform larger computations to check our results. These can take some time and for the reader which is not interested in such numerical business we just report the difference between our results and the ones present in literature , so (s)he can skip the next two subsections:

numerical checks

We perform higher order estimates, this take some time.

ω = 20

parity = 0; ω = 20;

{matM,matL,matN,matK}=N[getMatrices[parity,ω]];//Timing

(eneStates=Table[λ=findParameter[i]; Flatten[calcEneNR[i,λ]],{i,1,5}]);//Timing

thisComputation[ω]=Flatten[eneStates[[All,1]]];

ω = 25

parity = 0; ω = 25;

{matM,matL,matN,matK}=N[getMatrices[parity,ω]];//Timing

(eneStates=Table[λ=findParameter[i]; Flatten[calcEneNR[i,λ]],{i,1,5}]);//Timing

thisComputation[ω]=Flatten[eneStates[[All,1]]];

Dimensions[matM]

ωList = {15,20,25};

Here are the differences between our computation and the one cited above. We can see an agreement which clearly goes better with higher ω.

TableForm[Table[thisComputation[ωList[[i]]] - ene1SLiterature,{i,1,3}],
TableHeadings-> {Table[Row[{"ω = ",ωList[[i]]}],{i,1,3}],
Table[Row[{"  δE 1s",i,"s"}],{i,1,5}]
}]

Triplet

The same computations can be carried on for triplet states, we present only the result for ω = 15.

parity = 1; ω = 15;

{matM,matL,matN,matK}=N[getMatrices[parity,ω]];//Timing

(eneStates=Table[λ=findParameter[i]; Flatten[calcEneNR[i,λ]],{i,1,5}]);//Timing

NumberForm[TableForm[eneStates],10]

eneStates[[All,1]]-ene3SLiterature

References

1) A. Burgers, D. Wintgen and J.-M. Rost, J. Phys. B, At. Mol. Opt.  Phys.  28, 3163, (1995).