Decay chains as a linear system

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, book every branch to some member, and every column of A sums to zero, so 1ᵀ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.

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 its norm is below a threshold of order one (about 5.4 for the degree-13 approximant), 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, stable only for steps shorter than 2/λₘₐₓ — Ra-226 into Rn-222 is a rate ratio of 1.5 × 10⁵, so crossing one radium half-life takes at least 5.3 × 10⁴ steps, where the propagator needs none.

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 = 1 and Δt = 3.3333 × 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) with the bracket computed as −expm1(−Δt), since 1 − e^(−Δt) typed as written cancels just as badly; the second factor then tends smoothly to 1 — or call expm and stop thinking about it.

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 from a fresh elution 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.

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)ᵀ. 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: from N₂(0) = 0, 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.

dN/dt = A N is the mean of a Markov jump process, not a trajectory. With 10¹⁵ nuclei the relative fluctuation is about 3 × 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.