Theory and Computation
Solution to the decay chain differential equations
radioactivedecay calculates an analytical solution to the decay chain
differential equations using the method outlined by Amaku et al. [1].
Consider a system of \(n\) radionuclides, where the vector \(\mathbf{N}(t)\) has elements \(N_{i}(t)\) representing the number of atoms of radionuclide \(i\) at time \(t\). The radioactive decay chain differential equations can be expressed using matrix notation as
\(\varLambda\) is a lower triangular matrix whose elements are
where \(\lambda_{j}\) is the decay constant of radionuclide \(j\), and \(b_{ji}\) is the branching fraction from radionuclide \(j\) to \(i\).
\(\varLambda\) is a diagonalizable matrix so its eigendecomposition can be used to rewrite the matrix differential equation as
Here \(\varLambda_d\) is a diagonal matrix whose elements along the diagonal are the negative decay constants (\(\varLambda_{dii} = -\lambda_{i}\), these are the eigenvalues of \(\varLambda\)). Matrix \(C\) and its inverse \(C^{-1}\) are both lower triangular matrices. The column vectors of \(C\) are the eigenvectors of \(\varLambda\). Amaku et al. showed the elements of \(C\) can be calculated as
The elements of the inverse matrix \(C^{-1}\) can then be calculated as
The analytical solution to the matrix differential equation given the intial condition \(\mathbf{N}(t=0)=\mathbf{N}(0)\) is
This is the equation that is calculated by radioactivedecay upon each call
to decay(). \(e^{\varLambda_{d} t}\) is the matrix exponential of
\(\varLambda_{d} t\). It is a diagonal matrix with elements
\(e^{\varLambda_{d} t}_{ii} = e^{-\lambda_i t}\).
The final equations that are needed are for converting between various quantities. Converting between mass (\(\mathbf{M}\), in grams) and number of atoms (\(\mathbf{N}\)) uses the vector of atomic masses (\(\mathbf{M_u}\)) and the Avogadro constant (\(N_a\)):
Converting between activity (\(\mathbf{A}\)) and number of atoms (\(\mathbf{N}\)) uses the vector of decay constants (\(\mathbf{\lambda}\)):
Implementation in radioactivedecay
As matrices \(C\) and \(C^{-1}\) are independent of time, they can be
pre-calculated and imported into radioactivedecay as part of a radioactive
decay dataset. \(C\) and \(C^{-1}\) are pre-calculated for the
ICRP-107 [2] dataset in
this Jupyter notebook.
Note that \(C\), \(C^{-1}\) and \(e^{\varLambda_{d} t}\) are sparse matrices. Storing the matrices in a sparse matrix format greatly reduces the amount of memory that they occupy and the computational expense for the matrix multiplications in each decay calculation.
Notes on computation and numerical precision
Numerical issues can arise when using double-precision floating-point numbers to compute decays for some chains [3]. For example, the Es-254 decay chain contains U-238 (half-life 4.468 billion years) and Po-214 (half-life 164.3 microseconds), which is a 20 orders of magnitude difference in half-life. Round-off errors and loss of significance can occur causing unphysical results, e.g.
>>> inv = rd.Inventory({'Es-254': 1.0})
>>> inv.decay(0.0).activities()
{'At-218': -8.24439035981494e-30, 'Bi-210': 2.5308844932098316e-26,
'Bi-214': -4.256549745172888e-26, 'Bk-250': 0.0,
'Cf-250': 0.0, 'Cm-246': 8.802967479989175e-21,
'Es-254': 1.0, 'Fm-254': 0.0,
'Hg-206': -3.4696439711117526e-34, 'Pa-234': 2.330729590281097e-29,
'Pa-234m': -1.5696690930108473e-26, 'Pb-206': 0.0,
'Pb-210': 2.673060958594837e-26, 'Pb-214': -7.310828272597407e-27,
'Po-210': -1.048466176320909e-27, 'Po-214': 2.3260114484256133e-26,
'Po-218': -1.1433437709020225e-26, 'Pu-242': 1.3827905917787723e-22,
'Ra-226': -1.0811575068833228e-26, 'Rn-218': -1.618765025703667e-33,
'Rn-222': -1.581593359682259e-26, 'Th-230': -1.2628442466252288e-26,
'Th-234': -2.6140879622245746e-27, 'Tl-206': -4.332210492987691e-34,
'Tl-210': 2.2028710112960294e-31, 'U-234': -1.0389580591195201e-26,
'U-238': -8.466705440297454e-27}
All the progeny of Es-254 should have an activity of exactly zero for this calculation.
radioactivedecay thus offers a decay calculation mode using SymPy
[4] arbitrary precision computation routines for when high
numerical precision is needed:
>>> inv = rd.InventoryHP({'Es-254': 1.0})
>>> inv.activities()
{'At-218': 0.0, 'Bi-210': 0.0,
'Bi-214': 0.0, 'Bk-250': 0.0,
'Cf-250': 0.0, 'Cm-246': 0.0,
'Es-254': 1.0, 'Fm-254': 0.0,
'Hg-206': 0.0, 'Pa-234': 0.0,
'Pa-234m': 0.0, 'Pb-206': 0.0,
'Pb-210': 0.0, 'Pb-214': 0.0,
'Po-210': 0.0, 'Po-214': 0.0,
'Po-218': 0.0, 'Pu-242': 0.0,
'Ra-226': 0.0, 'Rn-218': 0.0,
'Rn-222': 0.0, 'Th-230': 0.0,
'Th-234': 0.0, 'Tl-206': 0.0,
'Tl-210': 0.0, 'U-234': 0.0,
'U-238': 0.0}
The InventoryHP class decay() method carries exact SymPy expressions
through decay calculations as far as is practicable. At the final step, the
decayed activity for each radionuclide is evaluated to high numerical precision
and cast to a double-precision float to return the decayed Inventory.
In practice using SymPy to exactly evaluate the exponential terms in the above
analytical solution to the radionuclide decay equations can be very time
consuming. Therefore by default the exponential terms are evaluated numerically
to 320 significant figures of precision mid-decay calculation. Empirically this
was found to give results for a range of test decay calculations, i.e. using a
higher number of significant figures offered no improvement in the numerical
accuracy of the results after the outputs are cast to double-precision floats.
You can also select your own number of significant figures for the calculation
by setting the InventoryHP.sig_fig attribute of the InventoryHP
instance.
References
M Amaku, PR Pascholati & VR Vanin, Comp. Phys. Comm. 181, 21-23 (2010). DOI: 10.1016/j.cpc.2009.08.011
ICRP Publication 107: Nuclear Decay Data for Dosimetric Calculations. Ann. ICRP 38 (3), 1-96 (2008). PDF
RI Balkin et al. Atomic Energy 123, 406-411 (2018). 10.1007/s10512-018-0360-2
A Meurer et al. PeerJ Comp. Sci. 3, e103 (2017). DOI: 10.7717/peerj-cs.103