Chapter 7

Decay Chains & Radioactive Series

Radioactive Decay Chains

When a radioactive nucleus decays, the daughter product is often itself unstable. This leads to a chain of successive decays until a stable nuclide is reached. The mathematical description of such chains was first solved by Harry Bateman in 1910 and forms the foundation for understanding natural radioactivity, nuclear waste management, and radiometric dating.

A general decay chain can be written as:

$$N_1 \xrightarrow{\lambda_1} N_2 \xrightarrow{\lambda_2} N_3 \xrightarrow{\lambda_3} \cdots \xrightarrow{\lambda_{n-1}} N_n \text{ (stable)}$$

where $\lambda_i$ is the decay constant for species $i$, related to its half-life by $\lambda_i = \ln 2 / t_{1/2,i}$.

The Bateman Equations

The coupled differential equations governing a decay chain are:

$$\frac{dN_1}{dt} = -\lambda_1 N_1$$

$$\frac{dN_i}{dt} = \lambda_{i-1} N_{i-1} - \lambda_i N_i \quad (i = 2, 3, \ldots, n-1)$$

$$\frac{dN_n}{dt} = \lambda_{n-1} N_{n-1}$$

For a two-member chain ($N_1 \to N_2$ stable), starting with $N_1(0) = N_1^0$ and $N_2(0) = 0$:

$$N_1(t) = N_1^0 \, e^{-\lambda_1 t}$$

$$N_2(t) = N_1^0 \frac{\lambda_1}{\lambda_2 - \lambda_1}\left(e^{-\lambda_1 t} - e^{-\lambda_2 t}\right)$$

The general Bateman solution for the $n$-th member, starting from pure parent, is:

$$N_n(t) = N_1^0 \left(\prod_{j=1}^{n-1} \lambda_j\right) \sum_{i=1}^{n} \frac{e^{-\lambda_i t}}{\prod_{\substack{j=1 \\ j \neq i}}^{n}(\lambda_j - \lambda_i)}$$

This elegant closed-form solution assumes all decay constants are distinct. When two decay constants are nearly equal, numerical methods are preferred to avoid cancellation errors.

Secular and Transient Equilibrium

Two important limiting cases arise depending on the relative half-lives of parent and daughter:

Secular Equilibrium ($\lambda_1 \ll \lambda_2$)

When the parent half-life is much longer than the daughter's ($t_{1/2,1} \gg t_{1/2,2}$), the parent activity remains essentially constant. After several daughter half-lives, the daughter activity equals the parent activity:

$$A_2 \approx A_1 \implies \lambda_2 N_2 \approx \lambda_1 N_1$$

Example: $^{226}$Ra (1600 yr) → $^{222}$Rn (3.82 d). After ~30 days, the radon activity equals the radium activity.

Transient Equilibrium ($\lambda_1 < \lambda_2$)

When the parent half-life is longer but not overwhelmingly so, the daughter activity initially grows, reaches a maximum, then tracks the parent decay with a constant ratio:

$$\frac{A_2}{A_1} = \frac{\lambda_2}{\lambda_2 - \lambda_1} \quad \text{(at equilibrium)}$$

The daughter activity maximum occurs at:

$$t_{\max} = \frac{\ln(\lambda_2/\lambda_1)}{\lambda_2 - \lambda_1}$$

If $\lambda_1 > \lambda_2$ (parent decays faster), no equilibrium is established. The daughter accumulates, reaches a maximum, and then decays with its own half-life.

The Four Natural Decay Series

Since alpha decay reduces the mass number by 4, all naturally occurring heavy nuclides fall into one of four families defined by $A \mod 4$:

Series	Form	Head	Stable End	$t_{1/2}$ (Head)
Thorium	$4n$	$^{232}$Th	$^{208}$Pb	1.41 × 10¹⁰ yr
Neptunium	$4n+1$	$^{237}$Np	$^{209}$Bi	2.14 × 10⁶ yr
Uranium	$4n+2$	$^{238}$U	$^{206}$Pb	4.47 × 10⁹ yr
Actinium	$4n+3$	$^{235}$U	$^{207}$Pb	7.04 × 10⁸ yr

The Neptunium series ($4n+1$) is extinct in nature because its longest-lived member ($^{237}$Np) has a half-life much shorter than the age of the Earth. The other three series are observed in natural uranium and thorium ores.

Branching Ratios

Many nuclides in these chains can decay by more than one mode (e.g., both alpha and beta decay). The branching ratio $b_i$ for mode $i$ is:

$$b_i = \frac{\lambda_i}{\lambda_{\text{total}}} = \frac{\lambda_i}{\sum_j \lambda_j}$$

For example, $^{212}$Bi decays by alpha emission (35.94%) to $^{208}$Tl, or by beta emission (64.06%) to $^{212}$Po.

Radioactive Dating Methods

Carbon-14 Dating

$^{14}$C is produced in the atmosphere by cosmic ray neutrons: $^{14}$N(n,p)$^{14}$C. Living organisms maintain equilibrium with atmospheric $^{14}$C/$^{12}$C ratio. After death,$^{14}$C decays with $t_{1/2} = 5730$ yr:

$$t = -\frac{t_{1/2}}{\ln 2}\ln\left(\frac{A(t)}{A_0}\right) = -\frac{1}{\lambda}\ln\left(\frac{A(t)}{A_0}\right)$$

Effective range: ~300 to ~50,000 years.

Uranium-Lead Dating

Two independent decay chains provide a consistency check. The age is determined from the ratios of lead isotopes to uranium:

$$\frac{^{206}\text{Pb}}{^{238}\text{U}} = e^{\lambda_{238} t} - 1 \qquad \frac{^{207}\text{Pb}}{^{235}\text{U}} = e^{\lambda_{235} t} - 1$$

These can be combined into a concordia diagram. Points lying on the concordia curve give concordant ages; discordant points indicate open-system behavior (lead loss).

Potassium-Argon Dating

$^{40}$K decays by electron capture to $^{40}$Ar (branching ratio 10.72%) and by $\beta^-$ to $^{40}$Ca (89.28%), with combined $t_{1/2} = 1.248 \times 10^9$ yr:

$$t = \frac{1}{\lambda_{\text{total}}}\ln\left(1 + \frac{\lambda_{\text{total}}}{\lambda_{\text{EC}}}\frac{^{40}\text{Ar}}{^{40}\text{K}}\right)$$

Applicable to rocks from ~100,000 years to billions of years old. The trapped argon in minerals provides the clock.

Activity and Units

The activity of a radioactive sample is defined as the number of disintegrations per unit time:

$$A = \lambda N = \frac{\ln 2}{t_{1/2}} N$$

Units: 1 Becquerel (Bq) = 1 disintegration/s; 1 Curie (Ci) = 3.7 × 10¹⁰ Bq. For a chain in secular equilibrium, the total activity of each member is equal:

$$\lambda_1 N_1 = \lambda_2 N_2 = \lambda_3 N_3 = \cdots$$

This means that in an old uranium ore where secular equilibrium has been established, all 14 members of the $^{238}$U chain have the same activity. The number of atoms of each member scales inversely with its decay constant.

Simulation: Bateman Decay Chain

This simulation solves the Bateman equations for a three-member decay chain A → B → C (stable), showing both transient equilibrium behavior and activity evolution.

Python

script.py96 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ──────────────────────────────────────────────
# Bateman Equations: Three-Member Decay Chain
# Parent -> Daughter -> Granddaughter
# ──────────────────────────────────────────────

# Decay constants (per year)
# Example: Ra-226 -> Rn-222 -> Po-218
lambda1 = np.log(2) / 1600.0      # Ra-226: t1/2 = 1600 yr
lambda2 = np.log(2) / (3.8235/365.25)  # Rn-222: t1/2 = 3.8235 days
lambda3 = np.log(2) / (3.1e-3/365.25/24/3600)  # Po-218: t1/2 = 3.1 min

# Initial number of parent atoms
N1_0 = 1.0e6
N2_0 = 0.0
N3_0 = 0.0

# Time array (in years)
t = np.linspace(0, 8000, 5000)

# Bateman solution for 3-member chain
N1 = N1_0 * np.exp(-lambda1 * t)

# Daughter (exact Bateman)
N2 = N1_0 * lambda1 / (lambda2 - lambda1) * (np.exp(-lambda1 * t) - np.exp(-lambda2 * t))

# Activities (A = lambda * N)
A1 = lambda1 * N1
A2 = lambda2 * N2

# Secular equilibrium ratio
print(f"Ra-226 half-life: 1600 years")
print(f"Rn-222 half-life: 3.82 days")
print(f"Ratio lambda2/lambda1: {lambda2/lambda1:.2e}")
print(f"Secular equilibrium reached when A_daughter ~ A_parent")

# For plotting, use a simpler example with comparable half-lives
# A -> B -> C (stable)
t1_A = 10.0  # hours
t1_B = 3.0   # hours
lam_A = np.log(2) / t1_A
lam_B = np.log(2) / t1_B

N_A0 = 1000.0
t_hr = np.linspace(0, 50, 1000)

N_A = N_A0 * np.exp(-lam_A * t_hr)
N_B = N_A0 * lam_A / (lam_B - lam_A) * (np.exp(-lam_A * t_hr) - np.exp(-lam_B * t_hr))
N_C = N_A0 - N_A - N_B

A_A = lam_A * N_A
A_B = lam_B * N_B

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0a0a1a')

for ax in axes:
    ax.set_facecolor('#0a0a1a')
    for spine in ax.spines.values():
        spine.set_edgecolor('#334155')

# Left: Number of atoms
axes[0].plot(t_hr, N_A, color='#ef4444', linewidth=2.5, label='N_A (parent)')
axes[0].plot(t_hr, N_B, color='#f59e0b', linewidth=2.5, label='N_B (daughter)')
axes[0].plot(t_hr, N_C, color='#22c55e', linewidth=2.5, label='N_C (stable)')
axes[0].set_xlabel('Time [hours]', color='#94a3b8', fontsize=12)
axes[0].set_ylabel('Number of Atoms', color='#94a3b8', fontsize=12)
axes[0].set_title('Decay Chain: A -> B -> C (stable)', color='white', fontsize=14)
axes[0].tick_params(colors='#94a3b8')
axes[0].grid(True, alpha=0.2, color='#334155')
axes[0].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')

# Right: Activities
axes[1].plot(t_hr, A_A, color='#ef4444', linewidth=2.5, label='Activity A')
axes[1].plot(t_hr, A_B, color='#f59e0b', linewidth=2.5, label='Activity B')
axes[1].set_xlabel('Time [hours]', color='#94a3b8', fontsize=12)
axes[1].set_ylabel('Activity [decays/hr]', color='#94a3b8', fontsize=12)
axes[1].set_title('Transient Equilibrium (t1/2_A > t1/2_B)', color='white', fontsize=14)
axes[1].tick_params(colors='#94a3b8')
axes[1].grid(True, alpha=0.2, color='#334155')

# Mark transient equilibrium time
t_max = np.log(lam_B / lam_A) / (lam_B - lam_A)
axes[1].axvline(x=t_max, color='#64748b', linestyle='--', alpha=0.6, label=f't_max = {t_max:.1f} hr')
axes[1].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')

plt.tight_layout()
plt.savefig('decay_chain.png', dpi=130, bbox_inches='tight',
            facecolor=fig.get_facecolor())
print(f"Time of max daughter activity: {t_max:.2f} hours")
print(f"Transient equilibrium ratio A_B/A_A at late times: {lam_B/(lam_B - lam_A):.3f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Simulation: Four-Member Chain

This Fortran program integrates a four-member decay chain using 4th-order Runge-Kutta, demonstrating numerical solution of the coupled differential equations and verification of atom conservation.

Fortran

program.f9075 lines

program decay_chain_bateman
  implicit none
  ! ──────────────────────────────────────────
  ! Bateman equations: numerical integration
  ! of a 4-member decay chain A->B->C->D(stable)
  ! using 4th-order Runge-Kutta
  ! ──────────────────────────────────────────

integer, parameter :: nsteps = 2000
  real(8) :: dt, t, t_max_sim
  real(8) :: N(4), k1(4), k2(4), k3(4), k4(4), Ntmp(4)
  real(8) :: lam(3)  ! decay constants
  real(8) :: t_half(3)
  integer :: i

! Half-lives in hours
  t_half(1) = 20.0d0   ! Parent A
  t_half(2) = 5.0d0    ! Daughter B
  t_half(3) = 1.0d0    ! Granddaughter C
  ! D is stable

do i = 1, 3
    lam(i) = log(2.0d0) / t_half(i)
  end do

! Initial conditions
  N(1) = 10000.0d0  ! Parent
  N(2) = 0.0d0      ! Daughter
  N(3) = 0.0d0      ! Granddaughter
  N(4) = 0.0d0      ! Stable end product

t_max_sim = 80.0d0  ! hours
  dt = t_max_sim / dble(nsteps)
  t = 0.0d0

write(*,'(A)') '# Bateman Decay Chain (Runge-Kutta Integration)'
  write(*,'(A)') '# t(hr)    N_A        N_B        N_C        N_D'
  write(*,'(F8.2, 4F12.2)') t, N(1), N(2), N(3), N(4)

do i = 1, nsteps
    ! RK4 step
    call derivs(N, lam, k1)
    Ntmp = N + 0.5d0 * dt * k1
    call derivs(Ntmp, lam, k2)
    Ntmp = N + 0.5d0 * dt * k2
    call derivs(Ntmp, lam, k3)
    Ntmp = N + dt * k3
    call derivs(Ntmp, lam, k4)

N = N + dt / 6.0d0 * (k1 + 2.0d0*k2 + 2.0d0*k3 + k4)
    t = t + dt

if (mod(i, 200) == 0) then
      write(*,'(F8.2, 4F12.2)') t, N(1), N(2), N(3), N(4)
    end if
  end do

write(*,*)
  write(*,'(A)') '# Final state:'
  write(*,'(A,F10.2)') '# Total atoms conserved: ', N(1)+N(2)+N(3)+N(4)
  write(*,'(A,F10.4)') '# Activity ratio A_B/A_A: ', (lam(2)*N(2))/(lam(1)*N(1)+1d-30)
  write(*,'(A,F10.4)') '# Activity ratio A_C/A_A: ', (lam(3)*N(3))/(lam(1)*N(1)+1d-30)

contains

subroutine derivs(Nin, lamin, dNdt)
    real(8), intent(in) :: Nin(4), lamin(3)
    real(8), intent(out) :: dNdt(4)
    dNdt(1) = -lamin(1) * Nin(1)
    dNdt(2) =  lamin(1) * Nin(1) - lamin(2) * Nin(2)
    dNdt(3) =  lamin(2) * Nin(2) - lamin(3) * Nin(3)
    dNdt(4) =  lamin(3) * Nin(3)
  end subroutine derivs

end program decay_chain_bateman

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Gamma Decay Nuclear Isomers →