Helium

Quantum Mechanics
 K. Konishi, G. Paffuti

Study numerically the ground state of Helium using the Hylleraas coordinates.

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

The method

The general strategy adopted for a variational calculation of the energy of Helium states has been sketched in notebook [16], where some explicit examples are also carried out.

In this notebook  we shall limit our analysis to S - states, mainly the ground state of the atom. We will adopt  natural units (m = ℏ = e = 1), conversion factors to the standard units in spectroscopic analysis are given below.

In a first approximation the nucleus is taken fixed (infinite mass), the first order corrections in m/M (electron mass/ nucleus mass) are computed perturbatively.

The S states are completely described by a function of only three variables, the two electron - nucleus distances, (Subscript[r, 1])and (Subscript[r, 2])and their relative distance, (Subscript[r, 12]). A linear combination of these variables form the Hylleraas variables, s, t, u:

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.

The variational principle states that the energy eigenstates can be computed my minimizing the energy on a unit sphere of the Hilbert space (i.e. on vectors of unit norm):

δ ((⟨ψ|H|ψ⟩)/(⟨ψ|ψ⟩))= 0

⟨ψ | H | ψ⟩  and ⟨ψ | ψ⟩ are quadratic form in ψ, and in Hylleraas coordinates.  The ratio () can be expressed as

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

where

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

M is the kinetic energy, L the potential (or better the negative of potential energy) and N the norm.

In the variational approach used here ψ is a finite linear combination of basis functions (Subscript[f, i])[λ x]. The  parameter λ can be reabsorbed by a scale transformation x → x/λ and riappears as a coefficient of homogeneous functionals, M and L. Then the variational principicle requires stationarity both with respect to λ and to the coefficients of the linear combination of

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

The linear combination, with λ reabsorbed, is of the form

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

Stationarity with respect to c[i] imply the usual matrix form of Schrödinger equation (in a non orthogonal basis):

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

Stationarity with respect to λ requires minimization of the eigenvalue ℰ[λ] of this equation.

(Subscript[M, ij]) etc. are the quadratic forms computed on basis elements, i.e. matrix elements of the Hamiltonian.

In this notebook we follow Hylleraas' pioneering work and we will use as base functions:

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

The matrix elements in this basis have been computed in notebook [15] and are collected in the function matrices given below.  The function requires as input Z, the charge of the nucleus (i.e. we can do computations also for ions of lithium etc.) and the maximum value for indexes i, j, k in (), maxs, maxt, maxu.  To allow an easier manipulation with the usual matrix multiplication rules we order the wave functions with an index written as

I = i + (maxs+1)j + (maxs+1)(maxt+1)k + 1;

The range of the indices is

I = 1, 2, ... dim;      dim = (maxs+1)(maxt+1)(maxu+1);

as the reader can easily show.

Matrix elements depend on a global variable, shT which must be fixed before starting the computation. shT = 0 is for symmetric states, shT = 1 for antisymmetric ones (orthoHelium). The varable appears as a shift in t exponent, 2j+shT. In effect the exchange of the particles corresponds to t → - t, then even powers of t give symmetric states, odd powers antisymmetric ones.

The output of the function matrices gives the three matrices M, L, N and a fourth one, K, which gives matrix elements of the operator

K = (Subscript[p, 1]).(Subscript[p, 2])  .

In the text and in notebook [15] it has been shown that up to order m/M the Hamiltonian for Helium atom is

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 |ψ⟩)/(⟨ψ|ψ⟩)).

For some checks on the matrix elements written below see notebook [15].

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

Eigenvalues

N.B. For a basis with 120-150 elements the computation starts to show some instabilities, for larger dimensions it is better to switch to other basis functions, as described in notebook [15].

Initialization

We choose the atom, here Helium, Z = 2, and maximum values for powers of s, t, u in  Hylleraas basis. The starting value for λ is chosen as the one obtained in lowest order variational method with hydrogenoid wave functions, any other parameter of the same order of magnitude works the same.

Z=2;
maxs = 3;maxt = 4;maxu=5; 
λstart = N[2*(Z-(5/16))];

Here we fix the global variable shT, i.e. we choose between paraHelium (shT = 0, S=0) and orthoHelium (shT=1, S=1).

The reader must avoid the use of names shT and shT2 in modifications of the notebook.

shT = 0; 
shT2 = shT/2;

Here we compute matrices needed for eigenvalues. The time of computation depends heavily 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]];//Timing

functions

We need a function to compute eigenvalues and one to find minimum with respect to λ. The first one is simply Eigenvalues

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

For searching  the  value of  λ minimum  we could use Mathematica function FindMinimum, as in the following lines

(λ = xVar /. FindMinimum[energy[xVar, state], {xVar, 1, 5}, WorkingPrecision -> 20][[2]]) // Timing

We take instead another option, we start from λstart and make a binary search, i.e. we compute the function for λ ± dλ and look for a decrease in energy. If we found a smaller value we repeat the procedure with dλ → dλ/2. It is not very smart but works quite well. With FindMinimum we find, at least for Helium a slight different result, within eight digits which is in any case our precision. The reader can use FindMinimum if he want a simpler code.

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 results

As a general warm up we compute the energy for the starting value of the parameter. We choose to compute energy for ground state, then we put state = 1

state = 1;

NumberForm[ene0= energy[state,λstart],10]

Now we use all the above machinery and optimize with respect to λ :

NumberForm[λ = findParameter[state,λstart],6]

It is seen that λ has changed for a small amount, with this relatively large basis. This change is in any case essential if we want precise results.

With this new value for λ we can compute the energy and eigenvectors. In fact the energy could be extracted from findParameter, but we prefer to do again the computation both for clarity and for computing other eigenvalues and the eigenvectors.

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

NumberForm[ene = eigH[[state]],10]

The difference with respect the original value is only of  (10^(-6))a.u. , this is what we mean by a small change in λ.

To later comparison with experimental data we present the results in two ways, the energy is written in atomic units while the ionization potential, defined as the energy needed to extract one electron, is written in (cm^(-1)). The conversion factor is given below

unitRy = 109737.3156; unitcm = 2*unitRy;
J = (-((Z^2)/2)- ene)*unitcm;

The experimental value for J (including the case Z = 1, i.e. the negative ion (SuperMinus[H])) are, in (cm^(-1)):

Jexp={6100,198310.6691,610079., 1241242.};

Here is the comparison with experimental data

NumberForm[TableForm[
{{ene," ","a.u."},
{J,Jexp[[Z]],"(cm^(-1))"}},
TableHeadings-> {{"E","J"},{" calculated ", "   exp"}}
],12,NumberSigns-> {"-"," "}]

The difference in computed and experimental value is of the order  δJ/J ∼ 3×(10^(-5)), δJ ∼ 34 (cm^(-1)) for Helium.  In the next section we improve the result computing the main corrections to our result, which in any case is quite impressive.

Mass corrections and relativistic corrections

Mass corrections

The main correction to the previous result is quite trivial : we have to substitute the electron mass with the reduced mass. In problems with more than two bodies the reduced mass is not the only correction, a polarization term is also present, see eq.(). We first write ratios  m/M  for atoms up to Beryllium. For Helium we use a particularly accurate value, i.e. tha mass of an α particle (nucleus of Helium) because the experimental value for J is very accurate and we have to avoid any possible source of errors. For other atoms the atomic mass is used, which in fact differs from nucleus mass for higher orders in m/M. We note that also the energy of the ionized atom changes by the same factor, then the whole J must by ultiplied by μ/m - 1.

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);];
{dEmass ,dJmass}= {ene,J}*(massRatio - 1)

We see that the error in J in mainly due to this reduced mass effect.

Using the computed matix K and the known eigenvector we can compute also the polarization effect

dEpol = λ^2*((eigV[[state]].(matK.eigV[[state]]))/((eigV[[state]].(matN.eigV[[state]]))))(mElectronAtomicUnits/(mAtom[[Z]]));
dJpol = - dEpol*unitcm;
{dEpol,dJpol};

An updated comparison of this non relativistic (N.R.) calculation with experimental data is

ENR = ene+dEmass+dEpol;
JNR = J+dJmass+dJpol;
NumberForm[TableForm[
{{ene,dEmass,dEpol,ENR," ","a.u."},
{J,dJmass,dJpol,JNR,Jexp[[Z]],"(cm^(-1))"}},
TableHeadings-> {{"E","J"},{" variational ", " δ-mass ", " δ-polar."," tot. N.R."," exp"}}
],{10,6},
NumberPadding-> {"",""},
NumberSigns-> {"-"," "}]

The error now is only 1 part per million, i.e. is even too good (for Helium) we have not included relativistic corrections, which are expected of order (((Z α))^2).

Relativistic corrections

For S - states there are no L-S coupling but there are always relativistic corrections due to corrections to kinetic energy,  etc. These can be computed in our scheme, and have been computed. As the calculation is a bit long and has some interest on its own we have dedicated to it a notebook, [19]. We report here only the result for the ground state, if the reader wants to compute relativistic corrections to other states, she (he) is encouraged to look into notebook [19] and obtain the results to be used here.

Jrelat = -{0.304,0.544,7.83,xxx};

The peculiarity of this result is that relativistic correction is much smaller than expected.  This is due to a cancellation of terms in the computation.

The main correction to the result comes from another sourse, the radiative Lamb shift, i.e. corrections induced to the levels by vacuum fluctuations of the quantized electromagnetic field. The calculation of these corrections is beyond the range of what can be reasonably requested in an elementary course in Quantum Mechanics.  We report just the lowest order result, the source is the book [3]

JLambshift = - {.0037,1.341,-16.69,-75 - xxx}

The sum of these two corrections we call Jrad, and is listed below for atoms up to Beryllium,  the last number is just an estimate.

Jrad=-{.0037+0.304,1.341+0.544,-16.69+7.83,-75}

Our final comparison with experimental data  is then :

NumberForm[TableForm[{JNR,Jrad[[Z]],JNR+Jrad[[Z]],Jexp[[Z]]},TableDirections-> Row,
TableHeadings-> {{" J_NR ", " J_rad ", "  J_th"," J_exp"}}],{10,4}]

We have reproduced experimental results with a precision of (10^(-6))with our simple methods, and with a precision of (10^(-8))using the computed Lamb shift effect. The remaining small discrepancy is due to higher order radiative corrections, which have been computed, and to our knowledge today theoretical and experimental results have the same precision and they agree, and this is one of the great achievements of Quantum Electrodynamics.

N.B. To simplify Mathematica  instructions we have written all numbers in the same form, but we have to add that our variational computation can be believed correct only to 8 digits, i.e. J_th for Helium is really 198310.65 where the error is on the last digit. More precise results from numerical point of view will be obtained with other base functions, but the main error is not numerical, but, as we just commented, we have to add higher order relativistic corrections.

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