Complex scaling 2

Quantum Mechanics
K. Konishi, G. Paffuti

Use the complex scaling methods to study resonances in the potential   (x^2)/2 - g (x^4)/2.

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

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 clockwise sense. H[θ] show 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 and WKB

The choice of the basis is quite arbitrary. In the following we will use two different basis: a Fourier basis and an harmonic oscillator basis

Fourier basis

We will 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].

Harmonic oscillator basis

The computational details are discussed in notebook [NB-10.10]. The function hermiteFunctions for fixed parity compute the matrices (p^2), (x^2), (x^4) from which the Hamiltonian can be written.

WKB calculations

To check our results we report here a slightly modified version of the routines in [NB-11.4]. We compute E-i Γ/2 at first order in WKB. The calling sequence is

getWKBEK[g, state]

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]

hermiteMatrices[dim_,parity_Integer:-1]:=
Module[{x2,p2,x4,x1}, 
Block[{i,n,shift},
If[parity == 0 || parity == 1,n =2i-2+parity; shift = 1;,n=i-1;shift =2;];x2=SparseArray[{{i_,i_}->(1/4)(2n+1),{i_,k_}/;k-i==shift->(1/2)(Sqrt[(n+1)(n+2)])},{dim,dim}];x2 =(x2+Transpose[x2]);
p2 =SparseArray[{{i_,i_}->(1/4)(2n+1),{i_,k_}/;k-i==shift->-(1/2)(Sqrt[(n+1)(n+2)])},{dim,dim}];
p2 = (p2+Transpose[p2]);
x4=SparseArray[{{i_,i_}->(3/8)(2(n^2)+2n+1),{i_,k_}/;k-i==shift->((4n+6)/4)(Sqrt[(n+1)(n+2)]),
{i_,k_}/;k-i==2shift->(1/4)(Sqrt[(n+1)(n+2)(n+3)(n+4)])},{dim,dim}];
x4 =  (x4+Transpose[x4]);
If[parity == 0 || parity == 1,Return[{p2,x2,x4}];];
x1=SparseArray[{{i_,k_}/;k-i==1->(Sqrt[(n+1)/2])},{dim,dim}];
x1 = (x1+Transpose[x1]);
{p2,x2,x4,x1}]];

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θ]}}]]]

WKB functions

funzEK[ene_,g_]:= Module[{aL,aR,Jval},
aL = (Sqrt[ ((1 -(Sqrt[ 1-8 ene  g]))/(2g))]);aR = (Sqrt[ ((1 +(Sqrt[ 1-8 ene  g]))/(2g))]);
Jval =2(Sqrt[g]) (aR/(3π))(((aR^2)+(aL^2)) EllipticE[((aL^2)/(aR^2))] -  ((aR^2)-(aL^2)) EllipticK[((aL^2)/(aR^2))]); Jval]

corrEK[ene_,g_]:=((Sqrt[2]) g ((-1+24 ene g) (1-4 ene g+(Sqrt[1-8 ene g])) EllipticE[((1-(Sqrt[1-8 ene g]))/(1+(Sqrt[1-8 ene g])))]+(1+(Sqrt[1-8 ene g])+4 ene g (-5+24 ene g-4 (Sqrt[1-8 ene g]))) EllipticK[((1-(Sqrt[1-8 ene g]))/(1+(Sqrt[1-8 ene g])))]))/(3 (-1+8 ene g) (-1+(Sqrt[1-8 ene g])) (((1+(Sqrt[1-8 ene g])))^(5/2)) π);

fzeroEK[x_?NumberQ,g_,state_]:=(funzEK[x,g] +corrEK[x,g]- state-(1/2) );

getWKBEK[g_,state_]:=Module[{Vmax,TperEK,KbarEK,EneEK,Γ2EK},
Vmax = ((μ^2)/(8g));
EneEK = Re[x/.FindRoot[fzeroEK[x,g,state]==0,{x,Vmax/2,$MachineEpsilon,Vmax-$MachineEpsilon}]];
aR =(Sqrt[ ((1 + Sqrt[1-8EneEK g])/(2g))]);aL = (Sqrt[ ((μ - Sqrt[1-8EneEK g])/(2g))]);
TperEK = Re[NIntegrate[2./(Sqrt[2 EneEK -  (x^2)+ g (x^4)]),{x,-aL,aL}]];
KbarEK =Re[Exp[-2 NIntegrate[(Sqrt[Abs[2 EneEK -  (x^2)+ g (x^4)]]),{x,aL,aR}]]];
Γ2EK = (1/TperEK)KbarEK;  (* Note! this is Γ/2 *)
Return[EneEK - I Γ2EK ]]

Metastable states

The model

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

μ  = 1 ; g = 0.03;parita=0;   (* g è una variabile globale *)
V[x_]:= (1/2)μ x^2 - (1/2)g x^4
 Plot[V[x],{x,-1.1(Sqrt[(μ/g)]),1.1(Sqrt[(μ/g)])},AxesLabel-> {"x","V[x]"}]

Fixed resonances

This potential is unbounded below. In the text it has been shown that in this case one do not expect a continuous rotation of eigenvalues with θ, but instead only fixed resonances.

Let us check this point by using the harmonic oscillator basis. Only states with fixed (even) parity are selected for efficiency. For details on the procedure and the functions see the notebook [NB-10.12].

parity = 0;dim = 100; 
{p2,x2,x4} = hermiteMatrices[dim,parity];

θ = (π/48)Range[11];
Φ =Exp[2 I θ];

eigHTheta=Table[Sort[Eigenvalues[(p2/(2Φ[[i]]))+ ((Φ[[i]] x2)/2) - ((g ((Φ[[i]])^2)x4)/2)],Abs[#1]<Abs[#2]&],{i,1,Length[θ]}];//Timing

The following animation shows clearly that the first eigenvalues are fixed with θ. Large eigenvalues behavior depend mainly on the finite size of the basis, as the reader can check by performing a computation with larger dimensions.

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

The degree of stability is best appreciated by a table of the first four eigenvalues as θ varies

TableForm[eigHTheta[[All,Range[4]]],TableSpacing-> {3,3}]

For comparison these are the WKB results

getWKBEK[g,#1]&/@{0,2,4}

Fourier basis

As a final check we repeat the computation by using a Fourier basis, as in notebook [NB-13.1].

This is the grid

dim = 100; Npt=1025;
L=35; x0 = -L/2;

and this is the computation

eigHFourier=Table[
Sort[Eigenvalues[calcH[V,x0,L,dim,Npt,θ[[i]]]]
,Abs[#1]<Abs[#2]&],{i,1,Length[θ]}];//Timing

To compare we select eigenvalues 1,3,5... as in the previous case we computed only even parity states:

TableForm[eigHFourier[[All,Range[1,7,2]]],TableSpacing-> {3,3}]

There is a perfect agreement.

Exercise

The reader can perform an analysis with higher dimensional basis and varying g.