Coordinate systems for Helium

Quantum Mechanics
 K. Konishi, G. Paffuti

Study the Hylleraas and Pekeris coordinate systems for Helium.

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

Introduction

The variables

In a three body problem, like in the Helium atom, the positions of the three particles (nucleus + 2 electrons), in the center of mass, is fixed by 3 global rotations and by three lengths: (Subscript[r, 1]), (Subscript[r, 2]),(Subscript[r, 12]), the distances between electrons and nucleus and the distance between the two electrons respectively. We assume that the center of mass coincide with the position of the nucleus.  The corrections of order m/M (= electron mass/nucleus mass) due to such a simplification will be considered separately.

The geometry of the problem is clear: the invariant object constructed with the three particles is the triangle with vertices in their positions, the triangle is uniquely specified by the lengths of its three edges, (Subscript[r, 1]), (Subscript[r, 2]),(Subscript[r, 12]). The distance (Subscript[r, 12])can also be expressed through the angle θ formed by the two lines connecting nucleus to electrons:

(Subscript[r, 12]^2) = ((((Subscript[x, 1]) - (Subscript[x, 2])))^2)= (Subscript[r, 1]^2) + (Subscript[r, 2]^2) - 2 (Subscript[r, 1]) (Subscript[r, 2]) Cos[θ]

A triangle defines a reference frame, in which (x,y) plane coincides with the plane of the triangle, with the x axis, for example, being directed along the line nucleus-electron 1; the z axis is the axis perpendicular to this plane. The three global rotations specify the orientation of this frame with respect to a fixed reference frame, i.e. we need three Euler angle to pass from "comoving frame" to the fixed frame. (The Mathematica code for the figure is at the end of this section).

A wave function of the system will be of the form ψ[(Subscript[x, 1]), (Subscript[x, 2])] in Cartesian coordinates, and ψ[(Subscript[r, 1]), (Subscript[r, 2]),(Subscript[r, 12]),α, β, γ ] as a function of intrinsic coordinates and of the  Euler angles. In this notebook we will be interested only in S-states, which do not depend on orientation of the reference frame.  The dependence on the Euler angles will drop. Our problem is to write explicit formulae for this coordinate change and in particular to write quadratic forms for the kinetic and potential energies in the new coordinate system. The quadratic forms we are interested in are

M =(1/2) ∫(((((Subscript[∇, 1])Ψ))^2)+ ((((Subscript[∇, 2])Ψ))^2)) ;    L = ∫((Z/r1)+(Z/r2) -(1/r12) )(Ψ^2) ;    N = ∫(Ψ^2);    K = ∫(((Subscript[∇, 1])Ψ) ((Subscript[∇, 2])Ψ)) ;

The last form, K, is needed to compute m/M corrections to the spectrum.

The reference frame

Let us call ex, ey, ez the three basis vectors on the fixed frame and fx, fy, fz the corresponding basis in the moving frame. The two sets of vectors are related by a rotation matrix ℛ which can be expressed through the Euler angles, see Chapter 4 of the main text.

Here we define the three basic rotations around the x,y,z axes:

Rx[α_]:={{1,0,0},{0,Cos[α],Sin[α]},{0,-Sin[α],Cos[α]}};
Ry[α_]:= {{Cos[α],0,-Sin[α]},{0,1,0},{Sin[α],0,Cos[α]}};
Rz[α_]:={{Cos[α],Sin[α],0},{-Sin[α],Cos[α],0},{0,0,1}};

The rotation matrix is defined by Rz[γ].Ry[β].Rz[α].

myRotationMatrix3D[α_,β_,γ_]:=Rz[γ].Ry[β].Rz[α];
rotEul =Transpose[ myRotationMatrix3D[α,β,γ]];

The basis vectors transform as the inverse matrix (or transposed matrix for orthogonal matrices):

ex = {1,0,0}; ey = {0,1,0}; ez = {0,0,1};
fx = rotEul.ex; fy = rotEul.ey; fz = rotEul.ez;

Vectors fx, fy, fz are the  rows of the following table

TableForm[{fx,fy,fz}]

The coordinates p1 and p2 of the two electrons are now simple expressed by:

p1 = Simplify[r1*Cos[(θ/2)] fx - r1* Sin[(θ/2)]fy];  p2 = Simplify[r2*Cos[(θ/2)] fx + r2*Sin[(θ/2)]fy];

Instructions for the figure

pt2 = {0, 0, 0};
pt1 = {Sin[π/4] Cos[π/3], Sin[π/4] Sin[π/3], Cos[π/4]};
pt3 = {Sin[π/4] Cos[2 π/3], Sin[π/4] Sin[2 π/3], Cos[π/4]};
Graphics3D[{
  {Opacity[0.6], Sphere[pt3, r], Sphere[pt1, r]}, 
  {Opacity[0.6], Green, Sphere[pt2, 1.5 r]},
  {Blue, Opacity[0.4], Polygon[{pt1, pt2, pt3}]},
  Line[{pt2, pt2 + {1, 0, 0}}],
  Line[{pt2, pt2 + {0, 0, 1}}],
  Line[{pt2, pt2 + {0, 1, 0}}]
  }, ViewPoint -> {3.37, -0.9, 2.1},
 Boxed -> False]

Change of variables

Hylleraas variables

Our aim is to change the variables from Cartesian coordinates to Hylleraas coordinates, defined as

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

Let as name q1 and q2 the two positions in Cartesian coordinates

q1={x1,y1,z1};
q2 = {x2,y2,z2};

A trivial computation shows that

(Subscript[∇, 1])f[s,u,t]=((Subscript[x, 1])/r1)(∂_s) f + ((Subscript[x, 1])/(Subscript[r, 1]))(∂_t) f + (((Subscript[x, 1])-(Subscript[x, 2]))/(Subscript[r, 12]))(∂_u) f;
(Subscript[∇, 2])f[s,u,t]=((Subscript[x, 2])/(Subscript[r, 2]))(∂_s) f - ((Subscript[x, 2])/(Subscript[r, 2]))(∂_t) f - (((Subscript[x, 1])-(Subscript[x, 2]))/(Subscript[r, 12]))(∂_u) f;

Our differential operators will act on functions of s, t, u and, at most, simple powers of Cartesian coordinates. Then on a function f the operator (∂_i) acts as:

grad1[f_]:= Table[((∂_(q1[[i]])) f) ,{i,1,3}]+( (q1/r1)(∂_s) f)+( (q1/r1)(∂_t) f)+( ((q1-q2)/u)(∂_u) f);
grad2[f_]:= Table[((∂_(q2[[i]])) f) ,{i,1,3}]+( (q2/r2)(∂_s) f)-( (q2/r2)(∂_t) f)-( ((q1-q2)/u)(∂_u) f);

The suffixes 1 and 2 referring to the two particles. Expressing r1 and r2 as a function of s and t:

grad1[f_]:= Table[((∂_(q1[[i]])) f) ,{i,1,3}]+( ((2 q1)/((s+t)))(∂_s) f)+( ((2q1)/((s+t)))(∂_t) f)+( ((q1-q2)/u)(∂_u) f);
grad2[f_]:= Table[((∂_(q2[[i]])) f) ,{i,1,3}]+( ((2q2)/((s-t)))(∂_s) f)-( ((2q2)/((s-t)))(∂_t) f)-( ((q1-q2)/u)(∂_u) f);

The rest of the calculation is simply an application of the usual substitution rules.

Cartesian coordinates will be expressed through intrinsic coordinates and the Euler angles :

ruleQ = Join[Thread[Rule[q1,p1]],Thread[Rule[q2,p2]]];

We have to pass from the radial coordinates r1, r2, θ to Hylleraas coordinates. We use the trivial relations

(u^2)= (Subscript[r, 1]^2)+(Subscript[r, 2]^2) -2(Subscript[r, 1])(Subscript[r, 2])Cos[θ];        (Subscript[r, 1]) = ((s+t)/2);       (Subscript[r, 2])=(( s-t)/2)

ruleH = {Cos[θ]-> (r1^2 + r2^2 - u^2)/(2 r1 r2), 
Sin[θ]-> (Sqrt[1 - ((r1^2 + r2^2 - u^2)/(2 r1 r2))^2]),
r1 -> (s+t)/2, r2-> (s-t)/2};

Note on the function Thread

The function Rule[a,b] is just the explicit form of the substitution a -> b . The function Thread apply this rule to each couple of elements of the two lists. As an example:

temp1 = {a1,b1}; temp2 = {a2,b2};
Thread[Rule[temp1,temp2]]

Clear[a1,b1,a2,b2,temp1,temp2]

Quadratic forms

First we have to write the jacobian from Cartesian coordinates to Hylleraas coordinates. It is easy to check (see also below) that

Jacobian = u ((s^2)-(t^2)) Sin[β] ( dα dβ dγ ds dt du ).

We are interested only in quantities invariant under rotations, so we neglect α, β and γ and write simply

dJ = u (s^2-t^2);

We will use repeatedly the following function to collect powers of derivatives :

collectDerivatives[espr_]:=Collect[espr,{((Derivative[1,0,0][ψ][_,_,_])^2),((Derivative[0,0,1][ψ][_,_,_])^2),((Derivative[0,1,0][ψ][_,_,_])^2)}]

The only thing we have now to do is to write the quadratic forms we are interested in and apply the substitution rules :

quadM= Simplify[TrigReduce[(1/2)dJ*(grad1[ψ[s,u,t]].grad1[ψ[s,u,t]]+
grad2[ψ[s,u,t]].grad2[ψ[s,u,t]])/.ruleQ]//.ruleH,Assumptions-> {0<t<u<s}];
quadM=collectDerivatives[quadM]

quadK= Simplify[TrigReduce[dJ*(grad1[ψ[s,u,t]].grad2[ψ[s,u,t]])/.ruleQ]//.ruleH,Assumptions-> {0<t<u<s}];
quadK=collectDerivatives[quadK]

quadL = Simplify[TrigReduce[(dJ*(-1/u+ Z/r1 + Z/r2)((ψ[s,u,t])^2))/.ruleQ]//.ruleH,Assumptions-> {0<t<u<s}]

quadN = Simplify[TrigReduce[(dJ*((ψ[s,u,t])^2))/.ruleQ]//.ruleH,Assumptions-> {0<t<u<s}]

All dependence on the Euler angles disappear, as expected.

Text form

The quadratic forms written above are used in all notebooks on the Helium. We note for completeness that the quadratic form L is defined as -V where V is the potential, the sign is easily checked looking at Z dependent terms. The variational principle is written as

δ( M + V - λ N)=0    ⟺   δ( M - L - λ N)=0

To help the reader we rewrite above expressions in simpler notation using

shortForm={Derivative[0,0,1][ψ][s,u,t]-> "(((∂ψ)/(∂t)))",Derivative[0,1,0][ψ][s,u,t]-> "(((∂ψ)/(∂u)))" ,
Derivative[1,0,0][ψ][s,u,t]-> "(((∂ψ)/(∂s)))"};

The results have to be compared with

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) ; 
K=(∫_0^∞)ds (∫_0^s)du(∫_0^u)dt (-2s ((u^2)-(t^2)) (((∂Ψ)/(∂u)) ((∂Ψ)/(∂s)))+u ((s^2)+(t^2)-2 (u^2)) ((((((∂Ψ)/(∂s))))^2) -(((((∂Ψ)/(∂t))))^2) )-2t ((s^2)-(u^2)) (((∂Ψ)/(∂u)) ((∂Ψ)/(∂t)))-u ((s^2)-(t^2)) (((((∂Ψ)/(∂u))))^2)   );

Collect[quadM//.shortForm,((s^2)-(t^2)) u]

quadL//.shortForm

quadN//.shortForm

Collect[quadK//.shortForm,u ((s^2)+(t^2)-2 (u^2))]

Operators for relativistic corrections

With the same technique it is easy to write Laplacian for single particle variables. As an example consider the two operators which appear in the computation of relatvistic corrections to S-states:

Op1 = (Subscript[P, 1]^4) + (Subscript[p, 2]^4);    Op2 =(1/(Subscript[r, 12]))( (Subscript[p, 1])(Subscript[p, 2])+ (Subscript[n, 12])((Subscript[n, 12])·(Subscript[p, 1]))(Subscript[p, 2])).

(Subscript[n, 12]) is the direction (unit vector) going from electron 1 to electron 2.

Operator Op1

Operator Op1 is reduced, integrating by part to computation of

∫ ((Subscript[Δ, 1])ψ)((Subscript[Δ, 1])ψ)+ ∫ ((Subscript[Δ, 2])ψ)((Subscript[Δ, 2])ψ)

so we need the expression of the laplacian of each particle in Hylleraas coordinates. The definition is:

lap1[f_]:= Sum[grad1[grad1[f][[i]]][[i]],{i,1,3}];
lap2[f_]:= Sum[grad2[grad2[f][[i]]][[i]],{i,1,3}];

Explicitly:

collect2[espr_]:=Collect[espr,{Derivative[2,0,0][ψ][s,u,t],Derivative[0,0,2][ψ][s,u,t],Derivative[0,2,0][ψ][s,u,t],
Derivative[1,1,0][ψ][s,u,t],Derivative[1,0,1][ψ][s,u,t],Derivative[0,1,1][ψ][s,u,t]},Simplify]

laplaciano1 = collect2[Expand[
Simplify[TrigReduce[lap1[ψ[s,u,t]]/.ruleQ]//.ruleH,Assumptions-> {0<t<u<s}]]]

laplaciano2 = collect2[Expand[
Simplify[TrigReduce[lap2[ψ[s,u,t]]/.ruleQ]//.ruleH,Assumptions-> {0<t<u<s}]]]

An equivalent way of writing these operators is :

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

{Simplify[Delta1[s,u,t]== laplaciano1],Simplify[Delta2[s,u,t]== laplaciano2]}

Operator Op2

With the same technique we can give an expression for operator Op2.  We translate the derivative part, for calculation of matrix element it must be multiplied by the prefactor 1/u and by the Jacobian.

p1p2=Simplify[(1/2)TrigReduce[(Sum[grad1[grad2[ψ[s,u,t]][[i]]][[i]],{i,1,3}]
+Sum[grad2[grad1[ψ[s,u,t]][[i]]][[i]],{i,1,3}])/.ruleQ]//.ruleH,Assumptions-> {0<t<u<s}];
p1p2=collect2[p1p2];

ruleCos = {Cos[2θ]-> 2 ((Cos[θ])^2)-1};

n12n12p1p2 = Simplify[TrigReduce[
Sum[Sum[
(q1-q2)[[i]](q1-q2)[[j]]/(u^2)*(1/2)
(grad1[grad2[ψ[s,u,t]][[i]]][[j]]+
grad2[grad1[ψ[s,u,t]][[i]]][[j]]),{i,1,3}]
,{j,1,3}]/.ruleQ]//.ruleCos//.ruleH,Assumptions-> {0<t<u<s}];
n12n12p1p2=collect2[n12n12p1p2];

derivativeOp2=collect2[Simplify[p1p2+n12n12p1p2,Assumptions-> {0<t<u<s}]]

In textual form

shortForm={Derivative[0,0,1][ψ][s,u,t]-> "(∂_t)ψ",Derivative[0,0,2][ψ][s,u,t]-> "(∂_tt)ψ",Derivative[0,1,0][ψ][s,u,t]-> "(∂_u)ψ" ,Derivative[0,2,0][ψ][s,u,t]-> "(∂_uu)ψ" ,
Derivative[1,0,0][ψ][s,u,t]-> "(∂_s)ψ",Derivative[2,0,0][ψ][s,u,t]-> "(∂_ss)ψ",
Derivative[1,1,0][ψ][s,u,t]->"(∂_su)ψ",Derivative[0,1,1][ψ][s,u,t]->"(∂_tu)ψ",Derivative[1,0,1][ψ][s,u,t]->"(∂_st)ψ"
};

derivativeOp2/.shortForm

An equivalent expression for this operator is :

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

Simplify[derivOp2[s,u,t]== derivativeOp2]

The reader can check all expressions presented so far by a manual calculation of course.

Triangular (Pekeris) variables

A second set of variables very useful for the Helium atom is the set of triangular variables. These variables have been used by Pekeris in the variational study of  the Helium and will be used in notebook [NB-10.20-Helium-Pekeris.nb]. Their definition is

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

The great advantage over Hylleraas variables is that their variation is not constrained :

0 ≤ x,y,z ≤ ∞

Inverse transformation formulae are :

r1=((x+y)/2);       r2=((x+z)/2);        r12=((y+z)/2)

We can  repeat the previous steps, starting from definition of derivatives, from:

(Subscript[∇, 1])f[s,u,t]=(((Subscript[x, 1])/r1)-(((Subscript[x, 1])-(Subscript[x, 2]))/(Subscript[r, 12])) )(∂_x) f + (((Subscript[x, 1])/r1)+(((Subscript[x, 1])-(Subscript[x, 2]))/(Subscript[r, 12])) )(∂_y) f + (-((Subscript[x, 1])/r1)+(((Subscript[x, 1])-(Subscript[x, 2]))/(Subscript[r, 12])) )(∂_z) f;
(Subscript[∇, 2])f[s,u,t]=(((Subscript[x, 2])/r2)+(((Subscript[x, 1])-(Subscript[x, 2]))/(Subscript[r, 12])) )(∂_x) f + (-((Subscript[x, 2])/r2)-(((Subscript[x, 1])-(Subscript[x, 2]))/(Subscript[r, 12])) )(∂_y) f + (((Subscript[x, 2])/r2)-(((Subscript[x, 1])-(Subscript[x, 2]))/(Subscript[r, 12])) )(∂_z) f;

grad1[f_]:= Table[((∂_(q1[[i]])) f) ,{i,1,3}]+
2( ((( q1)/((x+y)))-((q1-q2)/(y+z)))(∂_x) f)+2( ((( q1)/((x+y)))+((q1-q2)/(y+z)))(∂_y) f)+2( (-(( q1)/((x+y)))+((q1-q2)/(y+z)))(∂_z) f);
grad2[f_]:= Table[((∂_(q2[[i]])) f) ,{i,1,3}]+
2( ((( q2)/((x+z)))+((q1-q2)/(y+z)))(∂_x) f)+2( (-(( q2)/((x+z)))-((q1-q2)/(y+z)))(∂_y) f)+2( ((q2/((x+z)))-((q1-q2)/(y+z)))(∂_z) f);

ruleQ = Join[Thread[Rule[q1,p1]],Thread[Rule[q2,p2]]];

We have to pass from radial coordinates r1, r2, θ to triangular coordinates. From previous formulas

r1=((x+y)/2);r2=((x+z)/2); r12=((y+z)/2) ⇒ Cos[θ] = (((x^2)-y z+x (y+z))/((x+y) (x+z)));

rulePek = {Cos[θ]-> (((x^2)-y z+x (y+z))/((x+y) (x+z))), 
Sin[θ]-> (Sqrt[1 - ((((x^2)-y z+x (y+z))/((x+y) (x+z))))^2]),
r1 -> ((x+y)/2), r2-> ((x+z)/2),r12-> ((y+z)/2)};

New variables are a simple linear transformation with respect Hylleraas variables

s= (1/2) (2 x+y+z),        u= ((y+z)/2),         t= (1/2) (-y+z)

and the jacobian from Cartesian coordinates can be found by a simple substitution

dJPek =Simplify[ dJ/.{s-> (1/2) (2 x+y+z),u-> ((y+z)/2),t-> (1/2) (-y+z)}]

The quadratic forms are now computed just by substitution, as before :

quadMPek= Simplify[TrigReduce[(1/2)dJPek*(grad1[ψ[x,y,z]].grad1[ψ[x,y,z]]+
grad2[ψ[x,y,z]].grad2[ψ[x,y,z]])/.ruleQ]//.rulePek,Assumptions-> {0<x && y>0 && z>0}];
quadMPek=collectDerivatives[quadMPek]

quadLPek = Simplify[TrigReduce[(dJPek*(-1/r12+ Z/r1 + Z/r2)((ψ[x,y,z])^2))/.ruleQ]//.rulePek,Assumptions-> {0<x && y>0 && z>0}]

quadNPek = Simplify[TrigReduce[(dJPek*((ψ[x,y,z])^2))/.ruleQ]//.rulePek,Assumptions-> {0<x && y>0 && z>0}]

quadKPek= Simplify[TrigReduce[dJPek*(grad1[ψ[x,y,z]].grad2[ψ[x,y,z]])/.ruleQ]//.rulePek,Assumptions-> {0<x && y>0 && z>0}];
quadKPek=collectDerivatives[quadKPek]

Text form

The quadratic forms written above are used in all notebooks on Helium. We note for completeness that the quadratic form L is defined as -V where V is the potential, the sign is easily checked looking at Z dependent terms. The variational principle is written as

δ( M + V - λ N)=0    ⟺   δ( M - L - λ N)=0

To help the reader we rewrite above expressions in simpler notation using

shortFormPek={Derivative[0,0,1][ψ][x,y,z]-> "(((∂ψ)/(∂z)))",Derivative[0,1,0][ψ][x,y,z]-> "(((∂ψ)/(∂y)))" ,
Derivative[1,0,0][ψ][x,y,z]-> "(((∂ψ)/(∂x)))"};

The results have to be compared with

M=(∫_0^∞)d𝓍 (∫_0^∞)d𝓎(∫_0^∞)d𝓏 (𝓍 ((𝓎^2)+(𝓏^2)+𝓍 (𝓎+𝓏)) (((((∂Ψ)/(∂𝓍))))^2) -2𝓍 𝓎 (𝓍+𝓎) (((∂Ψ)/(∂𝓎)) ((∂Ψ)/(∂𝓍)))- 2𝓍 𝓏 (𝓍+𝓏) (((∂Ψ)/(∂𝓏)) ((∂Ψ)/(∂𝓍)))+
𝓎 (2 (𝓍^2)+2 (𝓎+𝓏) 𝓍+𝓏 (𝓎+𝓏)) (((((∂Ψ)/(∂𝓎))))^2) +𝓏 (2 (𝓍^2)+2 (𝓎+𝓏) 𝓍+𝓎 (𝓎+𝓏)) (((((∂Ψ)/(∂𝓏))))^2) );
L=(∫_0^∞)d𝓍 (∫_0^∞)d𝓎(∫_0^∞)d𝓏  ((1/4) ((((𝓎-𝓏))^2)-(((2 𝓍+𝓎+𝓏))^2)+4 Z (𝓎+𝓏) (2 𝓍+𝓎+𝓏)) ((Ψ )^2));
N=(∫_0^∞)d𝓍 (∫_0^∞)d𝓎(∫_0^∞)d𝓏 ((1/2) ((𝓍+𝓎) (𝓍+𝓏) (𝓎+𝓏))(Ψ^2)) ;
K=(∫_0^∞)d𝓍 (∫_0^∞)d𝓎(∫_0^∞)d𝓏 ( 2 𝓍 𝓎 𝓏 (((((∂Ψ)/(∂𝓍))))^2) +2𝓍 𝓎 (𝓍+𝓎) (((∂Ψ)/(∂𝓎)) ((∂Ψ)/(∂𝓍)))+
2𝓍 𝓏 (𝓍+𝓏) (((∂Ψ)/(∂𝓏)) ((∂Ψ)/(∂𝓍)))-2 𝓍 𝓎 (𝓍+𝓎+𝓏) (((((∂Ψ)/(∂𝓎))))^2) -2𝓎 𝓏 (𝓎+𝓏) (((∂Ψ)/(∂𝓏)) ((∂Ψ)/(∂𝓎)))-2 𝓍 𝓏 (𝓍+𝓎+𝓏) (((((∂Ψ)/(∂𝓏))))^2) );

quadMPek//.shortFormPek

quadLPek//.shortFormPek

quadNPek//.shortFormPek

Simplify[Expand[quadKPek//.shortFormPek]]

Operators for relativistic corrections

For completeness let us now compute the operators needed for relativistic corrections in terms of Pekeris variables. For S-states the operators are

Op1 = (Subscript[P, 1]^4) + (Subscript[p, 2]^4);    Op2 =(1/(Subscript[r, 12]))( (Subscript[p, 1])(Subscript[p, 2])+ (Subscript[n, 12])((Subscript[n, 12])·(Subscript[p, 1]))(Subscript[p, 2])).

(Subscript[n, 12]) is the direction (unit vector) going from electron 1 to electron 2.

Operator Op1

Operator Op1 is reduced, integrating by part to computation of

∫ ((Subscript[Δ, 1])ψ)((Subscript[Δ, 1])ψ)+ ∫ ((Subscript[Δ, 2])ψ)((Subscript[Δ, 2])ψ)

so we need the expression of the laplacian of each particle in Pekeris coordinates. The definition is:

lap1Pk[f_]:= Sum[grad1[grad1[f][[i]]][[i]],{i,1,3}];
lap2Pk[f_]:= Sum[grad2[grad2[f][[i]]][[i]],{i,1,3}];

Explicitly:

collect2Pk[espr_]:=Collect[espr,{Derivative[2,0,0][ψ][x,y,z],Derivative[0,0,2][ψ][x,y,z],Derivative[0,2,0][ψ][x,y,z],
Derivative[1,1,0][ψ][x,y,z],Derivative[1,0,1][ψ][x,y,z],Derivative[0,1,1][ψ][x,y,z]},Simplify]

laplaciano1Pk = collect2Pk[Expand[
Simplify[TrigReduce[lap1Pk[ψ[x,y,z]]/.ruleQ]//.rulePek,Assumptions-> {0<x&& y.0 && z>0}]]]

laplaciano2Pk = collect2Pk[Expand[
Simplify[TrigReduce[lap2Pk[ψ[x,y,z]]/.ruleQ]//.rulePek,Assumptions->  {0<x&& y.0 && z>0}]]]

These operators can also be written as

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

{Simplify[Delta1[x,y,z]== laplaciano1Pk],Simplify[Delta2[x,y,z]== laplaciano2Pk]}

Operator Op2

With the same technique we can give an expression for operator Op2.  We translate the derivative part, for calculation of matrix element it must be multiplied by the prefactor 1/(Subscript[r, 12]) and by the Jacobian.

p1p2=Simplify[(1/2)TrigReduce[(Sum[grad1[grad2[ψ[x,y,z]][[i]]][[i]],{i,1,3}]
+Sum[grad2[grad1[ψ[x,y,z]][[i]]][[i]],{i,1,3}])/.ruleQ]//.rulePek,Assumptions-> {0<x&& y.0 && z>0}];
p1p2=collect2Pk[p1p2];

ruleCos = {Cos[2θ]-> 2 ((Cos[θ])^2)-1};

n12n12p1p2 = Simplify[TrigReduce[
Sum[Sum[
(4/(((y+z))^2))(q1-q2)[[i]](q1-q2)[[j]]*(1/2)
(grad1[grad2[ψ[x,y,z]][[i]]][[j]]+
grad2[grad1[ψ[x,y,z]][[i]]][[j]]),{i,1,3}]
,{j,1,3}]/.ruleQ]//.ruleCos//.rulePek,Assumptions-> {0<x&& y.0 && z>0}];
n12n12p1p2=collect2Pk[n12n12p1p2];

derivativeOp2=collect2Pk[Simplify[p1p2+n12n12p1p2,Assumptions-> {0<x&& y.0 && z>0}]]

In textual form

shortForm={Derivative[0,0,1][ψ][x,y,z]-> "(∂_z)ψ",Derivative[0,0,2][ψ][x,y,z]-> "(∂_zz)ψ",Derivative[0,1,0][ψ][x,y,z]-> "(∂_y)ψ" ,Derivative[0,2,0][ψ][x,y,z]-> "(∂_yy)ψ" ,
Derivative[1,0,0][ψ][x,y,z]-> "(∂_x)ψ",Derivative[2,0,0][ψ][x,y,z]-> "(∂_xx)ψ",
Derivative[1,1,0][ψ][x,y,z]->"(∂_xy)ψ",Derivative[0,1,1][ψ][x,y,z]->"(∂_yz)ψ",Derivative[1,0,1][ψ][x,y,z]->"(∂_xz)ψ"
};

derivativeOp2/.shortForm

An equivalent expression for this operator is :

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

Simplify[derivOp2[x,y,z]== derivativeOp2]