Relativistic corrections-Pekeris variables

Quantum Mechanics
 K. Konishi, G. Paffuti

Compute the relativistic corrections to Helium atom using triangular coordinates.

Unprotect[dimIndR,indR,segnoParita,αKP,βKP,γKP,formamatM,formamatL,formamatN,formamatK];
Unprotect[x,y,z,u,s,t,i1,i2,i3,j1,j2,j3,fromInd];

ClearAll["Global`*"];
Off[General::spell, General::spell1,InterpolatingFunction::"dmval"];
Protect[x,y,z,u,s,t,i1,i2,i3,j1,j2,j3,fromInd];

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 in ± 1 depending on exchange symmetry of the wave function, indR and dimIndR are a list of indexes 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 extend results of notebook [20] to the computation of relativistic corrections, as we have done in notebook [19] for Hylleraas variables. Final results, with added the Lamb shift effect, are compared to experimental data.

Short review on triangular valiables

The triangular variables are defined by

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

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 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 chosen as λ (2,1,1), i.e. we allow only an overall scale factor.

The basis functions are

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

The variational problem

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

The relevant matrix elements are:

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.

To use these matrices two function call are required :

The section Description of functions give some details on numerical aspects.

The relativistic corrections

In notebook [19] we have seen that relativistic corrections for S-states are described by four effective potentials

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

Operators in triangular variables

(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 triangular coordinates, see notebook [15]

E-V = E + (Z/(Subscript[r, 1]))+ (Z/(Subscript[r, 2]))- (1/(Subscript[r, 12]))= E + 2Z((1/(x+y))+(1/(x+z)))-(2/(y+z));

(Subscript[Δ, 1])ψ = ((4 x z (∂_(z,z))ψ)/((x+y) (y+z)))+((4 y (x+y+z) (∂_(y,y))ψ)/((x+y) (y+z)))+4 ((((x-z)(∂_z)ψ+(x+2 y+z) (∂_y)ψ+(-x+z) (∂_x)ψ))/((x+y) (y+z)))-((8 x z (∂_(x,z))ψ)/((x+y) (y+z)))+((4 x z (∂_(x,x))ψ)/((x+y) (y+z)));

(Subscript[Δ, 2])ψ = ((4 z (x+y+z) (∂_(z,z))ψ)/((x+z) (y+z)))+((4 x y  (∂_(y,y))ψ)/((x+z) (y+z)))+((4 ((x+y+2 z) (∂_z)ψ+(x-y) ((∂_y)ψ-(∂_x)ψ)))/((x+z) (y+z)))-((8 x y (∂_(x,y))ψ)/((x+z) (y+z)))+((4 x y (∂_(x,x))ψ)/((x+z) (y+z)));

In the computation of integrals one must add the Jacobian between Cartesian and Pekeris coordinates, i.e.

dJ = (1/2)(x+y)(x+z)(y+z).

(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]))ψ = -((4 x z (x+y+z) (y+2 z) (∂_(z,z))ψ)/((x+y) (x+z) (((y+z))^2)))-((8 y z ((x^2)+x (y+z)+(((y+z))^2)) (∂_(y,z))ψ)/((x+y) (x+z) (((y+z))^2)))-((4 x y (x+y+z) (2 y+z) (∂_(y,y))ψ)/((x+y) (x+z) (((y+z))^2)))-((4 ((∂_z)ψ+(∂_y)ψ-(∂_x)ψ))/(y+z))+((8 x z (∂_(x,z))ψ)/((x+y) (y+z)))+((8 x y (∂_(x,y))ψ)/((x+z) (y+z)))+((4 x y z (-x+y+z) (∂_(x,x))ψ)/((x+y) (x+z) (((y+z))^2)))

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

(1/2)(x+y)(x+z)(y+z)(2/(y+z)) = (x+y)(x+z)

(Subscript[V, 4])

For (Subscript[r, 12]) = 0 we have y = z = 0 → r1 = x/2 (x, y, z must be non negative). Then the operator is reduced to an x integration. To fix the correct jacobian consider 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)  = 64 (π^2)

The conversion factor is (π^2)/2 in fact

((Pi^2)/2)Integrate[(1/2)(x+y)(x+z)(y+z) Exp[-x - y/2 - z/2],
{x,0,∞},{y,0,∞},{z,0,∞}]

In particular, in terms of our quadratic form N:

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

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^∞)𝕕x (x^2)(ψ^2)[x,0,0] = (1/N)(∫_0^∞)𝕕x (x^2)(ψ^2)[x,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]= 2Z (π^2) (∫_0^∞)(Subscript[𝕕r, 1])(Subscript[r, 1]^2)(ψ^2)[(Subscript[r, 1])]

In Pekeris coordinates (Subscript[x, 2])= 0 imply x=z=0; r1= y/2 while (Subscript[x, 1])=0imply  x=y=0; r2=z/2  than the previus result for both elecrons is written as

Z 2(π^2) ((1/8)(∫_0^∞)𝕕y (y^2)(ψ^2)[0,y,0]+(1/8)(∫_0^∞)𝕕z (z^2)(ψ^2)[0,0,z])

Dividing by the norm and taking account both terms :

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

Description of functions

Functions for 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. With respect the function used in notebook [20] we have added the computation of a global vector, fromInd which allow to pass form index to components i1,i2,i3.

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

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 construct 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}];

Functions for relativistic corrections

We just copy the above definitions

Delta1[x_,y_,z_]:=((4 x z (∂_(z,z))ψ[x,y,z])/((x+y) (y+z)))+((4 y (x+y+z) (∂_(y,y))ψ[x,y,z])/((x+y) (y+z)))+(1/((x+y) (y+z)))4 ((x-z)(∂_z)ψ[x,y,z]+(x+2 y+z) (∂_y)ψ[x,y,z]+(-x+z) (∂_x)ψ[x,y,z])-((8 x z (∂_(x,z))ψ[x,y,z])/((x+y) (y+z)))+((4 x z (∂_(x,x))ψ[x,y,z])/((x+y) (y+z)))

Delta2[x_,y_,z_]:=((4 z (x+y+z) (∂_(z,z))ψ[x,y,z])/((x+z) (y+z)))+((4 x y  (∂_(y,y))ψ[x,y,z])/((x+z) (y+z)))+((4 ((x+y+2 z) (∂_z)ψ[x,y,z]+(x-y) ((∂_y)ψ[x,y,z]-(∂_x)ψ[x,y,z])))/((x+z) (y+z)))-((8 x y (∂_(x,y))ψ[x,y,z])/((x+z) (y+z)))+((4 x y (∂_(x,x))ψ[x,y,z])/((x+z) (y+z)))

derivOp2[x_,y_,z_]:= -((4 x z (x+y+z) (y+2 z) (∂_(z,z))ψ[x,y,z])/((x+y) (x+z) (((y+z))^2)))-((8 y z ((x^2)+x (y+z)+(((y+z))^2)) (∂_(y,z))ψ[x,y,z])/((x+y) (x+z) (((y+z))^2)))-((4 x y (x+y+z) (2 y+z) (∂_(y,y))ψ[x,y,z])/((x+y) (x+z) (((y+z))^2)))-((4 ((∂_z)ψ[x,y,z]+(∂_y)ψ[x,y,z]-(∂_x)ψ[x,y,z]))/(y+z))+((8 x z (∂_(x,z))ψ[x,y,z])/((x+y) (y+z)))+((8 x y (∂_(x,y))ψ[x,y,z])/((x+z) (y+z)))+((4 x y z (-x+y+z) (∂_(x,x))ψ[x,y,z])/((x+y) (x+z) (((y+z))^2)))

dJ[x_,y_,z_]:= (1/2)(x+y)(x+z)(y+z)

Utilities and units

We collect here some functions used several times in other notebooks.

At the end we have collected the various conversion factors and masses for atomic nuclei.

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

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

calcEneψNR[state_,λ_]:= Module[{de1,de2,ene,eigH,eigV},
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[[state]];
(* these are mass corrections *)
de1=ene*(massRatio[[Z]] - 1);
de2 = λ^2*((eigV[[state]].(matK.eigV[[state]]))/((eigV[[state]].(matN.eigV[[state]]))))*(mElectronAtomicUnits/(mAtom[[Z]]));
{{de1,de2},ene,eigV[[state]]}]

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

Computation

Orders of magnitude and errors

The energy to be computed are of the order ε ∼ (Z^2)/2 in atomic units. For Helium ground state ε ∼2.9 a.u. With the parameter ω = 15 used below energies for ground state and first excited state have a numerical precision of (10^(-12)) and (10^(-10)) respectively. The first order mass corrections (order m/M ϵ) have the same precision.  Relativistic corrections are computed in the form δV (α^2), where α is the fine structure constant. δV are computed with an accuracy of 5-6 digits. A pessimistic estimate of numerical error is then (for the ground state of Helium)

Row[{"(Subscript[δ, num])E ∼ ",(10^(-5))(alpha^2)," au  =",  (10^(-5))(alpha^2) unitcm , " (cm^(-1))"}]

This is much less than  corrections di to neglected physical effects.  We put "by hand" the correction due to Lamb shift, (some unit in (cm^(-1))) but we neglect order (((m/M))^2)corrections, higher order relativistic corrections, of order ε (((Zα))^4)etc. These effects contribute, as an order of magnitude for Helium:

ε0 = 2; TableForm[{{"(m/M)^2ε ∼ ",δϵ1=((((mElectronAtomicUnits/4)))^2)ε0," au  =",  δϵ1* unitcm , " (cm^(-1))"},
{"(Z α)^4ε ∼ ",δε2 = (((2 alpha))^4)ε0 ," au  =",  δε2* unitcm , " (cm^(-1))"} }]

This means that we have to expect agreement with experimental data up to 0.02 - 0.03 (cm^(-1))and that our precision in numerical calculation is sufficient for this level of comparison.

The integrals appearing in the computation of relativistic corrections can be computed analytically, expanding in powers products of Laguerre polynomials, but the procedure would be too long to present here and would not improve our precision, so we limit to a slow and relatively inefficient computation using explicit integrals computed with NIntegrate.

The atom

We take as an example a computation of the ground state of the Helium atom, the reader can generalize to other states and other two electrons ions

Z=2;numMaxEig=5;

This compute algebric forms for matrices:

computeAlgebricMatrices[Z];

The calculation is slow, if necessary the reader can choose a smaller value for ω.

We fix the exchange parity and the parameter ω, and compute matrices:

parity = 0; ω = 15;

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

Here we fix the state and compute the relative optimized scale parameter :

state = 1;

λ = findParameter[state];
{{dEnemass1,dEpol},ene,cfψ} =calcEneψNR[state,λ];
ene

We write the wave function using the list of coefficient cfψ just coputed and the base functions. The indexes of relevant Laguerre polynomial are extracted from the global vector fromInd, computed within getMatrices with a call to funzIndice.

ψ[x_,y_,z_]:= Exp[- λ x - λ y/2 - λ z/2]*
Sum[cfψ[[m]]*LaguerreL[fromInd[[m,1]],2 λ x]*(
LaguerreL[fromInd[[m,2]],λ y]LaguerreL[fromInd[[m,3]],λ z]+segnoParita[parity]*LaguerreL[fromInd[[m,2]],λ z]LaguerreL[fromInd[[m,3]],λ y]),{m,1,dimIndR}]

Now we start to compute the four terms for relativistic corrections. We leave in the inputs below the function Timing to give an idea of execution time to the reader:

(normψ = NIntegrate[dJ[x,y,z]((ψ[x,y,z])^2),{x,0,∞},{y,0,∞},{z,0,∞},Method->{"MultiDimensionalRule","Generators"-> 9}])//Timing

(meanV1=Module[{ris1,ris2},
ris2=(1/(4normψ))NIntegrate[dJ[x,y,z] ( Delta1[x,y,z]*Delta2[x,y,z]),
{x,0,∞},{y,0,∞},{z,0,∞},Method->{"MultiDimensionalRule","Generators"-> 9},AccuracyGoal-> 6];
ris1= ((- 1)/(2normψ))NIntegrate[dJ[x,y,z] (((ene+2Z((1/(x+y))+(1/(x+z)))-(2/(y+z))))^2)((ψ[x,y,z])^2),{x,0,∞},{y,0,∞},{z,0,∞},Method->{"MultiDimensionalRule","Generators"-> 9},AccuracyGoal-> 6];
ris1+ris2])//Timing

(meanV2 = (1/(normψ))NIntegrate[
(x+y)(x+z)*(1/2)ψ[x,y,z]*derivOp2[x,y,z],{x,0,∞},{y,0,∞},{z,0,∞},
AccuracyGoal-> 6,Method->{"MultiDimensionalRule","Generators"->9}])//Timing

meanV3=(Z/(2normψ))(NIntegrate[(y^2) ((ψ[0,y,0])^2),{y,0,∞}]+NIntegrate[(z^2) ((ψ[0,0,z])^2),{z,0,∞}])

meanV4 =  (1/(normψ))NIntegrate[ (x^2) ((ψ[x,0,0])^2),{x,0,∞}]

{meanV1,meanV2,meanV3,meanV4}

Comparison with experimental data

Jexp={6100,198310.6691,610079., 1241242.};
JLambshift = - {.0037,1.341,-16.69,-75 };

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

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

{dEmass ,dJmass}= {ene,J}*(massRatio[[Z]] - 1);
dJpol = - dEpol*unitcm;

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,7},
NumberPadding-> {"",""},
NumberSigns-> {"-"," "}]

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

δJrel=(δEion-(alpha^2)*Total[{meanV1,meanV2,meanV3,meanV4}])unitcm

Jrad = δJrel + JLambshift[[Z]];

NumberForm[TableForm[{JNR,JNR+δJrel,JLambshift[[Z]],JNR+Jrad,Jexp[[Z]]},TableDirections-> Row,
TableHeadings-> {{" J_NR ", " J_th0 "," J_LS ", "  J_th"," J_exp"}}],{10,4}]

The agreement with experimental data is within our theoretical approximations and the reader can check that a computation with ω = 20 gives exactly the same results, i.e. the computation is stable from numerical point of view.