University Physics V · Radioactive Decay · 15.3

Alpha Decay as a WKB Barrier Integral

Alpha half-lives run from a fraction of a microsecond to longer than the age of the Earth, and nothing else in nuclear physics spans that. This topic shows you where the range comes from: one integral across a Coulomb barrier, evaluated numerically, with everything the integral cannot know pushed into a single fitted number out front.

Open the lab ↓Practice & worked examples ↓

Build the model

Connect the measurement to the mechanism.

An α emitter is not a stationary state, and no bound-state calculation will ever hand you a half-life. What produces one is a boundary condition: impose the Siegert condition, purely outgoing waves at large r, on the radial equation for α–daughter relative motion and the problem stops being Hermitian. Its eigenvalue moves off the real axis to E = E₀ − iΓ/2, the interior norm decays as exp(−Γt/ħ), and λ = Γ/ħ is read straight off the pole.

Solving that non-Hermitian problem exactly is hard, so Gamow factorises it: an α is assumed already assembled inside the well, striking the wall about 10²¹ times a second, each strike leaking with probability T = exp(−2∫κ dr), where κ = √(2μ(V − Q))/ħ between the classical turning points. Outside the strong-force radius the α sees only the daughter's Coulomb repulsion plus a centrifugal term, both known exactly, and Q comes from a mass table, so the exponent is a quadrature with no adjustable parameters. What it costs is the prefactor.

Whether four nucleons are really clustered at the surface, what the interior wavefunction looks like, where exactly R is — all of it collapses into a preformation factor P_α the model cannot compute and that varies by a factor of forty across emitters. The exponent is physics; the number in front of it is a fit.

Simple definition: The WKB barrier integral reduces α decay to one number: the exponent 2G that κ = √(2μ(V − Q))/ħ accumulates between the nuclear surface and the radius where the Coulomb barrier falls back to Q.
Example: For ²³⁸U → ²³⁴Th, Q = 4.270 MeV puts the outer turning point at b = 60.7 fm against a surface at R = 9.30 fm, giving 2G = 88.0 and a transmission of 10^-38.2 per strike.

Gamow state: the outgoing-wave boundary conditionu(r) → C e^(ikr) as r → ∞, E = E₀ − iΓ/2

Names what is being solved. λ = Γ/ħ is read off the pole, so a 1 keV width is a 6.6 × 10⁻¹⁹ s lifetime.

Γ is the width in MeV; ħ = 6.582 × 10⁻²² MeV s

The barrier beyond the strong-force radiusV(r) = 2Z_d e²/(4πε₀ r) + l(l+1)ħ²/(2μr²)

Everything outside R is known exactly, with no nuclear model in it. That is the only reason the estimate is worth making.

2Z_d e²/4πε₀ = 2Z_d × 1.440 MeV fm; μ = m_α M_d/(m_α + M_d) ≈ 3665 MeV/c²

WKB transmission per attemptT = exp(−2G), G = ∫_R^b κ(r) dr, κ = √(2μ(V(r) − Q))/ħ

One dimensionless number. Between ²¹²Po and ²³⁸U, 2G runs from 32.3 to 88.0 and carries the whole 10²⁴ spread.

κ in fm⁻¹ with ħc = 197.33 MeV fm; the outer limit b solves V(b) = Q

The s-wave integral in closed formG = (b√(2μQ)/ħ)[arccos√x − √(x(1−x))], x = R/b

Checks the quadrature to six figures, and is where the Geiger–Nuttall line comes from when x is small.

l = 0 and a pure Coulomb tail only; Q in MeV, b and R in fm

From transmission to decay constantλ = P_α f T, f = v/2R, v = c√(2Q/μc²), T_half = ln2/λ

For ²³⁸U: f = 7.8 × 10²⁰ s⁻¹, T = 6 × 10⁻³⁹, λ = 4.7 × 10⁻¹⁸ s⁻¹, so T_half = 4.6 Gyr against 4.468 measured.

f ≈ 10²¹ s⁻¹ and moves by only 1.5 across all emitters; P_α is fitted, never derived

The Geiger–Nuttall linelog₁₀ T_half = a Z_d/√Q + c, a = 1.47, c = −46.7

The small-x limit predicts slope 3.925/ln10 = 1.70; finite R softens it to 1.47. It puts ²²²Rn at 2.8 d against 3.82 d.

T_half in s, Q in MeV, Z_d the daughter's proton number; fitted to ²³⁸U and ²¹²Po alone

The boundary condition, not the potential, makes it decay

A bound state has a real eigenvalue and a norm that never changes, so no bound-state calculation can produce a half-life. What produces one is the boundary condition. Impose the Siegert condition on the radial equation for α–daughter relative motion — u(r) → C e^(ikr) as r → ∞, purely outgoing, nothing coming back in — and the operator stops being Hermitian on that domain. Its eigenvalues move off the real axis to E = E₀ − iΓ/2, the time factor e^(−iEt/ħ) becomes e^{(−iE₀t/ħ)} e^(−Γt/2ħ), and the probability of still finding the α inside falls as e^(−Γt/ħ). The decay constant is read straight off the pole: λ = Γ/ħ, with ħ = 6.582 × 10⁻²² MeV s, so a width of 1 keV is a lifetime of 6.6 × 10⁻¹⁹ s. Everything after this is a scheme for getting Γ without solving the non-Hermitian eigenvalue problem exactly.

Build the barrier out of quantities you can look up

Beyond the strong-force radius the potential contains no nuclear model, and that is what makes the estimate worth making. Q comes from a mass table as Q = [M(A, Z) − M(A−4, Z−2) − M(⁴He)]c²; because Z = (Z−2) + 2 the electron masses cancel, so atomic masses go in unaltered. For ²³⁸U that is 0.004584 u × 931.494 MeV/u = 4.270 MeV. The relative motion uses the reduced mass μ = m_α M_d/(m_α + M_d), which for any heavy daughter is 3665 MeV/c² — just under m_α c² = 3727 MeV, and near enough constant across the actinides that you can treat it as a fixed number. Outside R the α sees V(r) = 2Z_d e²/(4πε₀ r) + l(l+1)ħ²/(2μr²), with 2Z_d e²/4πε₀ = 2Z_d × 1.440 MeV fm. Take R = 1.2 fm × (A_d¹⁄³ + 4¹⁄³): for ²¹²Po → ²⁰⁸Pb that is 9.02 fm, where V(R) = 26.2 MeV — nearly three times the 8.95 MeV the α actually carries.

Do the integral, and remove the turning-point singularity

The limits are the nuclear surface R, where the potential is cut off, and the radius b at which the barrier drops back to the α's energy: V(b) = Q, so b = 2Z_d e²/(4πε₀ Q) when l = 0. Between them κ(r) = √(2μ(V(r) − Q))/ħ is real, and G = ∫_R^b κ dr. Work in MeV and fm with ħc = 197.33 MeV fm and κ comes out in fm⁻¹ directly. The integrand vanishes as √(b − r) at the upper limit, so its derivative is infinite there and plain Simpson quadrature loses order; substitute r = R + (b − R)(1 − w²) and the Jacobian 2w(b − R) cancels the square root, restoring clean convergence. For ²³⁸U with l = 0 the closed form checks it: b = 60.70 fm, x = R/b = 0.1532, prefactor b√(2μQ)/ħ = 54.42, bracket arccos√x − √(x(1−x)) = 0.8085, so G = 44.0, 2G = 88.0, and T = e^(−88.0) = 6 × 10^-39.

The prefactor is a frequency, and nearly a constant

T is a probability per attempt, not per second, so a rate has to be supplied. Gamow's estimate is an assault frequency f = v/2R, the reciprocal of the time an α of speed v takes to cross the well and return. Using the asymptotic speed v = c√(2Q/μc²) gives 0.048c for ²³⁸U and f = 7.8 × 10²⁰ s⁻¹, and then λ = P_α f T with T_half = ln2/λ. Now weigh the two factors against each other. Between ²¹²Po at Q = 8.95 MeV and ²³⁸U at Q = 4.27 MeV, f moves from 1.16 × 10²¹ to 7.8 × 10²⁰ s⁻¹ — a factor of 1.5. Over the same pair 2G moves from 32.3 to 88.0, which is 10^24.2 in transmission, and the measured half-lives differ by 10^23.7. Every decade of the spread sits in the exponent; the frequency is bookkeeping.

Geiger and Nuttall's line falls out of the small-x limit

Let x = R/b → 0. Then arccos√x − √(x(1−x)) → π/2 and 2G → 2πη, with η = 2Z_d e²/(4πε₀ ħv) the Sommerfeld parameter. In numbers 2G → 3.925 Z_d/√Q with Q in MeV, so log₁₀ T_half is linear in Z_d/√Q with slope 3.925/ln10 = 1.70. Real emitters have x between 0.15 and 0.34, not zero, and that finite-R correction pulls the bracket down to 0.8085 for ²³⁸U — cutting the exponent almost in half and flattening the slope with it. Fit the line to two measured points, ²³⁸U (Z_d = 90, Q = 4.270 MeV, 4.468 Gyr) and ²¹²Po (Z_d = 82, Q = 8.954 MeV, 0.299 μs), and you get log₁₀ T_half = 1.47 Z_d/√Q − 46.7 with T_half in seconds. That two-parameter line then puts ²²⁶Ra at 850 yr against a measured 1600 yr and ²²²Rn at 2.8 d against 3.82 d, across 24 decades of rate.

What the exponent leaves for the preformation factor

Set P_α = 1 and compare. ²³⁸U comes out at 4.6 Gyr against 4.468 Gyr and ²³²Th at 1.46 × 10¹⁰ yr against 1.405 × 10¹⁰ yr, which is far better than this crude R and assault frequency deserve; ²¹²Po lands 4.8 times too fast and ²¹⁰Po 41 times too fast. So the fitted P_α runs from about 1.0 down to 0.025 — a factor of forty, and that spread is precisely the many-body physics the barrier integral discarded. The cleanest reading is the two polonium isotopes, where R, Z_d and the barrier shape are nearly identical: ²¹²Po sits two protons and two neutrons outside doubly magic ²⁰⁸Pb, so its α is effectively pre-assembled, while ²¹⁰Po has N = 126 closed and must break that shell to make one. Their fitted preformation factors differ by 8.5. Nothing in κ(r) knows about shell closures, and nothing in it ever will.

Change one variable at a time

Make the relationship visible.

Interactive model

Decay energy Q4.3 MeV

Daughter proton number Z_d90

Nuclear surface R9.3 fm

Push Q from 4.3 to 5.3 MeV at Z_d = 90: b walks in from 60 to 49 fm, the shaded region loses a fifth of its width, and the half-life exponent drops seven decades. Then nudge R by half a femtometre and watch a surface nobody measured move the answer by one.

TURNING POINT b60.3 fm

BARRIER EXPONENT 2G87.4

log₁₀ T per strike-38.0

log₁₀ T_half (s)16.9

Live interpretationTURNING POINT b: 60.3 fm. BARRIER EXPONENT 2G: 87.4. log₁₀ T per strike: −38.0. log₁₀ T_half (s): 16.9

Catch the common trap

Explain before calculating.

Misconception check

The energy that goes into the barrier integral is Q, not the α's measured kinetic energy. The two differ by the daughter's recoil, E_α = Q M_d/(M_d + m_α), which for ²³⁸U means a detector reading 4.198 MeV where Q = 4.270 MeV — a gap of 1.7%. That looks negligible until you see what the exponent does with it: feeding 4.198 MeV through the same calculation gives 2G = 89.41 instead of 87.99 and a half-life of 1.93 × 10¹⁰ yr instead of 4.63 × 10⁹, wrong by a factor of 4.2 from an error smaller than most people's rounding. The reason Q is the right energy is structural. The barrier problem is the two-body relative motion written with the reduced mass μ, and the energy in that frame is the full release Q. E_α is what is left after the recoil has been subtracted; it belongs on a spectrum, not inside κ(r).

Alpha half-lives span about 10²⁴, from 0.299 μs for ²¹²Po to 4.468 Gyr for ²³⁸U. Working from λ = P_α f exp(−2G), which factor actually produces that range?

Choose an answer to test the model.

Practice & worked examples

Reason from the model, then test the result.

EasyFor ²¹²Po → ²⁰⁸Pb + α, the atomic masses are 211.988868 u, 207.976652 u and 4.002603 u. Find Q, the α kinetic energy a detector would record, the outer turning point b, and the height of the barrier at the nuclear surface. Use 2Z_d e²/4πε₀ = 2Z_d × 1.440 MeV fm and R = 1.2 fm × (A_d¹⁄³ + 4¹⁄³).

Q from the mass difference: the parent has 84 electrons and the products have 82 + 2 = 84, so the electron masses cancel and atomic masses go straight in. 211.988868 − 207.976652 − 4.002603 = 0.009613 u, and 0.009613 × 931.494 = 8.954 MeV.
Split Q by momentum conservation: the α takes E_α = Q M_d/(M_d + m_α) ≈ Q(A − 4)/A = 8.954 × 208/212 = 8.786 MeV, and the recoiling ²⁰⁸Pb carries the remaining 0.169 MeV. The barrier integral needs Q, not this 8.786 MeV.
Outer turning point, where the Coulomb barrier falls back to Q: b = 2Z_d × 1.440/Q = 2 × 82 × 1.440/8.954 = 236.2/8.954 = 26.4 fm.
Nuclear surface: R = 1.2 × (208¹⁄³ + 4¹⁄³) = 1.2 × (5.925 + 1.587) = 9.02 fm, so V(R) = 236.2/9.02 = 26.2 MeV.
Compare: the α holds 8.95 MeV and faces a barrier peaking at 26.2 MeV, 2.9 times its own energy, with 26.4 − 9.0 = 17.4 fm of classically forbidden ground between the two turning points.

AnswerQ = 8.954 MeV, E_α = 8.786 MeV, b = 26.4 fm, V(R) = 26.2 MeV = 2.9Q, and a forbidden region 17.4 fm wide.

MediumEstimate the half-life of ²³⁸U against α decay to the ground state of ²³⁴Th. Atomic masses: 238.050788 u, 234.043601 u, 4.002603 u. Take μc² = 3665 MeV, R = 9.30 fm, l = 0, and P_α = 1, and use the closed form G = (b√(2μQ)/ħ)[arccos√x − √(x(1−x))] with x = R/b.

Q = (238.050788 − 234.043601 − 4.002603) u × 931.494 MeV/u = 0.004584 × 931.494 = 4.270 MeV.
Turning point and shape parameter: b = 2 × 90 × 1.440/4.270 = 259.2/4.270 = 60.70 fm, so x = R/b = 9.30/60.70 = 0.1532.
Prefactor: √(2μQ) = √(2 × 3665 × 4.270) = √31299 = 176.9 MeV, so b√(2μQ)/ħ = 60.70 × 176.9/197.33 = 54.42, using ħc = 197.33 MeV fm.
Bracket: arccos√0.1532 − √(0.1532 × 0.8468) = arccos(0.3914) − 0.3602 = 1.1686 − 0.3602 = 0.8085. Hence G = 54.42 × 0.8085 = 44.0, 2G = 88.0, and T = e^(−88.0) = 6.1 × 10^-39.
Assault frequency: v/c = √(2Q/μc²) = √(8.540/3665) = 0.0483, so v = 1.447 × 10²² fm/s and f = v/2R = 1.447 × 10²²/18.6 = 7.78 × 10²⁰ s^-1.
Rate: λ = f T = 7.78 × 10²⁰ × 6.1 × 10⁻³⁹ = 4.7 × 10⁻¹⁸ s⁻¹, so T_half = ln2/λ = 1.46 × 10¹⁷ s = 4.6 × 10⁹ yr.

AnswerT_half ≈ 4.6 Gyr, within 4% of the measured 4.468 Gyr — closer than the model earns, since P_α was set to 1 and moving R from 9.30 to 9.70 fm alone would divide the answer by five.

Hard²³⁸U (0⁺) also decays to the first excited 2⁺ state of ²³⁴Th at 49.55 keV. Assign the orbital angular momentum, then estimate the intensity of that α₁ branch relative to the ground-state branch, taking the same R = 9.30 fm and μc² = 3665 MeV. The ground-state values are Q₀ = 4.2699 MeV and 2G₀ = 87.993.

Angular momentum and parity: Jᵢ = 0 and J_f = 2 force l = 2 exactly, and parity (−1)^l = +1 matches the 0⁺ → 2⁺ pair, so l = 2 is allowed and the centrifugal term is switched on.
Energy cost first: Q₁ = 4.2699 − 0.04955 = 4.2204 MeV. Redoing the l = 0 closed form at that Q gives b = 61.42 fm and 2G = 88.968, up 0.974 on the ground-state branch from the energy loss alone.
Centrifugal cost second: l(l+1)ħ²c²/(2μc²r²) = 6 × 197.33²/(2 × 3665 × 9.30²) = 0.369 MeV at r = R, falling to only 8 keV at r = b. It also pushes b out to 61.54 fm. Re-integrating numerically gives 2G = 89.486, a further 0.518.
Total penalty Δ(2G) = 89.486 − 87.993 = 1.493, so the transmission ratio is exp(−1.493) = 0.2249, and the assault frequencies differ only by √(Q₁/Q₀) = 0.9942.
Branch ratio: λ₁/λ₀ = 0.2249 × 0.9942 = 0.224, so the α₁ intensity is 0.224/(1 + 0.224) = 18.3% of all α emissions. The measured value is 20.9%.
Note what made this work: P_α, the assault-frequency convention and R all cancel in a ratio taken within one nuclide, which is why WKB predicts branch ratios far better than it predicts absolute half-lives.

Answerl = 2; Δ(2G) = 1.49, of which 0.97 is the 49.55 keV energy loss and 0.52 the centrifugal term; predicted α₁ intensity 18.3% against 20.9% measured.