Density of states and scattering phases 1

Quantum Mechanics
 K. Konishi, G. Paffuti

Compute the resonance parameters from scattering phases in a potential (Subscript[V, 0])(r^2)Exp[-r].

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

Introduction

Summary of the theory

In the text it has been pointed out that the spectral density (at fixed angular momentum L) is related to the scattering phases by

(Subscript[g, L])[E] = (Subscript[g, 0])[E] + (1/π)((Subscript[dδ, L])/dE).

Where (Subscript[g, 0]) is the free particle contribution. In particular near a resonance

(Subscript[g, L])[E]≃ (1/π) ((Γ/2)/((Γ^2)/4 + (((E-E0))^2))).

We used the computation of scattering phases to locate resonances in our study of WKB approximation, see notebook [NB-11.7]. Here we extend that study by considering higher angular momentum states.

The definition of phases

Phase scattering is defined (for a radial potential) by the asymptotic behavior

u(r) ∼ Sin[kr - L(π/2) + δ]

of the regular solution of Schrödinger equation.

More accurately if the potential V[r] goes to 0 sufficiently fast as r→∞ the asymptotic form of the Schrödinger equation is

0=(d/(d (r^(2 ))))u[r] + (k^2)u[r] - ((L(L+1))/(r^2))u[r] = (d/(d (r^(2 ))))u[r] + (k^2)u[r] - (((((L+(1/2)))^2)-(1/4))/(r^2))u[r].

The general solution of this equation is

u[r]= (Subscript[A, 1])(Sqrt[r]) (Subscript[J, L+1/2])[k r]+ (Subscript[A, 2])(Sqrt[r]) (Subscript[Y, L+1/2])[k r].

J and Y are Bessel functions of the first and second kind. From their asymptotic behavior (⋁ = L+1/2)

(Subscript[J, ν])(x) ∼ (Sqrt[(2/(π x))])Cos[x - ν(π/2)-(π/4)]  = (Sqrt[(2/(π x))])Sin[x - L(π/2)]   ;   
(Subscript[Y, ν])(x) ∼ (Sqrt[(2/(π x))])Sin[x - ν(π/2)-(π/4)]  = - (Sqrt[(2/(π x))])Cos[x - L(π/2)]  ;

We see that

(Subscript[A, 1])= A Cos[δ]; (Subscript[A, 2])= A Sin[δ];          Tan[δ] =- (Subscript[A, 2])/(Subscript[A, 1]).

Computation of δ

We will use the more accurate formula (). To compute (Subscript[A, 1])and (Subscript[A, 2])we solve the Schrödinger equation for the regular solution (see below) and consider two points (Subscript[r, 1])and (Subscript[r, 2]) in the asymptotic region, inserting the computed value of u[r] at this points in () we get the coefficients. The points (Subscript[r, 1]),(Subscript[r, 2])must lie in the classically allowed region and for sure r > (L+1/2)/k.

The regular solution

For L = 0 and a regular potential one can simply take u[0]=0, u'[0] = 1, this give the regular solution, the arbitrary normalization do not enter in the definition of δ as only the ratio of asymptotic coefficients  (Subscript[A, 1])and (Subscript[A, 2]) matters.

If L > 1 or if the potential is singular at the origin some care must be taken as numerically routines like NDSolve do not accept a singular point as a boundary. We assume now that the potential is regular, with a Taylor expansion at the origin

V[r] ≃ V0 + V1 r + V2 (1/2) (r^2) + ...

Near the origin a regular wave function behaves like u[r]∼ (r^(L+1)) (see text). Next order terms will be (the normalization being arbitrary)

u[r]= (r^(L+1))(1 + (Subscript[α, 1]) r+ (Subscript[α, 2])((r^2)/2)+...).

Inserting the Taylor expansion of the potential in the Schrödinger equation (we use m=ℏ=1)

(d/(d (r^(2 ))))u[r] + ((k^2)-2 V)u[r] - ((L(L+1))/(r^2))u[r]=0,

one easily obtains (we leave the exercise to the reader)

u[r] ∼ (r^(ℓ+1))  (1 - (((k^2)+2V0)/(2 (3+2 ℓ)))(r^2)   -(V1/(3 (2+ℓ))) (r^3)    + (((k^4)+4 (k^2) V0+4 (V0^2)-6 V2-4 ℓ V2)/(8 (3+2 ℓ) (5+2 ℓ))) ((r )^4)  +... )

This fix the behavior of u[r] near the origin and allow to impose the boundary condition at a point (Subscript[r, 0])near 0 but not coincident with the origin.

Functions used in this notebook

The function psi0

psi0[dV_,L, k, r]

write the expansion () for a given triplet dV = {V[0], V'[0], V''[0]}.

The function getPhase solve the linear 2×2 linear system for computing (Subscript[A, 1])and (Subscript[A, 2]) in (). As a result it give the phase δ, computed with ArcTan

getPhase[val1, val2, x1, x2, L, k]

To smoothen the discontinuity introduced by of ArcTan we use the function smoothPhase, introduced in notebook [NB-11.7].

smoothPhase[phase]

Code of the functions

psi0[dV_,ℓ_,k_,r_]:=Module[{V0,V1,V2},
{V0,V1,V2}=dV;
ℓ!(((k r))^(ℓ+1))(1-((((k^2)+2V0)*(r^2))/(2 (3+2 ℓ))) -(V1/(3 (2+ℓ))) (r^3)    + (((k^4)+4 (k^2) V0+4 (V0^2)-6 V2-4 ℓ V2)/(8 (3+2 ℓ) (5+2 ℓ))) ((r )^4) )];

getPhase[val1_,val2_,x1_,x2_,L_,q_]:= Module[{c11,c12,c21,c22,a1,a2,d1,d2,ph},
c11 = (Sqrt[x1]) BesselJ[L+(1/2),q x1];c21 = (Sqrt[x1]) BesselY[L+(1/2),q x1];
c12 = (Sqrt[x2]) BesselJ[L+(1/2),q x2];c22 = (Sqrt[x2]) BesselY[L+(1/2),q x2];
{d1,d2}={a1,a2}/.First[NSolve[{a1*c11 + a2 *c21 == val1,a1*c12 + a2 *c22 == val2},{a1,a2}]];
ph =  ArcTan[d1,-d2]]

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]

Computation of spectral density and phase shift

The model

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

We know from several sources, see notebooks [NB-13.1 , NB-13.4], that this model has the lowest resonance at

{Eth,Γth} = {3.42639,0.025549};

At the origin the first derivatives of V[r] are

dVnear0 = {V[0],V'[0],V''[0]};

To study the scattering problem we define the effective potential

Veff[L_,r_]:= 2*(V0 (r^2)(E^(-r)) +  ((L(L+1))/(2  (r^2))));

We will compute phase scattering for a set of energies running from a low value,  0.00005, to about 100. The grid is not equal spaced as we will be interested in energies around E ∼ 3

dErun[x_]:= Which[ x≤ 0.0015,0.00005,
 ( 3.35≤ x ≤ 3.5), 0.005,
x < 8, 0.05, x≥ 8,1.0];
Emin = 0.00005; erun = Emin;
ene = Prepend[Table[erun = erun + dErun[erun],{320}],Emin];

To solve the Schrödinger equation we fix the left extremum xL, where the regularity condition is imposed, and the right boundary xR ∼ 40. The asymptotic coefficients are computed from the value of the regular solution  at xR and xR-0.5.

k =  Sqrt[2*ene];
phase = 0*ene;
xL = 0.1;
xR = 40.+0.*k;

In the following a fixed value for xR will be used, an alternative choice is to use a k-dependent point as in the commented line below:

(* xR = Max[30.,1.5*(ℓ+1/2)/# ]&/@ k; *)

The equation to be solved is

eq[L_,k_]:= {u''[x] + (k^2) u[x]==Veff[L,x]u[x],u[xL]==psi0[dVnear0,L,k,xL],(u'[xL]==(∂_r) psi0[dVnear0,L,k,r]/.r-> xL)}

S wave

We start with the simple s - wave scattering

Lmom = 0;

The solution for all values of the energy parameter is stored in the vector sol:

sol = Table[ NDSolve[eq[Lmom,k[[i]]],u,{x,xL,xR[[i]]}],{i,1,Length[ene]}];//Timing

Here we compute the asymptotic phase with getPhase and eliminate discontinuities induced by ArcTan with smoothPhase. The result is stored in the vector δ

phase = Table[{x1,x2}= {xR[[i]],xR[[i]]-0.5};
{val1,val2} = { (u[x1]/.sol[[i]][[1]]), (u[x2]/.sol[[i]][[1]])};
getPhase[val1,val2,x1,x2,Lmom,k[[i]]],{i,1,Length[ene]}];
delta = smoothPhase[phase];

This show the function δ[E]

ListPlot[Transpose[{ene,delta}],Joined-> True,PlotRange-> {{0,25},{-12,0}},AxesLabel-> {"E","δ"}]

The sharp variation near E=3 is a signal of the resonance.

To study the resonance region in more detail we define an interpolating function fdelta

fdelta=Interpolation[Transpose[{ene,delta}]];

The derivative δ'[E] gives directly the resonance contribution to the spectral density. Here is an overall graph and an expansion of the zone with the relatively sharp resonance:

ymax =Ceiling[Max[ (2/π)fdelta'[ene]]]/2;
{xmin,xmax} = {3.3,3.6};
Row[ListPlot[Transpose[{ene,(1/π)fdelta'[ene]}],Joined-> True,PlotRange-> {#1,{-0.5,ymax}},AxesLabel-> {"E","δg(E)"},ImageSize-> 4 72]&/@{{0.1,Max[ene]/2},{xmin,xmax}}]

The resonance parameter can be evaluated by a fit on data obtained by the interpolation function. We choose a grid ptEne to generate data:

ptEne = Range[xmin,xmax,0.005];
ptDelta = fdelta'[ptEne]/π;
data = Transpose[{ptEne,ptDelta }];
fit=FindFit[data,b/2/(Γ^2/4+ (x-E0)^2),{b,E0,Γ},x];
{"from δ: ",{E0,Γ}/.fit," expected: ",{Eth,Γth}}

The agreement with the known values is excellent. Here we show the relative plot

Plot[b/2/(Γ^2/4+ (x-E0)^2)/.fit,{x,xmin,xmax},
PlotLabel-> "L = " <> ToString[Lmom],AxesLabel-> {"E","δg(E)"},
PlotRange-> {{xmin,xmax},{0,ymax}},
Epilog->{PointSize[0.015],Point/@Transpose[{ptEne,ptDelta}]}]

P wave

The computation is identical.

Lmom = 1;

The effective potential has the form

Plot[Veff[Lmom,r],{r,0,10},PlotRange->{All,{0,Ceiling[Max[Veff[Lmom,Range[0.5,4,0.1]]]]}},
AxesLabel-> {"r","(Subscript[V, eff])(r)"}]

The solution and the computation of δ proceed as before:

sol = Table[ NDSolve[eq[Lmom,k[[i]]],u,{x,xL,xR[[i]]}],{i,1,Length[ene]}];

phase = Table[{x1,x2}= {xR[[i]],xR[[i]]-0.5};
{val1,val2} = { (u[x1]/.sol[[i]][[1]]), (u[x2]/.sol[[i]][[1]])};
getPhase[val1,val2,x1,x2,Lmom,k[[i]]],{i,1,Length[ene]}];
delta = smoothPhase[phase];

This show the function δ[E]

ListPlot[Transpose[{ene,delta}],Joined-> True,PlotRange-> {{0,25},{-12,0}},AxesLabel-> {"E","δ"}]

fdelta=Interpolation[Transpose[{ene,delta}]];

ymax =1.5;
{xmin,xmax} = {4.2,5.5};
Row[ListPlot[Transpose[{ene,(1/π)fdelta'[ene]}],Joined-> True,PlotRange-> {#1,{-0.5,1.5}},AxesLabel-> {"E","δg(E)"},ImageSize-> 4 72]&/@{{0.1,Max[ene]/2},{xmin,xmax}}]

A quite broad structure is present near E ∼ 4.6 and can be analyzed in terms of a resonance.

ptEne = Range[xmin,xmax,0.05];
ptDelta = fdelta'[ptEne]/π;
data = Transpose[{ptEne,ptDelta - Last[ptDelta] }];
fit=FindFit[data, b/2/(Γ^2/4+ (x-E0)^2),{b,E0,Γ},x];
{"{E, Γ) from δ: ",{E0,Abs[Γ]}/.fit}

Plot[Last[ptDelta]+b/2/(Γ^2/4+ (x-E0)^2)/.fit,{x,xmin,xmax},
PlotLabel-> "L = " <> ToString[Lmom],AxesLabel-> {"E","δg(E)"},
PlotRange-> {{xmin,xmax},{0,ymax}},
Epilog->{PointSize[0.015],Point/@Transpose[{ptEne,ptDelta}]}]

As expected the resonance structure is much less pronounced than in the previous case.