The WKB approximation: tunneling and alpha decay

PHYS-372 solved the Schrödinger equation for potentials made of flat pieces, wells and barriers with sharp edges, and for a few smooth potentials with special mathematical structure. Most potentials are neither. When the potential changes slowly compared with the de Broglie wavelength, there is an approximation that captures the physics directly and connects it with classical mechanics. It was developed in 1926 by Gregor Wentzel, Hendrik Kramers and Léon Brillouin, and is known by their initials, WKB, although Harold Jeffreys had found the essential mathematics in 1923. This lesson derives it, uses it to quantize bound states, and applies it to the problem that first made tunneling famous: the alpha decay of heavy nuclei, solved by George Gamow, and independently by Ronald Gurney and Edward Condon, in 1928.

The local wavelength

In one dimension, write the time-independent Schrödinger equation as \(\psi'' = -(p^2/\hbar^2)\psi\), where \(p(x) = \sqrt{2m[E - V(x)]}\) is the classical momentum a particle of energy \(E\) would have at \(x\). If \(V\) were constant, the solutions would be plane waves \(e^{\pm ipx/\hbar}\). When \(V\) varies slowly, it is natural to guess that the solution is still locally a wave, with a wavelength \(\lambda(x) = h/p(x)\) that changes from place to place. Write \(\psi(x) = A(x)e^{i\phi(x)}\) with real amplitude and phase. Substituting and separating real and imaginary parts gives two equations:

The second says \(A^2\phi'\) is constant, so \(A \propto 1/\sqrt{\phi'}\). The first is exact but hard; the WKB approximation drops \(A''\), assuming the amplitude varies slowly compared with the phase. Then \(\phi' = \pm p/\hbar\), and

The amplitude factor has a classical meaning: \(|\psi|^2 \propto 1/p \propto 1/v\), so the particle is most likely to be found where it moves slowly, exactly as for a classical particle that spends more time where it is slow. The condition for the approximation is that the wavelength changes little over one wavelength, \(|d\lambda/dx| \ll 1\), or equivalently that the potential changes by much less than the kinetic energy over a wavelength.

Inside a barrier

In a classically forbidden region, where \(E < V\), the momentum is imaginary. Writing \(\kappa(x) = \sqrt{2m[V(x) - E]}/\hbar\), the same reasoning gives growing and decaying exponentials:

For a barrier wide enough that the decaying solution dominates, the probability of getting through is the square of the ratio of amplitudes on the two sides:

where \(x_1\) and \(x_2\) are the classical turning points, the edges of the forbidden region. This is the most used result of the WKB method. Its strength is the exponential: the transmission depends on the barrier's height and width only through the integral \(\gamma\), and because \(\gamma\) appears in an exponent, small changes in energy or width change \(T\) by large factors. The prefactor is not captured; for a rectangular barrier the exact result of PHYS-372 is \(T \approx [16E(V_0 - E)/V_0^2]\,e^{-2\kappa a}\) for thick barriers, and the program shows that WKB misses only that factor, 2.56 for the case computed, while correctly giving \(T\) over sixteen orders of magnitude.

Turning points and the connection formulas

At a turning point, \(p \to 0\) and the WKB amplitude \(1/\sqrt p\) blows up: the approximation fails exactly where the classical particle stops and turns around. Near such a point the potential is nearly linear, \(V(x) \approx E + V'(x_0)(x - x_0)\), and the Schrödinger equation with a linear potential can be solved exactly in terms of Airy functions. Matching the Airy solution to the WKB forms on both sides, a procedure worked out by Kramers, gives connection formulas. The most important says that a solution decaying into the forbidden region on the right of a turning point at \(x_2\) joins a standing wave on the left with a definite phase:

The extra phase of \(\pi/4\) at each smooth turning point is the quantum correction to the classical picture of a particle bouncing off a wall.

Bound states: Bohr and Sommerfeld revisited

For a particle trapped between two smooth turning points, the wave function must decay into both forbidden regions. Applying the connection formula at each end, the two standing-wave forms inside agree only if the phase accumulated between the turning points satisfies

or, written as the integral of \(p\) around a full classical orbit, \(\oint p\,dx = (n + \tfrac12)h\). This is the old quantum condition of Bohr and Arnold Sommerfeld (1913 and 1916) with the half-integer correction that the full theory requires. For the harmonic oscillator it gives \(E_n = (n + \tfrac12)\hbar\omega\) exactly, a coincidence of that potential. For other smooth potentials it is accurate for large \(n\), where the wavelength is short compared with the scale of the potential, and usually good even for the ground state. If one wall is vertical, the wave function must vanish there, and the phase correction at that end becomes \(\pi/2\) instead of \(\pi/4\); with two vertical walls the condition gives the exact infinite-well energies.

Alpha decay

Heavy nuclei such as uranium-238 emit alpha particles, helium-4 nuclei, with kinetic energies of 4 to 9 MeV. The puzzle in the 1920s was twofold. Rutherford had scattered alpha particles of about 8.8 MeV off uranium and found that they did not reach the nucleus: the Coulomb barrier there is higher than 8.8 MeV, yet uranium emits alpha particles of only 4.2 MeV. And the half-lives depended on energy in an extraordinary way: as the alpha energy rises from about 4 to 9 MeV, half-lives fall from billions of years to fractions of a microsecond, the empirical Geiger–Nuttall law of 1911. Gamow's explanation was that the alpha particle, preformed inside the nucleus, tunnels out through the Coulomb barrier.

Model the potential seen by the alpha particle as a well inside the nuclear radius \(R\) and the Coulomb repulsion \(V(r) = 2Z'e^2/4\pi\varepsilon_0r\) outside, where \(Z'\) is the charge of the daughter nucleus. The alpha particle of energy \(E\) leaves the barrier at the outer turning point \(r_c = 2Z'e^2/4\pi\varepsilon_0E\), tens of femtometres out. The Gamow integral can be done exactly:

where the approximation holds for \(R \ll r_c\). The leading term is proportional to \(Z'/\sqrt E\). The alpha particle rattles back and forth inside the nucleus at a speed \(v\), striking the barrier about \(v/2R\) times per second, roughly \(10^{21}\) per second, and escapes with probability \(e^{-2\gamma}\) at each attempt. The mean lifetime is therefore \(\tau = (2R/v)\,e^{2\gamma}\), and

which is the Geiger–Nuttall law, with constants that the theory predicts. The program evaluates the full Gamow expression for six alpha emitters, using a nuclear radius \(R = 1.2\,(A_d^{1/3} + 4^{1/3})\) fm, the distance at which the alpha particle and daughter nucleus touch.

import math
hbarc = 197.327          # MeV fm
mac2 = 3727.379          # alpha-particle rest energy, MeV
e2 = 1.439965            # e^2 / (4 pi eps0) in MeV fm
c = 2.99792458e23        # fm per second

# 1. WKB transmission through a rectangular barrier compared with the exact result (electron, eV and nm)
def exact_T(E, V0, a, mc2=510998.95, hc=197.327):
    k = math.sqrt(2 * mc2 * E) / hc; q = math.sqrt(2 * mc2 * (V0 - E)) / hc   # per nm when hc in eV nm
    return 1 / (1 + V0 ** 2 * math.sinh(q * a) ** 2 / (4 * E * (V0 - E)))
print("Electron transmission through a 5 eV barrier: exact T versus the WKB estimate exp(-2 q a)")
for a in (0.2, 0.5, 1.0, 2.0):
    E, V0 = 1.0, 5.0
    q = math.sqrt(2 * 510998.95 * (V0 - E)) / 197.327
    print(f"  width {a:3.1f} nm: exact {exact_T(E, V0, a):.3e}   WKB {math.exp(-2 * q * a):.3e}   ratio {exact_T(E, V0, a) / math.exp(-2 * q * a):.2f}")
print()

# 2. Gamow theory of alpha decay
def half_life(Z, A, Q, r0=1.2):
    Zd, Ad = Z - 2, A - 4
    E = Q * Ad / A                                    # alpha kinetic energy after recoil
    mu = mac2 * Ad * 931.494 / (mac2 + Ad * 931.494)  # reduced mass, MeV
    R = r0 * (Ad ** (1 / 3) + 4 ** (1 / 3))           # touching radius, fm
    rc = 2 * Zd * e2 / E                              # classical turning point, fm
    x = R / rc
    gamma = math.sqrt(2 * mu * E) / hbarc * rc * (math.acos(math.sqrt(x)) - math.sqrt(x * (1 - x)))
    v = c * math.sqrt(2 * E / mu)                     # inside speed, approximated by the outside speed
    tau = 2 * R / v * math.exp(2 * gamma)             # mean life = (time between wall hits) / T
    return math.log(2) * tau, gamma, rc

YEAR = 3.15576e7
data = [("U-238", 92, 238, 4.270, 4.468e9 * YEAR), ("Th-232", 90, 232, 4.082, 1.405e10 * YEAR),
        ("Ra-226", 88, 226, 4.871, 1600 * YEAR), ("Po-210", 84, 210, 5.407, 138.4 * 86400),
        ("Po-214", 84, 214, 7.833, 164.3e-6), ("Po-212", 84, 212, 8.954, 0.299e-6)]
print("Alpha decay: WKB (Gamow) estimate versus measured half-life")
print("  nucleus   Q (MeV)  turning point (fm)  2 gamma   log10 t(s) model  measured")
for name, Z, A, Q, t_meas in data:
    t, g, rc = half_life(Z, A, Q)
    print(f"  {name:7s}  {Q:6.3f}       {rc:6.1f}         {2 * g:6.1f}      {math.log10(t):7.2f}      {math.log10(t_meas):7.2f}")

Electron transmission through a 5 eV barrier: exact T versus the WKB estimate exp(-2 q a)
  width 0.2 nm: exact 4.209e-02   WKB 1.660e-02   ratio 2.54
  width 0.5 nm: exact 9.085e-05   WKB 3.549e-05   ratio 2.56
  width 1.0 nm: exact 3.224e-09   WKB 1.259e-09   ratio 2.56
  width 2.0 nm: exact 4.060e-18   WKB 1.586e-18   ratio 2.56

Alpha decay: WKB (Gamow) estimate versus measured half-life
  nucleus   Q (MeV)  turning point (fm)  2 gamma   log10 t(s) model  measured
  U-238     4.270         61.7           89.4        17.78        17.15
  Th-232    4.082         63.2           90.6        18.30        17.65
  Ra-226    4.871         51.8           74.1        11.10        10.70
  Po-210    5.407         44.5           62.5         6.03         7.08
  Po-214    7.833         30.7           40.0        -3.81        -3.78
  Po-212    8.954         26.9           33.2        -6.79        -6.52

With no adjustable parameters beyond the radius formula, the model reproduces measured half-lives spanning 25 orders of magnitude, from 0.3 microseconds for polonium-212 to 14 billion years for thorium-232, typically to within a factor of about 10. The exponent \(2\gamma\) ranges from 33 to 91; a change of 1 MeV in the alpha energy changes it by roughly 10, a factor of \(e^{10} \approx 20{,}000\) in half-life. This extreme sensitivity, a hallmark of tunneling, is the whole explanation of the Geiger–Nuttall law. The remaining discrepancies reflect the crude nuclear model and the probability that an alpha particle is preformed, which nuclear structure calculations supply.

Tunneling elsewhere

The same exponential governs many processes. In the scanning tunneling microscope of Gerd Binnig and Heinrich Rohrer (1981), the current between a sharp tip and a surface falls by roughly a factor of 10 for every 0.1 nm of extra gap, which lets the instrument resolve single atoms. In the Sun, protons fuse only by tunneling through their mutual Coulomb barrier, and the Gamow factor sets the temperature dependence of the fusion rate, the subject of Astrophysics, PHYS-430. Flash memory stores data by tunneling electrons through a thin oxide onto a floating gate. And the reverse process to alpha decay, a positive ion tunneling into a nucleus, sets the rates of the nuclear reactions in stars that build the elements. In each case, a WKB estimate of \(\gamma\) gives the dependence on energy and distance in a single line.