Metastable states 1

Quantum Mechanics
                       K. Konishi, G. Paffuti

Study in the WKB approximation the metastable states for the radial potential

V[r]= (Subscript[V, 0]) (r^2)Exp[-r];   (Subscript[V, 0])= 7.5 .

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

Introduction

In this notebook metastable s-wave states in a model radial potential

V[r] = V0 (r^2)Exp[-r]

will be studied.

In the first part we shall perform a WKB analysis. while in the second part the Gamow - Siegert approach is applied. We take V0 = 7.5. We will use natural units m = ℏ = 1.

V0 = 7.5;
V[r_]:= V0 (r^2)Exp[-r];

Plot[V[r],{r,0,10},AxesLabel-> {"r","V[r]"}]

WKB analysis

For radial s-states the correct form of semiclassical quantization condition is (see text)

J = (1/π)(∫_0^(Subscript[x, 2]))p[x]dx = n + (3/4)

The potential has a maximum at r =2, (Subscript[V, M]) = 4 (Subscript[V, 0])(e^(-2)). An upper bound for J is obtainted for E = (Subscript[V, M])

(Subscript[J, max])=(1/π)(∫_0^2)(Sqrt[2((Subscript[V, M])- (Subscript[V, 0])(r^2)(e^(-r)))])dr ≃ 0.365 (Sqrt[(Subscript[V, 0])]) ≃1.;   3/4 < 1 < 7/4

then we expect at most 1 WKB state.

The basic functions need to solve this kind of problem have been introduced in notebooks [*].  They are given below with small modifications ().

Vmax = 4 V0/(E^2);
rmax = 2.;

getTurningPoints[ene_,V_]:= Module[{aL,aR},
aL = r/.FindRoot[ene -  V[r],{r,0.1rmax,0,rmax}];
aR = r/.FindRoot[ene -  V[r],{r,1.1rmax,rmax,50}];
{aL,aR}];

calcJ[ene_,V_]:= Module[{aL,aR,Jval},
{aL,aR}= getTurningPoints[ene,V];
Jval =(1/π)Re[NIntegrate[(Sqrt[2(ene - V[x])]),{x,0,aL}]]];

quantize[ene_?NumberQ,V_,n_]:=calcJ[ene,V] - (3/4) -n;

getEneWKB[n_,V_]:=  y/. FindRoot[quantize[y,V,n] ==0,{y,0.5Vmax,0.1Vmax,0.9Vmax}]

The lowest order WKB energy is

eneWKB0=getEneWKB[0,V]

Let us note that in the quadratic approximation near r = 0

V[r]∼ (1/2) 2(Subscript[V, 0]) (r^2);

The ground state in this (radial) harmonic oscillator would have an energy

(Subscript[E, HO]) ∼ (3/2)(Sqrt[2(Subscript[V, 0])]) ∼ 5.81.

The factor 3/2 = 1 + 1/2 is due to the fact that in radial coordinates the states, vanishing at r= 0, correspond to odd states of the analogous 1-d problem, as the reader must know.

The width of the state is approximatively given by

Γ ≃ (1/T)Exp[-2 (∫_(Subscript[a, L])^(Subscript[a, R]))p[x]dx] = (1/T)Exp[-2 (∫_(Subscript[a, L])^(Subscript[a, R]))(Sqrt[2(V[r]-E)])dx]

where T is the classical period of motion :

T = 2(∫_0^(Subscript[a, L]))((dx)/(p[x]))

getGamma[ene_,V_]:= Module[{aL,aR,T,K},
{aL,aR}= getTurningPoints[ene,V];
T=Re[NIntegrate[2./(Sqrt[2(ene -  V[r])]),{r,0,aL}]];
K =Re[Exp[-2 NIntegrate[(Sqrt[Abs[2(ene -  V[r])]]),{r,aL,aR}]]];
K/T]

The lowest order result is

ΓWKB0 = getGamma[eneWKB0,V]

First order result

At next order in WKB approximation the quantization condition becomes

J + δJ = n + (3/4)

where, in general

J = (1/π)(∫_(Subscript[x, 1])^(Subscript[x, 2]))p[x]dx ;     δJ = - ((ℏ^2)/(24π))(∂/(∂E))(∫_(Subscript[x, 1])^(Subscript[x, 2]))((V''[z])/(((2m(E-V[x])))^(1/2))) .

To compute the correction δJ a possible strategy is to create an interpolating function for the integral. We will limit to a small interval around a fixed energy, in the present case eneWKB0.

getList[V_]:= Module[{ene,Linteg,aL,aR},
ene=Range[0.10 Vmax,0.98V[rmax],0.05];
Linteg =Table[0.,{i,1,Length[ene]}];
Linteg=Table[
{aL,aR}=Re[getTurningPoints[ene[[i]],V]];
- (1./(24π))Re[NIntegrate[((V''[x])/(Sqrt[2(ene[[i]]-V[x])])),{x,0,aL},
PrecisionGoal-> 3]],
{i,1,Length[ene]}];
Return[Transpose[{ene,Linteg}]]]

The interpolation function is

funcδJ=Interpolation[getList[V]];

and the quantization condition change as

δJ[x_] = funcδJ'[x];
quantize[ene_?NumberQ,V_,n_]:=calcJ[ene,V]+δJ[ene] - (3/4) -n;

The improved energy and width are :

eneWKB1=getEneWKB[0,V]

ΓWKB1 = getGamma[eneWKB1,V]

Comparing with lowest order results

{eneWKB0,ΓWKB0}

we see a substantial change in Γ, as expected due to the exponential dipendence on energy.

Gamow - Siegert approach

In this approach we look for solutions with a radial asymptotic behavior Exp[i k r]/r  or Exp[i k r] for the reduced wave function. We shall work with the reduced radial wave function u[r]. The boundary condition u[0] = 0 will consitute an eigenvalue equation for k, with complex solutions. The with is related to k, in natural units, by

((k^2)/2)= E - i (Γ/2)

With W[r] = 2 V[r] the equation to be solved is

u''[r] + ((k^2)- W[r])u[r] = 0;               u[0]=0;   u[r](Underscript[→, r→ ∞]) Exp[i k r]

We solve two equations for r < xM and for r > xM, with relative boundary conditions, obtaining, for each k, two solutions uL and uR. The matching condition uL'/uL= uR'/uR is solved with FindRoot. The matching condition does not depend on normalizations, which are arbitrarily set to 1.

xL = 0; xR = 30; xM =1;

eqL[k_]:= {uL''[x] + (k^2) uL[x]==2 V[x]uL[x],uL[xL]==0,uL'[xL]==1}
eqR[k_]:= {uR''[x] + (k^2) uR[x]==2 V[x]uR[x],uR[xR]==Exp[I k xR],uR'[xR]== I k  Exp[I k xR]}

This is the matching condition :

fzero[k_?NumberQ]:= Module[{solL,solR,fuL,fuR,fduL,fduR},
solL = NDSolve[eqL[k],uL,{x,xL,xM}];
solR = NDSolve[eqR[k],uR,{x,xR,xM}];
fuL = uL[xM]/.solL[[1]]; fduL =  uL'[xM]/.solL[[1]];
fuR = uR[xM]/.solR[[1]]; fduR =  uR'[xM]/.solR[[1]];
fduL/fuL - fduR/fuR]

and this is the solution :

ksol=k/.FindRoot[fzero[k],{k, 2.5-  0.05I}]

Energy and width are given by

{eneGS,ΓGS}= {Re[0.5 (ksol^2)],-2 Im[0.5 (ksol^2)]}

To be compared with WKB results :

NumberForm[TableForm[Transpose[{{eneWKB0,eneWKB1,eneGS},{ΓWKB0,ΓWKB1,ΓGS}}],
TableHeadings-> {{"WKB 0","WKB 1","G.S."},{"  E","   Γ"}},
TableSpacing-> {2,3}],6]

Resonances above threshold

As will be seen in more detail on the chapter on metastable states, resonances can exist even above the classical threshold, here (Subscript[V, M]). Their meaning in terms of phase scattering will be explained in next notebook. Here we look for them by solving Gamow - Siegert equation.

We make a "scan" of k and collect different solutions.

kList = Range[2.5,4.1,0.1]-0.05I;

Off[FindRoot::"lstol"];
kvalues=(k/.FindRoot[fzero[k],{k,#1}])&/@ kList;
On[FindRoot::"lstol"];

The built-in Mathematica function (in Mathematica 6.0) Tally collect identical elements of a list with the specified test, here the difference between absolute values must be less then (10^(-3)).

kvalues=Sort[Tally[kvalues,Abs[#1-#2]≤ (10^(-4))&][[All,1]],Re[#1]< Re[#2]&]

Energies and width are given by

Transpose[{Re[0.5 (kvalues^2)],-2 Im[0.5 (kvalues^2)]}]

As the widths are very large it is difficult to assign a significance of resonance to all these values, in particular it one can hardly distinguish almost degenerate values, like E = 7.91 and E = 7.94.