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

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

# 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)$

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

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

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

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$

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

Comparison with field-free wf: larger effects for smaller $\ell$:

Results: cross section

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();

Ramsauer-Townsend effect at lower $E$, nonclassical - interference of few $\ell$-waves

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/

Hard sphere at different energies