University Physics V · Radioactive Decay · 15.6

Decay Chains as a Linear System

Stop deriving Bateman formulae one nuclide at a time. Assemble the chain as a single lower-triangular matrix, exponentiate it, and the parent, every daughter, both equilibria and the awkward equal-λ case fall out of one call — provided you read the matrix's spectrum before trusting its output.

Open the lab ↓Practice & worked examples ↓

Build the model

Connect the measurement to the mechanism.

A decay chain is exactly linear, and that is not an approximation: a nucleus decays at rate λᵢ whatever else is in the sample, so the inventory obeys dN/dt = A N with A constant and members coupled only by feeding terms. Order them by generation and A is lower triangular — −λᵢ on the diagonal, bᵢ λᵢ one row below — which hands you the spectrum for nothing: the eigenvalues are the decay constants themselves, no characteristic polynomial required. The solution is the matrix exponential N(t) = e^(At) N(0), one object that carries any starting vector to any time and composes as e^(A(t+s)) = e^(At) e^(As), so an elution becomes a projection between two propagations rather than a fresh derivation.

The Bateman formulae are not a rival method: they are the entries of e^(At) written out, and their 1/(λₖ − λⱼ) coefficients are the partial-fraction form of the same eigendecomposition V diag(e^{(−λᵢ t)}) V⁽⁻¹⁾. That is also where the cost sits. The diagonal form exists only while the λᵢ are distinct; as two approach, their eigenvectors turn parallel and cond(V) grows like 1/|λ₂ − λ₁|, and at equality A is defective, one eigenvector short of a basis, so V⁽⁻¹⁾ does not exist even though the physics stays finite and merely acquires a Jordan term t e^(−λt).

Scaling and squaring never forms V, which is the whole reason to prefer it.

Simple definition: A decay chain is the constant-coefficient linear system dN/dt = A N, with A carrying −λᵢ on the diagonal and the feeding term bᵢ λᵢ one row below, whose solution is the matrix exponential acting on the starting inventory, N(t) = e^(At) N(0).
Example: For Mo-99 (λ₁ = 0.010512 h⁽⁻¹⁾) feeding Tc-99m (λ₂ = 0.115396 h⁽⁻¹⁾) through b = 0.8757, A = [[−0.010512, 0], [0.009205, −0.115396]] h⁽⁻¹⁾; the matrix is lower triangular, so its eigenvalues are just −0.010512 and −0.115396 h⁽⁻¹⁾, read off the diagonal.

The chain matrixdN/dt = A NAᵢᵢ = −λᵢA_(i+1, i) = bᵢ λᵢ

Triangular means the eigenvalues are the diagonal itself — the chain's time constants need no characteristic polynomial.

Nᵢ in nuclei, λᵢ in s⁽⁻¹⁾, bᵢ the fraction of i's decays feeding member i+1; ordered by generation, A is lower triangular

Matrix exponential solutionN(t) = e^(At) N(0)e^(At) = Σ (At)ⁿ/n!

One propagator moves any starting vector to any time, and e^(A(t+s)) = e^(At) e^(As) lets an elution be a projection applied between two steps.

A matrix function, not element-by-element: numpy.exp(A) is a different object; use scipy.linalg.expm(A*t). At is dimensionless.

Diagonalised route and its conditioninge^(At) = V diag(e^{(−λᵢ t)}) V⁽⁻¹⁾

Says why the answer is a sum of pure exponentials, and warns exactly where the numerical route dies: eigenvectors turning parallel.

Columns of V are eigenvectors of A; valid only for distinct λᵢ. For two members cond(V) ≈ √((λ₂−λ₁)² + b²λ₁²)/|λ₂ − λ₁|.

Bateman, two membersN₂(t) = b λ₁ N₁(0) (e^(−λ₁t) − e^(−λ₂t))/(λ₂ − λ₁) + N₂(0) e^(−λ₂t)

The closed form to test code against — and its (λ₂ − λ₁) denominator is the defect warning printed in plain sight.

λ in s⁽⁻¹⁾, N in nuclei; requires λ₁ ≠ λ₂. These are the (2,1) and (2,2) entries of e^(At) written out.

Transient and secular equilibriumA₂/A₁ → b λ₂/(λ₂ − λ₁)→ b when λ₁ ≪ λ₂

Mo-99 into Tc-99m: λ₂/(λ₂ − λ₁) = 1.100, so the daughter settles at 0.8757 × 1.100 = 0.963 of the parent activity, within 5% after 29 h.

Activity Aᵢ = λᵢ Nᵢ in becquerel; the limit is reached to 5% once (λ₂ − λ₁)t ≥ 3, whatever N₂(0) was.

The defective limitλ₂ → λ₁ = λ: N₂(t) = b λ N₁(0) t e^(−λt)

Peaks at t = 1/λ with N₂ = b N₁(0)/e = 0.368 b N₁(0). expm returns it; the eigendecomposition returns noise.

The 0/0 limit of the Bateman form; at exact equality A has only one eigenvector for a doubled eigenvalue.

Assemble A, then read its spectrum for free

Number the members by generation, one row and one column each. The diagonal carries the total loss, Aᵢᵢ = −λᵢ, where λᵢ is the full decay constant including any branch that leaves the chain; the entry directly below it, A_(i+1, i) = bᵢ λᵢ, carries the fraction bᵢ of those decays that lands on the next member. Include the stable end member with λ = 0 and every column of A sums to zero, so 1^T A = 0: the all-ones row vector is a left null vector and total nuclei are conserved. That identity is the cheapest audit available — build the matrix, sum its columns, and a mistyped branch shows up at once. For Mo-99 into Tc-99m into Tc-99, the 12.4% of Mo-99 decays that go straight to the ground state put an entry in row 3, column 1, so A stops being bidiagonal; with the members still ordered by generation it stays lower triangular, and its eigenvalues are still exactly −λ_Mo, −λ_Tc99m and 0.

The exponential is the solution, expm is the algorithm

Define e^(At) by the series Σ (At)ⁿ/n!, which converges for every square A. Differentiate term by term and d/dt e^(At) = A e^(At), while e^(A⋅0) = I, so N(t) = e^(At) N(0) satisfies the initial-value problem — and, the system being linear with a bounded constant matrix, it is the only thing that does. For two members the propagator is [[e^(−λ₁t), 0], [b λ₁(e^(−λ₁t) − e^(−λ₂t))/(λ₂ − λ₁), e^(−λ₂t)]]. Do not evaluate the series directly: its terms grow before they shrink and cancellation ruins them. scipy.linalg.expm scales At down by 2^s until the norm is below about 1/2, evaluates a Padé approximant with a few matrix products and one well-conditioned solve, then squares s times. Two traps: numpy.exp(A) exponentiates element by element and is a different matrix, and forward Euler is throttled by the fastest member — Ra-226 into Rn-222 is a rate ratio of 1.5 × 10⁵, needing about 1.5 × 10⁶ steps to cross one radium half-life, where the propagator needs none.

Bateman is a matrix element, and its denominators warn you

The Bateman formulae are not an alternative to the matrix; they are the entries of e^(At) for a lower-triangular A, written as a partial-fraction expansion. Fed only at the top, Nₙ(t) = N₁(0)(Π bᵢ λᵢ) Σⱼ e^{(−λⱼ t)}/Π_(k≠j)(λₖ − λⱼ), and every denominator is a difference of decay constants — the same differences that appear when V is inverted. That structure is clean on paper and treacherous in floating point. Write the two-member numerator as e^(−λ₁t)(1 − e^(−Δt)) with Δ = λ₂ − λ₁: at Δt = 3.3 × 10⁻⁶ the two exponentials are 0.3678794412 and 0.3678782149, agreeing to better than four parts in a million, so the subtraction throws away about five and a half of the sixteen digits a double carries, and the 1/Δ prefactor then magnifies whatever noise is left. Near degeneracy, evaluate it as t e^(−λ₁t) (1 − e^(−Δt))/(Δt), whose second factor tends smoothly to 1 — or call expm and stop thinking about it.

Both equilibria are one solution read in two limits

Divide the two-member solution by the parent and the transient separates cleanly: A₂/A₁ = b λ₂/(λ₂ − λ₁) [1 − e^{(−(λ₂−λ₁)}t)] + (λ₂N₂(0)/λ₁N₁(0)) e^{(−(λ₂−λ₁)}t). Whatever the daughter started at, the memory of it decays as e^{(−(λ₂−λ₁)}t) and the ratio settles on b λ₂/(λ₂ − λ₁): transient equilibrium, in which the daughter's activity is a fixed multiple of the parent's and both then fall at the slower rate λ₁. Take λ₁ ≪ λ₂ and that multiple collapses to b — secular equilibrium, the daughter simply repeating what an effectively constant parent feeds it. For Mo-99 into Tc-99m the multiple is 0.8757 × 1.100 = 0.963, reached to within 5% after 3/(λ₂ − λ₁) = 29 h. Reverse the inequality, λ₁ greater than λ₂, and there is no equilibrium at all: the long-time decay follows the smallest decay constant in the chain, the dominant eigenvalue of A, whichever member owns it.

Equal decay constants break the basis, not the physics

Set λ₂ = λ₁ = λ in the two-member matrix: A = [[−λ, 0], [bλ, −λ]]. The eigenvalue −λ now has algebraic multiplicity 2, but A + λI = [[0, 0], [bλ, 0]] has rank 1, so its null space is one-dimensional and there is a single eigenvector, (0, 1)^T. A is defective — a Jordan block, not a diagonalisable matrix — and no V with e^(At) = V diag(e^(−λt)) V⁽⁻¹⁾ exists. The differential equation is untroubled: N₂(t) = b λ N₁(0) t e^(−λt), peaking at t = 1/λ with b N₁(0)/e = 0.368 b N₁(0), and t e^(−λt) is precisely the generalised-eigenvector term a Jordan block contributes. The numerical failure is quiet rather than loud. numpy.linalg.eig hands back two near-parallel columns instead of refusing, V⁽⁻¹⁾ amplifies rounding by roughly the reciprocal of the angle between them, and the output still looks like a curve. Never build a chain solver on eig.

What the constant matrix leaves out

dN/dt = A N is the mean of a Markov jump process, not a trajectory. With 10¹⁵ nuclei the relative fluctuation is about 10⁻⁷ and the distinction is academic; with ten atoms of a superheavy element the population is a staircase and the process itself has to be simulated. A is also constant, which asserts that nothing outside the chain touches it. Irradiation in a reactor adds a source, dN/dt = A N + s, solved by variation of constants or by exponentiating the augmented matrix [[A, s], [0, 0]]. A chemical separation is not a change to A at all but a projection P applied between propagations, N(t₂) = e^{(A(t₂−t₁)}) P N(t₁), which is exactly what eluting a generator does to the daughter. And the stiffness that makes expm attractive also offers a shortcut: once a fast member sits in secular equilibrium, its differential equation can be replaced by the algebraic constraint λᵢ Nᵢ = b_(i−1) λ_(i−1) N_(i−1), shrinking the system.

Change one variable at a time

Make the relationship visible.

Interactive model

Rate ratio λ₂/λ₁4.00

Feeding branch b0.90

Initial daughter N₂(0)/N₁(0)0.00

Drive λ₂/λ₁ towards 1: the equilibrium ratio and the eigenvector condition number both run away while the dashed set-up line slides off the right edge — A going defective. Then move the initial daughter: all curves merge onto one line, since the fast mode dies and the slow eigenvector lives.

EQUILIBRIUM RATIO A₂/A₁1.20

SET-UP TIME, (λ₂−λ₁)t = 31.00 1/λ₁

A₂/A₁ AT t = 2/λ₁1.20

EIGENVECTOR CONDITION1.0

Live interpretationEQUILIBRIUM RATIO A₂/A₁: 1.20. SET-UP TIME, (λ₂−λ₁)t = 3: 1.00 1/λ₁. A₂/A₁ AT t = 2/λ₁: 1.20. EIGENVECTOR CONDITION: 1.0

Catch the common trap

Explain before calculating.

Misconception check

The Bateman denominator λ₂ − λ₁ going to zero is not the physics blowing up. The denominator and the numerator e^(−λ₁t) − e^(−λ₂t) vanish together, and their ratio tends to the perfectly finite t e^(−λt): dN/dt = A N is linear with a bounded constant matrix, so it has exactly one solution for every choice of λ, degenerate or not. What fails at λ₂ = λ₁ is the representation — A is defective, one eigenvector short of a basis, so "a sum of pure exponentials" is no longer a complete description and V⁽⁻¹⁾ does not exist. The usual patch, nudging λ₂ by 10⁻⁹, is worse than it looks: it trades a clean division by zero for a V with condition number near 10⁹, and the code then returns a plausible-looking answer carrying nine digits of silent damage. Compute e^(At) directly instead.

A two-member chain has λ₁ = λ₂ = λ exactly, with N₂(0) = 0 and feeding branch b. What does the daughter population N₂(t) do?

Choose an answer to test the model.

Practice & worked examples

Reason from the model, then test the result.

EasyThe Pb-210 → Bi-210 → Po-210 → Pb-206 tail of the uranium series has half-lives 22.20 y, 5.012 d and 138.4 d, every branching ratio equal to 1, and Pb-206 stable. Write the matrix A of dN/dt = A N in units of d⁽⁻¹⁾, state its eigenvalues without forming a characteristic polynomial, and give one arithmetic check on the matrix itself.

Convert with λ = ln2/T_half, working in days: 22.20 y = 8108.6 d, so λ(Pb-210) = 0.693147/8108.6 = 8.548 × 10⁻⁵ d⁽⁻¹⁾; λ(Bi-210) = 0.693147/5.012 = 0.13830 d⁽⁻¹⁾; λ(Po-210) = 0.693147/138.4 = 5.008 × 10⁻³ d⁽⁻¹⁾; λ(Pb-206) = 0.
One row and column per member, ordered by generation. The diagonal is the loss −λᵢ, and because every branch is 1 the entry directly below is the full +λᵢ, so A is 4 × 4 lower bidiagonal with diagonal (−8.548 × 10⁻⁵, −0.13830, −5.008 × 10⁻³, 0) and subdiagonal (8.548 × 10⁻⁵, 0.13830, 5.008 × 10⁻³).
A is lower triangular, so det(A − μI) is the product of its diagonal entries minus μ: the eigenvalues are the diagonal itself, (0, −8.548 × 10⁻⁵, −0.13830, −5.008 × 10⁻³) d⁽⁻¹⁾. No polynomial is solved, and each eigenvalue is one member's decay constant.
Audit it: each column sums to zero (−λᵢ on the diagonal, +λᵢ in the feed below) and the stable member's column is all zeros, so 1^T A = 0 and total nuclei are conserved. The zero eigenvalue's right eigenvector is (0, 0, 0, 1)^T, so e^(At) N(0) drifts to pure Pb-206.

AnswerA is 4 × 4 lower bidiagonal, diagonal (−8.548 × 10⁻⁵, −0.13830, −5.008 × 10⁻³, 0) d⁽⁻¹⁾ with the same magnitudes on the subdiagonal; its eigenvalues are exactly those four numbers, and zero column sums confirm nuclei are conserved.

MediumA Mo-99/Tc-99m generator (half-lives 65.94 h and 6.0067 h, feeding branch b = 0.8757) is eluted at t = 0, leaving 40.0 GBq of Mo-99 and no Tc-99m. Find when the Tc-99m activity peaks and how large the peak is, give the ratio A₂/A₁ the chain settles to, and check the peak independently against the condition dN₂/dt = 0.

λ₁ = ln2/65.94 h = 0.010512 h⁽⁻¹⁾ and λ₂ = ln2/6.0067 h = 0.115396 h⁽⁻¹⁾, so λ₂ − λ₁ = 0.104884 h⁽⁻¹⁾ and λ₂/λ₁ = 10.978.
With N₂(0) = 0 the daughter follows e^(−λ₁t) − e^(−λ₂t); setting its derivative to zero gives λ₁e^(−λ₁t) = λ₂e^(−λ₂t), so tₘₐₓ = ln(λ₂/λ₁)/(λ₂ − λ₁) = 2.3959/0.104884 = 22.8 h.
A₂ = [b λ₂/(λ₂ − λ₁)] A₁(0) (e^(−λ₁t) − e^(−λ₂t)) = 0.9635 × 40.0 GBq × (0.78653 − 0.07165) = 0.9635 × 40.0 × 0.71488 = 27.6 GBq.
Independent check: dN₂/dt = 0 means feed equals loss, b λ₁N₁ = λ₂N₂, so at the peak A₂ = b A₁ exactly. A₁(22.8 h) = 40.0 e^(−0.2401) = 31.46 GBq, and 0.8757 × 31.46 = 27.55 GBq — the same number.
Long-time behaviour: A₂/A₁ = 0.9635 [1 − e^{(−0.104884 t)}], which tends to 0.963 and is within 1% of it after ln(100)/0.104884 = 43.9 h. This is transient equilibrium, not secular: λ₂/λ₁ is only 11, so the ratio settles at 0.963 rather than collapsing onto b = 0.8757.

Answertₘₐₓ = 22.8 h; peak Tc-99m activity 27.6 GBq, equal to 0.8757 × the Mo-99 activity there; the ratio A₂/A₁ then settles at 0.963, reached to 1% after about 44 h.

HardA two-member chain has b = 1, N₁(0) = 1.000 × 10⁶, N₂(0) = 0 and λ₁ = 0.3000 h⁽⁻¹⁾, with λ₂ = λ₁ + δ. At t = 10/3 h: (a) show the Bateman expression has a finite δ → 0 limit and evaluate it; (b) evaluate it at δ = 1.0 × 10⁻⁶ h⁽⁻¹⁾ and say what the subtraction costs in double precision; (c) state what the eigendecomposition does at δ = 0 and why expm does not.

Factor the difference: e^(−λ₁t) − e^(−λ₂t) = e^(−λ₁t)(1 − e^(−δt)), so N₂ = λ₁N₁(0) e^(−λ₁t) (1 − e^(−δt))/δ. Expanding, (1 − e^(−δt))/δ = t(1 − δt/2 + …), so δ cancels out of the denominator and the limit exists.
δ → 0 gives N₂ = λ₁N₁(0) t e^(−λ₁t). Here λ₁t = 0.3000 × 10/3 = 1.000, so N₂ = 1.000 × 10⁶ × e⁽⁻¹⁾ = 3.6788 × 10⁵ — and t = 1/λ₁ is exactly where that Jordan term peaks, at N₁(0)/e.
At δ = 1.0 × 10⁻⁶ h⁽⁻¹⁾: δt = 3.333 × 10⁻⁶, so the correction factor is 1 − δt/2 = 1 − 1.667 × 10⁻⁶ and N₂ = 367879.4 − 0.61 = 3.67879 × 10⁵. The physical answer has moved by 1.7 parts per million — nothing.
The arithmetic fares worse. e^(−λ₁t) = 0.3678794412 and e^(−λ₂t) = 0.3678782149 agree to 3.3 parts per million, so subtracting them discards about five and a half of the sixteen digits a double carries, and the 1/δ = 10⁶ prefactor magnifies what survives. Halve δ repeatedly and the surviving digits go with it.
At δ = 0, A = [[−λ, 0], [λ, −λ]] and A + λI = [[0, 0], [λ, 0]] has rank 1, so the doubled eigenvalue −λ owns a single eigenvector, (0, 1)^T. A is defective, V is singular, and numpy.linalg.eig returns two near-parallel columns whose inverse is meaningless.
scipy.linalg.expm forms no eigenvectors: it scales At down by 2^s, evaluates a Padé approximant using matrix products and one well-conditioned solve, then squares s times. It returns λ t e^(−λt) = 3.6788 × 10⁵ to full precision.

AnswerN₂ → λ₁N₁(0) t e^(−λ₁t) = 3.6788 × 10⁵ at t = 1/λ₁; at δ = 10⁻⁶ h⁽⁻¹⁾ it is 3.67879 × 10⁵, lower by 0.61 nuclei. At δ = 0 the matrix is defective and V is singular, while expm still returns the t e^(−λt) term exactly.