Relativistic corrections for Helium

Quantum Mechanics
K. Konishi, G. Paffuti

Compute the relativistic corrections to the S-state levels of Helium atom.

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

Introduction

The effective potential for relativistic corrections in S - states of Helium is, (see problems):

(Subscript[V, rel])= (α^2)((Subscript[V, 1]) + (Subscript[V, 2])+ (Subscript[V, 3])+ (Subscript[V, 4]))

with

(Subscript[V, 1]) =- (1/8)((Subscript[p, 1]^4)+(Subscript[p, 2]^4));    (Subscript[V, 2])=- (1/2)(1/(Subscript[r, 12]))( (Subscript[p, 1])(Subscript[p, 2])+ (Subscript[n, 12])((Subscript[n, 12])·(Subscript[p, 1]))(Subscript[p, 2]));
(Subscript[V, 3])= Z (π/2)((δ^((3)))((Subscript[x, 1]))+ (δ^((3)))((Subscript[x, 2])));  (Subscript[V, 4])= (Subscript[c, S]) π (δ^((3)))((Subscript[x, 12])); (Subscript[c, S])= (Piecewise[{{1, singlet}, {-(5/3), triplet}}])

(Subscript[n, 12]) is the direction (unit vector) going from electron 1 to electron 2. Z is the nuclear charge. (Subscript[V, 4 ])mean value vanishes on antisymmetrical (triplet) states.

In this notebook we will compute relativistic corrections to S-states of Helium, first for simple hydrogenoid functions for pedagogical reasons, then for more accurate trial functions obtained with the variational principle in Hylleraas coordinates.

Operators in Hylleraas coordinates

The coordinates 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 ≤ ∞

Fro a discussion of these coordinates see notebook [15]. Here we collect some formulas derived in other notebooks which help to compute mean values of operators ().

(Subscript[V, 1])

The mean value can be computed by an integration by parts as

((⟨ψ|(Subscript[V, 1])|ψ⟩ )/(⟨ψ|ψ>))=- (1/8)(1/(⟨ψ|ψ>))∫{((((Subscript[Δ, 1])ψ))^2) +((((Subscript[Δ, 2])ψ))^2) }

From now on we omit to write the normalization factor. For numerical stability it is better to transform these integrals using Schrödinger equation at lowest order, If E is the eigenvalue and V the potential term we have

⟨ψ|(1/4)(Subscript[p, 1]^4)+(1/4)(Subscript[p, 2]^4)|ψ⟩= ⟨ψ|(1/2)(Subscript[p, 1]^2)(E - V - (1/2)(Subscript[p, 2]^2))|ψ⟩+⟨ψ|(1/2)(Subscript[p, 2]^2)(E-V -(1/2)(Subscript[p, 1]^2))|ψ⟩=
 ⟨ψ|((1/2)(Subscript[p, 1]^2)+ (1/2)(Subscript[p, 2]^2))(E - V )|ψ⟩- (1/2) ⟨ψ|(Subscript[p, 1]^2)(Subscript[p, 2]^2)|ψ⟩ = ⟨ψ|(((E - V ))^2)|ψ⟩- (1/2) ⟨ψ|(Subscript[p, 1]^2)(Subscript[p, 2]^2)|ψ⟩

Then we obtain a bit more stable form:

⟨ψ|(Subscript[V, 1])|ψ⟩ = - (1/2) ⟨ψ|(((E - V ))^2)|ψ⟩ + (1/4) ⟨ψ|(Subscript[p, 1]^2)(Subscript[p, 2]^2)|ψ⟩= - (1/2)∫(((E - V ))^2)(ψ^2) +(1/4)∫((Subscript[Δ, 1])ψ)((Subscript[Δ, 2])ψ)

Clearly the two expressions are equivalent only if ψ is an exact solution of the Schrödinger equation, in approximate case they can differ, indicating the approximation made in the computation.

In Hylleraas coordinates

E-V = E + (Z/(Subscript[r, 1]))+ (Z/(Subscript[r, 2]))- (1/(Subscript[r, 12]))= E + 2Z((1/(s+t))+(1/(s-t))) - (1/u);

(Subscript[Δ, 1])ψ = (4/(s+t))((∂_s) ψ+(∂_t) ψ)+(2/u)(∂_u) ψ+((∂_(s,s))ψ + (∂_(t,t))ψ+(∂_(u,u))ψ)+2 (∂_(s,t))ψ+2((s t + (u^2))/(u(s+t)))((∂_(s,u))ψ+(∂_(t,u))ψ);

(Subscript[Δ, 2])ψ = (4/(s-t))((∂_s) ψ-(∂_t) ψ)+(2/u)(∂_u) ψ+((∂_(s,s))ψ + (∂_(t,t))ψ+(∂_(u,u))ψ)-2 (∂_(s,t))ψ+2((-s t + (u^2))/(u(s-t)))((∂_(s,u))ψ-(∂_(t,u))ψ);

In the computation of integrals one must add the Jacobian between Cartesian and Hylleraas coordinates, i.e. u((s^2)- (t^2)).

(Subscript[V, 2])

This is the most complicated term. In notebook [] we have computed the expression for the derivative part of this operator:

-( (Subscript[p, 1])(Subscript[p, 2])+ (Subscript[n, 12])((Subscript[n, 12])·(Subscript[p, 1]))(Subscript[p, 2]))ψ = -((((t^2) (u^2)-3 (u^4)+(s^2) ((t^2)+(u^2))) (∂_(t,t))ψ)/(((s^2)-(t^2)) (u^2)))-((2 (∂_u) ψ)/u)+((4 t (-(s^2)+(u^2)) (∂_(u,t))ψ)/(((s^2)-(t^2)) u))-2 (∂_(u,u))ψ+((4 s ((t^2)-(u^2)) (∂_(s,u))ψ)/(((s^2)-(t^2)) u))+((((u^2) ((t^2)-3 (u^2))+(s^2) ((t^2)+(u^2)))(∂_(s,s))ψ)/(((s^2)-(t^2)) (u^2)))

In matrix elements this must be multiplied by the prefactor 1/u and by the Jacobian.

(Subscript[V, 4])

For (Subscript[x, 12])= 0 we have u = 0 and consequently t=0 and s = 2 (Subscript[r, 1]). The only problem is to fix the correct normalization. The overall factor for the jacobian in 2(π^2), i.e. for invariant functions

(d^3)(Subscript[x, 1]) (d^3)(Subscript[x, 2]) → 2(π^2) u((s^2)-(t^2))ds du dt;  0≤ t ≤ u ≤ s ≤ ∞

This can be checked in the following way, take the integral of Exp[-(Subscript[r, 1])- (Subscript[r, 2])]: In Cartesian coordinates evidently the integral is

(((4π Integrate[Exp[-r1](r1^2),{r1,0,∞}]))^2)

In Hylleraas coordinates

2 (π^2)Integrate[u((s^2)-(t^2))Exp[-s],{s,0,∞ },{u,0,s},{t,0,u}]

which proves the statement. In Cartesian coordinates, integrating the δ function

((⟨ψ|(Subscript[V, 4])|ψ⟩)/(⟨ψ|ψ⟩))= (π/(⟨ψ|ψ⟩)) ∫(ψ^2)[(Subscript[x, 1]),(Subscript[x, 1])](Subscript[𝕕x, 1]) =((4 (π^2))/(⟨ψ|ψ⟩))(∫_0^∞)(Subscript[𝕕r, 1])(Subscript[r, 1]^2)(ψ^2)[(Subscript[r, 1]),(Subscript[r, 1])]=((4 (π^2))/(⟨ψ|ψ⟩))(1/8)(∫_0^∞)𝕕s (s^2)(ψ^2)[s,0,0]

The norm N, in usual calculations, is defined without the factor 2(π^2), i.e.

⟨ψ|ψ⟩ = 2(π^2) N

Then

((⟨ψ|(Subscript[V, 4])|ψ⟩)/(⟨ψ|ψ⟩))= (1/N)(1/4)(∫_0^∞)𝕕s (s^2)(ψ^2)[s,0,0]

(Subscript[V, 3])

One proceeds in the same way for (Subscript[V, 3]). For example the term with (Subscript[x, 2])= 0 (second δ function) is, in Cartesian coordinates

Z(π/2)∫(Subscript[𝕕x, 1])(ψ^2)[(Subscript[x, 1]),0]= Z 2(π^2) (∫_0^∞)(Subscript[𝕕r, 1])(Subscript[r, 1]^2)(ψ^2)[(Subscript[r, 1])]

In Hylleraas coordinates (Subscript[x, 2])= 0 imply s=t=u= (Subscript[r, 1])than previus result is written as

Z 2(π^2) (∫_0^∞)𝕕s (s^2)(ψ^2)[s,s,s]

Dividing by the norm and taking account both terms :

((⟨ψ|(Subscript[V, 3])|ψ⟩)/(⟨ψ|ψ⟩)) = ((2 Z)/N) (∫_0^∞)𝕕s (s^2)(ψ^2)[s,s,s]

Matrices

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[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;];
λ]

Operators

Delta1[s_,u_,t_]:=(4/(s+t))((∂_s) ψ[s,u,t]+(∂_t) ψ[s,u,t])+(2/u)(∂_u) ψ[s,u,t]+
((∂_(s,s))ψ[s,u,t] + (∂_(t,t))ψ[s,u,t]+(∂_(u,u))ψ[s,u,t])+2 (∂_(s,t))ψ[s,u,t]+2((s t + (u^2))/(u(s+t)))((∂_(s,u))ψ[s,u,t]+(∂_(t,u))ψ[s,u,t])

Delta2[s_,u_,t_]:=  (4/(s-t))((∂_s) ψ[s,u,t]-(∂_t) ψ[s,u,t])+(2/u)(∂_u) ψ[s,u,t]+
((∂_(s,s))ψ[s,u,t] + (∂_(t,t))ψ[s,u,t]+(∂_(u,u))ψ[s,u,t])-2 (∂_(s,t))ψ[s,u,t]+2((-s t + (u^2))/(u(s-t)))((∂_(s,u))ψ[s,u,t]-(∂_(t,u))ψ[s,u,t])

derivOp2[s_,u_,t_]:= -((((t^2) (u^2)-3 (u^4)+(s^2) ((t^2)+(u^2))) (∂_(t,t))ψ[s,u,t])/(((s^2)-(t^2)) (u^2)))-((2 (∂_u) ψ[s,u,t])/u)+((4 t (-(s^2)+(u^2)) (∂_(u,t))ψ[s,u,t])/(((s^2)-(t^2)) u))-2 (∂_(u,u))ψ[s,u,t]+((4 s ((t^2)-(u^2)) (∂_(s,u))ψ[s,u,t])/(((s^2)-(t^2)) u))+((((u^2) ((t^2)-3 (u^2))+(s^2) ((t^2)+(u^2)))(∂_(s,s))ψ[s,u,t])/(((s^2)-(t^2)) (u^2)))

Computation

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

state = 1;
shT = 0; 
shT2 = shT/2;

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

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

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

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

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

The reader can verify that the order of indices in Table is correct, i.e. first is incremented i1, then j1 then k1, as in the definition of index in the function matrices.

ψ[s_,u_,t_]:= Exp[-λ s/2]*eigV[[state]].Flatten[Table[(((s^i1)(t^(2j1+shT))(u^k1) (λ^(i1+2j1+k1+shT)))/((i1+2*j1+k1)!)),
{k1,0,maxu},{j1,0,maxt},{i1,0,maxs}]];

Off[NIntegrate::"slwcon"]

The norm of the wave function is computed with

normψ = NIntegrate[u((s^2)-(t^2))*ψ[s,u,t]^2,{s,0,∞},{u,0,s},{t,0,u}];

Let us compute the mean values of different operators. For (Subscript[V, 1])it is convenient to perform two integrals and expand the factor (((E-V))^2)in order to cancel some of the apparent singularities. The reader should observe that in all the integrals the jacobian factor is inserted.

meanV1=Module[{ris1,ris2},
ris1= (1/(normψ))NIntegrate[((-1)/2)u((s^2)-(t^2)) (((ene+2Z(((2s)/((s^2)-(t^2))))-(1/u)))^2)((ψ[s,u,t])^2),{s,0,∞},{u,0,s},{t,0,u},Method->{"MultiDimensionalRule","Generators"-> 9}];
ris2= (1/(normψ))NIntegrate[
(1/4)u((s^2)-(t^2))Delta1[s,u,t]Delta2[s,u,t],{s,0,∞},{u,0,s},{t,0,u},Method->{"MultiDimensionalRule","Generators"-> 9}];
ris1+ris2];

For (Subscript[V, 2]) we insert also the prefactor 1/u

meanV2 = (1/(normψ))NIntegrate[
((s^2)-(t^2))*(1/2)ψ[s,u,t]*derivOp2[s,u,t],{s,0,∞},{u,0,s},{t,0,u},Method->{"MultiDimensionalRule","Generators"->9}];

The integrals (Subscript[V, 3]) and (Subscript[V, 4])are easy:

meanV3= ((2 Z)/(normψ))NIntegrate[(s^2)((ψ[s,s,s])^2),{s,0,∞}];

meanV4=If[shT== 1,0.,(1/(4normψ))NIntegrate[(s^2)((ψ[s,u,t])^2)/.{u->0,t-> 0},{s,0,∞}]];

The partial results are

{meanV1,meanV2,meanV3,meanV4}

On[NIntegrate::"slwcon"]

The reader has problably noticed that the computation is quite long. At the end of the notebook we have collected in a function all previous steps and we have added the options PrecisionGoal → 5 for integrals. The reader can try on his computer several options, as Accuracy → 6 etc. to compute the integrals. Another possibility is to change the dimension of the basis.

Tables

Here we report the results for the first three states of paraHelium (1(s^2),1s2s, 1s3s) and for the first two of orthoHelium (1s2s, 1s3s).

tableδErelativisticS0 = {{-13.5222,-0.139128,11.373,0.33491},{-10.2746,-0.00927233,8.21768,0.0273481},{-10.2746,-0.00927233,8.21768,0.0273481}};

Total/@tableδErelativisticS0

tableδErelativisticS1 = {{-10.45874,-0.00164,8.29546,0.},{-10.02746,-0.00058,7.94098,0.}}

Energy corrections

We convert now our numbers to atomic units by inserting the missing factor (α^2)

alpha = 1/137.035999765369;

TableForm[{Append[(alpha^2)*Total/@tableδErelativisticS0,"  a.u."],
Append[Join[{""},(alpha^2)*Total/@tableδErelativisticS1],"  a.u."]},
TableHeadings-> {{"S=0","S=1"},{"    1(s^2) ", "     1s2s ", "     1s3s "}}
]

To compute relativistic corrections to ionization potential we have to take into account also the analogous corrections to the ion, i.e. hydrogenic atom in 1s state. These are known to be

δEion = - (Z^4)(alpha^2)/8;

The correction to the ionization potential will be

δJ = δEion- δV

In units (cm^(-1))

unitRy = 109737.3156; unitcm = 2*unitRy;

TableForm[{Append[(δEion-(alpha^2)*Total/@tableδErelativisticS0)unitcm ,"  (cm^(-1))"],
Append[Join[{""},(δEion-(alpha^2)*Total/@tableδErelativisticS1)unitcm],"  (cm^(-1))"]},
TableHeadings-> {{"S=0","S=1"},{"    1(s^2) ", "  1s2s ", "  1s3s "}}
]

These results can be used as corrections in the notebook [17].

Function

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

relativisticδE[state_]:=
Module[{ris1,ris2,λstart},
λstart = N[2*(Z-5/16)];
λ = findParameter[state,λstart];
{eigH,eigV}=Re[Eigensystem[{λ^2 matM-λ matL,matN}]];
order = Ordering[eigH];
eigH = eigH[[order]];
eigV = eigV[[order]];
ene = eigH[[state]];
ψ[s_,u_,t_]:= Exp[-λ s/2]*eigV[[state]].Flatten[Table[(((s^i1)(t^(2j1+shT))(u^k1) (λ^(i1+2j1+k1+shT)))/((i1+2*j1+k1)!)),
{k1,0,maxu},{j1,0,maxt},{i1,0,maxs}]];normψ = NIntegrate[u((s^2)-(t^2))*ψ[s,u,t]^2,{s,0,∞},{u,0,s},{t,0,u},
Method->{"MultiDimensionalRule","Generators"-> 9},PrecisionGoal-> 5];
Off[NIntegrate::"slwcon"];
ris1= (1/(normψ))NIntegrate[((-1)/2)u((s^2)-(t^2)) (((ene+2Z(((2s)/((s^2)-(t^2))))-(1/u)))^2)((ψ[s,u,t])^2),{s,0,∞},{u,0,s},{t,0,u},Method->{"MultiDimensionalRule","Generators"-> 9},
PrecisionGoal-> 5];
ris2= (1/(normψ))NIntegrate[
(1/4)u((s^2)-(t^2))Delta1[s,u,t]Delta2[s,u,t],{s,0,∞},{u,0,s},{t,0,u},Method->{"MultiDimensionalRule","Generators"-> 9},PrecisionGoal-> 5];
meanV1=ris1+ris2;
meanV2 = (1/(normψ))NIntegrate[
((s^2)-(t^2))*(1/2)ψ[s,u,t]*derivOp2[s,u,t],{s,0,∞},{u,0,s},{t,0,u},Method->{"MultiDimensionalRule","Generators"->9},PrecisionGoal-> 5];
meanV3= ((2 Z)/(normψ))NIntegrate[(s^2)((ψ[s,s,s])^2),{s,0,∞}];
meanV4=(1/(4normψ))NIntegrate[(s^2)((ψ[s,u,t])^2)/.{u->0,t-> 0},{s,0,∞}];
{meanV1,meanV2,meanV3,meanV4}]

Example of use of the function ( is supposed that the reader has fixed the parameter shT)

(withfunct=relativisticδE[#1]&/@{1,2,3})//Timing