Complex scaling and Helium

Quantum Mechanics
 K. Konishi, G. Paffuti

Use the complex scaling methods to give a qualitative study of metastable Helium states.

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

Introduction

Taking away one electron from the Helium atom leaves an ion with an electron in one of the hydrogen-type level 1s, 2s,   with energy

(Subscript[E, ion])=- ((Z^2)/(2(n^2))) .

From these points a cut starts in the complex energy plane, corresponding to a branch of the continuous spectrum of the resolvent. Under complex scaling the cut rotates exposing metastable states with Im[E]<0.

In the following we shall compute eigenvalues for a complex scaling of the coordinates and extract the parameters of the resonances. The computation uses the function matrices[] discussed in [NB-10.17] which computes the matrices needed to write the Hamiltonian in the Hylleraas basis.

In the variational approach one has to solve the equation

((k^2) ℳ - k ℒ )ψ = 𝒩 ψ ;

and optimize eigenvalues with respect to k . All scaling dependence is in k. In a complex scaling transformation (k transforms as the inverse of  a coordinate):

k→ (Subscript[k, 0])Exp[-i θ].

A careful computation would require optimization with respect to (Subscript[k, 0]), in our qualitative study we will keep (Subscript[k, 0])fixed.

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

Computation

We fix some parameters of the variational computation. maxs, maxt, maxu are the maximum power of the variables s, t, u used in the Hylleraas basis. The parameter shT has to be set to 0 for singlet states.

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

These are the expected branching points

bPoints = - (Z^2)/(2. ((Range[3])^2));

Here we compute the matrices

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

We choose a list of θ values and write a function to plot eigenvalues together with lines indicating the expected position of rotated spectrum (a rotation of 2θ in clockwise direction).

ntheta=8;
θ=  Pi/96.Range[0,ntheta];

plotEigs[i_]:= ListPlot[Transpose[{Re[eigctot[[i]]],Im[eigctot[[i]]]}],PlotStyle->PointSize[0.015],PlotRange-> {{-3,1},{-0.5,0.1}},
Epilog-> Table[Line[{{bPoints[[k]],0},{bPoints[[k]]+5 Cos[2θ[[i]]],-5 Sin[2θ[[i]]]}}],
{k,1,3}]]

A few trials show that a reasonable value of (Subscript[k, 0])for the first branch point is the one used as a starting point in the computation of helium spectrum, see [NB-10.17]. Here we plot the result for the last θ value

k0= N[2*(Z-(5/16))];;
eigctot = Table[k=k0*Exp[-I θ[[i]]];Eigenvalues[{k^2 matM-k matL,matN}],{i,1,Length[θ]}];

plotEigs[9]

and an animation to show the rotation of the spectrum.

ListAnimate[Table[plotEigs[k],{k,1,Length[θ]}],AnimationRunning->False]

The region around the first branch point does not contain resonant states.

A more reasonable value for the region around the second branch point, corresponding to an ion in the 2s state, is given by (Subscript[k, 0])= 2. We repeat the computation and the plot:

k0=2.;
eigctot = Table[k=k0*Exp[-I θ[[i]]];Eigenvalues[{k^2 matM-k matL,matN}],{i,1,Length[θ]}];
plotEigs[9]

Here we list the first few "stable" states under scaling, i.e. the resonances

(tab = Table[pos=Flatten[
Position[(-.79 < Re[#] && Re[#] <-0.55&&-0.1< Im[#]<0.001)&/@eigctot[[itheta]],True]];eigctot[[itheta]][[pos]],{itheta,1,Length[θ]}])//TableForm

A much  more accurate computation is performed in ref.[1] where the following values for the low lying resonances are found

Resonances = {-0.777867 -I 0.0022706,-0.621927 - I 0.0001078, -0.589894 -I 0.000681}

We see a reasonable agreement, even with our simple treatment of the optimization parameter.

References

1) A. Bürgers, D. Wintgen, J.-M. Rost, J. Phys. B, At. Mol. Opt. 28, 3163 (1995).