Table of Contents
Radioactive Decay
Radioactive decay is a random process. For one unstable atom, it is completely random when it decays. For a large quantity of atoms, though, radioactive decay follows certain rules. The most important rule is that a radioactive sample, a collection of atoms of the same type, decay with a constant half-life. After one half-life, half of the original amount of atoms will have decayed. We can write this mathematically as $$ N(t)=N(0)\left(\frac{1}{2}\right)^{\frac{t}{t_{1/2}}} $$ where $N(t)$ is the number of atoms at time $t$, and $t_{1/2}$ is the half-life. To prove to ourselves that this works, we can plug in $t=0$, which yields $$ N(0)=N(0)\left(\frac{1}{2}\right)^0=N(0) $$ and plugging in $t=t_{1/2}$ gives $$ N(t_{1/2})=N(0)\left(\frac{1}{2}\right)^1=\frac{N(0)}{2} $$
We can rearrange this equation to make our lives easier. We start by writing that $$ N(t)=N(0)\left(\frac{1}{2}\right)^{\frac{t}{t_{1/2}}}=N(0)2^{-\frac{t}{t_{1/2}}} $$ By definition, $\exp(\log(x))=x$. We use this to say that $$ \begin{align} 2^{-\frac{t}{t_{1/2}}}&= \exp\left(\log\left(2^{-\frac{t}{t_{1/2}}}\right)\right) \\ &=\exp\left(-\frac{t}{t_{1/2}}\log(2)\right) \\ &=\exp\left(-\lambda t\right) \end{align} $$ where we define $\lambda=\frac{\log(2)}{t_{1/2}}$. Plugging this into the decay equation gives an equivalent but more ubiquitous form: $$ N(t)=N(0)e^{-\lambda t} $$ $\lambda$ is called the decay constant, and represents a rate. A more natural quantity to think about, however, is the number of atoms that decay per unit of time in a sample, called the activity. The SI units for activity are the Becquerel (Bq), where 1 Bq is 1 decay per second. The number of atoms that decay per unit of time in a sample is the time rate-of-change of the number of atoms in a sample, $\frac{dN}{dt}$. We can take this derivative to get $$ \begin{align} \frac{dN(t)}{dt}&=N(0)\frac{d}{dt}\left(e^{-\lambda t}\right) \\ &=N(0)\left(-\lambda e^{-\lambda t}\right) \\ &=-\lambda\left(N(0)e^{-\lambda t}\right) \\ &=-\lambda N(t) \end{align} $$ Denoting the activity of the sample as $A$, we write $$ A(t)=-\lambda N(t) $$
Decay Chains
When atoms decay, they emit some of their mass as radiation. What remains is a different atom, either a new element or a different isotope of the old one. This new isotope may itself be unstable and decay with a different half-life into yet another atom. Depending on what isotope you start from, you can end up with long decay chains. One such decay chain is called the “thorium series”, which starts at $^{232}$Th and ends on stable $^{208}$Pb. The thorium series is shown below.
Th232
└── Ra228
└── Ac228
└── Th228
└── Ra224
└── Rn220
└── Po216
└── Pb212
└── Bi212
├── Po212
│ └── Pb208
└── Tl208
└── Pb208
Each isotope in the chain has its own half-life (and therefore its own decay constant). The decay of $^{232}$Th results in the production of $^{228}$Ra, which itself decays, producing $^{228}$Ac, and so on.
Looking at these chains, we can ask the question, “how much of each isotope exists at any given time?” We can formulate an answer to this question using the equations we derived above. To simplify notation, we denote $^{232}$Th as isotope 1 and $^{228}$Ra as isotope 2. Let’s say we start with $N_1(0)$ atoms of $^{232}$Th, and $N_2(0)=0$ atoms (no atoms) of $^{228}$Ra. We can immediately write the time rate-of-change of the thorium atoms: $$ \frac{dN_1(t)}{dt}=-\lambda_1N_1(t) $$ If we are losing thorium atoms at a rate of $\lambda_1N_1(t)$, then we must be gaining radium atoms at a rate of $\lambda_1N_1(t)$. For radium then, we have $$ \frac{dN_2(t)}{dt}=\lambda_1N_1(t) $$ In addition to gaining atoms from the decay of the thorium, we are also losing the radium atoms through its own decay at a rate of $\lambda_2N_2(t)$. We therefore need to subtract out this term to get $$ \frac{dN_2(t)}{dt}=\lambda_1N_1(t)-\lambda_2N_2(t) $$ We therefore get a system of two coupled differential equations to describe the number of atoms of thorium and of radium at any point in time. We can expand our analysis to the next isotope in the chain, actinium. The number of actinium atoms depends on the rate of decay of radium and the rate of decay of itself: $$ \frac{dN_3(t)}{dt}=\lambda_2N_2(t)-\lambda_3N_3(t) $$ If we keep up this numbering scheme, we can easily write an equation for each isotope in a decay chain that is $n$ isotopes long. The last element in the decay chain (lead-212 in our example) is stable, so it does not have the loss term from its own decay, only the gain term from the decay of the isotope above it in the chain. $$ \frac{dN_n(t)}{dt}=\lambda_{n-1}N_{n-1}(t) $$
Solving the Differential Equations
Numerical Integration
The goal of this method is to turn the scary derivatives into sums without doing any difficult math. We start by defining a time-step, $\delta t$, that is small compared to the half-lives of all the isotopes in the decay chain. If we really squint, the decay equations $$ \begin{align} \frac{dN_1(t)}{dt}&=-\lambda_1N_1(t) & \frac{dN_2(t)}{dt}&=\lambda_1N_1(t)-\lambda_2N_2(t) \end{align} $$ look like $$ \begin{align} \frac{\delta N_1(t)}{\delta t}&=-\lambda_1N_1(t) & \frac{\delta N_2(t)}{\delta t}&=\lambda_1N_1(t)-\lambda_2N_2(t) \end{align} $$ This does actually work. The derivative is the slope of the curve at a point. We can approximate the slope of the curve at a point by the average slope of the curve around the point. This is an approximation, but as long as the interval ($\delta t$) is short compared to how quickly the slope of the curve is changing (the half-life of the isotope), then it is not a bad one. We can multiply both sides of the equations by $\delta t$ to get $$ \begin{align} \delta N_1(t)&=-\lambda_1N_1(t)\delta t & \delta N_2(t)&=(\lambda_1N_1(t)-\lambda_2N_2(t))\delta t \end{align} $$ If we keep track of amount of atoms for each isotope at each time-step, we simply add $\delta N$ to the amount of atoms at the current time-step to get the number of atoms at the next time-step. In order to calculate this $\delta N$, we need the number of atoms of the isotope of interest and the number of atoms of the isotope one position above it in the decay chain at the current time-step. Starting at $t=0$ where we know the initial number of atoms of all isotopes, we can calculate in time-steps going forward in time the number of atoms for each isotope in the chain.
This method is very easy to implement in a computer program, and is often accurate enough. Care does have to be taken, however, that the time-step is short enough. If there are very short-lived isotopes in the decay chain, the required time-step may be so small that the method is not computationally feasible.
We can write a simple program in Python which calculates the decay of $^{232}$Th into $^{228}$Ra. First, we define the half-lives of the isotopes. The half-lives are given in years. We will also need the numpy library, which we import.
import numpy as np
half_life_Th = 1.405e10
half_life_Ra = 5.73
Next, we calculate the decay constants from the half-lives.
lambda_Th = np.log(2) / half_life_Th
lambda_Ra = np.log(2) / half_life_Ra
Since both the half-lives are multiple years, we can set our time-step as half of a year and still get reasonable results. We also need to set the initial number of atoms of thorium and radium. We put the initial numbers as the first elements in arrays, which we will append the new number of atoms to at each time-step.
dt = 0.5
N_Th = [1e15]
N_Ra = [0.]
Finally, we make a loop that calculates the number of atoms of each isotope every year for 20 years. In the loop, we first calculate the change in the number of atoms of each isotope from decay over the year. Next, we subtract the amount of atoms lost from the previous value in the array, and for radium, we add the number of atoms lost from thorium.
for i in range(20 * 2):
dN_Th = lambda_Th * N_Th[i] * dt
dN_Ra = lambda_Ra * N_Ra[i] * dt
N_Th.append(N_Th[i] - dN_Th)
N_Ra.append(N_Ra[i] + dN_Th - dN_Ra)
This leaves us with two arrays containing the number of atoms of thorium and radium at half year increments. We can plot the data using the following code:
import matplotlib.pyplot as plt
times = np.linspace(0, 20, 41)
plt.plot(times, N_Th, label="Thorium Atoms")
plt.plot(times, N_Ra, label="Radium Atoms")
plt.xlabel("Time [years]")
plt.ylabel("Number of Atoms")
plt.yscale("log")
plt.legend()
plt.show()
The Bateman Equation
The Bateman Equation is an exact solution to the system of coupled differential equations. It is given by $$ N_i(t)=N_1(0)\left(\prod\limits^{i-1}_{k=1}\lambda_k\right) \sum\limits^i_{k=1}\frac{e^{-\lambda_kt}} {\prod\limits^i_{j=1,j\neq k}(\lambda_j-\lambda_k)} $$ The Bateman Equation solves the decay chain exactly, and is fairly easy to implement in a computer program. It can cause issues, however, if two isotopes in the chain have decay constants $\lambda_j\approx\lambda_k$. Floating-point error can cause the difference to go to zero causing a divide by zero.
The equation as written does not allow for isotopes down the chain to have non-zero starting numbers. Only the first isotope’s starting number is reflected in the solution. To treat this case, the chain can be split into components where each component is a chain with only one of the isotopes having a starting number of atoms. The component chains can be solved individually, and the results superimposed on one another. For example, given the chain
A ───► B ───► C
we can solve three separate chains
A ───► 0 ───► 0
0 ───► B ───► 0
0 ───► 0 ───► C
and sum the results to get $C(t)$.
Branching
If the decay chain has branching, a similar strategy can be employed. Consider the following chain:
A
├───── B ───┐
│ 25% │
│ │
└───── C ───┴─ D
75%
The goal is to solve for the number of atoms of isotope D at some time $t$. In order to do this, all possible paths to the isotope D must be found. For this simple chain, this task is fairly simple. For more complicated chains, it is beneficial to understand the structure of the chain. In a decay chain, there are no cycles. Following the path of decay, there will not be any closed loops (incidentally, this structure is known as a directed acyclic graph). This property allows us to employ a depth-first search algorithm to traverse the chain, finding all paths to the nuclide of interest.
For the chain given above, the two paths are
A ───► B ───► D
A ───► C ───► D
Each of these paths can be treated individually with the Bateman equation. This yields two solutions for the number of atoms of isotope D, D$_1$ and D$_2$ for the solution of the first and second paths, respectively. In this case, it is evident that 25% of the number of atoms of D come from the first path, and 75% come from the second path. The solution for the number of atoms of D is then $$ \mathrm{D=0.25D_1 + 0.75D_2} $$
For more complicated chains, with multiple non-trivial paths to the nuclide of interest, the contribution of each path can be represented as a path weight. The weight of a path is found as follows. Assume a path $p$ comprising nuclides $N_1$ to $N_n$. The branching ratio from $N_1$ to $N_2$ is represented as $b_{1\rightarrow2}$. The weight $w_p$ of the path $p$ is $$ w_p=\prod\limits_{i=2}^{n}b_{(i-1)\rightarrow i} $$ Assuming the solution to the Bateman equation for nuclide $n$ on the path $p$ is $N_{n,p}$, the total number of atoms of nuclide $n$ is $$ N_n=\sum\limits_{p}w_pN_{n,p} $$