Scattering 3

Quantum Mechanics
                      K. Konishi, G. Paffuti

Solve numerically the scattering problem for a  spherical potential well of radius a.  Study the cross section and the possible rainbow.

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] = {({{(Subscript[U, 0]), r<a}, {0, r>0}}) .

Show[
Graphics[{Thickness[.02],Line[{{0,-3},{1,-3},{1,-0.04}}],Text["a",{1,0.4}],
Text["U0",{1.3,-2.2}]}],PlotRange-> {{0,4},{-4,2}},Axes->Automatic,AxesLabel-> {"r","U[r]"},Ticks-> {{2,3,4},Automatic},ImageSize-> 2 72]

We will consider to fix the notations the attractive case, (Subscript[U, 0])< 0 and put

(β^2)= - 2 ((m (Subscript[U, 0]))/(ℏ^2));      (q^2)= (k^2)+ (β^2).

The repulsive case is reached with the substitution β → i β.

The Schrödinger equation reads

r< a:   Δ ψ + (q^2)ψ = 0 ;         r>a:  Δ ψ + (k^2)ψ = 0 ;

The asymptotic form of the solution, which is identical to the true solution as U = 0 there, is (apart a multiplicative factor

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

For r < a one has again the same equation with k→q. The regular solution can contain only the non diverging term sJ[q r] then the solution has the form

(Subscript[R, int]) = (Subscript[C, L]) sJ[q r].

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

Form of the solution

If R[L, k, r] is the regular radial solution, with the overall normalization () the complete solution is

ψ = (∑_(L=0)^∞) (Subscript[A, L]) R[L,k,r] (Subscript[P, L])[Cos[θ]];        (Subscript[A, L])= ((i^L)/(2k))(2L+1)Exp[i (Subscript[δ, L])].

The scattered asymptotic wave is conveniently written as

(Subscript[ψ, sc]) = (∑_(L=0)^∞) (1/2)(i^L)(2L+1)((e^(2i (Subscript[δ, L])))-1)(Subscript[P, L])[Cos[θ]] (sJ[L,k r] + i sY[L, k r]).

The combination in the rhs is the spherical Hankel function of the first kind,  sH[L,k r] in our notations.

In our particular problem the asymptotic solution written in terms of spherical Bessel functions coincides with the exact solution, as U=0 for r > a. Then the whole solution can be alternatively written as

r> a:  ψ = (e^(i k z)) + (Subscript[ψ, sc]) = (e^(i k z))+ (∑_(L=0)^∞) (1/2)(i^L)(2L+1)((e^(2i (Subscript[δ, L])))-1)(Subscript[P, L])[Cos[θ]] (sJ[L,k r] + i sY[L, k r])
r<a: ψ = (∑_(L=0)^∞) (Subscript[A, L]) (Subscript[C, L]) sJ[q r]  (Subscript[P, L])[Cos[θ]];

The coefficients (Subscript[C, L])are fixed below.

Boundary conditions and  phase shifts

The boundary condition are

(Subscript[ψ, out])=(Subscript[ψ, in]);  (Subscript[ψ, out]^′) = (Subscript[ψ, in]^′);     at    r=a.

The continuity of ψ implies from () and ():

(Subscript[C, L])=2k(((Cos[δ] sJ[L,k a]- Sin[δ]sY[L,k a]))/(sJ[q a])) .

The phases (Subscript[δ, L]) can be most easily find with the continuity of ψ'/ψ

q((sJ'[λ, qa])/(sJ[λ,q a]))= k (( sJ'[λ,k a]-Tan[δ] sY'[λ,k a])/( sJ[λ,k a]-Tan[δ] sY[λ,k a]));

which gives

Tan[δ]=(( k  sJ[λ, q a]sJ'[λ,k a] - q  sJ'[λ, q a]sJ[λ,k a])/(-q sJ'[λ, q a] sY[λ,k a]+ k sJ[λ,q a] sY'[λ,k  a]))

This formally solves our problem.

Numerical computation

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,θ];)

makeδAS[maxL_,kIn_,{a_,β_}]:= Module[{δ,𝒜,𝒮,mom,djLqa,djLka,dyLka,q,k,valE},
k = N[kIn]; 
q = (Sqrt[k^2 + β^2]);
djLqa = ((- 1)/(2q a)) jL[q a][[Range[2,maxL]]]+(1/2)(jL[q a][[Range[1,maxL-1]]] - 
jL[q a][[Range[3,maxL+1]]]);
valE =( ((- 1)/(2q a))(Sqrt[(π/(2q a))])BesselJ[#+(1/2),q a] +
 (Sqrt[(π/(2q a))])(1/2) (BesselJ[-(1/2)+#,q a]-BesselJ[(3/2)+#,q a]))&/@{0,maxL};
djLqa = Append[Prepend[djLqa,valE[[1]]],valE[[2]]];


djLka = ((- 1)/(2k a)) jL[k a][[Range[2,maxL]]]+(1/2)(jL[k a][[Range[1,maxL-1]]] - jL[k a][[Range[3,maxL+1]]]);
valE =( ((- 1)/(2k a))(Sqrt[(π/(2k a))])BesselJ[#+(1/2),k a] +
 (Sqrt[(π/(2k a))])(1/2) (BesselJ[-(1/2)+#,k a]-BesselJ[(3/2)+#,k a]))&/@{0,maxL};
djLka = Append[Prepend[djLka,valE[[1]]],valE[[2]]];

dyLka = ((- 1)/(2k a)) yL[k a][[Range[2,maxL]]]+(1/2)(yL[k a][[Range[1,maxL-1]]] - yL[k a][[Range[3,maxL+1]]]);
valE =( ((- 1)/(2k a))(Sqrt[(π/(2k a))])BesselY[#+(1/2),k a] +
 (Sqrt[(π/(2k a))])(1/2) (BesselY[-(1/2)+#,k a]-BesselY[(3/2)+#,k a]))&/@{0,maxL};
dyLka = Append[Prepend[dyLka,valE[[1]]],valE[[2]]];

δ = Re[ArcTan[(k jL[q a]*djLka - q djLqa*jL[k a])/
(- q  djLqa  yL[k a] + k jL[q a] dyLka)]];
 mom = Table[I^j(2 j +1),{j,0,maxL}];
(* 𝒜 are connection coefficients *)
 𝒜 = mom*Exp[I δ]*(Cos[δ] jL[k a] - Sin[δ] yL[k a])/jL[q a];
(* these factors are essentially the S matrix *)
𝒮 = (1/2)mom*(Exp[2I δ]-1); 
{δ,𝒜,𝒮}];

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

Checks

Lists

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

maxLCheck = 20; aCheck=1.; rCheck = 1.7;θCheck = 0.65;

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

Derivatives

This is a check for the definition of the derivative

{dsJ[λ,x] == (∂_x)sJ[λ,x],dsY[λ,x] == (∂_x)sY[λ,x]}

Algebraic checks

We check that out expression for δ is correct

Simplify[q((dsJ[λ, q a])/(sJ[λ,q a]))== ( k (( dsJ[λ,k a]-Tan[δ] dsY[λ,k a])/( sJ[λ,k a]-Tan[δ] sY[λ,k a])))/.Tan[δ]-> (( k  sJ[λ, q a]dsJ[λ,k a] - q  dsJ[λ, q a]sJ[λ,k a])/(-q dsJ[λ, q a] sY[λ,k a]+ k sJ[λ,q a] dsY[λ,k  a])),
Assumptions-> {k>0,r>0,λ≥0,q>0,a>0,λ∈ Integers};]

Checks for the phases

We check the expression for the phases generated by makeδAS with that computed directly from its definition

maxL = 20; k =0.5; a=1.; β = 10.;

makeList[maxL];
{δ,𝒜,𝒮} = makeδAS[maxL,k,{a,β}];

This is the analytical definition given above, the agreement is excellent

δDirect= Module[{q },
q = (Sqrt[(k^2)+(β^2)]);
Table[ArcTan[(( k  sJ[λ, q a]dsJ[λ,k a] - q  dsJ[λ, q a]sJ[λ,k a])/(-q dsJ[λ, q a] sY[λ,k a]+ k sJ[λ,q a] dsY[λ,k  a]))],{λ,0,maxL}]];
Max[Abs[δ-δDirect]]

The same is done for a repulsive potential

maxL = 20; k =0.5; a=1.; β = I 10.;
makeList[maxL];
{δ,𝒜,𝒮} = makeδAS[maxL,k,{a,β}];
δDirect= Module[{q },
q = (Sqrt[(k^2)+(β^2)]);
Table[ArcTan[(( k  sJ[λ, q a]dsJ[λ,k a] - q  dsJ[λ, q a]sJ[λ,k a])/(-q dsJ[λ, q a] sY[λ,k a]+ k sJ[λ,q a] dsY[λ,k  a]))],{λ,0,maxL}]];
Max[Abs[δ-δDirect]]

Clear[a,β,q,k,maxL];

The scattering solution

Low energy and scattering lengths

The analytical expression for δ phases is

δfunct[λ_,k_,β_,a_]:=Module[{q },
q = (Sqrt[(k^2)+(β^2)]);
ArcTan[(( k  sJ[λ, q a]dsJ[λ,k a] - q  dsJ[λ, q a]sJ[λ,k a])/(-q dsJ[λ, q a] sY[λ,k a]+ k sJ[λ,q a] dsY[λ,k  a]))]]

Attractive potentials

The low energy behavior for (Subscript[δ, 0])is

Collect[Simplify[Normal[Series[δfunct[0,k,β,a],{k,0,3}]],Assumptions->{ β>0,a>0}],k,Simplify]

In agreement with the general law (Subscript[δ, L])∝ (k^(2L+1)). The scattering length is

αLength = a-((Tan[a β])/β)

This formula breaks down at  a β=(2 n + 1) π/2. As it is known from our study of radial well at these values of the depth a zero energy bound state appear. As Tan[x]/x > 1 for x < π/2 we see that the scattering length is negative for small depths, then becomes positive when a bound state appears as β increases.

Plot[1- Tan[x π/2]/(π x/2),{x,0,4.1},AxesLabel-> {"2βa / π","α/a"}]

We can easily write the radial solutions, for  r > a

Rext[λ_,k_,r_,β_,a_]:= 2 k( Cos[δfunct[λ,k,β,a]] sJ[λ,k r] - Sin[δfunct[λ,k,β,a]] sY[λ,k r]);

and for all r > 0

R[λ_,k_,r_,β_,a_]:=
Module[{q},q = (Sqrt[(k^2)+(β^2)]);
 If[ r ≥ a, Rext[λ,k,r,β,a], Rext[λ,k,a,β,a]/sJ[λ,q a] sJ[λ,q r] ]];

The following graph shows the geometrical meaning of α, the scattering length, and show the shape of the wave function when bound states are absent or present.

{β1,β2}={1.1,2.1};
Module[{valF, valdF,β,k,a,f,fext,sign,pl1,pl2},
a=1;
 k = 0.5;
β := β1/a;
sign:=Sign[R[0,k,0.1 a,β,a]];
f[r_]:=  R[0,k,r,β,a]sign;
fext[r_]:= Rext[0,k,r,β,a]sign;
valF :=a* fext[a];
valdF := fext[a]+a*fext'[a];
pl1=Plot[{If[r>0,r*f[r],0],If[r<a,valF + valdF(r-a),r fext[0,k,r,β]]},{r,-1.5,2},PlotStyle-> {Automatic,{Red,Dashed}},
AxesLabel-> {"r", "χ[r]"},
 Epilog->{Text["α",{(-valF + a valdF)/valdF,0.12}],
Text["V < 0",{1.5,1}],Text["without bound states",{1.3,0.7}]}];
β := β2/a;
pl2=Plot[{If[r>0,r*f[r],0],If[r>a,valF + valdF(r-a),r fext[0,k,r,β]]},{r,0,2},PlotStyle-> {Automatic,{Red,Dashed}},
AxesLabel-> {"r", "χ[r]"},
 Epilog->{Text["α",{(-valF + a valdF)/valdF,0.12}],
Text["V < 0",{0.6,0.5}],Text["with bound state",{0.7,0.3}]}];
GraphicsRow[{pl1,pl2},ImageSize-> Large]]

Repulsive potential

For repulsive potential it is sufficient to make the substitution β → i β.

β1 = 3 I ;
Module[{valF, valdF,β,k,f,a,fext,sign,pl1,pl2},
a=1;
 k = 0.5;
β := β1/a;
sign:=Sign[R[0,k,0.1 a,β,a]];
f[r_]:=  R[0,k,r,β,a]sign;
fext[r_]:= Rext[0,k,r,β,a]sign;
valF :=a* fext[a];
valdF := fext[a]+a*fext'[a];
pl1=Plot[{If[r>0,r*f[r],0],If[r<a,valF + valdF(r-a),r fext[0,k,r,β]]},{r,0,2},PlotStyle-> {Automatic,{Red,Dashed}},
AxesLabel-> {"r", "χ[r]"},
 Epilog->{Text["α",{(-valF + a valdF)/valdF,-0.12}],
Text["V > 0",{1.0,1}]}]]

Zero energy bound states

As β increases the well becomes deeper and more bound states appear. The values of β for which the L-wave states becomes bounded are easily found: they correspond to zero energy normalizable solutions. These solutions can be approached from the "scattering side", i.e. from E>0. At zero energy the normalizable solution for r > a behaves (r^(-(L+1)))(it is a free particle equation), for S-waves the result is to be understood as a limit L→0. Its logarithmic derivative at r=a=1 is -(L+1).

This is confirmed by a direct computation: for bound states with E = -(p^2)/2 the normalizable solution is (Subscript[h, L]^((1)))[i p r]and its logarithmic derivative at r=a=1 is, as p→0:

Limit[(Table[(∂_r)SphericalHankelH1[L,p r]/SphericalHankelH1[L,p r],{L,0,2}])/.r-> 1,p-> 0]

The connection formula with internal solution becomes

q((sJ'[λ, q a])/(sJ[λ,q a]))= - (L+1) ;   q = (Sqrt[(k^2)+(β^2)]) → β

This is an equation for β and its successive solutions gives the depth at which L-wave bound states appear. The first few equations are

zeroEnergyEquations=Simplify[FunctionExpand[(β dsJ[L,β]+ (L+1)sJ[L,β])/.L-> Range[0,3]]]

The relative depths are

zeroEnergy=β/.Table[FindRoot[zeroEnergyEquations[[i]] == 0 ,{β,1.0+i}],{i,1,4}]

Complete solution

The definition of the solution

NOTE : From now on the variables k, a, β, 𝒮, 𝒜 will be global variables. We clear all possible stored value.

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

We choose units such that a = 1.

a=1.;

This define ψ. For r > a we choose to write the solution as the incident wave plus the scattered wave.

ψ[r_,θ_]:= Module[{q =(Sqrt[k^2 + β^2])}, If[ Re[r] < Re[a],Tr[𝒜* jL[N[q r]]*legP[N[θ]]],
Exp[I k r Cos[θ]] +Tr[𝒮*( jL[N[k r]] +I yL[N[k r]])*legP[N[θ]]]]];

Example

β = 10;a=1;k=0.5; maxL=20;

makeList[maxL];
{δ,𝒜,𝒮} = makeδAS[maxL,k,{a,β}];
Plot[Abs[ψ[r,0.1]],{r,0,1.2 a}]

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 angular momenta higher than k a are suppressed, we choose as a typical  cutoff  maxL= 2 k a.

Density plot for the solution

This function make a density plot for ψ

makeContour[numPtGraph_,kIn_,βIn_,{xmin_,xmax_,ymax_}]:=
Module[{maxL},
{k,β}=N[{kIn,βIn}]; (* we fix the global variables *)
maxL = Max[10,Floor[2 k a]];
makeList[maxL];
{δ,𝒜,𝒮} = makeδAS[maxL,k,{a,β}];
βPrint = If[Im[β]== 0,ToString[NumberForm[β,{3,2}]],
ToString[NumberForm[Im[β],{3,2}]]<>" I"];
 DensityPlot[((Abs[ψ[(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,{2,1}]]<> 
" β= "<>βPrint]];

a=1.;k = 10.0;
β = 20. I;

makeContour[100,k,β,{-5,10,4}]

Low energy resonant scattering.

At low energy one expect a quasi - resonant scattering if bound states exist with E ∼ 0. We have computed above the critical values of β for the appearance of bound states: {1.5708, 3.14159, 4.49341, 5.76346}, for L = {0,1,2,3}. These value are stored in the list zeroEnergy.

It is interesting to show how change the shape of the scattering solution as β cross one of this critical values:

k=0.1;

GraphicsGrid[{makeContour[70,k,#1,{-3,3,3}]&/@{0.5π,π,4.46,4.495}},ImageSize-> Large]

The appearance of the p - wave (second plot) and of the d-wave is quite evident.

Cross section

The amplitude, the differential cross section and the total cross section are:

f[θ] = (∑_(λ=0)^∞) ((((e^(2 i δ[λ]))-1) (1+2 λ))/(2 i k))LegendreP[Cos[θ]] ;  
d σ = |f[θ]((|)^2)  d Ω;
σ= ∫dΩ|f[θ]((|)^2)   = ((4 π)/(k^2))(∑_(λ=0)^∞)(1+2 λ) Sin[δ[λ]]^2

Low energy and high energy

Low energy

s-wave is dominant as  k-> 0. In this regime f = - α, where α is the scattering length computed above. Then

σ → 4π (α^2);        α =  a-((Tan[a β])/β).

High energy

For k >> 1/a has a forward peak and the eikonal approximation can be used (see text):

f = - i k (∫_0^∞)( Exp[2 i Δ[b]] - 1) BesselJ[0,k θ b]b db;
Δ[b]= - (m/(2(ℏ^2)k))(∫_(-∞)^∞)U[(Sqrt[b^2 + z^2])]dz;

The range of validity of this approximation is

Abs[U] ≪ ((ℏ^2)/(m (a^2)))(((k a))^2).     In this problem: |β((|)^2)≪ 2(k^2).

At even higher energies the scattering is described by the Born approximation (see text):

f → fB= - (( m)/(2 π (ℏ^2)))∫U Exp[- i q r]dV;       q = k - k';    q = 2 k Sin[(θ/2)];
Abs[U] ≪ ((ℏ^2)/(m (a^2)))(k a).     In this problem: |β((|)^2)≪ 2k/a.

For a central potential the angular integral can be carried out giving:

fB = - ((2 m)/( (ℏ^2)))(∫_0^∞)U[r] ((Sin[q r])/q) r d r.

Analytical results

Born approximation

In our notations - 2 m U/(ℏ^2) = (β^2)  and the Born approximation for a radial well takes the form

fB = (β^2)(∫_0^a) ((Sin[q r])/q) r d r = ((β^2)/(q^3))(∫_0^(a q)) Sin[x] x d x = ((β^2)/(q^3))(Sin[q a] - q a Cos[q a]).

To compute the total cross section one has to perform the angular integral. As fB do not depends on the azimuthal angle this last integration gives a factor 2π. The integration over the polar angle θ can be transformed to an integral over q using  q = 2 k Sin[θ/2]:

dΩ → ((2π)/(k^2)) q dq .

The total cross section is

σB = (∫_0^(2k))(fB^2)((2π q)/(k^2)) dq = (β^4)π (-((a^2)/(8 (k^4)))+((a^4)/(2 (k^2)))-(((Sin[2 a k])^2)/(32 (k^6)))+((a Sin[4 a k])/(16 (k^5)))) (Underscript[→, k → ∞]) (β^4)π((a^4)/(2 (𝓀^2))) .

Eiconal approximation

The total cross section can be computed by the optical theorem

σ = ((4π)/k)Im[f[0]].

As BesselJ[0,0] = 1, in the eikonal approximation one gets

f[0] = - i k (∫_0^∞)( Exp[2 i Δ[b]] - 1) b db;         Δ[b]= - (m/(2(ℏ^2)k))(∫_(-∞)^∞)U[(Sqrt[b^2 + z^2])]dz.

In our case the integral inn z is bounded as (b^2)+(z^2)≤ (a^2), then

Δ[b]= - (m/(2(ℏ^2)k))(∫_(-∞)^∞)U[(Sqrt[b^2 + z^2])]dz = - (m/((ℏ^2)k))(∫_0^∞)U[(Sqrt[b^2 + z^2])]dz 
          = - (m/((ℏ^2)k))U (Sqrt[(a^2)- (b^2)]) = ((β^2)/(2 k))(Sqrt[(a^2)- (b^2)]).

The impact parameter b is bounded by b ≤ a. With the change of variables

x = ((Sqrt[(a^2)- (b^2)])/a);  (b^2) =(a^2) - (a^2) (x^2);    b db=-(a^2)x dx ;

f[0] can be written as

f[0]= - ik (a^2)(∫_0^1) (Exp[ ((((iβ)^2)a)/k)x]-1)x dx = -i (a^2) k (-(1/2)+((k (-k+(e^((i a (β^2))/k)) (k-i a (β^2))))/((a^2) (β^4)))).

The optical theorem gives for the cross section:

σ= 2 ((a^2) π+((2 (k^2) π)/(β^4))-((2 (k^2) π Cos[((a (β^2))/k)])/(β^4))-((2 a k π Sin[((a (β^2))/k)])/(β^2))).

Let us note that the result is even in (β^2) then it is the same both for attractive ( (β^2) > 0) and repulsive ( (β^2) < 0) potentials.

In the limit k → ∞ the Born approximation is reproduced.

σ -> σB = (1/2) (a^4) π (β^4) ((((1/k)))^2).

In the formal limit  β → ∞,  k/β fixed one gets σ = 2 π (a^2), the cross section for an infinite well studied in the notebook [NB-16.2].

The check of the above manipulations is left to the reader.

Numerical results

a=1.0;

Using the function makeδAS it is easy to write a function getFandSigma which computes phase shifts and angular distribution. The sum over angular momenta has an upper limit maxL.  The function needs as input (k, θ, β), k and θ can be vectors. We define a global variable deltaGlobal which stores the last set of phase shifts.

getFandSigma[kvIn_,θvIn_,β_]:= Module[{f2,σ,maxL,δ,kv,θv},
kv = Flatten[{N[kvIn]}];
θv = Flatten[{θvIn}];
f2 = Table[0,{Length[kv]},{Length[θv]}];
σ = Table[0,{Length[kv]}];
Do[
k = kv[[j]];
maxL = Max[15,Floor[2 a 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]}];
deltaGlobal = δ;
{σ,f2}]

To fix the problem let us consider a repulsive potential.

β = 8.1 I;

Low energy

The scattering length is

αScatt =a-((Tan[a β])/β);

We define a vector of k values:

kv =Range[$MachineEpsilon,1.,0.05];

The total cross section can be computed with getFandSigma giving as input a single angle.

{σ,f2}=getFandSigma[kv,0.1,β];

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

High energy

kv =Union[Range[1.,20.,0.5],Range[20.,100.,1.5]];

{σ,f2}=getFandSigma[kv,0.1,β];

The Born approximation and the eikonal approximation are

σB[p_,β_]:=(β^4)π (-((a^2)/(8 (p^4)))+((a^4)/(2 (p^2)))-(((Sin[2 a p])^2)/(32 (p^6)))+((a Sin[4 a p])/(16 (p^5))));
σEIK[p_,β_]:= 2 ((a^2) π+((2 (p^2) π)/(β^4))-((2 (p^2) π Cos[((a (β^2))/p)])/(β^4))-((2 a p π Sin[((a (β^2))/p)])/(β^2)))

Show[ListPlot[Transpose[{kv,σ}],Joined-> True,PlotRange-> {{0,Max[kv]},All},
AxesLabel-> {"k","σ"}],
Plot[{σEIK[p,β],σB[p,β]},{p,0.5 Abs[β],Max[kv]},
PlotStyle->{ {Red,Dashed},{Blue,Dashed}}],PlotRange-> {Automatic,{0,2Max[σ]}}]

For  k≥|β((|)^2) the cross section is described by the Born approximation, the eikonal approximation works for k greater than 3 β.

A similar analysis can be done for the angular distribution. The eikonal and Born approximations are:

f2Eik[θ_,k_,β_]:= ((Abs[I k (a^2)NIntegrate[(Exp[I (((β^2)a)/k)(Sqrt[1-(x^2)])]-1)BesselJ[0,k a x Sin[θ]]x,{x,0,1}]])^2)

f2Born[θ_,k_,β_]:= Module[{q},
q= 2 k Sin[(θ/2)];
If[θ == 0.,((((1/3)(β^2)(a^3)))^2),
(((((β^2)/(q^3))(Sin[q a] - q a Cos[q a])))^2)]]

Here is an example :

k=30; β = 8.0 I;
θvalues = Range[0,0.1 Pi,0.001 Pi];

{σ,f2}=getFandSigma[k,θvalues,β];

ListPlot[{Transpose[{θvalues,f2[[1]]}],
Transpose[{θvalues,f2Eik[θvalues,k,β]}],
Transpose[{θvalues,f2Born[#1,k,β]&/@θvalues}]
}
,Joined-> True,PlotStyle->{Automatic,{Red,Dashed},{Blue,Dashed}}, 
AxesLabel->{"θ","dσ/dΩ"},PlotLabel-> Row[{"k=",k,"  β=",Abs[β]}],
PlotRange-> All]

Rainbow

The wave equation inside the potential well

Δ ψ + ((k^2)+(β^2))ψ = 0;

is formally identical to the wave equation for the light moving in a medium with refractive index

(n^2)= (((k^2)+ (β^2))/(k^2));      (β^2) = ((n^2)-1)(k^2).

The problem studied in this notebook is then equivalent to the scattering of light from a spherical droplet.

A well known phenomenon in the scattering of sun light from water droplet is the rainbow. In ray optics language the phenomenon is due to a double refraction of light rays inside a droplet. The geometry is shown in the figure below (graphics commands are in a closed cell for easier reading)

Block[{i1,r1,p1,p2,p3,p4,n1=4./3,angp2,angp3},
i1=π/4; 
 r1=ArcSin[Sin[i1]/n1];
p1={-Cos[i1],Sin[i1]};
angp2=π - i1 - (π - 2 r1);
p2={Cos[angp2],Sin[angp2]};
angp3= (π - 2 r1)-angp2;
p3={Cos[angp3],-Sin[angp3]};
def =-( π +2 i1 - 4 r1);
p4= p3 + {Cos[def],Sin[def]};
x0 =  p3[[1]]-p3[[2]]/Tan[def];
Show[Graphics[{
Circle[{0,0},1],Line[{{-2,Sin[i1]},p1,p2,p3,p4}],
Line[{{-2,0},{2,0}}],
{Dashed,Line[{{0,0},1.5p1}],Line[{{0,0},p2}],Line[{{0,0},1.5p3}],
Line[{p3,{x0,0}}]},Text["i",#1]&/@{1.4p1 -{0,0.15},1.5p3+{-0.1,0.15}},
Text["r",#1]&/@ {0.8p1+{0.12,0}, 0.8p2-{0.12,-0.05} ,0.8p2-{0.12,0.20},
   0.8p3+{0.12,0.20}},Text["θ",{x0,0}+ {0.1,-0.1}]
}]]]

The incident and refraction angles i and r are related by Snell law

Sin[r] = (1/n)Sin[i];

The angle of the incident ray change by i-r in the first refraction, then by π-2r in the internal reflection and finally by i-r when leaves the droplet.  The total scattering angle is

θ = π + 2 i - 4 r.

The incident beam consists of parallel rays with different incident angle with respect to the droplet. This family of rays form a caustic when θ is stationary with respect to variation of i (in words more rays concur to form the same θ). The equation ∂θ/∂i=0 gives easily

Cos[i]= (Sqrt[((n^2)-1)/3]);

which fixes the scattering angle for the rainbow via eq.(3.). This formula is due to Descartes. It can be generalized to higher order reflections and it is easy to show that for k internal reflections

Cos[i]= (Sqrt[((n^2)-1)/(2k+(k^2))]).

k is the order or the rainbow, below we will limit to k=1. π-θ is the angle with respect the sun direction at which an observer see the rainbow arc

If we want to reproduce this physics we have to use very short wavelengths, and for fixed index of refraction β has to change as in(3.). Our dimensionless variables are ka and βa and we see that keeping the index of refraction fixed as k→∞ means to consider bigger droplets, this is indeed the physical situation for sunlight and water droplet in the atmosphere.

The computation is not completely trivial as the backward scattering distribution is highly fluctuating. We choose as an example an index of refraction n∼1.33 which is similar to the index of refraction relevant for sunlight rainbows. Here are the relevant angles:

nInd = 4./3; 
order = 1;
incAngle = ArcCos[(((((nInd^2)-1)/(2order + order^2)))^(1/2))];
refAngle = ArcSin[Sin[incAngle]/nInd];
defAngle = Mod[2 incAngle - 2 refAngle +order(π -2refAngle),2π];
{incAngle,refAngle,defAngle}

In degrees

(π - defAngle)180/π

which match  the first bow in the sky for sunlight.

We choose a set of k and a set of angles in the backward direction

kvalues=N[{100.,500.,1000.}];
θvalues = Range[Pi/2., Pi,0.0005 Pi];

and compute the angular distribution of scattering

f2Rainbow=Table[k0= kvalues[[j]];β0=(Sqrt[(nInd^2)-1])k0;
{σ,f2}=getFandSigma[k0,θvalues,β0];f2[[1]],{j,1,Length[kvalues]}];

Here we plot the results, we see that as k grows the distribution has a growing peak moving toward the predicted value (the dashed line)

GraphicsRow[Table[ListPlot[Transpose[{θvalues,f2Rainbow[[j]]}],Joined-> True,
PlotLabel->"k = "<>ToString[N[ kvalues[[j]],3]],
PlotRange-> {{π/2,0.99π},{0,1}},AxesLabel-> {"θ",""},ImageSize-> 2.85 72,
Epilog-> {Red,Dashed,Line[{{defAngle,0},{defAngle,2}}]}],{j,1,Length[kvalues]}],Scaled[0.01]]

Airy approximation

A very simplified discussion of this problem has been given in the text (see complements to the chapter 16). For more details we refer the interested reader to the vast literature on the subject, see as an example ref.[2] and references quoted there.

In the text it has been shown that the differential cross section near a caustic has the form

((dσ)/(dΩ))= ((2π (Subscript[L, 0]))/((α^(2/3)) (k^2)))((Ai[((- (θ - (Subscript[θ, 0])))/(α^(1/3)))])^2).

Where (Subscript[θ, 0])is the classical scattering angle stationary with respect variation in L, i.e. ∂θ/∂L = 0.  The parameter α is the coefficient of the Taylor expansion at L=(Subscript[L, 0]):

θ ∼ (Subscript[θ, 0])+ α (((L-(Subscript[L, 0])))^2).

In our case the classical momentum is related to the incidence angle i by L = k a Sin[i], then we see that the Decartes construction ∂θ/∂i = 0 is exactly the construction of the caustic ∂θ/∂L = 0.

Parameters

In our notation

Sin[i] = (L/(k a));   θ[L] = π + 2 ArcSin[(L/(k a))] - 4 ArcSin[(L/(n k a)) ]

Clear[k,a,n];

theta[L_]:=  π + 2 ArcSin[(L/(k a))] - 4 ArcSin[(L/(n k a)) ];

Solve[theta'[L]== 0,L]

For (Subscript[L, 0])we have

(Subscript[L, 0])=((a k (Sqrt[4-(n^2)]))/(Sqrt[3])) = a k Sin[(Subscript[i, c])]  ; as   Cos[(Subscript[i, c])] =(Sqrt[((((n^2)-1))/3)]) ⇒ Sin[(Subscript[i, c])]= (Sqrt[(1/3)(4-(n^2))])

L0 = ((a k (Sqrt[4-(n^2)]))/(Sqrt[3]));

The Taylor expansion gives

Normal[Simplify[Series[theta[L],{L,L0,2}],Assumptions-> {n>0,a>0,k>0}]]

which gives for the parameter α

α = ((9 (Sqrt[-4+5 (n^2)-(n^4)]) )/(4 (a^2) (k^2) (((-1+(n^2)))^2)));

The Airy function has a maximum for an argument -1.02 and this explain the small shift in the position of the peak in the angular distribution:

(Subscript[θ, max])= (Subscript[θ, c])+ 1.02 (α^(1/3)).

Reflection

The formula (3.) cannot be directly used in our case, as in fact the classical motion does not have a reflection. To take account the reflection let us remember (see problems in Chapter 3) that a plane wave impinging on plane of separation between to half spaces with relative refractive index (Subscript[n, R]) has for the reflected and transmitted amplitude (Fresnel formulas):

(Subscript[A, R])= ((Cos[α]-(Subscript[n, R]) Cos[β])/(Cos[α]+(Subscript[n, R]) Cos[β])) ;  (Subscript[A, T])= ((2 Cos[α])/(Cos[α]+(Subscript[n, R]) Cos[β]));

where α and β are the incidence angle and the refraction angle respectively. In our case we have two refractions and one reflection with corresponding Fresnel coefficients

AFresnel[α_,n_]:= Module[{β},β = ArcSin[Sin[α]/n];
{((Cos[α]-n  Cos[β])/(Cos[α]+n  Cos[β])),((2 Cos[α])/(Cos[α]+n  Cos[β]))}]

coeff={AFresnel[incAngle,nInd],AFresnel[refAngle,1/nInd],AFresnel[refAngle,1/nInd]}

We take as an overall coefficient for the intensity

cReflection = (((coeff[[1,2]]*coeff[[2,1]]*coeff[[3,2]]))^2)

With our value for n, n= 4/3, the final amplitude is 8/27. This correction has to be done as the semiclassical approximation breaks down both in the reflection point and in the refraction, when the particle hit the sphere the potential has a discontinuity. In the framework of tunneling processes the first and third correction take into account the ``reflection over the barrier'' which lower the transmitted wave, the second correction is itself a reflection over  the barrier, i.e. in classical mechanics the particle cannot be reflected.

Airy approximation

a=1.;

Here we write the basic Airy approximation for the scattering cross section

f2Airy[θ_,k1_,n_]:= Module[{incAngle,refAngle,defAngle,a1},
incAngle = ArcCos[(((((n^2)-1)/(2 +1.)))^(1/2))];
refAngle = ArcSin[Sin[incAngle]/n];
defAngle = Mod[π + 2 incAngle - 4 refAngle,2π];
a1=1.;
L0 = k1  a1((Sqrt[4-(n^2)])/(Sqrt[3]));
α = (((9 (Sqrt[-4+5 (n^2)-(n^4)]) )/(4 (((-1+(n^2)))^2))))((((1/(k1 a1))))^2);
(( 2π L0)/((α^(2/3))(k1^2)))(1/(Sin[defAngle]))((AiryAi[- (((θ - defAngle))/(α^(1/3)))])^2)]

We compute the exact numerical solution and make a comparison with this approximation.

Warning: We need a very high k , so the computation takes some time, if the reader in interested only in the qualitative aspects of the problem can change to 1000 or 2000 the value of k defined below, in kvalues.

kvalues=N[{5000.}];
θvalues = Range[2.35, 2.55,1./(10^3)];

This is the angular distribution of scattering

f2Rainbow2=Table[k0= kvalues[[j]];β0=(Sqrt[(nInd^2)-1])k0;
{σ,f2}=getFandSigma[k0,θvalues,β0];f2[[1]],{j,1,Length[kvalues]}];

{mt1,mt2}={Min[θvalues],Max[θvalues]};

ListPlot[
{Transpose[{θvalues,f2Rainbow2[[1]]}],
Transpose[{θvalues,cReflection f2Airy[θvalues,kvalues[[1]],nInd]}]},Joined-> True,
PlotRange-> {{mt1,mt2}, All},PlotStyle-> {Automatic,{Red,Dashed}}]

We can explicitly check that the (small) discrepancy is due to interference with neglected terms in the saddle point approximation. We can sum only angular momenta near (Subscript[L, 0])in the numerical solution:

saddlePointSum[θv_,L_,k_]:= Module[{L1,L2,mom,phase1},
L1 = Floor[0.8L]+1; L2 = Floor[1.2L]+1;
phase1=( 2Range[L1,L2]+1)*(Exp[2I deltaGlobal[[L1;;L2]]]-1);
Abs[(1/(2I k))Tr[phase1*legP[#][[L1;;L2]]]]^2&/@θv]

This function use the global variable deltaGlobal defined in getFandSigma, it stores the values of the phase shifts for the last value of k (here we have only one k).

f2gauss = saddlePointSum[θvalues,L0,kvalues[[1]]];

ListPlot[
{Transpose[{θvalues,f2gauss}],
Transpose[{θvalues,cReflection f2Airy[θvalues,kvalues[[1]],nInd]}]},Joined-> True,
PlotRange-> {{mt1,mt2}, All},PlotStyle-> {Automatic,{Red,Dashed}}]

The agreement with Airy approximation in the peak region is better, this means that the Airy approximation, for this value of k, miss a small contribution from corrections to saddle point.

References