Complex scaling

Quantum Mechanics
K. Konishi, G. Paffuti

Use the complex scaling method to study resonances in a radial potential (Subscript[V, 0])(r^2)Exp[-r] and in a one-dimensional potential  (x^2)/(2 (1+g (x^4))).

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

Introduction

Summary of the theory

In the text it has been shown that the complex scaling transformation

r→ r (e^(i θ))  ⇒  H[θ] = (1/(2 (e^(2i θ))))(p^2) + V[ r (e^(i θ))],

allows an analytic continuation of the resolvent for Im[E] < 0. In this transformation momenta rotate as (e^(-iθ)), consequently the continuous spectrum of H has a rotation of 2θ in the clockwise sense. H[θ] shows also discrete eigenvalues in the half plane Im[E] < 0 and these correspond to resonant states of H. Intuitively the cut corresponding to the continuous spectrum rotates and expose poles in the second sheet of the Riemann surface on which the resolvent is defined. A resonance is a pole at E = (Subscript[E, 0])-iΓ/2 then it appears for 2θ > ArcTan[Γ/2 E].

The above transformation works also for complex parameters, and writing α Exp[i θ] for the transformation parameter

r→ r α (e^(i θ))  ⇒  H[α,θ] = (1/(2 (α^2)(e^(2i θ))))(p^2) + V[ r α (e^(i θ))],

The transformation can also be considered as a basis transformation. If (Subscript[f, n])[x]is a basis in the Hilbert space, the matrix elements of the transformed Hamiltonian are the same as those of the original Hamiltonian with respect to the basis (Subscript[f, n])[x/α (e^(i θ))]

The computation of eigenvalues of H[α,θ] can be done using standard techniques, and in the following will be approximated by the use of a finite size basis (Subscript[f, n])[x]. As it is well known from variational computations the finite size of the basis gives a factitious dependence on the parameters for the eigenvalues, the dependence disappear for infinite complete basis. We will look for ``stability regions'' where the eigenvalues are stable, as explained in the notebooks on variational computations.

Basis for the variational computation

The choice of the basis is quite arbitrary. In the following we shall use a Fourier basis.

We discretize a "large" volume L in Npt points starting with x0. We will take the first dim Fourier components to write the Hamiltonian. The function

calcH[V,x0, L, dim, Npt, θ]

computes with a given function V for the potential the Hamiltonian matrix, rotated by an angle θ. The computational details are discussed in notebook [NB-10.12].

Code for functions

calcV[Vx_,dim_]:= Module[{mat,fV,nV,fV1,p2},
mat =ConstantArray[0.,{dim,dim}];
nV = Length[Vx]-1;
fV = FourierDCT[Re[Vx],1] + I  FourierDCT[Im[Vx],1];
fV1 = Drop[fV,-1];
fV = Join[fV,Reverse[fV1]];
Table[mat[[i,j]] = fV[[j-i+1]] - fV[[i+j+1]],{i,1,dim-1},{j,i+1,dim}];
mat = mat+Transpose[mat];
Table[mat[[i,i]] = fV[[1]] - fV[[2 i+1]],{i,1,dim}];
mat =(1/(Sqrt[2 nV]))mat]

calcH[V_,x0_,L_,dim_,Npt_,θ_]:= Module[{p2,Φ,xpt,Vx,matV,ham},
Φ = Exp[I θ];
p2 =(1/(Φ^2))N[DiagonalMatrix[((( π/L*Range[1,dim]))^2)]];
xpt = Φ (x0+(L/Npt)Range[0,Npt]);
Vx = V[xpt];
matV = calcV[Vx,dim];
ham = 0.5*p2 + matV]

Metastable states

Radial potential model

We study a model for which we know that there are some metastable states, computed both with semiclassical methods and with Gamow Siegert approximation, see [NB-11.6]

V0 = 7.5;rmax = 2; 
V[r_]:= V0 (r^2)(E^(-r));
Plot[V0 (r^2)(E^(-r)),{r,0,20},PlotRange->All,AxesLabel-> {"r","V(r)"}]

Example with fixed θ

First we define a grid and choose a dimension for our matrices. We perform the computation for several values of L, the maximum distance allowed. In this example the rotation angle θ is kept fixed and we check the stability in the choice of L

dim = 50; θ=π/12; Npt=513;
x0 = 0;
Lmin = 3.5; Lmax = 15; dL = 0.5;
L = Range[Lmin,Lmax,dL];

The eigenvalues are computed using the output of the function calcH, which gives the Hamiltonian

eigH=Table[Sort[Eigenvalues[calcH[V,x0,L[[i]],dim,Npt,θ]]],{i,1,Length[L]}];//Timing

In the following we shall repeatedly do a plot of eigenvalues together with a line showing the hypothetical rotation 2θ of the cut of continuous eigenvalues. We will use a function

plotEig[k, eigH, θ, px]

which perform the plot of the k-th row of eigenvalues of eigH, i.e. for the k-th value of L. px is a parameter which fixes the right end on the plot, allowing an easy zoom of the interesting zone.

plotEig[k_,eigH_,θ_,px_]:=Module[{pl1},
pl1=ListPlot[Transpose[{Re[eigH[[k]]],Im[eigH[[k]]]}],
PlotRange-> {{0,px},{-(px/2),0.2}},PlotStyle-> PointSize[0.018],
AxesLabel-> {"Re[E]","Im[E]"}, 
PlotLabel-> Row[{"θ = ",InputForm[θ]}],
Epilog-> Line[{{0,0},{px,-px Tan[2θ]}}]]]

Here we see that a couple of eigenvalues are left behind the cut rotation, the zoomed image can help the reader. These values correspond to resonances.

GraphicsRow[{plotEig[20,eigH,θ,40],plotEig[20,eigH,θ,12]},ImageSize-> {9,3} 72]

Once it is known how many eigenvalues are interesting we can select them by a cut on their argument, |Arg[λ]| < 2θ. The function selectEig do the work for all L values:

selectEig[eigH_,θ_,numeig_]:=
Module[{numval,lis,val},
val = {};
Do[
lis = Select[eigH[[j]],Arg[#1]> -0.9*2θ &];
lis = PadRight[lis,numeig];
val = Append[val,lis];,
{j,1,Length[eigH]}];
val]

These are all selected resonances

resonances = selectEig[eigH,θ,2];

TableForm[resonances,
TableHeadings-> {L,{ "    E0 - i Γ/2", "    E1 - i Γ/2"}},TableSpacing-> {2,3}]

We see that for large enough L the two resonances are clearly stable. Their real part and the associated width are

eigdef=resonances[[-1]]; 
E0 = Re[eigdef[[1]]]; Γ0 = -2 Im[eigdef[[1]]];
E1 = Re[eigdef[[2]]]; Γ1 = -2 Im[eigdef[[2]]];
Print["(E0,Γ0) = ",{E0,Γ0},"  (E1,Γ1) = ",{E1,Γ1}]

This result is in agreement with the one obtained in notebook [NB-11.6] by the Gamow Siegert method

(E0,Γ0) = {3.42639,  0.025549}  (E1,Γ1) = {4.85242, 2.1188}   ⟵ Gamow - Siegert

Let us note that the second resonant state is higher than potential maximum, then it cannot be interpreted as a simple barrier transmission effect.

Changing θ and new resonances

The flexibility of the method allows to give an explicit example of the mechanism which produces the appearance of new resonances as θ grows. Let us do the same calculation for θ = π/6, it is clear that new eigenvalues have been left behind by the cut rotation:

θ=π/6;
eigH=Table[Sort[Eigenvalues[calcH[V,x0,L[[i]],dim,Npt,θ]]],{i,1,Length[L]}];

GraphicsRow[{plotEig[20,eigH,θ,30],plotEig[20,eigH,θ,8]},ImageSize-> {9,3} 72]

Stability in θ

A crucial check is that our results do not depend on θ, or better there is  a stability region for the results.

Let us fix L and vary θ:

L = 15;
θ = Range[12]*π/48;

eigH=Table[Sort[Eigenvalues[calcH[V,x0,L,dim,Npt,θ[[i]]]]],{i,1,Length[θ]}];

Here as usual we select eigenvalues which pass the cut test

resonances = {};
Do[{
ris = Select[eigH[[j]],Arg[#1]> -0.8*2θ[[j]] &];
ris = PadRight[ris,3];
resonances=Append[resonances,ris];
},{j,1,Length[θ]}];

The results show the existence of a third broad resonance

TableForm[Join[Transpose[{θ}],resonances,2],
TableHeadings-> {None,{"θ", "    E0 - i Γ/2", "    E1 - i Γ/2","    E2 - i Γ/2"}},
TableSpacing-> {3,2}]

It is apparent a stabilization up to 2θ = π/2, as expected.

The third resonance appear starting from θ = π/8. We can compare Arg[λ]/2 with the θ values

-1/2Table[Arg[resonances[[6,k]]],{k,1,3}]

N[θ]

and we see a qualitative agreement with the appearance of the resonances only for θ> Abs[0.5 Arg[Γ/2 E]].

One dimensional model

As a second exampe we consider a one dimentional model:

g = 0.1;xmax = (g^(-1/4));   (* g is a global variable in the following *)
V[x_]:= (1/2)((x^2)/(1+g (x^4)))
 Plot[V[x],{x,-5xmax,5xmax},AxesLabel-> {"x","V[x]"}]

The system can have metastable states, a particle trapped in the well can escape from both sides of the potential. We will use again a Fourier basis.

This is the grid. The parameter x0 is now fixed to - L/2 in order to have a symmetric grid.

dim = 100; Npt=1025;
L=35;
x0 = -L/2;
θ = (π/48)Range[1,12];θ[[12]]=(97/100)*(π/4);

The eigenvalues are computed as in the previous case

eigHTheta=Table[Eigenvalues[calcH[V,x0,L,dim,Npt,θ[[i]]]],{i,1,Length[θ]}];//Timing

An animation can clearly show the evolution of the cut and the appearance of the resonances

ListAnimate[Table[plotEig[k,eigHTheta,θ[[k]],5],{k,1,Length[θ]}],AnimationRunning->False]

These are some numerical values which can be used to judge the stability

tabH=Table[PadRight[Sort[Select[eigHTheta[[k]],Arg[#1]>-1.6θ[[k]]&]],3],
{k,2,Length[eigHTheta]-1}];
TableForm[tabH]

Exercise

Write a program to compute the WKB approximation for the eigenvalues in the second example above.