radioactivedecay: A Python package for radioactive decay calculations

radioactivedecay is a Python package for radioactive decay modelling. It contains functions to fetch decay data, define inventories of nuclides and perform decay calculations. The default nuclear decay dataset supplied with radioactivedecay is based on ICRP Publication 107, which covers 1252 radioisotopes of 97 elements. The code calculates an analytical solution to a matrix form of the decay chain differential equations using double or higher precision numerical operations. There are visualization functions for drawing decay chain diagrams and plotting activity decay curves.

Summary radioactivedecay is a Python package for radioactive decay modelling. It contains functions to fetch decay data, define inventories of nuclides and perform decay calculations. The default nuclear decay dataset supplied with radioactivedecay is based on ICRP Publication 107, which covers 1252 radioisotopes of 97 elements. The code calculates an analytical solution to a matrix form of the decay chain differential equations using double or higher precision numerical operations. There are visualization functions for drawing decay chain diagrams and plotting activity decay curves.

Statement of Need
Calculations for the decay of radioactivity and the ingrowth of progeny underpin the use of radioisotopes in a wide range of research and industrial fields, spanning from nuclear engineering, medical physics, radiation protection, environmental science and archaeology to non-destructive testing, mineral prospecting, food preservation, homeland security and defence. radioactivedecay is an open source, cross-platform package for decay calculations and visualization. It supports decay chains with branching decays and metastable nuclear isomers. It includes a high numerical precision decay calculation mode, which resolves numerical problems with using double-precision floating-point numbers to calculate decay chains involving radionuclides with disparate half-lives (Bakin et al., 2018). This set of features distinguishes radioactivedecay from other commonly-used decay packages, such as Radiological Toolbox (Hertel et al., 2015) and PyNE (Scopatz et al., 2012). Radiological Toolbox is a closed-source Windows application, so it is not easily scriptable and its use of double-precision arithmetic makes it susceptible to numerical round-off errors. PyNE uses approximations to help mitigate numerical issues, however these may potentially affect accuracy. Moreover as of v0.7.5, PyNE does not correctly model metastable nuclear isomers within decay chains, which means, for example, it cannot simulate the production of 99m Tc from 99 Mo for medical imaging applications.
Theory and Implementation radioactivedecay implements the solution to the decay differential equations outlined by Amaku et al. (2010). If vector N contains the number of atoms of each radionuclide in a system, its elements N i can be ordered such that no progeny (either first or subsequent generation) of radionuclide i has itself an index lower than i. Ordering in this manner is possible because natural radioactive decay processes do not increase the mass number of the decaying radionuclide, and there are no cyclic decay chains where radionuclide i can decay to other radionuclides then reform itself (Ladshaw et al., 2020). Note metastable nuclear isomers have distinct indices from their ground states in N.
The radioactive decay chain differential equations expressed in matrix form are: Λ is a lower triangular matrix with elements: (2) λ j is the decay constant of radionuclide j, and b ji is the branching fraction from radionuclide j to i. Λ is diagonalizable so its eigendecomposition can be used to rewrite Equation 1 as: Λ d is a diagonal matrix whose elements are the negative decay constants, i.e. Λ dii = −λ i . Matrix C and its inverse C −1 are both lower triangular matrices that are calculated as: The analytical solution to Equation 3 given an initial condition of N(0) at t = 0 is: radioactivedecay evaluates Equation 5 upon each call for a decay calculation.
Matrices C and C −1 are independent of time so they are pre-calculated and imported from files into radioactivedecay. C and C −1 are stored in sparse matrix data structures to minimize memory use and maximize efficiency when computing the matrix multiplications in Equation 5. For decay calculations with double-precision floating-point operations, C and C −1 are stored in SciPy (Virtanen et al., 2020) Compressed Sparse Row (CSR) matrix data structures. Conversely, they are stored in SymPy (Meurer et al., 2017) SparseMatrix data structures for high numerical precision calculations.
The high numerical precision decay calculation mode resolves numerical issues arising from using double-precision floating-point numbers for decay calculations for chains containing nuclides with disparate half-lives. One example is the decay chain for 254 Es, which contains 238 U (4.468 billion year half-life) and 214 Po (t 1/2 is 164.3 µs half-life). This a 20 orders of magnitude difference in half-life. Loss of numerical precision inevitably occurs when evaluating the off-diagonal elements of C and C −1 in Equation 4 with double-precision floating-point numbers (which hold approximately 15 decimal places of numerical precision). Note loss of precision also occurs in the converse scenario, i.e. when a decay chain contains radionuclides with similar half-lives. However this scenario does not occur in the ICRP Publication 107 decay dataset, as the relative difference between half-lives of any two radionuclides in the same decay chain is always greater than 0.1%.
The default operation of the high precision decay mode is to evaluate Equation 5 using floating-point numbers with 320 significant figures of precision. This is sufficient precision to ensure accurate results for any physically relevant decay calculation users may wish to perform. Moreover, computations in the high precision mode are still fast, taking less than one second on a notebook equipped with an Intel Core i5-8250U processor.

Decay & Atomic Mass Datasets
The default dataset supplied with radioactivedecay uses decay data from ICRP Publication 107 (Eckerman & Endo, 2008) and atomic masses from the Atomic Mass Data Center (AMDC) Kondev et al., 2021;Wang et al., 2021). Endo et al. (2005) and Endo & Eckerman (2007) describe the development of the ICRP Publication 107 decay dataset. Raw data from ICRP 107 and AMDC were converted into dataset files suitable for radioactivedecay in a Jupyter notebook. Along with SciPy and SymPy versions of the sparse matrices C and C −1 , the dataset files contain radionuclide half-lives, decay constants, progeny, branching fractions, decay modes and atomic masses. Although there is a default dataset, radioactivedecay allows the import and use other decay data. The main functionality of radioactivedecay is based around Nuclide, Inventory and InventoryHP classes. The Nuclide class is used for fetching atomic and decay data about a single nuclide, such as its atomic mass, half-life, decay modes, progeny and branching fractions. It creates diagrams of the nuclide's decay chain (ex. Figure 1(a)) using the NetworkX library (Hagberg et al., 2008).

Main Functionality
An Inventory can contain multiple nuclides, each with an associated quantity (the number of atoms of the nuclide). Nuclides can be stable or radioactive. The decay() method calculates the decay of the radioactive nuclides in an Inventory, adding any ingrown progeny automatically. The numbers(), activities(), masses(), and moles() methods output the inventory of nuclides as different quantities using the atomic data stored in the decay dataset. Additional activity_fractions(), mass_fractions(), and mole_fractions() methods provide the relative amounts of each nuclide in the inventory with respect to different quantities. Plots can be made of the variation of nuclide activities, masses and moles over time (ex. Figure 1(b)) using Matplotlib (Hunter, 2007).
The InventoryHP class is the high numerical precision complement of the Inventory class. It has the same API as the Inventory class, but uses SymPy high numerical precision routines for all calculations.

Validation
Decay calculations with radioactivedecay v0.4.2 were cross-checked against Radiolog ical Toolbox v3.0.0 (Hertel et al., 2015) and PyNE v0.7.5 (Scopatz et al., 2012) (see Jupyter notebooks in the comparisons repository). Radiological Toolbox employs the ICRP Publication 107 decay data. Fifty radionuclides were randomly selected and a decay calculation was performed for 1 Bq of each for a random decay time within a factor of 10 −3 to 10 3 of the half-life. Differences between decayed activities reported by each code were within 1% of each other in 64% of cases. Discrepancies greater than 1% were attributed to rounding differences, erroneous results from Radiological Toolbox, or numerical issues relating to decay chains containing radionuclides with disparate half-lives.
A dataset was prepared for radioactivedecay with the same Evaluated Nuclear Structure Data File (ENSDF, 2019) decay data as used by PyNE v0.7.5. Bugs in PyNE v0.7.5 cause incorrect decay calculation results for chains containing metastable nuclear isomers, 183 Pt, 172 Ir or 152 Lu. Thus the affected chains were not used for the comparisons. The decay of 1 Bq of every radionuclide was calculated for multiple decay times varying from 0 to 10 6 times the radionuclide's half-life. The absolute difference between the decayed activities reported by each code was less than 10 −13 Bq. Relative differences depended on the magnitude of the activity. Relative errors of greater than 0.1% only occurred when the calculated activity was less than 2.5 × 10 −11 Bq, i.e. 10 orders of magnitude smaller than the initial activity of the parent radionuclide. The discrepancies between the two codes were attributed to methodological differences for computing decay chains with radionuclides with large disparities between half-lives, and numerical issues arising from double-precision floating-point operations.
Limitations radioactivedecay does not model neutronics, so cannot evaluate radioactivity produced from activations or induced fission. It does not support external sources of radioactivity input or removal from an inventory over time. Caution is required if decaying backwards in time, as this can cause floating-point overflows when computing the exponential terms in Equation 5.
There are also some limitations associated with the ICRP Publication 107 decay dataset. It does not contain data for the radioactivity produced from spontaneous fission decay pathways and the minor decay pathways of some radionuclides. More details on limitations are available in the documentation, Endo et al. (2005), and Endo & Eckerman (2007).