Scattering 4

Quantum Mechanics
                    K. Konishi, G. Paffuti

Solve numerically the scattering problem for a  spherical potential V[r]=(Subscript[V, 0])(r^2)Exp[-r].  In particular make a density plot for the wave function for resonant scattering. Verify the eikonal approximation.

NOTE 1: To avoid interference with other notebooks it is recommended to quit and start again the Kernel.

Global variables and initialization

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

Theory and method of solution

Notations

The notations and normalization are those of [NB-16.1]:

ψ = (∑_L) (Subscript[A, L]) LegendreP[Cos[θ]] R[L,k,r].

(Subscript[A, L]) = (2L+1)(i^L) (e^(i (Subscript[δ, L]))) ;

Ras[L,k,r] = (Cos[δ] sJ[L,k r]- Sin[δ]sY[L,k r]).

Note: In this notebook we avoid the factor 2k in (Subscript[A, L]) and Ras to simplify  some formulas below. The factor was used in notebooks [NB-16.2] [NB-16.3],

ψ ∼ Exp[i k z] + f[θ] ((e^(i k r))/r);           f[θ]= (1/(2i k))(∑_L) (2L+1) ((e^(2i (Subscript[δ, L])))-1)(Subscript[P, L])[Cos[θ]]

The problem

The potential in our problem is

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

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

Phase shifts

We follow the general strategy outlined in [NB-16.1]. This particular problem has been studied with the same techniques in [NB-13.5].

The radial equation for the reduced L - wave function is (we use units ℏ = m = 1):

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

As V[r] goes rapidly to zero as r grows the u[r] asymptotic behavior is

(Subscript[A, 1])(Subscript[J, ℓ+1/2])(k r)+(Subscript[A, 2]) (Sqrt[r]) (Subscript[Y, ℓ+1/2])(kr).

From the asymptotic behavior of Bessel functions it is easy to see ( [NB-16.1]):

Tan[δ] = ((Subscript[A, 2])/(Subscript[A, 1])).

A numerical solution of () allows a computation of δ. The numerical solution is obtained with ().

Numerical computation

Solution of the radial equations

The equations will be solved between xL, xR, with xL∼0 and xR in the classical allowed region, xR > L/k . Two points near xR are used to compute the coefficients (Subscript[A, 1]),(Subscript[A, 2])in () and compute the phase shifts. To avoid discontinuities introduced by ArcTan we smooth the results with a function smoothPhase, the procedure is explained and used in notebook [NB-13.5] as an example. The phases are computed by getPhase.

The point xL is chosen greater than zero to avoid singularities in the effective potential. As explained in notebook [NB-13.5] the regular solution at small r has the behavior

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)  +... )

where

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

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. The small r function is constructed by the function psi0[L, k, r].

Lists of spherical Bessel functions and Legendre polynomials

To write a finite numerical expression we put a limit on our angular momenta, maxL, and verify the consistency of the procedure.

We will have to compute sequences of Legendre polynomials and spherical Bessel functions. This can be done by using Mathematica built in functions, but it turn out to be more efficient to use recursion relation to produce the whole list. The function

makeList[maxL]

define a sequence of functions sJ, sY, LegendreP, sH in the form of vectors  jL, yL, legP and hL.

For spherical Bessel function the recursion relation is

f[n-1, x] + f[n+1, x] = (2n+1)/ x   f[n,x];

while for Legendre polynomials

(n + 1) P[n + 1, x] - (2 n + 1) x P[n, x] + n P[n - 1, x] = 0 .

From numerical analysis it is known that stability suggest that for sJ functions the sequence has to be used in a decreasing order, see ref.[1].

The function makeδAS[maxL,k,a]

{δ,𝒜,S} = makeδAS[maxL,k,{a,β}]

compute the vectors of phase scattering and of factors

(Subscript[S, L])= (i^L)(((2L+1))/2)*(Exp[2i δ]-1);      (Subscript[𝒜, L])= (2L+1)((Exp[i δ]*(Cos[δ] jL[k a] - Sin[δ] yL[k a]))/(jL[q a])).

The coefficients 𝒜 are the coefficients connecting external and internal solution at r=a, i.e. are (Subscript[A, L])(Subscript[C, L]) in the above notation. It uses explicitly the recursive definition of derivatives to build a list for  δ's.

Note :

The names jL, yL, legP and hL  are used as global names. The function makeδAS must be used after a call to makeList. sJ and sY are used as alias for BesselFunctionJ etc. and have not to be changed by the user.

The names sJ, sY, sJ', sY'  are used as global names are not to be changed by the user.

Derivative and spherical Bessel functions

We will use the following relations for the derivatives of sj and sY:

sJ'[L,x]= ((- 1)/(2x))sJ[L,x] + (1/2) (sJ[L-1,x]-sJ[L+1,x]);
sY'[L,x]=- ((- 1)/(2x)) sY[L,x]+ (1/2) (sY[L-1,x]-sY[L+1,x]);

The relation is checked below. The derivative will be named dsJ, dsY.

Code for lists and functions

Lists

makeListJ[maxL_,x_]:= If[ x== 0.,Table[If[m==0,1.,0.],{m,0,maxL}],
Module[{JList=Table[1.,{maxL+1}]},
JList[[maxL+1]]= SphericalBesselJ[maxL,x];  (* (Sqrt[(π/(2x))])BesselJ[maxL+(1/2),x]; *)
JList[[maxL]]= SphericalBesselJ[maxL-1,x];(* (Sqrt[(π/(2x))])BesselJ[maxL-(1/2),x]; *)
(* i = L+1 *)
Do[ JList[[i]] = -JList[[i+2]]+(2 i + 1)JList[[i+1]]/x;,{i,maxL-1,1,-1}];JList]];

makeListY[maxL_,x_]:= 
Module[{YList=Table[1.,{maxL+1}]},
YList[[1]] =SphericalBesselY[0,x]; (* (Sqrt[(π/(2x))])BesselY[(1/2),x];*)
YList[[2]] =SphericalBesselY[1,x]; (* (Sqrt[(π/(2x))])BesselY[(3/2),x];*)
(* i = λ+1 *)
Do[ YList[[i]] = -YList[[i-2]]+(2 i -3)YList[[i-1]]/x;,{i,3,maxL+1}];YList];

makeListlegP[maxL_,θ_]:= Module[{t,legP = Table[1,{maxL+1}]},
 t = Cos[θ]; 
legP[[2]]=t;
Do[legP[[λ+2]] = (2λ+1) t/(λ+1) legP[[λ+1]] - λ /(λ+1)legP[[λ]];,{λ,1,maxL-1}];
legP]

makeList[maxL_]:= (
jL[x_]:= makeListJ[maxL,x];
yL[x_]:= makeListY[maxL,x];
hL[x_]:= jL[x]+ I yL[x];
legP[θ_]:= makeListlegP[maxL,θ];)

Radial equation

psi0[ℓ_,k_,r_]:=ℓ!(((k r))^(ℓ+1))(1-((((k^2)+2V[0])*(r^2))/(2 (3+2 ℓ))) -((V'[0])/(3 (2+ℓ))) (r^3)    + (((k^4)+4 (k^2) V[0]+4 ((V[0])^2)-6 V''[0]-4 ℓ V''[0])/(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]

Spherical functions

sJ = SphericalBesselJ;
sY = SphericalBesselY;
sH = SphericalHankelH1;
dsJ[L_,x_]:= ((- 1)/(2x))sJ[L,x] + (1/2) (sJ[L-1,x]-sJ[L+1,x]);
dsY[L_,x_]:= ((- 1)/(2x))sY[L,x] + (1/2) (sY[L-1,x]-sY[L+1,x])

Low energy scattering

For completeness and to check our procedure we repeat briefly the computations of notebook [NB-13.5].

An effective potential is defined by (the factor 2 is to simplify equations below)

Veff[ℓ_,r_]:= 2*(V0 (r^2)(E^(-r)) +  ((ℓ(ℓ+1))/(2  (r^2))));
xL:= 0.01;

The radial equation to be solved is

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

The boundary point for small r is xL (= 0.1) and there the boundary conditions are imposed using the analytic form ψ0 computed above, eq.(1.).

We will solve the equation for a set of energies:

Clear[L];
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];
k =  Sqrt[2*ene];

The right bound for the computation is defined by (L has been cleared above, so this definition will work for all L)

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

L = 0

Here we solve the radial equation for all k :

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

We compute the phase shifts

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,L,k[[i]]],{i,1,Length[ene]}];
delta = smoothPhase[phase];

and plot them

ListPlot[Transpose[{ene,delta}],Joined-> True,PlotRange-> {{0,25},{-12,4}}]

We can check the (k^(2L+1)) behavior of (Subscript[δ, L])by fitting the first few values of the computed shifts

data = {k[[#]],delta[[#]]}&/@ Range[25];

coeff=slope/.FindFit[data,slope x^(2L+1),slope,x];

Plot[coeff x^(2L+1),{x,0,(Last[data])[[1]]},Epilog-> {PointSize[0.02],Point/@ data}]

The scattering length is

αScattering = - coeff;

Born approximation

The high k behavior must match Born approximation:

(Subscript[δ, L])= ∫dr r V[r] (((Subscript[J, L+(1/2)])[k r])^2).

kBorn =Take[k,{Length[k],Length[k]-100,-5}];

These are phases in the Born approximation

deltaBorn = Table[-Pi NIntegrate[V[r] r*BesselJ[L+1/2,kBorn[[i]]*r]^2,{r,0,Infinity},PrecisionGoal-> 5],
{i,1,Length[kBorn]}];

The comparison with numerical results show the agreement

ListPlot[Transpose[{k,delta}],Joined-> True,PlotRange-> All,Epilog-> {PointSize[0.01],Point/@ 
Transpose[{kBorn,deltaBorn}]}]

Resonances

We know form the text and from the notebooks of chapter 11 and 13 that dδ/dE has a peak for E∼(Subscript[E, R]), for a metastable state of energy (Subscript[E, R])

δ'[E] = ((Γ/2)/((Γ^2)/4 + (((E-E0))^2))).

Here we check again this result for this model. We first define an interpolating function for the data

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

and make an overall plot to locate possible resonant states

ymax =Ceiling[Max[ (2/π)fdelta'[ene]]]/2;
ListPlot[Transpose[{ene,(1/π)fdelta'[ene]}],Joined-> True,
AxesOrigin-> {0,0},PlotRange-> {{-0.1,Max[ene]/4},{-0.5,ymax}},AxesLabel-> {"E","δg(E)"}]

A resonance is there at E ∼ 3.4 . Its exact location and parameters can be obtained with a fit

{xmin,xmax} = {3.3,3.6};
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]

Show[Plot[b/2/(Γ^2/4+ (x-E0)^2)/.fit,{x,xmin,xmax},PlotRange-> All],
 ListPlot[Transpose[{ptEne,ptDelta}],PlotStyle-> PointSize[0.012]]]

The result is consistent with that found in previous chapters for this model.

We leave to the reader the study of p and d waves.

Scattering solution

Clear[k,δ,maxL,sol];

To compute the complete scattering solution we have to solve the radial equations (up to a cutoff angular momentum) and sum them with formula (1.).

The function

{δ, normCoeff, sol} = getRadialSolution[k, L]

complete this task for a given momentum k and a fixed L. The solution is computed with NDSolve acting on the expression eq[k,L] defined above. A distant point xR is selected as right boundary. At points xR-0.5, xR-0.7 (two others points also work) we compute the value of the solution and by solving the linear system (1.) we compute the phase shifts δ and the normalization of the solution, stored in the variable normCoeff.

getRadialSolution[k_,λ_]:= Module[{x1,x2,c11,c12,c21,c22,d1,d2,δ,val1,val2,sol,normCoeff},
(* xR = Max[20.,1.5(λ+1/2)k]; *)
xL = 0.1;
xR = 40.;
sol = NDSolve[eq[k,λ],u,{x,xL,xR}][[1]];
{x1,x2}= {xR-05,xR- 0.7};
{val1,val2} = { (u[x1]/.sol), (u[x2]/.sol)};
{c11 ,c21}= x1 {jL[k x1][[λ+1]],yL[k x1][[λ+1]]};
{c12,c22} = x2 {jL[k x2][[λ+1]],yL[k x2][[λ+1]]};
{d1,d2}={a1,a2}/.First[NSolve[{a1*c11 - a2 *c21 == val1,a1*c12 - a2 *c22 == val2},{a1,a2}]];  
δ = ArcTan[d1,d2];
normCoeff = 1./Sqrt[(d1^2 + d2^2)];
{δ,normCoeff,sol}];

The function

{𝒜0, 𝒜, normPsi, sol} = getSolution[k, maxL]

assemble together the single solutions obtained with getRadialSolution to form a list of solutions, and computing the vectors of coefficients

𝒜0 = (i^L)(2L+1);    𝒜 = (i^L)(2L+1)(e^(2i (Subscript[δ, L])));

needed to write the ψ function. The function psi0List just make a list of the function ψ0 for all needed values of L.

getSolution[k_,maxL_]:= Module[{delta,normCoeff,sol,mom,𝒜0,𝒜},
{δ,normCoeff,sol}= Transpose[getRadialSolution[k,#]&/@Range[0,maxL]];
mom = Range[0,maxL];
𝒜0 = (I^mom) (2 mom + 1);
𝒜 = (I^mom) (2 mom + 1)*Exp[I δ];
{𝒜0,𝒜,normCoeff,sol}
]

Warning: the phases δ computed by this function are global, the reader must not use this variable.

psi0List[k_,r_,maxL_]:=Table[psi0[λ,k,r],{λ,0,maxL}]

Finally we define

ψ[k, r, θ]

For x<xL the solution ψ0 is taken while for x > xL we write the function in the form

ψ = Exp[i k x] + (∑_L)(i^L)(2L+1)(e^(2i (Subscript[δ, L]))) (1/r) (Subscript[u, L])[r]

The coefficient normPsi computed with getSolution takes care of the normalization. This coefficient is global and getSolution has to be called before any use of ψ.

ψ[k_,r_,θ_]:= Module[{},
Exp[I k r Cos[θ]] +Tr[
(𝒜*normPsi/r*If[r ≥ xL,( u[r]/.sol) ,psi0List[k,r,Length[normPsi]-1]]
- 𝒜0* jL[N[k r]])*legP[N[θ]]]];

Example

k1=1.5; maxL=30;

makeList[maxL];
{𝒜0,𝒜,normPsi,sol} = getSolution[k1,maxL];

GraphicsRow[{Plot[Abs[ψ[k1,r,0.]],{r,0,2xL}],Plot[Abs[ψ[k1,r,0.2]],{r,0,10}]},ImageSize-> Large]

Angular momentum cutoff

We have to verify how the angular momentum cutoff works. We plot the contribution of scattered solution, obtained by subtracting from the exact solution the free L-wave radial solution, for each L at several random radii and angles

maxL = 30;
makeList[maxL];
kv = {1.5,3.,4.5};

GraphicsRow[Table[{𝒜0,𝒜,normPsi,sol} = getSolution[kv[[j]],maxL];
Show[Table[r1 = RandomReal[{0.1,3}];θ1 = RandomReal[{0,Pi}];
ListPlot[Log[10,Abs[(𝒜*normPsi/r1*( u[r1]/.sol) - 𝒜0* jL[N[kv[[j]] r1]])*legP[N[θ1]]]]
,Joined-> True],{50}]],{j,1,Length[kv]}],ImageSize-> Large]

We see that up to k ∼ 5 a cutoff of 30 for angular momentum suppress higher L by 8-10 orders of magnitude.

To have an overall picture of ψ we can use a density plot, as in the function below.

makeContour[numPtGraph_,k_,{xmin_,xmax_,ymax_}]:=
Module[{maxL},
maxL = 30;
makeList[maxL];
{𝒜0,𝒜,normPsi,sol} = getSolution[k,maxL];
 DensityPlot[((Abs[ψ[k,(Sqrt[x^2+y^2]),ArcTan[x,y]]])^2),{x,xmin,xmax},{y,-ymax,ymax}, 
  PlotPoints-> numPtGraph,MaxRecursion-> 1, 
PerformanceGoal-> "Quality",ColorFunction->GrayLevel,
AspectRatio-> Automatic,
PlotLabel-> "k = "<>ToString[NumberForm[k,{3,2}]]]];

Let us use this function to show how fast and dramatic is the change in the density when E cross a resonant value

kRis = (Sqrt[2*3.42647])

GraphicsGrid[{makeContour[50,#1,{-10,10,10}]&/@{2.59,kRis,2.7}},ImageSize-> Large]

We see that at resonance the ψ is "trapped" in the potential.

Cross sections

The function getFandSigma[k, θ] compute the total cross section and the differential cross section for a list of k and θ.

getFandSigma[kvIn_,θvIn_]:= Module[{f2,σ,kv,θv},
kv = Flatten[{N[kvIn]}];
θv = Flatten[{θvIn}];
f2 = Table[0,{Length[kv]},{Length[θv]}];
σ = Table[0,{Length[kv]}];
Do[
k = kv[[j]];
makeList[maxL];
{𝒜0,𝒜,normPsi,sol} = getSolution[k,maxL];
(* we remember that 𝒜 - 𝒜0 = (Exp[2 I δ]-1) (I^L) (2L+1)  *)
f2[[j]] =Abs[ Tr[(𝒜 - 𝒜0)/(2 I k)/(I^Range[0,maxL]) * legP[#]]]^2&/@θv;
σ[[j]] = 4π /k^2 Tr[(2 Range[0,maxL]+1)* Sin[δ]^2];
,
{j,1,Length[kv]}];
{σ,f2}]

Here we do a short calculation for σ and angular distribution. Some warning messages are produced by Mathematica but this do not affect the results

maxL = 50;
kv =Union[Range[0.05,0.4,0.05], Range[0.4,3.,0.2],Range[2.58,2.64,0.01],Range[3.,8.,0.5]];
{σ,f2}=getFandSigma[kv,0.1];//Timing

This is the behavior of σ as a function of k. The low energy limit 4π (α^2) is shown

ListPlot[Transpose[{kv,σ}],Joined-> True,PlotRange-> {{0,Max[kv]},{0,1000}},
AxesLabel-> {"k","σ"},
Epilog-> {
Dashing[{0.015,0.015}],Line[{{0,4π αScattering^2},{Max[kv],4π αScattering^2}}]}]

Born approximation

For the low momenta considered the Born approximation is surely inaccurate, but let us compute it anyway.

Clear[V0,q]

With q=2kSin[θ/2]the Born amplitude is

fB = (- 2/q Integrate[ V[r] Sin[q r] r ,{r,0,Infinity},Assumptions -> q > 0])/.q-> 2 p Sin[θ/2]

The total cross section is computed by integrating on the solid angle

σBp=2π Integrate[fB^2 Sin[θ],{θ,0,π},Assumptions -> p>0]

Its numerical values are in complete disagreement with the numerical data

V0 = 7.5;
σB = σBp/.{p-> kv}

Eikonal approximation

Clear[V0]

In the text it has been shown that in the eikonal approximation

f[0] = - i k (∫_0^∞)(Exp[2 i Δ[b]]-1) bdb ;     Δ[b] = ((-m)/(2 k (ℏ^2)))2(∫_0^∞)V[(Sqrt[(b^2)+(t^2)])]dt  .

By the optical theorem :

σ = ((4π)/k)Im[f[0]] = 8 π (∫_0^∞)Sin[Δ[b]]^2 bdb.

To compute Δ first let us note that

V[r] = (∂_(α,α)) V0 Exp[-α r];    α→ 1

With (r^2)=(b^2)+(t^2) by changing integration variable to r we have

(∫_0^∞)Exp[-α (Sqrt[(b^2)+(t^2)])]dt = (∫_b^∞)Exp[-α r](r/(Sqrt[(r^2)- (b^2)]))dr = (Subscript[(Exp[-α r](Sqrt[(r^2)- (b^2)])), b,∞])+ α (∫_b^∞)Exp[-α r](Sqrt[(r^2)- (b^2)]).

The first term vanishes while the integral is

α Integrate[Exp[-α r] (Sqrt[r^2-b^2]),{r,b,Infinity},Assumptions->{ b>0,α > 0}]

Then Δ can be written (in units ℏ = m = 1):

Δb =  (1/p)Simplify[(- V0 (∂_(α,α)) (b BesselK[1,b α]))/.α-> 1]

To compute the final integral over b to get σ we can use the fact that Bessel k functions go to zero exponentially at infinity, then a cutoff integral is sufficient

V0=7.5;
σEIK = 8 π NIntegrate[((Sin[(((b^3) V0 (3 BesselK[1,b]+BesselK[3,b]))/(4 #1))])^2) b ,{b,0,50}]&/@kv;

Here we show the numerical data together with the eikonal approximation :

ListPlot[{Transpose[{kv,σ}],Transpose[{kv,σEIK}]},
PlotStyle-> {Black,{Red,Dashed}},
AxesLabel-> {"k","σ"},
Joined-> True,PlotRange-> {{0,Max[kv]},{0,1000}}]

The agreement is quite good even if the k values are not particularly high.

References