Helium, excited states

Quantum Mechanics
 K. Konishi, G. Paffuti

Extend the variational calculation to excited states of Helium atom.

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

Introduction

This notebook is an extension of notebooks [16] and [17]. We do not repeat theoretical basis and we assume that the reader has used functions constructed in those notebooks.

We neglect completely relativistic corrections and our effective Hamiltonian has the form

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

μ is the reduced mass electron-nucleus. At lowest order in the ratio m/M diagonalization of the Hamiltonian is converted in the variational formalism (approximatively) to the solution of a finite matrix problem

((λ^2)(Subscript[M, ij]) - (Subscript[λL, ij]))(Subscript[c, j]) = ℰ (Subscript[N, ij])(Subscript[c, j])

and to the minimization with respect to λ of the eigenvalue. The coefficients (Subscript[c, i])form the expansion of eigenstates in the basis of Hylleraas

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

the parameter shT has value 0, 1 for even and odd states respectively, while variables s, t, u are defined by

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

Corrections to eigenvalues found with eq () are mainly given by the mass orrections, the only ones considered in this notebook. They arise from the change from electron mass to the reduced mass in  () and from the last term in (), which will be represented by a matrix K. The four matrices M,N,L,K are computed by the function matrices, which needs as input Z, the nuclear charge, and maxs, maxt, maxu, the maximum allowed exponents i, j, k in the expansion of ψ.

The energy will be computed by the function energy[state, λ] (for a givenvalue of the scale parameter) and the function FindParameter[state,λ] will minimize the result with respect to λ.

In this notebook our aim is to give an idea of λ dependence of the results and a rough estimate of excited energy levels in parahelium and orthohelium.

Functions

These funtions are copied from notebook  [17].

Matrix Elements

matrices[Z_,maxs_,maxt_,maxu_]:=
(
maxinds = (maxs+1); maxindt = (maxt+1);maxindu = (maxu+1);
dim = maxinds*maxindt*maxindu;

matN = matL = matM = matN = matK = 0*IdentityMatrix[dim];

Do[  Do[ Do[Do[Do[Do[
i=i1+maxinds*j1+maxinds*maxindt*k1+1;
j=i2+maxinds*j2+maxinds*maxindt*k2+1;
a=i1+i2;b=2*(j1+j2 + shT);c=k1+k2;
 norma=(i1+2*j1+k1)!*(i2+2*j2+k2)!;
fatt3=(3+a+b+c)!;
fatt4=fatt3*(4+a+b+c);

fattN=fatt4/((1+b)*(3+b)*(3+b+c)*(5+b+c));
matN[[i,j]]=2*(5+a+b+c)*(6+2*b+c)*fattN/norma;

fattM=fatt3/(2*(1+b)*(3+b+c));
x1n=((6+2*b+c)*(20+a-a^2+b^2+9*c+c^2+b*(9+2*c)+4*i1*i2));
x1d=(3+b)*(5+b+c);x1=x1n/x1d;
x2=(16*(2+2*b+c)*(j1+shT2)*(j2+shT2))/((-1+b)*(1+b+c));
x3=-(4*(4+a+b+c-2*i1)*k2)/(3+b);
x4=(16*(j1+shT2)*k2)/(1+b+c);
x5=(4*(4+2*b+c)*k1*k2)/((3+b)*(1+b+c));
matM[[i,j]]=(x1+x2+x3+x4+x5)*fattM/norma;

matL[[i,j]]=fatt4/((1+b)*(3+b))*(2*(5+2*b+c)/(-(2+b+c)*(4+b+c))+4*(3+b)*Z/(3+b+c))/norma;

fattK=fatt3/(2*(3+b+c));
kin1=-((c*(20+a-a^2+b^2+9*c+c^2+b*(9+2*c)+4*i1*i2))/((1+b)*(3+b)*(5+b+c)));
kin2=(16*c*(j1+shT2)*(j2+shT2))/((-1+b^2)*(1+b+c));
kin3=(4*(4+a+b+c-2*i2)*k1)/((1+b)*(3+b));
kin4=(-16*(j2+shT2)*k1)/((1+b)*(1+b+c));
kin5=(-4*(4+2*b+c)*k1*k2)/((1+b)*(3+b)*(1+b+c));
matK[[i,j]]=fattK*(kin1+kin2+kin3+kin4+kin5)/norma;,
{k2,0,maxu}],{j2,0,maxt}],{i2,0,maxs}],
{k1,0,maxu}],{j1,0,maxt}],{i1,0,maxs}];
    matL =(1/2)(matL+Transpose[matL]);
matM =(1/2)(matM+Transpose[matM]);
matN =(1/2)(matN+Transpose[matN]);
matK =(1/2)(matK+Transpose[matK]);
{matL,matM,matN,matK})

energy

energy[state_?IntegerQ,λ_?NumberQ]:=Sort[Re[Eigenvalues[{λ^2 matM-λ matL,matN}]]][[state]];

findParameter[state_,λstart_]:= Module[
{dE=1,maxTry=50,dEprec=(10^(-10)),dλprec =(10^(-3)),ene,dλ,λ},
λ = λ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;];
λ]

The dependence on λ

Z = 2;
shT = 0; 
shT2 = shT/2;
maxs = 3;maxt = 5;maxu=5;

Here we compute the matrices needed for eigenvalues. The time of computation depends eavily on the computer. On very slow systems the reader can take smaller values for parameters maxs, maxt, maxu.

{matL,matM,matN,matK} = N[matrices[Z,maxs,maxt,maxu]];

Let us see how different eigenvalues depend on λ

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

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

From the previous graph it is clear that with this basis there is no range in λ for which all eigenvalues are stable and, as expected, higher eigenvalues show a sharper λ dependence.

To compare with experimental data let us collect in a single function the computation of both variational energy and mass corrections. We can specify as third argument a list of wanted states, otherwise the list will consist of a single element, the state choosen as first argument. The list can help to compare results obtained minimizing λ for each state (the correct procedure) with eigenvalues computed for a single λ, the same for all states.

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);];

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}]];
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}];
ene+de1+de2]

First we adopt the correct procedure, i.e. we minimize each eigenvalue with respect to λ :

λstart =  N[2*(Z-(5/16))];

eneStates=Flatten[Table[λ=findParameter[i,λstart];
If[i==1,λGround = λ;]; calcEneNR[i,λ],{i,1,5}]]

As a second option we compute eigenvalues with the same value of λ, for example the one obtained for the ground state

λGround = findParameter[1,λstart];
eneλGroundState =calcEneNR[1,λGround,Range[1,5]]

The ionization energy of the system is about - (Z^2)/2 = -2 a.u. than the two last estimates are surely wrong, having an energy greater than this. This stress the importance of using the correct λ, especially for relatively small basis. The first two excited levels are more consistent, the stablity would be increasd by taking a larger basis. Unfortunately with our choices the problem becomes quite unstable if we adopt a larger basis, we will cure this problem by switching to another base system, based on the Pekeris triangular variables in notebook [].

As far as experimental data is concerned, these are the experimental energy differences between excited states and ground state in parahelium, expressed in (cm^(-1))

If[Z==2 && shT== 0,diffExp = {166277.54,184864.93,190940.33,193663.61},diffExp=Table["",{j,1,4}]];

If the reader wishes to check other ions or orthohelium (s)he can write relevant numbers in the list diffExp.

unitRy = 109737.3156; unitcm = 2*unitRy;

These are estimated energy differences

diffEne1 = (eneStates[[2;;5]]-eneStates[[1]])unitcm;

and these one are expected to be a worse estimate

diffEne2 = (eneλGroundState[[2;;5]]-eneStates[[1]])unitcm;

string = {"1s2s-1(s^2)","1s3s-1(s^2)","1s4s-1(s^2)","1s5s-1(s^2)"};

NumberForm[TableForm[Transpose[{diffExp,diffEne1,diffEne2}],
TableHeadings-> {string,{"  DATA","λ ∀ state"," λ-fixed",}}
],8]