Simple numerical techniques for quantum scattering¶

Jay Wang and Owen Tower, UMass Dartmouth, www.faculty.umassd.edu/j.wang/¶

• #### Scattering: important but challenging

• common in nature, key to discovery of fundamental forces
• rich physics, continuum, asymptotics (not normalizable, even poles)
• requiring considerable math phys incl. special funcs, multipoles
• #### Computation: breaking down the barrier

• makes scattering accessible & broadens class scope
• lets us explorer wave function, phase shift, real effects
• offers in-depth hands-on activities, projects/asynchronous learning
• #### Scattering in central field potentials
• Numerical partial wave analysis with Numerov's method
• Example: Yukawa, ie, screened Coulomb, potential
• Wave function, phase shift, cross sections, nonclassical effect
• #### talk available at PICUP virtual conf site (.html and zipped .ipynb)

Partial wave analysis in a nutshell¶

• #### Radial Schrödinger Eq., $u_l'' + qu=0, \quad q= 2\left(E-V - \frac{l (l+1)\hbar^2}{2mr^2}\right), \quad u_l(E,r)=r R_l(E,r)$
• #### Incident and scattered waves:
• #### Scattering amplitude: $f(\theta)= \frac{1}{k} \sum_{l=0}^{\infty} (2l+1) e^{i\delta_l} \sin \delta_l \; P_l(\cos\theta)$
• #### Phase shift: $\tan \delta_l = \frac{(\gamma a -1) j_l(ka) - k a \, j_l'(ka)} {(\gamma a -1) n_l(ka) - k a \, n_l'(ka)}, \quad \gamma = \frac{u'_l(k,a)}{u_l(k,a)}$
• Requires wave function and its derivative in asymptotic region
• $u_l$ and $u'_l$ to be computed numerically
• Cross section: $\sigma(\theta) = |f(\theta)|^2$

Numerov's method: highly accurate, simple, efficient¶

• Given $u''+q(x)\, u=0$, $h=$ grid size, $c=\frac{h^2}{12}$, 3-term recursions for $u$ & $u'$:
• $u_{i+1} = \left \{ 2 \left(1 - 5c\, q_i \right) u_i - \left(1 + c\, q_{i-1} \right) u_{i-1} \right\}/ \left( 1+ c\, q_{i+1} \right ) + O(h^6)$
• $u_i' = \frac{1}{2h} \left \{ \left( 1+ 2c\, q_{i+1}\right ) u_{i+1} - \left( 1+ 2c\, q_{i-1}\right ) u_{i-1} \right \} + O(h^4)$
• Iteration: start with $u(0)=0, \; u(h)=1$, obtain $u(2h) \cdots u(nh)$
In [1]:
def numerov(q, u, n, x, h):     # Numerov integrator for u''+q(x)u=0
nodes, c = 0, h*h/12.       # given [u_0,u_1], return [u_0,u_1,...,u_{n+1}]
f0, f1 = 0., q(x+h)
for i in range(n):
x += h
f2 = q(x+h)             # Numerov method below
u.append((2*(1-5*c*f1)*u[i+1] - (1+c*f0)*u[i])/(1+c*f2))
f0, f1 = f1, f2
if (u[-1]*u[-2] < 0.0): nodes += 1
return u, nodes

In [2]:
# needed libraries
import scipy
if scipy.__version__[0] < '1':  # sph. Bessel func changed after ver 1.0.0
from scipy.special import sph_jn, sph_yn, lpn
else:
from scipy.special import spherical_jn as sph_jn, spherical_yn as sph_yn, lpn
import matplotlib.pyplot as plt, numpy as np
#%matplotlib notebook


Yukawa, ie screened Coulomb, potential, $V = -\frac{Z}{r} \exp(-r/s)$¶

In [3]:
def V(r, Z=2.0, s=1.0):                       # Yukawa potential
return -Z*np.exp(-r/s)/r

def q(r):                       # Sch eqn in Numerov form
return 2*(E - V(r)) - L*(L+1)/(r*r)

def wf(M, xm):                  # find w.f. and deriv at xm
c = (h*h)/6.
wfup, nup = numerov(q, [0,.1], M, 0., h)    # 1 step past xm
dup = ((1+c*q(xm+h))*wfup[-1] - (1+c*q(xm-h))*wfup[-3])/(h+h)
return wfup, dup/wfup[-2]


Yukawa potential, larger $s =$ more Coulomb like¶

In [4]:
plt.figure()
r = np.linspace(.01, 5,100)
plt.plot(r, V(r, Z=2., s=1.0))
plt.plot(r, V(r, Z=2., s=15.0))
plt.xlabel('r (a.u.)'), plt.ylabel('V (a.u.)')
plt.ylim(-5,0);


Computation¶

In [5]:
a, M, uL = 10., 200, []                 # matching point, steps
h, Lmax, E =a/M, 15, 2.                 # step size, max L, energy
k, ps = np.sqrt(2*E), np.zeros(Lmax+1)  # wave vector, phase shift
if scipy.__version__[0] < '1':
jl, dj = sph_jn(Lmax, k*a)          # (j_l, j_l') tuple
nl, dn = sph_yn(Lmax, k*a)          # (n_l, n_l')
else:
Lrange = range(Lmax + 1)
jl, dj = sph_jn(Lrange, k*a, False), sph_jn(Lrange, k*a, True) # (j_l, j_l')
nl, dn = sph_yn(Lrange, k*a, False), sph_yn(Lrange, k*a, True) # (n_l, n_l')

for L in range(Lmax+1):
u, g = wf(M, a)                     # g= u'/u
uL.append(u)
x = np.arctan(((g*a-1)*jl[L] - k*a*dj[L])/    # phase shift
((g*a-1)*nl[L] - k*a*dn[L]))
while (x < 0.0): x += np.pi         # handle jumps by pi
ps[L] = x
theta, sigma = np.linspace(0., np.pi, 100), []
cos, La = np.cos(theta), np.arange(1,2*Lmax+2,2)
for x in cos:                               # calc cross section
pl = lpn(Lmax, x)[0]                    # Legendre polynomial
fl = La*np.exp(1j*ps)*np.sin(ps)*pl     # amplitude
sigma.append(np.abs(np.sum(fl))**2/(k*k))


Results: wave function¶

In [6]:
plt.figure()                                # plot scaled w.f.
uL = np.array(uL)
for L in range(0, 5, 2):
plt.plot(h*np.arange(len(uL[L])), uL[L]/max(abs(uL[L])))
plt.xlabel('r (a.u.)'), plt.ylabel('$u$', fontsize=16);


Results: phase shift: positive $\delta_l$; larger $\delta_l$ for smaller $\ell$¶

In [7]:
plt.figure()                                # plot phase shift vs L
plt.plot(range(Lmax+1), ps, '-o')
plt.xlabel('$\ell$', fontsize=20), plt.ylabel(r'$\delta_\ell$', fontsize=20);


Results: cross section¶

In [8]:
plt.figure()
plt.plot(theta*180/np.pi, sigma)            # plot cross sections
xtck = [0, 30, 60, 90, 120, 150, 180]
plt.xticks(xtck, [repr(i) for i in xtck])   # custom ticks
plt.xlabel(r'$\theta$ (deg)'), plt.ylabel(r'$\sigma(\theta)$ (a.u.)'), plt.semilogy();


Summary¶

• #### Scattering with Numerov's method: simple & efficient
• #### Accessible & enriched class with computation
• #### Concepts of phase shifts, nonclassical Ramsauer-Townsend
• #### Active engagement, suitable to projects & online learning
• #### Other potentials, eg hard/fuzzy spheres, www.faculty.umassd.edu/j.wang/