Phase Shift Extraction Notebook

This notebook extracts elastic partial-wave phase shifts by solving the radial Schrödinger equation and matching the numerical solution to known free asymptotic functions. It turns the formal definition of $\delta_\ell$ into a reproducible scattering calculation.

Radial Scattering from a Gaussian Potential

Use a short-range central potential, for example

$$V(r)=V_0e^{-r^2/a^2}.$$

For reduced mass $\mu$ and energy

$$E=\frac{\hbar^2k^2}{2\mu},$$

the reduced radial wavefunction obeys

$$- \frac{\hbar^2}{2\mu} \frac{d^2u_\ell}{dr^2} + \left[ V(r) + \frac{\hbar^2\ell(\ell+1)}{2\mu r^2} \right] u_\ell = Eu_\ell.$$

The regular origin behavior is

$$u_\ell(r)\propto r^{\ell+1}$$

as $r\to0$.

Dimensionless Form

Set the range $a$ as the unit of length:

$$\rho=\frac{r}{a}.$$

Use the dimensionless energy and strength

$$\epsilon=k^2a^2, \qquad \nu=\frac{2\mu V_0a^2}{\hbar^2}.$$

Then the radial equation becomes

$$- \frac{d^2u_\ell}{d\rho^2} + \left[ \nu e^{-\rho^2} + \frac{\ell(\ell+1)}{\rho^2} \right] u_\ell = \epsilon u_\ell.$$

This form makes parameter sweeps independent of units.

Numerical Setup

For each $\ell$ and $k$:

1. choose a small starting radius $\rho_{\min}$;
2. initialize the regular solution with $u_\ell\propto\rho^{\ell+1}$;
3. integrate outward to a matching radius $\rho_m$ where $V(\rho_m)$ is negligible;
4. compute the logarithmic derivative

$$L = \frac{u_\ell'(r_m)}{u_\ell(r_m)};$$

5. match to free spherical solutions to extract $\delta_\ell$.

Numerov integration is a common choice for this second-order equation. A high-order adaptive ODE solver can also be used if the normalization and phase are handled consistently.

Matching Formula

Outside the potential range, write the reduced radial solution as

$$u_\ell(r) = C \left[ \hat j_\ell(kr)\cos\delta_\ell - \hat n_\ell(kr)\sin\delta_\ell \right],$$

where

$$\hat j_\ell(z)=zj_\ell(z), \qquad \hat n_\ell(z)=zn_\ell(z)$$

are Riccati–Bessel functions.

The conventions for $j_\ell$, $n_\ell$, and their Riccati forms are summarized in Bessel Functions.

Let

$$z=kr_m, \qquad J=\hat j_\ell(z), \qquad N=\hat n_\ell(z),$$

and let $J'$ and $N'$ denote derivatives with respect to $z$. The logarithmic derivative matching gives

$$\tan\delta_\ell = \frac{kJ'-LJ}{kN'-LN}.$$

This formula is insensitive to the arbitrary normalization of the numerical solution. Avoid dividing by $u_\ell$ at a node: instead compute

$$\begin{aligned} A&=ku_\ell(r_m)J'-u_\ell'(r_m)J,\\ B&=ku_\ell(r_m)N'-u_\ell'(r_m)N, \end{aligned}$$

and use $\delta_\ell=\operatorname{atan2}(A,B)$ modulo $\pi$. If the ODE was integrated in $\rho$, its returned derivative is $du/d\rho$, so use $ka$ in place of $k$ and $du/d\rho$ in place of $du/dr$ in both expressions. Mixing these derivatives gives a dimensionally inconsistent phase.

Run an s-Wave Example

This standard-library Python calculation uses fourth-order Runge–Kutta integration and the regular origin series. It computes one phase modulo $\pi$; an energy sweep crossing a bound-state branch or resonance needs phase unwrapping. The variable K is $ka$ and nu is the dimensionless strength above.

from math import atan2, ceil, cos, exp, sin

def phase_s(K, nu, step=0.002, match=6.0, start=1e-4):
    """Gaussian s-wave phase modulo pi; K=ka, nu=2*mu*V0*a**2/hbar**2."""
    n = ceil((match-start)/step)
    h = (match-start)/n
    r = start
    # Regular series u=r+(nu-K**2)*r**3/6, with one arbitrary common scale.
    c = nu-K*K
    u, v = r+c*r**3/6, 1+c*r*r/2
    def rhs(r, u, v):
        return v, (nu*exp(-r*r)-K*K)*u
    for _ in range(n):
        a,b = rhs(r,u,v)
        c,d = rhs(r+h/2,u+h*a/2,v+h*b/2)
        e,f = rhs(r+h/2,u+h*c/2,v+h*d/2)
        g,j = rhs(r+h,u+h*e,v+h*f)
        u += h*(a+2*c+2*e+g)/6
        v += h*(b+2*d+2*f+j)/6
        r += h
    z = K*match
    # Match u=C*sin(z+delta), v=K*C*cos(z+delta), without dividing by u.
    numerator = K*u*cos(z)-v*sin(z)
    denominator = K*u*sin(z)+v*cos(z)
    return atan2(numerator, denominator)

print(phase_s(0.3, -0.5))

At $ka=0.3$, $\rho_m=6$, and $\rho_{\min}=10^{-4}$, the result is:

Strength $\nu$	$h=0.02$	$h=0.01$	$h=0.005$
$0$	$0$	$0$	$0$
$+0.5$	$-0.0546827127$	$-0.0546827127$	$-0.0546827127$
$-0.5$	$+0.0759139869$	$+0.0759139870$	$+0.0759139870$

Phases are in radians; the free case is zero to the displayed precision. Repulsion gives a negative phase and weak attraction a positive one. Their unequal magnitudes show that the exact answer contains more than the first Born term, which changes sign without changing magnitude. Agreement under step refinement tests integration error. Separately vary match and start to check the neglected Gaussian tail and the truncated origin series.

Parameter Sweep

A useful first sweep is:

$$\ell=0,1,2, \qquad ka=0.1,0.3,1,3, \qquad \nu=-0.5,-2,-5.$$

At weak coupling, compare the numerical phase shifts with the Born-estimate sign and scale. At low $k$, verify $s$-wave dominance unless the potential is tuned near a threshold state.

Plots That Reveal the Physics

Use the calculated phases to plot:

• $\delta_\ell(k)$ versus $k$ for several $\ell$;
• $k\cot\delta_0(k)$ versus $k^2$ for extracting the scattering length and effective range;
• convergence of $\delta_\ell$ with matching radius;
• convergence with radial step size;
• optionally, the numerical radial wavefunction and its matched asymptotic form.

Phase shifts should be unwrapped or plotted with branch choices stated clearly.

Validation Checks

For $V=0$, the extracted phase shifts must be

$$\delta_\ell=0$$

up to numerical error.

For a sufficiently short-range potential away from an exact threshold state, choose the phase branch near an integer multiple of $\pi$ and check

$$\delta_\ell(k)-n_\ell\pi\propto k^{2\ell+1}.$$

For $\ell=0$, extract the scattering length from

$$a = - \lim_{k\to0} \frac{\tan\delta_0(k)}{k}.$$

The total elastic cross section computed from phase shifts,

$$\sigma_{\mathrm{el}} = \frac{4\pi}{k^2} \sum_{\ell=0}^{\ell_{\max}} (2\ell+1)\sin^2\delta_\ell,$$

should stabilize as $\ell_{\max}$ increases.

Before using a smooth potential, validate the implementation against Square-Well Scattering and Hard-Sphere Scattering. The first checks interface matching, scattering-length extrapolation, and resolution of a narrow $p$-wave feature. The second isolates the exterior boundary condition and tests partial-wave cutoff convergence from threshold to $kR\gg1$. Low-Energy S-Wave Scattering then tests extraction of $k\cot\delta_0$, the effective range, the unitarity window, and a shallow pole. Resonance from a Square Well tests phase unwrapping, background-aware fitting, and continuation to an outgoing $p$-wave pole.

Born Approximation Numerical Test specializes the same matching idea to a repulsive Gaussian and supplies a complete convergence study before using the phases to judge first Born accuracy.

Exercises

1. Why does logarithmic-derivative matching avoid dependence on the arbitrary normalization of $u_\ell$?

Solution

If $u_\ell$ is multiplied by a constant $C$, then both $u_\ell'$ and $u_\ell$ are multiplied by $C$. Their ratio

$$L=\frac{u_\ell'}{u_\ell}$$

is unchanged. Since the matching formula uses $L$, the extracted phase shift is independent of the normalization chosen during integration.

2. Why should the free-potential test give $\delta_\ell=0$?

Solution

When $V=0$, the regular radial solution is proportional to the free Riccati–Bessel function $\hat j_\ell(kr)$. In the matching form, this corresponds to $\cos\delta_\ell=1$ and $\sin\delta_\ell=0$, hence $\delta_\ell=0$ modulo $\pi$. The branch connected continuously to zero potential is $\delta_\ell=0$.

References

• B. R. Johnson, “The multichannel log-derivative method for scattering calculations,” Journal of Computational Physics 13, 445–449 (1973), doi:10.1016/0021-9991(73)90049-1.
• R. G. Newton, Scattering Theory of Waves and Particles, 2nd ed., Dover, 2002.
• J. R. Taylor, Scattering Theory: The Quantum Theory of Nonrelativistic Collisions, Dover, 2006.
• J. M. Thijssen, Computational Physics, 2nd ed., Cambridge University Press, 2007.