This notebook reproduces some results in the papers:
[1] C. L. Pekeris: "Ground State of Two-Electron Atoms", Phys. Rev. 112(5) (1958), 1649-1658. 
[2] C. L. Pekeris: "1 1S and 2 3S States of Helium", Phys. Rev. 115(5) (1959), 1216-1221.

Some formulas from the paper [1]

Some steps make use of our package HolonomicFunctions, that can be downloaded from:
http://www.risc.uni-linz.ac.at/research/combinat/software/

<<HolonomicFunctions.m

HolonomicFunctions package by Christoph Koutschan, RISC-Linz, Version 1.3 (25.01.2010)
⟶ Type  ?HolonomicFunctions  for help

The Schrödinger differential equation, formula (1). Note that E here is not Euler's constant.

pek1[n_]:=ToOrePolynomial[Sum[Der[x[i]]^2+Der[y[i]]^2,{i,n}]+2*(E+Z/Sqrt[Sum[x[i]^2,{i,n}]]+Z/Sqrt[Sum[y[i]^2,{i,n}]]-1/Sqrt[Sum[(y[i]-x[i])^2,{i,n}]]),OreAlgebra@@Join[Table[Der[x[i]],{i,n}],Table[Der[y[i]],{i,n}]]]

Formula (5):

pek5=ToOrePolynomial[
Der[r1]^2+(2/r1)**Der[r1]+Der[r2]^2+(2/r2)**Der[r2]+2*Der[r12]^2+(4/r12)**Der[r12]
+((r1^2-r2^2+r12^2)/(r1*r12))**Der[r1]**Der[r12]
+((r2^2-r1^2+r12^2)/(r2*r12))**Der[r2]**Der[r12]
+2*(E+Z/r1+Z/r2-1/r12),
OreAlgebra[Der[r1],Der[r2],Der[r12]]];

Formula (14)  (ϵ is represented by e):

pek14=ToOrePolynomial[e*((4u^2*v+4u*v^2+4u^2*w+4u*v*w+2*u*w^2)**Der[u]^2+(4u^2*v+4u*v^2+4v^2*w+4u*v*w+2v*w^2)**Der[v]^2+(8u^2*w+8v^2*w+4u*w^2+4v*w^2)**Der[w]^2-(4u*w*(2u+w))**Der[u]**Der[w]-(4v*w*(2v+w))**Der[v]**Der[w]+(-4u^2+4v^2+2w^2+4u*w+4v*w+8u*v-4u*v*w-4u*v^2-4u^2*v)**Der[u]+(4u^2-4v^2+2w^2+4u*w+4v*w+8u*v-4u*v*w-4u*v^2-4u^2*v)**Der[v]+(8u^2+8v^2-4w^2-2u*w^2-2v*w^2-4u^2*w-4v^2*w)**Der[w]-4*(u+v)*(u+v+w))+4*Z*(u+v)*(u+v+w)-(2u+w)*(2v+w),
OreAlgebra[Der[u],Der[v],Der[w]]];

Derive recurrence (22) in [1]

The starting point is the time-independent Schrödinger equation

ApplyOreOperator[pek1[3],ψ[x[1],x[2],x[3],y[1],y[2],y[3]]]//TraditionalForm

We write a little procedure that performs algebraic substitution in a differential equation.
Unfortunately, the command DFiniteSubstitute from our package cannot be applied, since the differential equation that we start with does not have polynomial coefficients.

DifferentialSubstitute[op1_,subs:{__Equal}]:=Module[{op=Flatten[{op1}],a1,a2,eqs,psi,xv,uv,vars,aes,gb},
op=ToOrePolynomial[op];
xv=First/@First[OreAlgebra[op]]; (* the original variables *)
uv=First/@subs;                                      (* the variables after substitution *)
vars=Join[xv,uv];

(* Get the algebraic equations for the substitutions. *)aes=HolonomicFunctions`Private`AlgebraicEquation[Last[#],xv,First[#]]&/@subs;
gb=GroebnerBasis[aes,vars,MonomialOrder->Lexicographic];

(* Do some rewriting with the equations, replacing old variables by new ones. *)eqs=ApplyOreOperator[op,psi@@(Last/@subs)];
eqs=eqs/.(If[Head[#[[2]]]===Power,#[[2,1]]^a1_.->#[[1]]^(a1/#[[2,2]]),#[[2]]->#[[1]]]&/@subs);
eqs=NormalizeCoefficients/@Together[ToOrePolynomial[eqs,psi@@uv,OreAlgebra@@(Der/@uv)]];

(* Reduce all coefficients with the aes Groebner basis. *)eqs=Function[opoly,OrePolynomial[Transpose[{Together[Last[PolynomialReduce[#,gb,vars,MonomialOrder->Lexicographic]]&/@OrePolynomialListCoefficients[opoly]],Last[Transpose[First[opoly]]]}],opoly[[2]],opoly[[3]]]]/@eqs;
Return[If[Head[op1]===List,#,First[#]]&[Expand[NormalizeCoefficients/@eqs]]];
];

With help of this procedure, the change of variables {x1, ..., y3} -> {r1, r2, r12} can be done:

pek=DifferentialSubstitute[pek1[3],
{r1==Sqrt[Sum[x[i]^2,{i,3}]],
r2==Sqrt[Sum[y[i]^2,{i,3}]],
r12==Sqrt[Sum[(y[i]-x[i])^2,{i,3}]]}]

The result agrees with formula (5) in the paper [1].

pek-NormalizeCoefficients[pek5]

Now the energy parameter E is replaced.

pek=pek/.E->-e^2;

The change of variables to perimetric coordinates can also be obtained by the above procedure DifferentialSubstitute:

pek=DifferentialSubstitute[pek,{u==e*(r2+r12-r1),v==e*(r1+r12-r2),w==2*e*(r1+r2-r12)}]

The package HolonomicFunctions can be used to perform the transformation of the differential equation in order to incorporate the exponential part of the solution. The above differential operator annihilates ψ and using the closure property product, we get an operator for exp((u+v+w)/2)*ψ:

pek=First[DFiniteTimes[{pek},Annihilator[Exp[(u+v+w)/2],{Der[u],Der[v],Der[w]}]]]

Again, our result agrees with formula (14) in [1]:

pek-pek14

Now the differential equation has to be transformed into a recurrence for the coefficients of the power series expansion with respect to Laguerre polynomials. This can be done by reduction with Gröbner bases. For this purpose, first an annihilating ideal for the product of Laguerre polynomials is computed

Timing[ann=Annihilator[LaguerreL[l,u]*LaguerreL[m,v]*LaguerreL[n,w],{Der[u],Der[v],Der[w],S[l],S[m],S[n]},MonomialOrder->EliminationOrder[3]]]

(* a faster variant, giving the same result *)
Timing[ann=SortBy[ToOrePolynomial[Join@@(Annihilator[LaguerreL@@#,{Der[#[[2]]],S[#[[1]]]}]&/@{{l,u},{m,v},{n,w}}),OreAlgebra[Der[u],Der[v],Der[w],S[l],S[m],S[n]],MonomialOrder->EliminationOrder[3]],LeadingExponent]]

The following reduction rewrites all D's in terms of shifts:

pek=OreReduce[pek,ann]

Next, the variables u,v,w have to be eliminated. Again, this can be obtained by means of Gröbner bases. An annihilating ideal w.r.t. the shift operators is computed, including negative shifts:

ann=Annihilator[LaguerreL[l,u]*LaguerreL[m,v]*LaguerreL[n,w],{S[l],S[m],S[n]}]

ann=OreGroebnerBasis[Join[ann,{S[l]*S[-1][l]-1,S[m]*S[-1][m]-1,S[n]*S[-1][n]-1}],OreAlgebra[u,v,w,S[l],S[m],S[n],S[-1][l],S[-1][m],S[-1][n]],MonomialOrder->EliminationOrder[3]]

op22=OreReduce[pek,ann];

This operator now represents the 33-term recurrence that appears as formula (22) in [1].

Length[Support[op22]]

Or, in a more human readable format:

rel22=Collect[ApplyOreOperator[op22,A[l,m,n]],A[__],Simplify]

The recurrence is self-adjoint and therefore we are done. This is checked by verifying Eq. (24) in [1]:

Union[Flatten[Expand[Table[Coefficient[op22,S[-1][l]^a*S[-1][m]^b*S[-1][n]^c]-(Coefficient[op22,S[l]^a*S[m]^b*S[n]^c]/.l->l-a/.m->m-b/.n->n-c),{a,0,2},{b,0,2},{c,0,2}]]]]

Matrix in the symmetric case (para)

We define the scheme that is used for ordering the indices l,m,n  (formula (28), corrected).

SchemePara[l_Integer/;l≥0,m_Integer/;m≥0,n_Integer/;n≥0]:=With[{w=l+m+n},(1/24)*w*(w+2)*(2w+5)+(1/16)*(1-(-1)^w)+(1/4)*(l+m)^2+(1/8)*(1-(-1)^(l+m))+(l+m)/2+l+1];

The following procedure generates the matrix for the symmetric ("para") case: A(l,m,n) = A(m,l,n).

MatrixPara[rel_,A_[l_,m_,n_],omega_Integer]:=Module[{a1,a2,a3,a,c,t,l1,lm,w,dim,rel2,row,matrix={}},dim=Ceiling[(omega+1)*(omega+3)*(2*omega+7)/24];
Do[
rel2=rel/.{l->l1,m->lm-l1,n->w-lm}
/.A[a1_,a2_,a3_]/;a1>a2->A[a2,a1,a3]
 /.A[a1_,a2_,a3_]/;a1+a2+a3>omega->0;
rel2=If[Head[rel2]===Plus,List@@rel2,{rel2}];
row=Table[0,{dim}];
Function[t,
a=First[Cases[t,A[__],Infinity]];
c=t/.a->1;
row[[SchemePara@@a]]+=c;]/@rel2;
AppendTo[matrix,row];
,{w,0,omega},{lm,0,w},{l1,0,lm/2}];
Return[matrix];
];

TableForm[MatrixPara[rel22,A[l,m,n],1]]

We now can reproduce Table II in [1]  (p. 1651). The rows are normalized as to make the matrix symmetric.

TableForm[Expand[MatrixPara[-rel22,A[l,m,n],2]/{2,2,1,2,1,1,2}]]

Matrix in the antisymmetric case (ortho)

We define the scheme that is used for ordering the indices l,m,n  (formula (28), corrected).

SchemeOrtho[l_Integer/;l≥0,m_Integer/;m≥0,n_Integer/;n≥0]:=With[{w=l+m+n},(1/24)*w*(w+2)*(2w-1)-(1/16)*(1-(-1)^w)+l*(m+n)+m];

The following procedure generates the matrix for the antisymmetric ("ortho") case: A(l,m,n) = -A(m,l,n).

MatrixOrtho[rel_,A[l_,m_,n_],omega_Integer]:=Module[{a1,a2,a3,a,c,t,l1,m1,w,dim,rel2,row,matrix={}},dim=Floor[(omega+1)*(omega+3)*(2*omega+1)/24];
Do[
rel2=rel/.{l->l1,m->m1,n->w-l1-m1}
/.A[a1_,a2_,a3_]/;a1>a2->-A[a2,a1,a3]
 /.A[a1_,a1_,a2_]->0
   /.A[a1_,a2_,a3_]/;a1+a2+a3>omega->0;
rel2=If[Head[rel2]===Plus,List@@rel2,{rel2}];
row=Table[0,{dim}];
Function[t,
a=First[Cases[t,A[__],Infinity]];
c=t/.a->1;
row[[SchemeOrtho@@a]]+=c;]/@rel2;
AppendTo[matrix,row];
,{w,0,omega},{l1,0,(w-1)/2},{m1,l1+1,w-l1}];
Return[matrix];
];

We now can reproduce Table IV in [2]  (p. 1219):

TableForm[Expand[MatrixOrtho[-rel22,A[l,m,n],3]/2]]

For relatively small dimensions, the determinant can be evaluated symbolically:

Det[MatrixOrtho[rel22/.Z->2,A[l,m,n],4]]

NSolve[%==0,e]

A simple procedure using chinese remaindering and polynomial interpolation can symbolically evaluate such determinants much faster:

(* it is assumed that the coefficients are integers. *)
DetUniv[mat_]:=
Module[{deg,var,det,det1,i,p,pp,vals,len,mats,xs},
var=Variables[mat];
If[Length[var]===1,var=First[var],Throw["matrix entries must be univariate polynomials."]];
deg=Total[Function[row,Max[Exponent[#,var]&/@row]]/@mat];
mats=Table[mat/.var->i,{i,Ceiling[-deg/2],Ceiling[deg/2]}];
xs=Range[Ceiling[-deg/2],Ceiling[deg/2]];
det={};det1=Null;
p=Developer`$MaxMachineInteger;pp=1;
While[det=!=det1,
p=NextPrime[p,-1];
vals={};
det1=det;
vals=Transpose[{xs,Det[#,Modulus->p]&/@mats}];
det=Mod[CoefficientList[InterpolatingPolynomial[vals,var,Modulus->p],var],p];
If[det1=!={},
len=Max[Length/@{det,det1}];
{det,det1}=PadRight[#,len]&/@{det,det1};
det=ChineseRemainder[#,{pp,p}]&/@Transpose[{det1,det}];
];
pp*=p;
det=det/.a_Integer/;a>pp/2:>a-pp;
];
Return[det.Table[var^i,{i,0,Length[det]-1}]];
];

Timing[Last[NSolve[Det[MatrixOrtho[rel22/.Z->2,A[l,m,n],10]]==0,e,WorkingPrecision->20]]]

Timing[Last[NSolve[DetUniv[MatrixOrtho[rel22/.Z->2,A[l,m,n],10]]==0,e,WorkingPrecision->20]]]

Generate Tables

We can generate (with symbolic determinants!) the first two columns of  Table III on page 1652 in [1].

Timing[tab=Table[Last[NSolve[DetUniv[MatrixPara[rel22,A[l,m,n],om]]==0,e,WorkingPrecision->12]][[1,2]],{Z,10},{om,8,9}];]

TableForm[tab]

For generating the values in the tables, the symbolic evaluation of the determinants is too slow. Therefore we use numeric methods which can be applied when the problem is reformulated as to find the generalized eigenvalues of the following two matrices:

mat=MatrixPara[rel22/.Z->2,A[l,m,n],8];
mA=mat/.e->0;
mB=-(mat-mA)/.e->1;

Eigenvalues[N[{mA,mB}],1]

Unfortunately, Mathematica does not allow for more digits:

Eigenvalues[N[{mA,mB},20],1];

Therefore, we will perform these numeric computations in Matlab, too (with higher precision and much faster than below).

PekerisTable[rel_,A_[l_,m_,n_],case_String,oms_List,zs_List]:=Module[{mat,mat1,matf,omi,om,zi,matA,matB,ev,col,table},matf=Switch[case,"para",MatrixPara,"ortho",MatrixOrtho,_,Throw["Unknown case. Give either \"para\" or \"ortho\"."]];
table=Join[{{"\\ ω"},{"Z\\n"},{"---"}},Transpose[{zs}]];
Do[
om=oms[[omi]];
mat=If[Length[zs]===1,matf[rel/.Z->First[zs],A[l,m,n],om],matf[rel,A[l,m,n],om]];
col={om,First[Dimensions[mat]],"--------"};
Do[
mat1=If[Length[zs]>1,mat/.Z->zs[[zi]],mat];
matA=mat1/.e->0;
matB=-(mat1-matA)/.e->1;
ev=Eigenvalues[N[{matA,matB}],1];
AppendTo[col,ev];
,{zi,Length[zs]}];
table=MapThread[Append,{table,col}];
,{omi,Length[oms]}];
Return[table];
];

Table III on page 1652 in [1].

Timing[tab1=PekerisTable[rel22,A[l,m,n],"para",Range[8,11],Range[10]];]

TableForm[tab1]

Table I on page 1217 in [2].

Timing[tab2=PekerisTable[rel22,A[l,m,n],"para",Range[12,21,3],{2}];]

TableForm[tab2]

Table V on page 1220 in [2].

Timing[tab3=PekerisTable[rel22,A[l,m,n],"ortho",Range[10,19,3],{2}];]

TableForm[tab3]