Scattering 2

Quantum Mechanics
                     K. Konishi, G. Paffuti

Solve numerically the scattering problem for an infinite spherical potential well of radius a.

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)/(2k))(i^L) (e^(i (Subscript[δ, L]))) ;

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

ψ ∼ 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 asymptotic scattered L-wave solution is easily written as

(i^L)(((2L+1))/2)((e^(2i (Subscript[δ, L])))-1)(Subscript[h, L]^((1)))[k r]

The problem

The potential in our problem is

U[r] = {({{∞, r<a}, {0, r>0}}) .

In the region r>a the Schrödinger equation describes a free particle

Δ ψ + (k^2)ψ = 0.

The asymptotic form of the solution

R = 2k (Cos[δ] sJ[L,k r]- Sin[δ]sY[L,k r])

coincides with the solution itself as the potential vanishes. The notation has been defined above.

Scattering phases

The boundary condition is   ψ = 0  for r=a. This implies from ():

Cos[δ] sJ[L,k a]- Sin[δ]sY[L,k a]= 0   ⇒  Tan[δ[L]] = ((sJ[L,k a])/(sY[L,k a])).

This formally solves our problem.

Numerical computation

To write a finite numerical expression we put a limit of 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 turns out to be more efficient to use recursion relation to produce the whole list. The function

makeList[maxL]

defines 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]

computes the vectors of phase scattering and of factors

(Subscript[S, L])= (i^L)(((2L+1))/2)*(Exp[2i δ]-1)

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.

Code for lists

makeListJ1[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]];

makeListY1[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];

makeListJ[maxL_,x_]:= If[ x== 0.,Table[If[m==0,1.,0.],{m,0,maxL}],
Module[{JList=Table[1.,{maxL+1}]},
JList[[maxL+1]]= (Sqrt[(π/(2x))])BesselJ[maxL+(1/2),x]; 
JList[[maxL]]=  (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]] = (Sqrt[(π/(2x))])BesselY[(1/2),x];
YList[[2]] = (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,θ];)

makeδAS[maxL_,q_,a_]:= Module[{δ,𝒜,𝒮,mom,k},
k = N[q]; 
δ = Re[ArcTan[ jL[k a]/ yL[k a]]];
 mom = Table[I^j(2 j +1),{j,0,maxL}];
(* these factors are essentially the S matrix *)
𝒮 = (1/2)mom*(Exp[2I δ]-1); 
{δ,𝒮}];

Checks

We perform some checks on our lists to avoid misprint errors in the code

maxLCheck = 20; aCheck=1.;

We define our lists

makeList[maxLCheck];

Here we compare the results between a straightforward computation and a computation with the lists (for Y and H1 the relative error is printed as the functions diverge at the origine.

errors=Transpose[Table[rCheck = Random[Real,{0.1,15.}]; θCheck =Random[Real,{0.01π,0.9π}];
{Max[Chop[Abs[jL[rCheck]-SphericalBesselJ[Range[0,maxLCheck],rCheck]]]],
Max[Chop[Abs[yL[rCheck]-SphericalBesselY[Range[0,maxLCheck],rCheck]]/Abs[yL[rCheck]]]],
Max[Chop[Abs[hL[rCheck]-SphericalHankelH1[Range[0,maxLCheck],rCheck]]/Abs[hL[rCheck]]]],
Max[Chop[Abs[legP[θCheck]-LegendreP[Range[0,maxLCheck],Cos[θCheck]]]]]},{i,1,10}]];
maxErrors = Max[#1]&/@errors

The scattering solution

The definition of the solution

Clear["jL*","yL*","hL*","legP*","𝒮","k*"]

We choose units such that a = 1.

a = 1.;

This is the definition of the scattering solution, as a sum of the plane wave and the asymptotic form, which in this problem is exact. We use Clear on the relevant variables in such a way that this definition will work below, not being affected by the previous use of the lists.

ψ[r_,θ_,k_,𝒮_]:=If[Re[r] ≤  Re[a],0, Exp[I k r Cos[θ]] +Tr[𝒮*( jL[N[k r]] +I yL[N[k r]])*legP[N[θ]]]];

The choice of  angular momentum cutoff

We choose a cutoff maxL = 2 k a and show that it is reasonable.

k = 50.;

Let us define our functions and compute phase shifts

maxL = Floor[2 k a];
makeList[maxL];
{δ, 𝒮} = makeδAS[maxL,k,a];

We plot at random values of r the contribution of each  momentum (in logarithmic scale). The horizontal axis is the index 1,2,3 .. enumerating angular momentum 0,1,2 ..

Show[Table[ListPlot[Log[10,Abs[𝒮*legP[Random[Real,{0,Pi}]]*hL[ k*Random[Real,{1.1,5}]]]],Joined-> True],{50}]]

We see that momenta of order 2 k a are depressed by at least 12 orders of magnitude, while angular momenta up to (k a) (half of the plot) have the same order of magnitude.

Some check for times

This is the time used in the previous computation

Show[Table[ListPlot[Log[10,Abs[𝒮*legP[Random[Real,{0,Pi}]]*hL[ k*Random[Real,{1.1,5}]]]],Joined-> True],{50}]]//Timing

Here we do the same graph (only random numbers change) without lists but using directly Mathematica functions

Show[Table[ListPlot[Log[10,Abs[
𝒮*LegendreP[Range[0,maxL],Cos[Random[Real,{0,Pi}]]]*SphericalHankelH1[Range[0,maxL], k*Random[Real,
{1.1,5}]]]],Joined-> True],{50}]]//Timing

It is evident that using lists is more convenient.

Density plot for the solution

This function make a density plot for ψ

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

An example of the total intensity distribution is given here

k=10.;
Show[makeContour[100,k,{-10,10,4}],ImageSize-> Medium]//Timing

The brighter color corresponds to higher values.  Note a forward luminosity, due to diffraction. The diffraction peak shrinks as k grows:

GraphicsRow[
Table[Show[makeContour[100,k,{-10,10,4}],ImageSize-> 4 72],{k,{5,10}}],Scaled[0.01]]

The shadow

As k increases a shadow region extends downstream. In the following plot we enhance the interesting zone for a higher value of k

k=40.;
Show[makeContour[100,k,{0,10,4}],ImageSize-> Medium]//Timing

The reader can use higher values for k, but the quality of the plots deteriorates: one has to increase the number of plot points or the parameter MaxRecursion in the body of the function.

The shadow can also be detected by plotting the intensity on a plane at fixed distance form the sphere

k=40;
x0=5;
maxL=Floor[2 k a];
makeList[maxL];
{δ,𝒮} = makeδAS[maxL,k,a];

Plot[((Abs[ψ[(Sqrt[(x0^2)+(y^2)]),ArcTan[x0,y],k,𝒮]])^2),{y,-3,3},
AxesLabel-> {"y","ψ"},PlotRange-> All]

or by plotting the intensity at a small transverse coordinate

kvalues = {10,30,50,100};
pl=Table[0,{Length[kvalues]}];

y0=0.1;

Table[
k = kvalues[[i]];
maxL=Floor[2 k a];
makeList[maxL];
{δ,𝒮} = makeδAS[maxL,k,a];
pl[[i]]=Plot[Tooltip[((Abs[ψ[(Sqrt[(x^2)+(y0^2)]),ArcTan[x,y0],k,𝒮]])^2),k],{x,1,10},
AxesLabel-> {"x","ψ"},PlotRange-> All];,
{i,1,Length[kvalues]}];

Show[pl]

A  simpler presence of the forward part in the scattering is obtained by plotting only the scattered wave intensity:

makeScatteringContour[numPtGraph_,k_,{xmin_,xmax_,ymax_}]:=
Module[{maxL},
maxL=Floor[2k a];
makeList[maxL];
{δ,𝒮} = makeδAS[maxL,k,a];
 DensityPlot[((Abs[ψ[(Sqrt[x^2+y^2]),ArcTan[x,y],k,𝒮]- Exp[I k x]])^2),{x,xmin,xmax},{y,-ymax,ymax},
ColorFunction->GrayLevel, PerformanceGoal-> "Quality",
  PlotPoints-> numPtGraph,MaxRecursion-> 1, 
AspectRatio-> Automatic,
PlotLabel-> "k = "<>ToString[NumberForm[k,{3,2}]],
Epilog->{{GrayLevel[0],Disk[{0,0},a]}}]];

GraphicsRow[
Table[Show[makeScatteringContour[100,k,{-5,20,4}],ImageSize-> 4 72],{k,{5,10}}],Scaled[0.01]]

Let us note that the high value of ψ - Exp[i k x] just after the sphere is the "shadowing": this part of the waves cancels the incoming wave and produce the shadow.

Low energy approximation

Let us verify the low k behavior of phase shift. We have found

Tan[(Subscript[δ, L])] = (((Subscript[j, L])[k a])/((Subscript[η, L])[k a]))

A simple Taylor expansion gives (we write b for the radius)

Block[{p,b},TableForm[Table[
Normal[Series[ArcTan[ SphericalBesselJ[j,p b]/ SphericalBesselY[j,p b]],{p,0,7}]],{j,0,3}],
TableHeadings-> {Range[0,3],{}}]]

In agreement with the expected behavior  (Subscript[δ, L])∼(k^(2L+1)). The scattering length is equal to the radius of the well.

Kirchhoff approximation

This scattering problem is identical to the diffraction problem for a scalar wave in electromagnetism or acoustics. Application of Kirchhoff principle give an approximate asymptotic solution in the form (see Problem 4 of this chapter)

ψScatteredKirchhoff[k_,r_,θ_]:= -I (a/2)(-I (E^(-2I k a Sin[θ/2]))- k a(((1+Cos[θ])/(k a Sin[θ])))BesselJ[1,k a Sin[θ]])((E^(I k r))/r)

The first term is a reflected wave of constant amplitude while the second term describes the diffractive forward peak. Both terms have a limited validity in angle, the first is not accurate around θ = 0, while the second one is valid in the forward region. Their interference is everywhere negligible.

The whole approximate solution take the form

ψKirchhoff[k_,r_,θ_]:= Exp[I k r Cos[θ]] + ψScatteredKirchhoff[k,r,θ];

Let us check our solution by comparison with this one.

k=40;
maxL=Floor[2 k a];
makeList[maxL];
{δ,𝒮} = makeδAS[maxL,k,a];

We plot the value at a fixed distance x0 on the positive and negative axis along the propagation of the plane wave, the well being at the origin.

x0=30;

GraphicsRow[{Plot[{Abs[ψ[(Sqrt[(x0^2)+((Tan[θ])^2)(x0^2)]),θ,k,𝒮]-Exp[I k x0]],Abs[ψScatteredKirchhoff[k,(Sqrt[(x0^2)+((Tan[θ])^2)(x0^2)]),θ]]},{θ,0.0001π,0.1π},
AxesLabel-> {"θ","ψScatt"},
PlotStyle-> {Black,Red},PlotRange-> All],Plot[{Abs[ψ[(Sqrt[(x0^2)+((Tan[θ])^2)(x0^2)]),θ,k,𝒮]-Exp[-I k x0]],Abs[ψScatteredKirchhoff[k,(Sqrt[(x0^2)+((Tan[θ])^2)(x0^2)]),θ]]},{θ,0.5001 π, π},
AxesLabel-> {"θ","ψScatt"},
PlotStyle-> {Black,Red},PlotRange-> All]},ImageSize-> {6 72,2 72}]

The accuracy is quite good, (Kirchhoff result is only approximate) and gives us confidence on the consistency of the computation.

We can also compute directly the scattering amplitude. First we note that from

(Subscript[J, 1])[z]∼ z/2

fKirchhoff[0,k_]:= ((- a)/2)+I ((k (a^2))/2);

In general

Clear[r0];

f[θ_,k_]:= (ψ[r0,θ,k,𝒮]-Exp[I k r0 Cos[θ]])Exp[-I k r0]r0;
fKirchhoff[θ_,k_]:= -I (a/2)(-I (E^(-2I k a Sin[θ/2]))- k a(((1+Cos[θ])/(k a Sin[θ])))BesselJ[1,k a Sin[θ]]);

As the asymptotic form can have (k a)/(r^2) we choose a large r0, with respect the k's studied:

r0=1000;
k = 100;
maxL=Floor[2 k a];
makeList[maxL];
{δ,𝒮} = makeδAS[maxL,k,a];

Here is the comparison between numerical and Kirchhoff approximation for real and imaginary part of the amplitude, respectively in backward and forward directions

Plot[{Re[f[θ,k]],Re[fKirchhoff[θ,k]]},{θ,0.7π,0.99π},
PlotStyle-> {Black,{Dashed,Red}}]

Plot[{Im[f[θ,k]],Im[fKirchhoff[θ,k]]},{θ,(10^(-8))π,0.1π},
PlotStyle-> {Black,{Dashed,Red}},PlotRange-> All]

Let us in particular verify the θ = 0 imaginary amplitude, directly related to the cross section

kValues = Range[40,400,10];

f1= Table[
k=k1;
maxL=Floor[2 k a];
makeList[maxL];
{δ,𝒮} = makeδAS[maxL,k,a];
Im[f[0,k]],{k1,kValues}];
fK= Im[fKirchhoff[10.^-9π,#]]&/@kValues;

ListPlot[Transpose[{kValues,((f1-fK)/f1)}],PlotRange-> {All,{0,Max[((f1-fK)/f1)]}},AxesLabel-> {"k","ΔIm[f]/Im[f]"}]

We see that the agreement is excellent in all cases.

Cross sections

The function getFandSigma computes the angular distribution |f[θ]((|)^2) and the total cross section for a predefined set of values for k and θ. The component f2[[i]] contains the angular distribution for the i-th k value.

getFandSigma[kvIn_,θvIn_]:= Module[{f2,σ,maxL,δ,kv,θv,k},
kv = Flatten[{N[kvIn]}];
θv = Flatten[{θvIn}];
f2 = Table[0,{Length[kv]},{Length[θv]}];
σ = Table[0,{Length[kv]}];
Do[
k = kv[[j]];
maxL = Max[10,Floor[3 k]];
makeList[maxL];
{δ,𝒮} = makeδAS[maxL,k,a];
(* 𝒮 was defined as 𝒮 = (1/2)(I^λ) (2 λ + 1) *(Exp[2I δ]-1) *)
f2[[j]] =Abs[ Tr[𝒮/(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}]

Let us have a look at the angular distribution for different k

θvalues = Range[0, Pi,0.005 Pi];
kvalues = {0.1, 6., 20.0};
{σ,f2}=getFandSigma[kvalues,θvalues];

Table[ListPlot[Transpose[{θvalues,f2[[j]]}],Joined-> True,
PlotLabel->"k = "<>ToString[N[ kvalues[[j]],3]],PlotRange-> All],{j,1,Length[f2]}]

Here instead we study the behavior of σ as a function of k.

At low energy only the S-wave gives a contribution and the cross section tends to  4 π (a^2) as it easy to show (see the above computation of the scattering length).

At high energy the cross section tends to 2 π (a^2), twice the geometrical section of the sphere. A proof of this will be given in the next notebook, see also below.

kvalues =Union[Range[0.1,1.,0.1],Range[1.,20.,0.5],Range[20.,50.,1.]];

To compute a list of σ with getFandSigma it is sufficient to fix a single value of the angle.

{σ,f2}=getFandSigma[kvalues,0.1];

ListPlot[Transpose[{kvalues,σ/π}],Joined-> True,PlotRange-> {All,{1,4}},
AxesLabel-> {"k","σ/π"},Epilog-> {
Dashing[{0.015,0.015}],Line[{{0,2},{Max[kvalues],2}}]}]

The high energy behavior

Ther are various ways to compute the short wavelength limit for σ. Here we show how it works with the Kirkoff result. We have quoted

f[θ] ≃ -i (a/2)(i (e^(-2i k a Sin[θ/2]))- k a(((1+Cos[θ])/(k a Sin[θ])))BesselJ[1,k a Sin[θ]]);

but we stressed the fact that the terms are valid in complementary angular regions.  In fact the cross section must be computed independently for the two and summed.

We can compute the total cross section by optical theorem, the forward amplitude being well described by the diffractive term. This term is purely imaginary and using

(Underscript[Lim, θ→0])(((1+Cos[θ])/(k a Sin[θ])))BesselJ[1,k a Sin[θ]] = 1 ;

we obtain

σ= ((4π)/k)Im[(Subscript[f, diffr])[0]] = ((4π)/k) ((k (a^2))/2)= 2π (a^2).

Let us note that the "reflection" term in the amplitude is almost always valid, except in the forward direction. It is easy then to compute the cross section due to reflection:

(Subscript[σ, reflection]) = ∫dΩ |(Subscript[f, refl])[θ]((|)^2) = 4π ((((a/2)))^2)= π (a^2)

In fact,  the reflected part of σ is equal to the geometrical cross section  as expected.

The integral of the diffractive term can in effect be evaluated with Mathematica. Putting z = k a

2π Integrate[((((((1+Cos[t])/(z  Sin[t])))BesselJ[1,z Sin[t]]))^2)Sin[t],{t,0,π},Assumptions-> z>0]

The last term decrease as 1/(z^3):

Series[HypergeometricPFQ[{(3/2)},{(5/2),3},-(z^2)],{z,Infinity,1}]

The Bessel function goes to zero at large z then the integral approach, for k→∞

I = 4(π/(z^2))= ((4π)/((k^2)(a^2))).

Inserting the pre-factor  (k (a^2)/2) of the amplitude we get a consistent result:

(Subscript[σ, diffr]) = ((4π)/((k^2)(a^2))) (((((k (a^2))/2)))^2)= π (a^2) = σ - (Subscript[σ, reflect])

References