Decay

Author

Anthony Scopatz

The Bateman equations governing radioactive decay are an important subexpression of generalized transmutation equations. In many cases, it is desirable to compute decay on its own, outside of the presence of an neutron or photon field. In this case radioactive decay is a function solely on intrinsic physical parameters, namely half-lives. This document recasts the Bateman equations into a form that is better suited for computation than the traditional expression.

Canonical Bateman Equations for Decay

The canonical expression of the Bateman equations for a decay chain proceeding from a nuclide \(A\) to a nuclide \(Z\) at time \(t\) following a specific path is as follows 1:

\[N_C(t) = \frac{N_1(0)}{\lambda_C} \cdot \gamma \cdot \sum_{i=1}^C \lambda_i c_{i} e^{-\lambda_i t}\]

The symbols in this expression have the following meaning:

symbol	meaning
\(C\)	length of the decay chain
\(i\)	index for ith species, on range [1, C]
\(j\)	index for jth species, on range [1, C]
\(t\)	time [seconds]
\(N_i(t)\)	number density of the ith species at time t
\(t_{1/2,i}\)	half-life of the ith species
\(\lambda_i\)	decay constant of ith species, \(ln(2)/t_{1/2,i}\)
\(\gamma\)	The total branch ratio for this chain

Additionally, \(c_{i}\) is defined as:

\[c_i = \prod_{j=1,i\ne j}^C \frac{\lambda_j}{\lambda_j - \lambda_i}\]

Furthermore, the total chain branch ratio is defined as the product of the branch ratio between any two species 2:

\[\gamma = \prod_{i=i}^{C-1} \gamma_{i \to i+1}\]

Minor modifications are needed for terminal species: the first nuclide of a decay chain and the ending stable species. By setting \(C=1\), the Bateman equations can be reduced to simply:

\[N_C(t) = N_1(0) e^{-\lambda_1 t}\]

For stable species, the appropriate equation is derived by taking the limit of when the decay constant of the stable nuclide (\(\lambda_C\)) goes to zero. Also notice that every \(c_i\) contains exactly one \(\lambda_C\) in the numerator which cancels with the \(\lambda_C\) in the denominator in front of the summation:

\[ \begin{align}\begin{aligned}\lim_{\lambda_C \to 0} N_C(t) = N_1(0) \gamma \left[e^{-0t} + \sum_{i=1}^{C-1} \lambda_i \left(\frac{1}{0 - \lambda_i} \prod_{j=1,i\ne j}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i} \right) e^{-\lambda_i t} \right]\\N_C(t) = N_1(0) \gamma \left[1.0 - \sum_{i=1}^{C-1} \left(\prod_{j=1,i\ne j}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i} \right) e^{-\lambda_i t} \right]\end{aligned}\end{align} \]

Now, certain chains have intermeadiate nuclides that are almost stable. For example, decaying from Es-254 to Po-210 goes through U-238, which is very close to stable relative to all of the other nuclides in the chain. This can trigger floating point precision issues, where certain terms will underflow or overflow or generate NaNs. Obviously this is a situation to be avoided, if at all possible. To handle this sitiuation, let’s call \(p\) the index of the nuclide that is almost stable. We can then note that the Bateman equations can be reduced by the observation that \(\lambda_p \ll \lambda_{i\ne p}\) after we separate out the p-term from the summation:

\[ \begin{align}\begin{aligned}\frac{N_C(t)}{N_1(0)} = \frac{\gamma}{\lambda_C}\sum_{i\ne p}^C \left[\lambda_i \frac{\lambda_p}{\lambda_p - \lambda_i} \left(\prod_{j\ne i,p}^C \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right] + \frac{\gamma}{\lambda_C} \lambda_p \left(\prod_{j\ne p}^C \frac{\lambda_j}{\lambda_j - \lambda_p} \right) e^{-\lambda_p t}\\\frac{N_C(t)}{N_1(0)} = \frac{\gamma}{\lambda_C}\sum_{i\ne p}^C \left[\lambda_i \frac{\lambda_p}{\lambda_p - \lambda_i} \left(\prod_{j\ne i,p}^C \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right] + \frac{\gamma}{\lambda_C} \lambda_p \left(\prod_{j\ne p}^C \frac{\lambda_j}{\lambda_j - \lambda_p}\right) e^{-\lambda_p t}\\\frac{N_C(t)}{N_1(0)} = \frac{\gamma}{\lambda_C}\sum_{i\ne p}^C \left[\lambda_i \frac{\lambda_p}{- \lambda_i} \left(\prod_{j\ne i,p}^C \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right] + \frac{\gamma}{\lambda_C} \lambda_p \left(\prod_{j\ne p}^C \frac{\lambda_j}{\lambda_j}\right) e^{-\lambda_p t}\\\frac{N_C(t)}{N_1(0)} = \frac{-\gamma\lambda_p}{\lambda_C}\sum_{i\ne p}^C \left[ \left(\prod_{j\ne i,p}^C \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right] + \frac{\gamma\lambda_p}{\lambda_C} e^{-\lambda_p t}\end{aligned}\end{align} \]

The above expression for intermediate nuclides that are almost stable is valid when the last nuclide in the chain is unstable. When the last nuclide is stable, both the pth (almost stable nuclide) and the Cth (last and stable nuclide) must be removed can be split off from the summation and handled separately. As previously, then take \(\lambda_C \to 0\) and \(\lambda_p \ll \lambda_{i\ne p,C}\).

\[ \begin{align}\begin{aligned}\frac{N_C(t)}{N_1(0)} = \frac{\gamma}{\lambda_C}\sum_{i\ne p}^{C-1} \left[\lambda_i \frac{\lambda_C}{\lambda_C - \lambda_i} \frac{\lambda_p}{\lambda_p - \lambda_i} \left(\prod_{j\ne i,p}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right] + \frac{\gamma}{\lambda_C} \lambda_p \frac{\lambda_C}{\lambda_C - \lambda_p} \left(\prod_{j\ne p}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_p} \right) e^{-\lambda_p t} + \frac{\gamma}{\lambda_C} \lambda_C \frac{\lambda_p}{\lambda_p - \lambda_C} \left(\prod_{j\ne p}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_C} \right) e^{-\lambda_C t}\\\frac{N_C(t)}{N_1(0)} = \gamma\sum_{i\ne p}^{C-1} \left[\frac{\lambda_i \lambda_p}{(\lambda_C - \lambda_i)(\lambda_p - \lambda_i)} \left(\prod_{j\ne i,p}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right] + \frac{\gamma\lambda_p}{\lambda_C - \lambda_p} \left(\prod_{j\ne p}^{C-1} \frac{\lambda_j}{\lambda_j} \right) e^{-\lambda_p t} + \frac{\gamma\lambda_p}{\lambda_p - \lambda_C} \left(\prod_{j\ne p}^{C-1} \frac{\lambda_j}{\lambda_j} \right) e^{-\lambda_C t}\\\frac{N_C(t)}{N_1(0)} = -\gamma\sum_{i\ne p}^{C-1} \left[\left(\prod_{j\ne i}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right] + \frac{\gamma\lambda_p}{\lambda_C - \lambda_p} e^{-\lambda_p t} + \frac{\gamma\lambda_p}{\lambda_p - \lambda_C} e^{-\lambda_C t}\\\frac{N_C(t)}{N_1(0)} = -\gamma\sum_{i\ne p}^{C-1} \left[\left(\prod_{j\ne i}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right] + \frac{\gamma\lambda_p}{\lambda_C - \lambda_p} \left(e^{-\lambda_p t} - e^{-\lambda_C t}\right)\\\frac{N_C(t)}{N_1(0)} = -\gamma\sum_{i\ne p}^{C-1} \left[\left(\prod_{j\ne i}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right] -\gamma e^{-\lambda_p t} + \gamma\end{aligned}\end{align} \]

Lastly, we must handle the degenerate case where two nuclides in a chain have the same exact half-lives. This unfortunate situation arrises out of the fundemental nuclear data. Let’s call these the pth and qth species. To prevent underflow, overflow, and NaNs, we must separate these nuclides out of the summation and then take the limit as \(\lambda_q \to \lambda_p\).

\[ \begin{align}\begin{aligned}\frac{N_C(t)}{N_1(0)} = \frac{\gamma}{\lambda_C}\sum_{i\ne p,q}^{C} \left[\lambda_i \left(\prod_{j\ne i}^{C} \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right] + \frac{\gamma}{\lambda_C} \lambda_p \frac{\lambda_q}{\lambda_q - \lambda_p} \left(\prod_{j\ne p,q}^{C} \frac{\lambda_j}{\lambda_j - \lambda_p} \right) e^{-\lambda_p t} + \frac{\gamma}{\lambda_C} \lambda_q \frac{\lambda_p}{\lambda_p - \lambda_q} \left(\prod_{j\ne p,q}^{C} \frac{\lambda_j}{\lambda_j - \lambda_q} \right) e^{-\lambda_q t}\\\frac{N_C(t)}{N_1(0)} = \frac{\gamma}{\lambda_C}\sum_{i\ne p,q}^{C} \left[\lambda_i \left(\prod_{j\ne i}^{C} \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right] + \frac{\gamma\lambda_p^2}{\lambda_C} \left(\prod_{j\ne p,q}^{C} \frac{\lambda_j}{\lambda_j - \lambda_p} \right) \lim_{\lambda_q\to\lambda_p}\frac{e^{-\lambda_p t} - e^{-\lambda_q t}}{\lambda_q - \lambda_p}\\\frac{N_C(t)}{N_1(0)} = \frac{\gamma}{\lambda_C}\sum_{i\ne p,q}^{C} \left[\lambda_i \left(\prod_{j\ne i}^{C} \frac{\lambda_j}{\lambda_j - \lambda_i}\right) e^{-\lambda_i t}\right] + \frac{\gamma\lambda_p^2}{\lambda_C} \left(\prod_{j\ne p,q}^{C} \frac{\lambda_j}{\lambda_j - \lambda_p} \right) t e^{-\lambda_p t}\end{aligned}\end{align} \]

Binary Reformulation of Bateman Equations

There are two main strategies can be used to construct a version of these equations that is better suited to computation, if not clarity.

First, lets aim for minimizing the number of operations that must be performed to achieve the same result. This can be done by grouping constants together and pre-calculating them. This saves the computer from having to perform the same operations at run time. It is possible to express the Bateman equations as a simple sum of exponentials

\[N_C(t) = N_1(0) \sum_{i=1}^C k_{i} e^{-\lambda_i t}\]

where the coefficients \(k_i\) are defined as:

Single Nuclide in Chain:

\[k_i = 1\]

Last Nuclide Unstable:

\[k_i = \frac{\gamma}{\lambda_C} \lambda_i \prod_{j\ne i}^C \frac{\lambda_j}{\lambda_j - \lambda_i}\]

Last Nuclide Stable:

\[ \begin{align}\begin{aligned}k_{i\ne C} = -\gamma \prod_{j=1,i\ne j}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i}\\k_C = \gamma\end{aligned}\end{align} \]

Last Nuclide Unstable and pth Almost Stable:

\[ \begin{align}\begin{aligned}k_{i\ne p} = -\frac{\gamma\lambda_p}{\lambda_C} \prod_{j\ne i,p}^C \frac{\lambda_j}{\lambda_j - \lambda_i}\\k_p = \frac{\gamma\lambda_p}{\lambda_C}\end{aligned}\end{align} \]

Last Nuclide Stable and pth Almost Stable:

\[ \begin{align}\begin{aligned}k_{i\ne p,C} = -\gamma \prod_{j\ne i}^{C-1} \frac{\lambda_j}{\lambda_j - \lambda_i}\\k_p = -\gamma\\k_C = \gamma\end{aligned}\end{align} \]

Half-life Degeneracy Between pth and qth:

\[ \begin{align}\begin{aligned}k_i = \frac{\gamma}{\lambda_C} \lambda_i \prod_{j\ne i}^C \frac{\lambda_j}{\lambda_j - \lambda_i}\\k_p = \frac{\gamma\lambda_p^2}{\lambda_C} t \prod_{j\ne p,q}^C \frac{\lambda_j}{\lambda_j - \lambda_p}\\k_q = 0\end{aligned}\end{align} \]

If \(k_i\) are computed at run time then the this expression results in much more computational effort that than the original Bateman equations since \(\gamma/\lambda_C\) are brought into the summation. However, when \(k_i\) are pre-caluclated, many floating point operations are saved by avoiding explicitly computing \(c_i\).

The second strategy is to note that computers are much better at dealing with powers of 2 then then any other base, even \(e\). Thus the exp2(x) function, or \(2^x\), is faster than the natural exponential function exp(x), \(e^x\). As proof of this the following are some simple timing results:

In [1]: import numpy as np

In [2]: r = np.random.random(1000) / np.random.random(1000)

In [3]: %timeit np.exp(r)
10000 loops, best of 3: 26.6 µs per loop

In [4]: %timeit np.exp2(r)
10000 loops, best of 3: 20.1 µs per loop

This is a savings of about 25%. Since the core of the Bateman equations are exponentials, it is worthwhile to squeeze this algorithm as much as possible. Luckily, the decay constant provides an intrinsic mechanism to convert to base-2:

\[ \begin{align}\begin{aligned}N_C(t) = N_1(0) \sum_{i=1}^C k_{i} e^{-\lambda_i t}\\N_C(t) = N_1(0) \sum_{i=1}^C k_{i} e^{\frac{-\ln(2)\cdot t}{t_{1/2,i}}}\\N_C(t) = N_1(0) \sum_{i=1}^C k_{i} 2^{\frac{-t}{t_{1/2,i}}}\end{aligned}\end{align} \]

This expression can be further collapsed by defining \(a\) to be the precomputed exponent values:

\[a_i = \frac{-1}{t_{1/2,i}}\]

Thus, the final form of the binary representation of the Bateman equations are as follows:

General Formulation:

\[N_C(t) = N_1(0) \sum_{i=1}^C k_{i} 2^{a_i t}\]

where the \(k_i\) are as listed above. However, for practical purposes, it is better to compute the \(k_i\) from half-lives rather than decay constants. This is because they provide less floating point error, fewer oppurtunities to underflow or overflow to NaN or infinity, and a better mechanism for detecting stability. Thus, alternatively, the \(k_i\) are computed as:

Single Nuclide in Chain:

\[k_i = 1\]

Last Nuclide Unstable:

\[k_i = \gamma t_{1/2,i}^{C-2} t_{1/2,C} \prod_{j\ne i}^{C} \frac{1}{t_{1/2,i} - t_{1/2,j}}\]

Last Nuclide Stable:

\[ \begin{align}\begin{aligned}k_i = -\gamma t_{1/2,i}^{C-2} \prod_{j\ne i}^{C-1} \frac{1}{t_{1/2,i} - t_{1/2,j}}\\k_C = \gamma\end{aligned}\end{align} \]

Last Nuclide Unstable and pth Almost Stable:

\[ \begin{align}\begin{aligned}k_{i\ne p} = -\frac{\gamma t_{1/2,C}}{t_{1/2,p}} t_{1/2,i}^{C-2} \prod_{j\ne i,p}^C \frac{1}{t_{1/2,i} - t_{1/2,j}}\\k_p = \frac{\gamma t_{1/2,C}}{t_{1/2,p}}\end{aligned}\end{align} \]

Last Nuclide Stable and pth Almost Stable:

\[ \begin{align}\begin{aligned}k_{i\ne p,C} = -\gamma t_{1/2,i}^{C-2} \prod_{j\ne i}^{C-1} \frac{1}{t_{1/2,i} - t_{1/2,j}}\\k_p = -\gamma\\k_C = \gamma\end{aligned}\end{align} \]

Half-life Degeneracy Between pth and qth:

\[ \begin{align}\begin{aligned}k_i = \gamma t_{1/2,i}^{C-2} t_{1/2,C} \prod_{j\ne i}^{C} \frac{1}{t_{1/2,i} - t_{1/2,j}}\\k_p = \gamma\ln(2) t_{1/2,p}^{C-4} t_{1/2,C} t \prod_{j\ne p,q}^C \frac{1}{t_{1/2,p} - t_{1/2,j}}\\k_q = 0\end{aligned}\end{align} \]

With completely precomputed \(k\), \(a\), and the exp2() function, this formulation minimizes the number of floating point operations while completely preserving physics. No assumptions were made aside from the Bateman equations themselves in this proof.

Note that it is not possible to reduce the number of operations further. This is because \(k\) and \(a\) cannot be combined without adding further operations.

Implementation Specific Approximations

The above formulation holds generally for any decay chain. However, certain approximations are used in practice to reduce the number of chains and terms that are calculated.

1. Decay chains coming from spontaneous fission are only optionally tallied as they lead to an explosion of the total number of chains while contributing to extraordinarily rare branches.

2. Decay alphas are not treated as He-4 production.

3. The \(k_i\) and \(a_i\) are filtered to reject terms where \(|k_i| / \max(|k_i|) < 10^{-16}\). This filtering prevents excessive calculation from species which do not significantly contribute to end atom fraction. The threshold \(10^{-16}\) was chosen as because it is a reasonable naive estimate of floating point error after many operations. Note that we may filter only on the \(k_i\) because \(2^{a_i t} \le 1\). That is, the exponentional component can only reduce the magnitude of a term, not increase it.

In principle, each of these statements is reasonable. However, they may preclude desired behavior by users. In such a situation, these assumptions should be revisited.

Additional Information

For further discussion, please see:

• the mailing list post,

• the pull request, #614, and

• the benchmark study.

Note that the benchmark study shows quite high agreement between this method and ORIGEN v2.2.

References

1

Jerzy Cetnar, General solution of Bateman equations for nuclear transmutations, Annals of Nuclear Energy, Volume 33, Issue 7, May 2006, Pages 640-645, http://dx.doi.org/10.1016/j.anucene.2006.02.004.

2

Logan J. Harr. Precise Calculation of Complex Radioactive Decay Chains. M.Sc thesis Air Force Institute of Technology. 2007. http://www.dtic.mil/dtic/tr/fulltext/u2/a469273.pdf