Metastable states and  scattering phases

Quantum Mechanics
                       K. Konishi, G. Paffuti

Compute the s-wave scattering phase as a function of  energy for the radial potential

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

Compare the results with the WKB approximation.

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

Introduction

The existence of metastable states can be investigated by several means: WKB analysis, Gamow-Siegert equation etc. One of the methods more closely connected to real experiments is the study of scattering phases.

Let us consider the s-wave scattering. The Schrödinger equation for reduced s-wave function u[r] is (in natural units ℏ = m = 1)

u''[r] +( (k^2) - 2 V[r]) u[r] = 0

k is related to energy of the incident particle by

E = ((k^2)/2)

The asymptotic form of a solution of () is, apart form an arbitrary normalization factor,

u[r] (Underscript[→, r → ∞]) Sin[kr + δ]

The phases δ are fixed by imposing regularity at the origin, u[0] = 0. They will be functions of k, δ[k], or of energy E.

In the text it has been shown that around a resonance state the derivative with respect to energy of δ[E] has a characteristic Breit-Wigner shape:

δ'[E] ≃ ((Γ/2)/((((E-(Subscript[E, 0])))^2)+ ((Γ^2)/4)))

(Subscript[E, 0]) and Γ beeing energy and width of the metastable state. Equation () allows a determination of metastable states parameters.

In this notebook we shall first compute numerically the phase shifts δ for a potential

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

Then we shall use the WKB approximation for phase shifts and compare the results. The results are finally compared with analogous results obtained by the WKB bound state techniques and Gamow Siegert approach in notebook [NB-11.6].

Solution of the Schrödinger equation

The equation to be solved is

u''[r] +( (k^2) - 2V[r]) u[r] = 0;  u[0]=0;  u[r] (Underscript[→, r → ∞]) Sin[kr + δ]

As the normalization is arbitray we can just solve this equation by initial boundary conditions u[0]=0, u'[0] = 1. Using the asymptotic condition we will have

Exp[i (k r + δ)] = Cos[k r + δ]+i Sin[k r+δ] = (1/(i k))u' + u

then

Exp[i δ] = ((1/(i k))u' + u )Exp[-i k r]

As the overall (real) normalization is arbitrary δ is directly extracted from the argument of the above combination.

δ = Arg[((1/(i k))u' + u )Exp[-i k r]]

This method has a difficulty: δ[E] is a continuous function of energy, while extracting its value with () will give values mod(2π). The problem is easily solved by noticing that on general grounds δ→0 as k→0, this fixes the origin of the phases and subsequent values can be recovered by continuity (see below).

Equations

We first write the potential, choose a large r value to impose asymptotic conditions, and write the differential equation for a generic momentum k

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

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

xL = 0; xR = 40;

eq[k_]:= {u''[x] + (k^2) u[x]==2 V[x]u[x],u[xL]==0,u'[xL]==1}

Energy

There are two delicate regions, k∼0, where initial condition must be fixed, and regions around resonances. We know from previous notebooks that the system has a resonance aroun E = 3.4, so we will build a grid of values finer in these regions.

This is a possibility, the reader can change parameters at his will

Emin = 0.00001;Emax = 120; dE = 0.005; dE2 = 1;
e1 = Range[Emin,0.5,dE];
e2 = Range[Max[e1]+dE,3-dE,0.1dE2];
e3 = Range[3,4,dE];
e4 = Range[4+dE,6,3 dE];
e5 = Range[Max[e4]+dE2,Emax,dE2];
ene = Join[e1,e2,e3,e4,e5];
k =  Sqrt[2*ene];

Here we compute phases

δ = Table[sol = NDSolve[eq[k[[i]]],u,{x,xL,xR}];
Arg[((1/k[[i]]*u'[xR]/.sol[[1]] )+
 I* (u[xR]/.sol[[1]]))Exp[-I k[[i]] xR]],{i,1,Length[ene]}];

As predicted the phases suffer from 2 π discontinuities

ListPlot[Transpose[{ene,δ}],PlotRange-> All]

The following simple function smoothPhase reabsorbs 2π jumps by comparing subsequents levels. As δ change rapidly around a resonance it is essential that the grid is sufficiently finer to avoid to miss discontinuities.

smoothPhase[phase_]:= Module[{disc,jump,outphase},
disc = 0*phase;
(* This detect the jumps *)
Do[ 
jump = phase[[i]] - phase[[i-1]];
disc[[i]]=Which[jump < -1.7*Pi,1,jump > 1.7*Pi,-1,True,0];,
{i,2,Length[phase]}];
(* This add 2π for each jump *)
outphase = Drop[FoldList[Plus,0,disc],1]*2Pi+phase]

The instruction FoldList add recursively the jumps.

δ = smoothPhase[δ];

Here is the plot of δ[k]

plot1=ListPlot[Transpose[{k,δ }],PlotRange-> All,Joined-> True,
AxesLabel-> {"k","δ(k)"}]

We see a sharp jump, this will correspond to the resonance.

Resonances

We write an interpolating function to make derivatives and explore the range around the jump:

rec=Interpolation[Transpose[{ene,δ}]];
Plot[rec'[x],{x,3.2,3.6},PlotRange-> All,AxesLabel-> {"E","δ'(E)"}]

We see a perfect Lorentzian shape. To extract the resonance parameters we do a fit. The theoretical curve

δ'[E] ≃ ((Γ/2)/((((E-(Subscript[E, 0])))^2)+ ((Γ^2)/4)))

has unit integral. We can perform the fit with or without this constraint, to check consistency.

We subtract the "background" with the first value of δ'[E] in the small energy range around the resonance.

ptEne = Range[3,4,dE];
ptDelta = rec'[ptEne];
points = Transpose[{ptEne,ptDelta-ptDelta[[1]]}];

This is the fit without the constraint

FindFit[points,pb/2/(pΓ^2/4+ (x-pE0)^2),{pb,pE0,pΓ},x]

This is the fit with the constraint

fit=FindFit[points,pΓ/2/(pΓ^2/4+ (x-pE0)^2),{pE0,pΓ},x]

{E0,Γ}= {pE0,pΓ}/.fit

These value have to be compared with results obtained with WKB bound state technique and Gamow Siegert approach :

{eneGS,ΓGS}={3.42639,0.02555};
{eneWKB1,ΓWKB1}= {3.42489,0.02647};

NumberForm[TableForm[Transpose[{{eneWKB1,eneGS,E0},{ΓWKB1,ΓGS,Γ}}],
TableHeadings-> {{"WKB 1","G.S.","δ"},{"  E","   Γ"}},
TableSpacing-> {2,3}],5]

The agreement is evident.

WKB phases

The second part of the notebook is devoted to a computation of the scattering phases within WKB approximation. We will use two different WKB approximations, the usual one and the Langer modification for radial problems.

The free s-wave solution obeying u[0] = 0, is Sin[k r]. δ is the phase difference, at point r, between the scattering solution and the free solution.

In the WKB approximation the solution in the classically allowed region has the form

u∼Cos[S[b,r] - (π/4)]

for a usual turning point, while, if the turning point is the barrier r=0

u∼Sin[S[0,r]]= Cos[S[0,r] - (π/2)]

where b is the classical turning point, where E  = V. For a free particle the phase for s -waves is always given by

Sin[k r] = Cos[k r - (π/2)] ⇒ (Subscript[φ, 0])[r] = kr - (π/2)

We can write at once δ by taking the phase difference. In our case

δ(E)= S[0,r] - k r = (∫_0^∞)((Sqrt[(k^2)- 2 V[r]]) - k)dr ;   E > (Subscript[V, max])
δ(E) = (S[b,r] - (π/4)) - (k r - (π/2)) =
                            (∫_b^∞)((Sqrt[(k^2)- 2 V[r]]) - k)dr  + (π/4)- k b;    E < (Subscript[V, max])

rmax = 2.;
k0 = (Sqrt[2 V[rmax]]);

This compute δ

getδWKB[k_]:= Module[{b},
b = If[k>k0,0,x/. FindRoot[k^2 - 2V[x] == 0,
{x,20,rmax+1./10^6,200}]];
If[k>k0,
Chop[NIntegrate[(Sqrt[k^2 - 2 V[x] ])  - k,{x,0,Infinity}]],
b= x/. FindRoot[k^2 - 2V[x] == 0,{x,20,rmax+(10^(-6)),200}];
Chop[NIntegrate[(Sqrt[k^2 - 2 V[x] ])  - k,{x,b,Infinity}]]+Pi/4- k b]]

kList=Union[Range[0.1,2,0.1],Range[2,12,0.5]];

here is the computation

δWKB=getδWKB[#1]&/@kList;

and the comparison with numerical results :

Show[ListPlot[Transpose[{kList,δWKB}],PlotStyle-> Red],plot1]

Langer correction

The Langer correction has the net effect of a substitution of centrifugal term in the potential

((L(L+1))/(2 (r^2))) → (1/(2(r^2)))(((L+(1/2)))^2)

In this version even s - waves, with L = 0, have a 1/(r^2) potential. In this case the barrier is never touched, even for free particles, and the phase shift is simply (for s-waves):

δ = (∫_b^∞)(Sqrt[(k^2)- 2 V[r]- (1/(4 (r^2)))])dr - (∫_(Subscript[b, 0])^∞)(Sqrt[(k^2)- (1/(4 (r^2)))])dr

b and (Subscript[b, 0])are the turning points for intercting and free particle. Clearly (Subscript[b, 0])= 1/(2k), b is found using FindRoot.

The effective potential in this case is a deformation of the original one:

Plot[V[r]+(1/(4 (r^2))),{r,0.1,10},AxesLabel-> {"r","V[r]"}]

For a repulsive potential surely b > (Subscript[b, 0])and using

(Sqrt[(k^2)- 2 V[r]- (1/(4 (r^2)))]) - (Sqrt[(k^2)- (1/(4 (r^2)))]) =  ((-2 V[r])/((Sqrt[(k^2)- 2 V[r]- (1/(4 (r^2)))])+(Sqrt[(k^2)- (1/(4 (r^2)))])))

the expression for δ can be written in the form

δ = (∫_b^∞)((-2 V[r])/((Sqrt[(k^2)- 2 V[r]- (1/(4 (r^2)))])+(Sqrt[(k^2)- (1/(4 (r^2)))])))dr -  (∫_(Subscript[b, 0])^b)(Sqrt[(k^2)- (1/(4 (r^2)))])dr

We can easily compute δ' s also in this case

getδWKBLanger[k_]:= Module[{b,b0,I1,I2},
b0 = 1/(2 k);
b = x/.FindRoot[k^2 - 2V[x] -(1/(4 (x^2)))== 0,{x,1/(2k),200}];

I1=NIntegrate[-2 ((V[x])/((Sqrt[k^2 - 2 V[x]-(1/(4(x^2))) ])+(Sqrt[k^2 -(1/(4(x^2))) ])))  ,{x,b,Infinity}];
 I2 = NIntegrate[ (Sqrt[k^2 -(1/(4(x^2))) ])  ,{x,(1/(2k)),b}]; 
Re[I1-I2]]

δWKBLanger=getδWKBLanger[#1]&/@kList;

Th ecomparison with numerical data is again satisfactory

Show[ListPlot[Transpose[{kList,δWKBLanger}],PlotStyle-> Red],plot1]