Original source: https://notebookarchive.org/scattering-cross-section--2021-02-1w06xy1/
James A Miller (Physics & Astronomy, University of Alabama in Huntsville)
Scattering experiments are critical for probing nuclear and atomic structure and interactions. The central diagnostic tool is the cross section, which can be measured experimentally as well as calculated theoretically on the basis of a physical model. In this notebook, we demonstrate the numerical calculation of the cross section using partial wave phase shifts.
Suppose particles of momentum ℏ k are directed along the z axis toward a stationary potential or scattering center, from which they can scatter elastically. Further assume a central potential V(r) and azimuthal symmetry in the scattered particles. The differential cross section dσ(θ)/dΩ is a measure of the number of particles scattered into the solid angle dΩ about θ, the scattering angle measured from the z axis. The total cross section σ is obtained by integrating the differential cross section over all angles, and is the effective interaction area. Specifically, an incident particle moving through the cross sectional area centered on the potential will suffer scattering. The incident particles are states of definite linear momentum. Decomposing this state into states of definite angular momentum, the radial or r-dependent part of each angular momentum state suffers a shift in phase when the potential is introduced. These quantities are of fundamental importance, and can be used to calculate cross sections. A full discussion of partial waves can be found in standard quantum mechanics texts\^(1)
The first step to calculating a cross section is then finding these phase shifts. While phase shifts are not restricted to this case, we suppose the scattering potential V(r) is nonzero inside some radius Subscript[r, a] and zero outside. Phase shifts Subscript[δ, ℓ] are obtained by matching the radial wavefunction inside the potential with that outside. The exterior radial wavefunction (where V = 0) is exact and given in terms of the spherical Bessel functions by u(r) = r [ Subscript[A, ℓ] Subscriptj, ℓ + Subscript[B, ℓ] Subscriptn, ℓ]. Here, the quantity ℓ = 0,1, 2, … is the angular momentum quantum number. A phase shift is defined by tan Subscript[δ, ℓ ]= -Subscript[B, ℓ]/Subscript[A, ℓ]. We could match at the single point Subscript[r, a] through the logarithmic derivative, for example. But to avoid taking derivatives of a function given at discrete points, we instead fit the numerical solution with the exact solution outside the potential. (The numerical solution inside can easily be continued into the exterior region by extending the potential.) Hence, solve the radial Schrodinger equation for u(r) over the interval [Subscript[r, 1] > 0, Subscript[r, 2] > Subscript[r, a]]. Since the potential is zero in the interval [Subscript[r, 2], Subscript[r, a]], we can match the numerical solution with u(r) = r [ Subscript[A, ℓ] Subscriptj, ℓ + Subscript[B, ℓ] Subscriptn, ℓ] to obtain Subscript[A, ℓ], Subscript[B, ℓ]. The normalization of u(r) is arbitrary since we just need the ratio Subscript[B, ℓ]/Subscript[A, ℓ]. Schematically:
We typically need Subscript[r, 1] >0 since the centrifugal, and possibly the physical, potential diverge at the origin. The numerical solution of the radial equation requires the starting values u(Subscript[r, 1]) and u'(Subscript[r, 1]), which can then be obtained by a small-r solution of the radial equation.
| ------ | ----- | ----- | ----- | … | ----- | ----- | ----- |
r=0 r_1 r_a r_2
Numerical u(r)
Analytical exterior solution -->
Fitting region
To illustrate the method, consider a finite square well potential V(r) = -Subscript[V, 0], for r< Subscript[r, a], and zero otherwise. Use natural units ℏ = m = Subscript[r, a] = 1. Find the S-wave (low-energy, ℓ = 0) cross section for k = 1.
Specify the potential using the parameters of Merzbacher, Section 13.6:
V0=19.22;
ra = 1.0;
V=If[r <=ra,-V0, 0];
r1 = 0.01;
r2 = 2.0;
ell = 0;
Give the small-r forms of u(r) and u'(r) from the radial equation:
us[r_,ell_]:= r^(ell+1);
usprime[r_,ell_]:= (ell+1)r^ell;
Solve the radial equation out to Subscript[r, 2] for the unnormalized wavefunction:
k=1.0;
energy=k^2/2;
uE= NDSolveValue[{-(1/2)u''[r] + ((ell(ell+1))/(2 r^2)+V)u[r]==energy u[r], u[r1]== us[r1,ell],u'[r1]==usprime[r1,ell]}, u,{r,r1,r2}]
Select 10 points from [Subscript[r, a],Subscript[r, 2]]:
rArray=Subdivide[ra,r2,9];
points=Transpose[{rArray,uE/@rArray}]
Fit these points with the exact wavefunction and get the phase shift:
δ=ArcTan[-B/A]/.FindFit[points,A r SphericalBesselJ[ell,k r]+B r SphericalBesselY[ell,k r],{A,B},r];
Find the partial wave cross section at this wavenumber:
(4 π)/k^2 (2 ell + 1) Sin[δ]^2
Divide by π to compare with Merzbacher:
%/π
This is correct, as is the cross section for other values of k.
To be most useful, cross sections need to be presented for a range of incident particle energies. It is thus convenient to incorporate the above steps into a function.
Function to accept the potential, wavenumber, angular momentum state, and integration interval:
sigma[potential_,ra_?NumberQ,r1_?NumberQ,r2_?NumberQ,u1_?NumberQ,u1prime_?NumberQ,ell_?IntegerQ,k_?NumberQ]:=
Module[{energy,u,uE,rArray,points,δ,A,B },
energy=k^2/2;
uE= NDSolveValue[{-(1/2)u''[r] + ((ell(ell+1))/(2 r^2)+potential)u[r]==energy u[r], u[r1]== u1,u'[r1]==u1prime}, u,{r,r1,r2}];
rArray=Subdivide[ra,r2,9];
points=Thread[{rArray,uE/@rArray}];
δ=ArcTan[-B/A]/.FindFit[points,A r SphericalBesselJ[ell,k r]+B r SphericalBesselY[ell,k r],{A,B},r];
(4 π)/k^2 (2 ell + 1) Sin[δ]^2
];
Consider the same finite square well example as above, but consider the S and P (ℓ = 1) wave contributions to the total cross section.
Specify the angular momentum quantum numbers and the starting u and u' values:
ell0=0;ell1=1;
u10 = us[r1,ell0]; u1prime0 = usprime[r1,ell0];
u11 = us[r1,ell1];u1prime1 = usprime[r1,ell1];
Plot the partial wave cross sections:
Plot[{sigma[V,ra,r1,r2,u10,u1prime0,ell0,k]/π,sigma[V,ra,r1,r2,u11,u1prime1,ell1,k]/π},{k,0.01,4}, PlotStyle -> Thick, PlotTheme -> "Scientific",
LabelStyle -> {FontFamily -> "Arial", 16, GrayLevel[0]}, Frame -> True, ImageSize -> Large,
PlotRange -> {0, 30}, FrameLabel -> {{"\!\(\*SubscriptBox[\(\[Sigma]\), \(\[ScriptL]\)]\)(k)/\[Pi]", ""},
{"k", "S and P Wave Cross Sections"}}, RotateLabel -> True, ClippingStyle -> Automatic,
GridLines -> Automatic, PlotLegends -> LineLegend[{"\[ScriptL]=0", "\[ScriptL]=1"}]]
The P-wave cross section exhibits a resonance at k ≈ 0.62. Resonances are important diagnostics of the scattering potential and result from the formation of a quasi-bound state during the interaction. They are evident as approximate standing waves in a region of the effective (centrifugal plus physical) potential.
Find the radial wavefunction at resonance:
k = 0.62;
ell = ell1;
energy=k^2/2;
uE= NDSolveValue[{-(1/2)u''[r] + ((ell(ell+1))/(2 r^2)+V)u[r]==energy u[r], u[r1]== u11,u'[r1]==u1prime1}, u,{r,r1,r2}]
Plot the potential and the radial wavefunction:
Plot[{((ell(ell+1))/(2 r^2)+V),80 uE[r]},{r,r1,r2},Epilog -> {Dashed, Line[{{ra, -20}, {ra, 20}}]},
PlotRange -> {-20, 20}, Frame -> True, PlotTheme -> "Scientific",
LabelStyle -> {FontFamily -> "Arial", 16, GrayLevel[0]},
FrameLabel -> {{"u(r),V(r)", ""}, {"r", "Radial Standing Wave"}}, RotateLabel -> True, ImageSize -> Large,
PlotLegends -> LineLegend[{"V(r)", "u(r)"}]]
We see that approximately 3/4 of a wavelength fit inside the well, forming a temporary standing wave analogous to a sound wave in an semi-open ended pipe.
Introduction to Quantum Mechanics, by David Griffiths (ISBN-13: 978-1107189638)
Quantum Mechanics, by Eugen Merzbacher (ISBN-13: 978-0471887027)
Modern Quantum Mechanics, by J. J. Sakurai & J. Napolitano (ISBN-13: 978-1108422413)