Algorithms and algebraic solutions of decay chain differential equations for stable and unstable nuclide fractionation

Computer Physics Communications xxx (xxxx) xxx

Contents lists available at ScienceDirect

Computer Physics Communications journal homepage: www.elsevier.com/locate/cpc

Algorithms and algebraic solutions of decay chain differential equations for stable and unstable nuclide fractionation ∗

Austin Ladshaw a , , Alexander I. Wiechert a , Yong-ha Kim a,b , Costas Tsouris a,c , ∗ Sotira Yiacoumi a , a

Georgia Institute of Technology, Atlanta, GA 30332, United States of America Louisiana State University, Baton Rouge, LA 70803, United States of America c Oak Ridge National Laboratory, Oak Ridge, TN 37831-6181, United States of America b

article

info

Article history: Received 15 May 2019 Received in revised form 11 July 2019 Accepted 27 August 2019 Available online xxxx Keywords: Radioactive Nuclide Decay Modeling

a b s t r a c t Radioactive decay processes, such as alpha decay, produce decay chains where the mass numbers of nuclides decrease as larger nuclides expel energetic particles to form smaller nuclides. Under these conditions, the coefficient matrix that describes the differential rate expressions for radioactive decay can be made lower triangular. With this special structure, formulating an algebraic solution to the decay chains can be done by first formulating the eigenvectors that make up the coefficient matrix, which can then be solved using forward substitution for a lower triangular matrix. This work details the derivation of algebraic solutions for decay chains of any number of stable and unstable nuclides with any number of branching based on this eigenvector analysis. A prototype computational code was developed to validate and compare this methodology against a number of other methods for solving similar systems. A two-phase sorting algorithm yielding the lower triangular matrix structure was established to apply the developed algebraic solutions for decay chains involving beta-emitting radionuclides transformed into daughter nuclides without change in their mass number. The methodologies produced in this work provide an efficient way to estimate nuclide fractions from natural decay processes. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Radionuclides refer to nuclides, either naturally or artificially produced, undergoing radioactive decay. The major source of naturally generated radionuclides includes the decay of uranium and cosmic rays, while that of artificially produced radionuclides is human activities (e.g., nuclear reactions at nuclear reactors). Naturally generated radionuclides include uranium and its progeny (e.g., radon), and are widely distributed in the atmosphere, hydrosphere, and geosphere. Artificially created radionuclides involve cesium-134, cesium-137, and many others, which can be released into the environment by nuclear reactor accidents and explosions, such as the Chernobyl and Fukushima incidents [1,2]. Both natural and artificial radionuclides can have tremendous implications for human society and the environment because they can cause severe radiation damage to living organisms, including humans. For example, radon is carcinogenic, and the fatality rate of radon exposure in stand-alone houses may reach about 20,000 fatalities per year [3]. For the Fukushima case, radiation ∗ Corresponding authors. E-mail addresses: [email protected] (A. Ladshaw), [email protected] (S. Yiacoumi).

exposure caused by the fallout is currently projected to cause up to 1300 cancer mortalities and 2500 cancer morbidities around the world, respectively [4]. The concentrations of radionuclides in the environment are influenced by various factors such as the atmospheric conditions and the size of nuclear events [5,6]. In recent modeling and observation studies [7–10], it has been demonstrated that radioactive aerosols are more likely to travel further through the atmosphere than non-radioactive particles. This is attributed to charging of the radioactive nuclide containing aerosols through emission of electrons and diffusion of ions, which in turn reduce the size growth of the aerosols and keep them suspended in the atmosphere for longer periods of time. Therefore, including radioactive decay mechanisms is necessary for understanding how radionuclides travel in the environment. Due to the complexity introduced from the inclusion of particle charging into atmospheric modeling of aerosol aggregation and transport, it is critical that solutions to decay mechanisms be computationally efficient. For this reason, prior work in this area [10] has focused on the utilization of the Bateman Equations for simple decay chains. These equations, however, usually assume that decay chain involves no branching and that there is only one non-zero initial nuclide [10,11]. There are ways to extend the Bateman Equations to include branching systems

https://doi.org/10.1016/j.cpc.2019.106907 0010-4655/© 2019 Elsevier B.V. All rights reserved.

Please cite this article as: A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al., Algorithms and algebraic solutions of decay chain differential equations for stable and unstable nuclide fractionation, Computer Physics Communications (2019) 106907, https://doi.org/10.1016/j.cpc.2019.106907.

2

A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al. / Computer Physics Communications xxx (xxxx) xxx

[12,13] and even some numerical methodologies developed for more complicated nuclide fractionations [14]; however, few studies have produced a closed-form, analytical solution to decay chains of any size and any amount of branching for both stable and unstable nuclide fractionation. Thus, it is paramount to devise an efficient algorithm to evaluate complex nuclide fractionations in atmospheric simulations to better understand and quantify how radioactivity influences transport of radioactive aerosols. Most recently, Amaku et al. [15], Levy [16], and Yuan and Kernan [17] have developed various solution methodologies to the decay chain problems based on algebraic manipulations and matrix analyses for unstable nuclide decay equations that result in lower triangular coefficient matrices. The present study is built upon those methodologies to demonstrate an exact, algebraic solution for the unstable nuclides and extend that solution to the stable nuclides given any problem size or initial condition. Combined with algorithms to sort the nuclides in a descending atomic mass and parent–daughter order, the coefficient matrix can always be constructed to be of lower triangular shape. Once in this form, the coefficient matrix has eigenvectors that can be efficiently computed following the procedure of Amaku et al. [15] and stored in memory to formulate the solution to the fractionation of any nuclide at any point in time. To validate this approach, solutions are compared, both algebraically and numerically, against the simpler Bateman Equations for special case decay chains and the generalized solution developed in this work. Comparisons are also made against other approximate solutions that utilize matrix polynomials to exponentiate the coefficient matrix. The mathematical equivalence between the algebraic solutions developed in this work and the matrix analysis solutions developed by Levy [16] is also demonstrated. Additionally, the pros and cons of this methodology are discussed, as well as the solution limitations and ways to further improve computational efficiency. The ultimate aim for this methodology will be to incorporate radionuclide decay chains into mesoscale nuclear fallout transport models to predict movement of debris clouds and their deposition over short ranges. 2. Methods 2.1. Solution to eigenvectors In nuclear reactions, nuclides are produced through radioactive decay, neutron absorption, spontaneous fission, and occasionally fusion [14]. In the environment, neutron absorption and fusion reactions can be neglected. Under these conditions, the mass number of an nuclide in a decay chain will always either remain unchanged, or will be reduced as protons and/or neutrons are ejected from the nucleus [10–12]. A decay process with branching decay chains can be described as: dzi dt

=

∑

bij λj zj − λi zi

(1)

j
where λi is the decay constant (s−1 ) for the ith nuclide, bij is the branch fraction of the jth nuclide that forms the ith nuclide, and zi is the number concentration of the ith nuclide. Note that this version of the equation assumes that the nuclides can be ordered such that only higher indexed nuclides form from the low index nuclides (i.e., ith nuclides are only formed by nuclides whose index is less than i). This ordering is important as it creates a lower triangular coefficient matrix (note that the

indexing notation starts from 0 for all vectors and matrices in this paper):

⎡

0

···

···

−λ1

0

b21 λ1

−λ2 .. .

··· .. . .. . bN ,N −1 λN −1

−λ0

⎢ ⎢ b10 λ0 ⎢ ⎢ A = ⎢b λ ⎢ 20 0 ⎢ . ⎣ .. bN0 λ0

..

. ···

bN ,N −2 λN −2

0

⎤

.. . .. .

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ 0 ⎦ −λN

(2)

For any linear system of ordinary differential equations, there is a known analytical solution involving matrix exponentiation [12,14]: zt = eAt z0

(3)

where zt is a vector of nuclide number concentrations at time t and z0 is the initial number concentrations of all unstable nuclides. The exponential of a matrix is typically solved using a polynomial approximation, but the approximation is computationally very costly [14]. Alternatively, the coefficient matrix, multiplied by the scalar time (At), can be factored into a matrix of all eigenvectors (V), its inverse matrix (V−1 ), and a diagonal matrix of the exponential terms with the corresponding eigenvalues (Λ):

v00 ⎢ = ⎣...

⎡

eAt = VΛV−1

··· .. . ··· ⎤ N

vN0

u00

⎡

⎢ × ⎣...

u0N

··· .. . ···

⎤ ⎡ −λ t v0N e 0 .. ⎥ ⎢.. . ⎦ ⎣. vNN

0

··· .. . ···

0

.. .

e−λN t

⎤ ⎥ ⎦

u0

.. ⎥ . ⎦

(4)

uNN

Note that the vki elements of V are the ith row element of the kth eigenvector and uki elements are the same for V−1 [12]. Utilizing this method would produce an exact solution, but would involve the determination of all eigenvectors and eigenvalues of A. Since A is lower triangular, all eigenvalues for A are simply the diagonal elements of A, indicating that –λi is the ith eigenvalue of A. Additionally, since the coefficient matrix comes from a linear system of nuclides with distinctive decay constants, each eigenvalue will be unique. Therefore, combining these two factors together creates a set of corresponding eigenvectors that must also be unique, while making matrix V a unit lower triangular matrix [12,15]:

⎡

1

0

⎢ 0 ⎢v ⎣..

1

1 V=⎢ ⎢.

vN0

..

. ···

··· .. . .. . vNN−1

0

⎤

.. ⎥ .⎥ ⎥ ⎥ 0⎦

(5)

1

Given that A and V are triangular, each eigenvector can be obtained efficiently using the standard eigenvalue equation:

(A + λk I) vk = 0

(6)

Eq. (6) can be solved using forward substitution of the linear system defined by the kth eigenvalue (i.e., –λk ). When this system of equations is applied to solve for each eigenvector of A, each element of V can be defined as [15]:

vik =

⎧ 0 ⎪ ⎪ ⎨1

1

⎪ ⎪ ⎩− λk − λi

∑

vjk bij λj

for for

i
for

i>k

(7)

j
Please cite this article as: A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al., Algorithms and algebraic solutions of decay chain differential equations for stable and unstable nuclide fractionation, Computer Physics Communications (2019) 106907, https://doi.org/10.1016/j.cpc.2019.106907.

A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al. / Computer Physics Communications xxx (xxxx) xxx

Similarly, elements of V−1 in Eq. (4) can be calculated by a recurrence relationship formulated from the definition of inverse matrices (Eq. (8)) applied iteratively for each kth column of the identity matrix (ek ) to formulate the kth column of V−1 , as represented by uk (Eq. (9)) [15]. Vuk = ek

uki =

(8)

⎧ 0 ⎪ ⎪ ⎨1 ∑

− ⎪ ⎪ ⎩

j ukj i

v

for for

i
for

i>k

(9)

j
2.2. Solution for unstable nuclides Since matrices V and V−1 are both lower triangular, and Λ is diagonal, their multiplication together for the matrix exponential (Eq. (4)) yields:

⎡

f00

0

⎢ ⎢ f10 At e =F=⎢ ⎢ . ⎣ ..

f11

fN0

··· .. . .. .

..

. ···

fN −1,N −1

0

⎤

.. ⎥ . ⎥ ⎥ ⎥ 0 ⎦

(10)

fNN

where F is the exponential matrix and fij is the ith-jth matrix element. Thus, the exponential matrix of Eq. (10) can also be represented by a lower triangular matrix [16] with elements fij formed through a recursion relationship with the eigenvectors as follows:

fij =

⎧ 0 ⎪ ⎪ ⎪ ⎨e−λi t

for for

∑

for

i
i

⎪ ⎪ ⎪ ⎩

vik ujk e−λk t

i>j

(11)

⎞ ⎛ i ∑ ∑ j k −λ t t ⎝ vi uk e k ⎠ zj0 zi = j≤i

dt

bij λj zj

.. .

0

···

0

0

− 0

| |

z0

.. ⎥ ⎢ .. ⎥ ⎢ ⎥ .⎥ ⎥⎢ . ⎥ 0 ⎥ ⎢ zN ⎥ ⎥⎢ ⎥ −⎥ ⎢ − ⎥ ⎢ ⎥ 0⎥ ⎥ ⎢ N0 ⎥ ⎢ . ⎥ .. ⎥ . ⎦ ⎣ .. ⎦

.. .

| | + |

0

··· .. . ··· − ··· .. .

⎤⎡

⎤

(14)

NM

Eq. (14) demonstrates the separation of the coefficient matrix A from coefficient matrix B representing the production of stable species. If the rate equations (13) for the stable nuclides were included in Eqs. (1) and (2), it would introduce zeros on the diagonal of A and the matrix would have duplicate eigenvalues, thereby complicating the evaluation of eigenvectors. As discussed in Section 2.1, the solution methodology for unstable nuclides is specific to a linear system with unique eigenvalues and eigenvectors. Therefore, for simplicity the solution to the stable nuclides was formulated by simply integrating equation (13) using the algebraic solution developed for the unstable nuclides from Eq. (12). This produces the result shown in Eqs. (15) through (17), where Eq. (17) represents the final algebraic solution for any stable nuclide. dNi dt

=

∑

bij λj

( j ∑ ∑

∀j∈i

Nit = Ni0 +

k≤j

∫ t∑ 0

Nit = Ni0 +

∑

∀j∈i

) v

l k −λl t j ul e

zk0

(15)

l=k

bij λj

( j ∑ ∑ k≤j

) vjl ukl e−λl t zk0 dt

⎡ ( j ( )) ⎤ −λl t ∑ ∑ 1 − e zk0 ⎦ bij λj ⎣ vjl ukl λl k≤j

(16)

l=k

(17)

l=k

2.4. Sorting of the unstable nuclides

In some fields, such as nuclear forensics [18,19], evaluating only the decay of unstable nuclides and isotopes may be insufficient. For these instances, it may be necessary to calculate the accumulation of stable isotopes (i.e., the by-products of decay) and their ratios to their unstable counterparts. The stable byproducts were not considered in deriving the decay solution since they would complicate the matrix analysis and do not impact the decay of unstable nuclides. Similarly to Eq. (1), accumulation of stable nuclides (Ni ) can be expressed as:

∑

0

|

(12)

k=j

2.3. Solution for stable nuclides

=

A

][ ] | 0 z − + − − B | 0 N ⎡ −λ0 0 0 ⎢ .. .. ⎢ . . 0 ⎢ ⎢bN0 λN · · · −λN ⎢ − − =⎢ − ⎢b λ ⎢ 00 0 · · · b0N λN ⎢ . .. .. ⎣ .. . . bM0 λ0 · · · bMN λN

[

∀j∈i

Note that the compact version of the solution shown in Eq. (12) j takes advantage of the fact that vki uk = 1 when i = j, which is −1 due to the fact that VV = I. With the entire set of the computed eigenvectors evaluated, this solution form can be used to query any unstable nuclide at a given point in time (zit ) based on the initial conditions of all unstable nuclides (zi0 ).

dNi

over all j from A that affects Ni by producing the stable nuclide i. The matrix representation for stable and unstable nuclides would be shown as:

k=j

Multiplying the F matrix by the initial condition vector (Eq. (3)) provides an algebraic solution of the original problem (Eq. (1)). This approach creates a closed form, analytical solution for a decay chain of nuclides with any number of branches:

3

(13)

∀j∈i

In Eq. (13), the decay term for the current nuclide was ignored since the stable nuclides do not decay. Note that here the sum is

In order to achieve the lower triangular coefficient matrix, as shown in Eq. (2), one must sort the unstable nuclides such that each nuclide of a higher index (i.e., daughter nuclides) is formed strictly by nuclides of a lower index (i.e., parent nuclides) (Eq. (1)). This restriction will only be possible during a natural decay process, in which parent nuclides (P) decrease in mass number (A), or whose mass number remains unchanged as in beta decay, to some daughter nuclide (D) through alpha decay, proton and neutron emission, or other processes (Table 1) [20–23]. Thus, this methodology cannot be employed for nuclide fractionations that result from neutron capture/absorption or fusion reactions [14]. For beta decay, the mass number (A) of radionuclides remains unchanged, but the atomic number (Z) can either increase or decrease (Table 1) [20–23]. Thus, sorting on mass number alone is insufficient and should be combined with some form of sub-sorting by parent–daughter relationships [17]. This requirement then brings forth another restriction to application of this method; there can be no cyclical relationships in the decay process. In other words, nuclide z cannot decay to another nuclide that could eventually reform nuclide z. Sorting the nuclides under these restrictions can be achieved with a two-step sorting algorithm. First, an insertion sort, based on mass number, is used to add unstable nuclides (zi ) to a growing

Please cite this article as: A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al., Algorithms and algebraic solutions of decay chain differential equations for stable and unstable nuclide fractionation, Computer Physics Communications (2019) 106907, https://doi.org/10.1016/j.cpc.2019.106907.

4

A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al. / Computer Physics Communications xxx (xxxx) xxx

Table 1 Summary of decay modes that do not result in increasing mass numbers. Decay mode

Decay reaction

Alpha

A ZP

→ AZ−−24 D + 42 He

Beta (Beta-)

A ZP

→ AZ+1 D + v + e−

Proton emission

A ZP

→ AZ−−11 D + 11 H

Neutron emission

→ AZ −1 D + 10 n

Electron capture (Beta+)

A ZP A ZP

Spontaneous fission

A ZP

→ AZ′ D + AZ−−ZA′ D

→ AZ−1 D + v ′

′

list of nuclides [24]. The insertion sort starts by adding nuclides to the list, beginning with the nuclides known to be initially present, from highest to lowest mass numbers (A). Each decaying nuclide then has its unstable daughters and/or emitted particles added to that list also following the highest to lowest mass numbers. If two or more nuclides have the same mass number, then the nuclide being added gets placed just before the next nuclide with the lower mass number. At the end of this initial sorting, all highest mass nuclides will be at the top of the list, but there will likely be groups of nuclides within the list that may have the same mass numbers. These sub-groups also require sorting. The final sorting algorithm must iterate through the partially sorted list to find all the groups of nuclides with the same mass numbers. Those groups must be sorted in place within that list such that the daughters of unstable nuclides occur below the parents that they correspond to (i.e., the sub-list must follow genealogical order [17]). This is accomplished by iterating through the group to find the nuclide that has no parent within that group, then moving that nuclide to the front of the group, and continuing from the next position down. Fig. 1 provides a visualization of this final sorting procedure. 2.5. Implementation details In order to implement the algebraic decay solution described in this work, a computer code was devised that can handle the sorting methods, eigenvector calculations, and decay chain relationships when given some starting set of nuclides. To properly utilize the decay chain relationships and nuclide data, a python script developed by Hykes [25] was employed to read and store decay data associated with nuclides registered in the ASCII formatted Nuclide Wallet Cards [23] distributed by Brookhaven National Laboratory. Once the data were read, a different python script, developed in this work, was used to reformat that information and produce a yaml [26] formatted version of the nuclide library that could be more easily and readily used by other computer programs or persons. The new yaml version of the nuclide library is then implemented to develop the decay chains that make up the stable and unstable nuclide lists. This was accomplished in C++ by creating objects for both single nuclides and a chain of nuclides. The chains were made up of a list of nuclides registered in the library, sorted using the algorithms described in Section 2.4, and then produced the eigenvectors and eigenvalues according to the new, sorted relationships between all other nuclides in the list (as shown in Sections 2.1 through 2.3). Source code for the algorithms associated with the nuclide and decay chain objects has been implemented into the Ecosystem Software [27] and is dubbed the Implicit Branching Isotope System or IBIS, named for the sacred Egyptian wading birds. This prototype software will be used to evaluate the validity of the methodologies described in this work, but does not represent a fully optimized software for nuclide fractionation analysis.

3. Results and discussion 3.1. Comparison to matrix analysis method Eq. (12) provides a simple algebraic expression for the solution to the matrix exponential of a lower triangular matrix and was produced using the recurrence relationship described by Amaku et al. [15]. A matrix analysis methodology was proposed by Levy [16], which uses the coefficient matrix’s commutativity properties (i.e., AG = GA). Based on this matrix property of A, Levy derived the following relationship for the exponential matrix:

gij =

⎧ 0 ⎪ ⎪ ⎪ ⎨e−λi t

( ) i−1 ∑ bij λj e−λi t − e−λj t fik bkj λj − bik λk fkj ⎪ ⎪ + ⎪ ⎩ λj − λi λj − λi

for for

i
for

i>j

k=j+1

(18) where eAt = G and gij are the ith row and jth column elements of G [16]. While Levy [16] used a completely different methodology to exponentiate the matrix A, the result was very similar in form to the result in Eq. (11) (eAt = G = F, where G and F are both lower triangular matrices). If both approaches were developed to solve the same problem, the final results should be equivalent. To show this, consider a 3×3 system for coefficient matrix A:

−λ0 A = b10 λ0 b20 λ0

[

0

−λ1

b21 λ1

0 0

] (19)

−λ2

By expanding out the terms in the F and G matrices, it will be demonstrated that they are mathematically equal to each other (i.e., fij = gij for all i and j). For the cases of i < j and i = j, the demonstration can be seen by comparing Eqs. (11) and (18). Thus, it only needs to be shown that fij = gij for i > j. The analysis starts by using Eqs. (7), (9), and (11) to produce f10 , f21 , and f20 . Eqs. (20) through (23) given in Box I show the evaluation of these terms. Then, Eq. (18) is used to calculate g10 , g21 , and g20 , yielding:

⎧ ⎪ ) b10 λ0 ( −λ1 t ⎪ ⎪ e − e−λ0 t g10 = ⎪ ⎪ ⎪ λ0 − λ1 ⎪ ⎨ ) b21 λ1 ( −λ2 t (24) g21 = e − e−λ1 t ⎪ λ1 − λ2 ⎪ ⎪ ⎪ ⎪ ) g b λ − b21 λ1 g10 ⎪ b λ ( ⎪ ⎩g20 = 20 0 e−λ2 t − e−λ0 t + 21 10 0 λ0 − λ2 λ0 − λ2 b λ b20 λ0 20 0 g20 = e−λ2 t − e−λ0 t (λ0 − λ2 ) (λ0 − λ2 ) b21 λ1 b10 λ0 + e−λ2 t (λ0 − λ2 ) (λ1 − λ2 ) b21 λ1 b10 λ0 b21 λ1 b10 λ0 − e−λ1 t − e−λ1 t (λ0 − λ2 ) (λ1 − λ2 ) (λ0 − λ2 ) (λ0 − λ1 ) b21 λ1 b10 λ0 + e−λ0 t (λ0 − λ2 ) (λ0 − λ1 ) (25) Note that in Eq. (18) the calculation order matters, because g20 depends on g10 and g21 . Immediately, it can be seen by simple inspection of Eqs. (22) and (24) that f10 and f21 are both equivalent to g10 and g21 . Showing the equivalence of f20 and g20 (Eqs. (23) and (25)) will take some additional analysis. While portions of f20 and g20 are obviously equivalent (i.e., the exp[−λ0 t] terms), equivalence of other terms is less clear. For instance, in order for the exponential coefficients of the λ1 terms to

Please cite this article as: A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al., Algorithms and algebraic solutions of decay chain differential equations for stable and unstable nuclide fractionation, Computer Physics Communications (2019) 106907, https://doi.org/10.1016/j.cpc.2019.106907.

A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al. / Computer Physics Communications xxx (xxxx) xxx

⎡

⎤

1 ⎢ ⎢ b10 λ0 ⎢ − V=⎢ (λ0 − λ1 ) ⎢ ⎣ b20 λ0 b21 λ1 b10 λ0 − + (λ0 − λ2 ) (λ0 − λ2 ) (λ0 − λ1 ) ⎡

−

0

0⎥

1

0⎥ ⎥

b21 λ1

(λ1 − λ2 )

⎥

(20)

⎥ ⎦

1

1

b10 λ0 ⎢ ⎢ (λ0 − λ1 ) ⎣ b20 λ0 b21 λ1 b10 λ0 b21 λ1 b10 λ0 − + (λ0 − λ2 ) (λ0 − λ2 ) (λ0 − λ1 ) (λ1 − λ2 ) (λ0 − λ1 ) b10 λ0 −λ0 t b10 λ0 −λ1 t =− e + e λ0 − λ1 λ0 − λ1 b21 λ1 −λ1 t b21 λ1 −λ2 t =− e + e λ1 − λ2 λ1 − λ2

V−1 = ⎢

⎧ ⎪ ⎪ f10 ⎪ ⎪ ⎪ ⎨ f21 ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

5

⎤

0

0

1

0⎥

b21 λ1

(λ1 − λ2 )

1

⎥ ⎥ ⎦

(21)

(22)

f20 = v20 u00 e−λ0 t + v21 u01 e−λ1 t + v22 u02 e−λ2 t

f20 = −

b20 λ0

e−λ0 t +

b21 λ1 b10 λ0

e−λ0 t −

b21 λ1 b10 λ0

e−λ1 t

(λ0 − λ2 ) (λ0 − λ2 ) (λ0 − λ1 ) (λ1 − λ2 ) (λ0 − λ1 ) b20 λ0 b21 λ1 b10 λ0 b21 λ1 b10 λ0 + e−λ2 t − e−λ2 t + e−λ2 t (λ0 − λ2 ) (λ0 − λ2 ) (λ0 − λ1 ) (λ1 − λ2 ) (λ0 − λ1 )

(23)

Box I.

be identical, the statement of Eq. (26) must be shown equal:

−

b21 λ1 b10 λ0 b21 λ1 b10 λ0 e−λ1 t = − e−λ1 t (λ1 − λ2 ) (λ0 − λ1 ) (λ0 − λ2 ) (λ1 − λ2 ) b21 λ1 b10 λ0 − e−λ1 t (λ0 − λ2 ) (λ0 − λ1 )

(26)

To show this, a common denominator must be chosen such that the left and right hand sides of Eq. (26) can be directly compared. By choosing (λ1 − λ2 )(λ0 − λ1 ) as the common denominator, the right hand side of Eq. (26) can be simplified as shown in Eq. (27): b21 λ1 b10 λ0 e−λ1 t (λ1 − λ2 ) (λ0 − λ1 ) [ b21 λ1 b10 λ0 (λ0 − λ1 ) = − (λ0 − λ2 ) (λ1 − λ2 ) (λ0 − λ1 ) ] b21 λ1 b10 λ0 (λ1 − λ2 ) −λ1 t − e (λ0 − λ2 ) (λ0 − λ1 ) (λ1 − λ2 ) [ b21 λ1 b10 λ0 (λ0 − λ1 ) (27) = − (λ0 − λ2 ) (λ1 − λ2 ) (λ0 − λ1 ) ] b21 λ1 b10 λ0 (λ1 − λ2 ) − e−λ1 t (λ0 − λ2 ) (λ0 − λ1 ) (λ1 − λ2 ) [ ] b21 λ1 b10 λ0 (λ0 − λ2 ) = − e−λ1 t (λ0 − λ2 ) (λ1 − λ2 ) (λ0 − λ1 ) b21 λ1 b10 λ0 =− e−λ1 t (λ1 − λ2 ) (λ0 − λ1 ) This demonstrates the equivalence of the λ1 exponential terms.

−

This same analysis can be applied to all terms in f20 and g20 to show that they are equal. While the example presented in Eqs. (19) through (27) shows that the eigenvector solution (Eq. (11)) is mathematically equivalent to the matrix communication result from Levy [16] (Eq. (18)), it is worth noting some major differences in the methodologies.

Although the Levy solution is algebraically simpler, as it does not involve the formation of the eigenvectors, it may actually be less efficient to implement computationally for the following reasons: (i) every time the need to query the decay products at a new point in time arises, Levy’s solution would demand that the gij matrix elements be recalculated and (ii) the solution uses the Parlett Recurrence [16] that needs previous matrix elements to be computed in a particular order before moving on to the next element. To take a look at these differences more closely, consider the following discussion. For the methodology presented in this work, the fij matrix elements never need to actually be formed because of the compact solution in Eq. (12) that already handles the matrix multiplication operation of eAt z0 . All operations to form the make up of the exponential have been handled ahead of time through the eigenvector analysis (Eqs. (6) through (9)), which only has to be performed once then stored in memory. With the solution of Levy, each time the decay problem needs to be solved at a different point in time, all matrix elements for G would have to be recalculated by Eq. (18). This would unnecessarily use computational resources. In addition, when forming the matrix G, the recurrence term used to simplify the solution forces the computations to be done in a particular order. For example, when the formation of matrix element g20 was examined (Eq. (24)), the elements g21 and g10 had to already be calculated. With the eigenvector solution, if all the eigenvectors are determined ahead of time, then there is no such restriction on the order of computation. This is especially important for high performance computing, as the eigenvector methodology presented in this work could be implemented with significant parallelism of the time-dependent solution for the decay chains (i.e., after forming the eigenvectors, Eq. (12) could be implemented with parallelism to calculate nuclide fractionation of N unstable nuclides on N computer processors). 3.2. Reproducing Bateman results Since the methodology presented in this work was for a generalized decay chain with any amount of branching, any type of

Please cite this article as: A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al., Algorithms and algebraic solutions of decay chain differential equations for stable and unstable nuclide fractionation, Computer Physics Communications (2019) 106907, https://doi.org/10.1016/j.cpc.2019.106907.

6

A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al. / Computer Physics Communications xxx (xxxx) xxx

Fig. 1. Visualization of the sorting algorithm to properly order a sub-group of nuclides of the same mass number by their parent–daughter relationships. In this example, a sub-group of nuclides of mass number 20 are being sorted in place by parent–daughter relationships. This sub-group contains isotopes of fluorine (F), oxygen (O), and nitrogen (N), which should decay in the following order: N-20 → O-20 → F-20. The sub-group reduces in size by 1 each time and the parent nuclide is moved to the front.

initial condition, and any number of nuclides, it should be able to produce other, less rigorous analytical solutions to specific decay chains. To show this, start by considering the Bateman Equations. Moral and Pacheco [12] have shown how the Bateman Equation [10,13] could be modified to account for branching paths for decay chains. In their demonstration, they considered the branching decay of Bi-212 to Pb-208, a stable isotope of lead. The decay constant of Bi-212 is 0.01145 s−1 (λ0 ), and it can decay into either Tl-208 (b10 = 0.3594) as an intermediate species or directly down to Pb-208 (b00 = 0.6406). Tl-208 then decays directly to Pb-208 (b01 = 1), and its decay constant is 0.2270 s−1 (λ1 ). The full coefficient matrix to describe this system, including accumulation of Pb-208, is shown as: A

| + |

[

− B

0

⎡ −λ0 ⎢b10 λ0 =⎣ − b00 λ0

z

][ ]

−

−

0

N

0

−λ1 − b01 λ1

| | + |

dz0 /dt ⎢dz1 /dt ⎥ =⎣ ⎦ − dN0 /dt

⎡

⎤⎡

⎤

z0 0 0 ⎥ ⎢z1 ⎥

−⎦ ⎣− ⎦

0

N0

⎤ (28)

In this representation, the unstable nuclide coefficient matrix (A) was partitioned from the portion that accounts for stable nuclides (B), as done in Eq. (14). Having previously shown that the solution to the unstable nuclides is mathematically equivalent to Levy’s matrix analysis of the unstable nuclides (recall Section 3.1), the intent here is to show that the solution for stable nuclides (Eq. (17)) will mirror the stable nuclide solution based on the generalized Bateman approach from Moral and Pacheco [12]. For this effort, the analysis starts by utilizing the eigenvector solution (Eqs. (6) through (9)) to formulate matrices V and V−1 from A:

[ V=

−

1 b10 λ0

λ0 − λ1

0

]

1

[ V−1 =

1 b10 λ0

λ0 − λ1

0

]

1

(29)

Then, by applying Eq. (17), we obtain an exact solution for Pb-208 as shown in Eqs. (30) through (32).

[ N0t

N00

=

+b00 λ0

z00

( j=0 ∑



l=0

0

v

l k j ul

1 (

λl

)] −λl t

1−e

)

⎡ ) ⎢ (∑ ⎢ 0 j=1 l k 1 ( ) −λl t ⎢ v j ul 1−e +b01 λ1 ⎢z0 λl ⎣ l=0

(30)

⎤ )⎥ ) ⎥ 1 ( +z10 vjl ukl 1 − e−λl t ⎥ ⎥ λl ⎦ l=1    0 [ ( )] ( ) 0 0 0 1 t N0 = b00 λ0 z0 v0 u0 1 − e−λ0 t [ ( λ0 )] ) ) 1 ( 1 ( +b01 λ1 z00 v10 u00 1 − e−λ0 t + v11 u01 1 − e−λ1 t λ0 λ1 ( j=1 ∑

(31)

[

N0t = b00 − b00 e−λ0 t +

b01 b10 λ0

b01 b10 λ0

− e−λ1 t (λ0 − λ1 ) (λ0 − λ1 ) ] b01 b10 λ0 λ1 b01 b10 λ0 λ1 −λ0 t 0 − + e z0 (λ0 − λ1 ) λ0 (λ0 − λ1 ) λ0

(32)

Note that, for this analysis, following the assumption of Moral and Pacheco [11], we assume that initially there is no Pb-208 or Tl208. Moral and Pacheco [12] gave the solution to Pb-208 as shown in Eq. (33). N0t = 1 − 1.019e−λ0 t + 0.019e−λ1 t z00

(

)

(33)

Their solution was not written in terms of the decay constants and branch fractions, so the solution shown in Eq. (32) must be evaluated under the same conditions to compare the two results.

Please cite this article as: A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al., Algorithms and algebraic solutions of decay chain differential equations for stable and unstable nuclide fractionation, Computer Physics Communications (2019) 106907, https://doi.org/10.1016/j.cpc.2019.106907.

A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al. / Computer Physics Communications xxx (xxxx) xxx

Fig. 2. Comparison between the Bateman Equations and the IBIS code implementation of the methodology described in this work for decay of Te-132. These results provide visual and numerical evidence that the IBIS method can reproduce the Bateman solution.

After applying the values for each branch fraction and decay constant to the symbolic solution, it can be seen that the result produced from the modified Bateman solution Eq. (33) is exactly equal to the eigenvector analysis solution: N0t =

[( b20 +

b01 b10 λ0

−

b01 b10 λ0 λ1

)

(λ0 − λ1 ) λ0 ) ] b01 b10 λ0 −λ1 t 0 e−λ0 t − × − b20 − e z0 (λ0 − λ1 ) λ0 (λ0 − λ1 ) [ ] = 1 − 1.019e−λ0 t + 0.019e−λ1 t z00 (λ0 − λ1 )

b01 b10 λ0 λ1

(

(34) This demonstrates that the methodology presented in this work can recreate more simplistic solutions for special case problems, including the evaluation of the build-up of stable nuclides. In addition to demonstrating the algebraic equivalence between the Bateman Equations and the methodology presented in this work, results of the coded implementation of the method (the IBIS code [27]) were compared to results from literature for Te132 decay. The decay chain for Te-132 is a linear, non-branching chain to I-132 (with λ0 = 2.5 × 106 s−1 ) that then decays to the stable Xeon isotope of Xe-132 (with λ1 = 8.39×10−5 s−1 ). According to the Bateman Equations, the analytical solution to this non-branching system, where only Te-132 is initially present, would be as shown in Eqs. (35) through (37) [10,12]. Te − 132: I − 132: Xe − 132:

z0t = z00 e−λ0 t z1t

=λ

0 0 z0

(

(35)

e−λ0 t −e−λ1 t

)

(36)

λ1 −λ0

(

z2t = z00 1 −

λ1 e−λ0 t −λ0 e−λ1 λ1 −λ0

) t

(37)

When the results are compared between this solution and the IBIS solution developed in this work, it is shown that they are equivalent in both the decay of Te-132 and the formation of the other two nuclides over time (Fig. 2). 3.3. Comparison to numerical solutions The results shown previously demonstrate that these algebraic solutions (Eqs. (12) and (17)) can recreate the solutions from

7

Fig. 3. Comparison of nuclide fractionation between the CRAM solver in SCALE and the IBIS algorithms. The average Euclidean difference between the simulation cases is at most 10−4 with an overall average difference of approximately 10−5 .

other methodologies, such as the matrix analysis method [16] and the Bateman Equations [10,12]. Thus, the next logical comparison would be with other numerical methodologies or software. SCALE [14,28] is a computer program that can perform numerical approximations to decay chains. For this purpose, it utilizes either a numerical approximation to the matrix exponentiation Eq. (3) called MATREX [28] or employs the Chebyshev rational approximation method (CRAM) for the sparse A matrix [28,29]. Either of these methods can be chosen at runtime of the SCALE computer program. Comparisons between the IBIS solutions are compared to results from SCALE (version 6.2.3) for the same initial nuclides to evaluate accuracy and efficiency [28]. To evaluate the algebraic solution accuracy for stable and unstable nuclides (Eqs. (12) and (17)), simulations were run using IBIS for a decay chain of 56 unstable nuclides (ranging from V54 to H-3) and 42 stable nuclides (ranging from Cr-54 to He-3). Each unstable nuclide was given an initial condition of 1 mole, while the concentrations of the stable nuclides were assumed to be initially zero. Then, the fractionation produced from this system was simulated over 10 years of decay and evaluated at 10 different points within that 10-year period. The fractionation after each period in time was tabulated, and then the same simulation case was run with SCALE using the CRAM algorithm, which is regarded as the most accurate of the available numerical methods [28]. The difference between the CRAM and IBIS results was taken in terms of the average Euclidean norm across all time points for a single nuclide (i.e., average root-mean-square error). Fig. 3 shows the average norms across all nuclides, both stable and unstable. These results show that the IBIS methods compare very well against those from the CRAM method employed in SCALE, with the maximum Euclidean difference being around 10−4 and a total average Euclidean difference across all nuclides around 10−5 . There was no observable bias in the differences that would lead to higher errors in stable or unstable nuclides as the differences appear to be randomly distributed (Fig. 3). The observed Euclidean differences between CRAM and IBIS could be the result of a couple of factors. First, CRAM is an approximation method; therefore, the order of accuracy of CRAM is dependent on the number of residue parameters computed as part of the Chebyshev expansion of the matrix [29,30]. In SCALE, CRAM defaults to a 16th order expansion, which correlates to an average error of less than 0.01% for most nuclides [28]. With

Please cite this article as: A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al., Algorithms and algebraic solutions of decay chain differential equations for stable and unstable nuclide fractionation, Computer Physics Communications (2019) 106907, https://doi.org/10.1016/j.cpc.2019.106907.

8

A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al. / Computer Physics Communications xxx (xxxx) xxx

IBIS, there is no error produced from approximation, because the algebraic solutions derived in Eqs. (12) and (17) are exact representations of the nuclide solutions. The other source of difference between them could be minor changes in the branch fractions or decay constants for this set of nuclides. SCALE comes packaged with its own library of nuclide data, which is not necessarily the same as the library of nuclides that IBIS uses from the Nuclear Wallet Cards [23]. These two factors combined could account for the small observed differences in the nuclide fractionation for this simulation case. 4. Conclusions Decay chains that come from natural decay processes, such as alpha and beta decay, can be efficiently calculated based on a thorough analysis of the coefficient matrix of branch fractions and decay constants. In the environment where fusion and neutron absorption rarely occur, a decay chain will decrease in mass number, except in the case of beta decay, until a stable nuclide is formed. Because of this behavior, the coefficient matrix of the nuclide relationships can be made lower triangular in shape. This structural property of the matrix can be used to develop simple, algebraic solutions for the unstable nuclides. Then, the solution to the stable nuclides can be formulated by integrating their respective formation rate equations based on the solutions to the unstable nuclides. In order to achieve this solution methodology for any given set of unstable nuclides, a sorting algorithm was developed to construct the decay chain from a given set of starting nuclides. This algorithm starts by doing an insertion sort on the initial nuclides and their daughters, thus creating a pseudo ordered list by mass number. Several sub-groups of nuclides in that list will have the same mass number. Those sub-groups are then ordered based on parent–daughter relationships that exist between them. After the sub-groups are ordered, the list is finalized such that the resulting coefficient matrix has the necessary lower triangular structure. This solution methodology was demonstrated both analytically and numerically by comparing against known solutions to specific problems. Algebraic analyses showed the equivalence between the eigenvector solution for the lower triangular system and the matrix analysis solution proposed by Levy [16] as well as the relationship to the Bateman Equations [10,12] for various problems. Numerical simulations for larger systems of nuclides were compared against software solutions for established methods, such as CRAM [29,30], and demonstrated that the eigenvector solution (IBIS) proposed in this work was capable of reproducing results from more complex systems. While these results are very promising, the algorithms employed, as done in IBIS, are not yet fully optimized. IBIS was the prototype used to evaluate the validity of the algebraic eigenvector solutions produced in this work. As such, there are gains that can be made through more optimized programming methods in sorting and through the utilization of message passing interfaces (MPI) for parallel execution of the evaluation of nuclide fractionation. Once the eigenvectors have been formed, each nuclide solution can be solved completely independently of each other. As a result, there is significant potential for speed up if each nuclide, or groups of nuclides, were to be computed on separate nodes or processors, thereby vastly improving the overall computational efficiency of the method presented in this work. The efficiency of calculation for nuclide fractionation is especially important when integrating decay processes into other natural phenomena or environmental systems. One such end goal for these computations may be to incorporate decay processes in atmospheric models for post-detonation events or nuclear

accidents, as was done by Kim et al. [10] for the decay chain of Te-132. In realistic scenarios, nuclear events will produce various concentrations of nuclides, where each may decay and branch to form other unstable nuclides over a variety of time scales. The algebraic solution produced in this work provides a relatively simple and efficient way to incorporate these complex decay chains into other software meant to simulate natural decay processes and their impact on aerosol transport in environmental systems. To this end, the models elucidated here will be introduced into the Multiphysics Object-Oriented Simulation Environment (MOOSE) [31] as part of a larger effort to developed short-range simulations of radioactive aerosol microphysics and the impact of radioactive decay on aerosol aggregation and fallout. Acknowledgment This work was supported by the Defense Threat Reduction Agency under grant number HDTRA11810023. The project or effort depicted was or is sponsored by the Department of Defense, Defense Threat Reduction Agency. The content of the information does not necessarily reflect the position or the policy of the federal government, and no official endorsement should be inferred. Notice: This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the US Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paidup, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http: //energy.gov/downloads/doe-public-access-plan). References [1] K. Adachi, M. Kajino, Y. Zaizen, Y. Igarashi, Sci. Rep. 3 (2013) 2554, http://dx.doi.org/10.1038/srep02554. [2] U. Baltensperger, H.W. Gäggeler, D.T. Jost, B. Zinder, P. Haller, J. Aerosol Sci. 18 (1987) 685–688, http://dx.doi.org/10.1016/0021-8502(87)90097-8. [3] J.S. Puskin, C.B. Nelson, J. Air Waste Manag. Assoc. 39 (1989) 915–950, http://dx.doi.org/10.1080/08940630.1989.10466577. [4] J.E. Ten Hoeve, M.Z. Jacobson, Energy Environ. Sci. 5 (2012) 8743–8757, http://dx.doi.org/10.1039/c2ee22019a. [5] J.T. McGahan, E.J. Kownacki, Sensitivity of fallout predictions to initial conditions and model assumptions, Defense Nuclear Agency, DNA-3439F, 1974. [6] H.G. Norment, DELFIC: Department of Defense Fallout Prediction System: Vol I – Fundamentals, Defense Nuclear Agency, DNA-5159F-1, 1979. [7] G. Lujaniene, V. Aninkevicius, V. Lujanas, J. Environ. Radioact. 100 (2009) 108–119, http://dx.doi.org/10.1016/j.jenvrad.2007.07.015. [8] Y. Kim, S. Yiacoumi, C. Tsouris, J. Environ. Radioact. 143 (2015) 91–99, http://dx.doi.org/10.1016/j.jenvrad.2015.02.017. [9] Y. Kim, S. Yiacoumi, A. Nenes, C. Tsouris, Atmos. Chem. Phys. 16 (2016) 3449–3462, http://dx.doi.org/10.5194/acp-16-3449-2016. [10] Y. Kim, S. Yiacoumi, A. Nenes, C. Tsouris, J. Aerosol Sci. 114 (2017) 283–300, http://dx.doi.org/10.1016/j.jaerosci.2017.09.024. [11] R.C. Tompkins, Department of Defense Land Fallout Prediction System: Vol V – Particle Activity, US Army Nuclear Defense Laboratory, NDL-TR-102, 1968. [12] L. Moral, A. Pacheco, Amer. J. Phys. 71 (2003) 684–686, http://dx.doi.org/ 10.1119/1.1571834. [13] D.R. Vondy, Development of a General Method of Explicit Solution to the Nuclide Chain Equations for Digital Machine Calculations, Oak Ridge National Laboratory, ORNL/TM-361, 1962. [14] I.C. Gauld, G. Radulescu, G. Ilas, B.D. Murphy, M.L. Williams, D. Wiarda, Nucl. Tech. 174 (2011) 169–195, http://dx.doi.org/10.13182/NT11-3. [15] M. Amaku, P.R. Pascholati, V.R. Vanin, Comput. Phys. Comm. 181 (2010) 21–23, http://dx.doi.org/10.1016/j.cpc.2009.08.011. [16] E. Levy, Comput. Phys. Comm. 234 (2019) 188–194, http://dx.doi.org/10. 1016/j.cpc.2018.07.011.

Please cite this article as: A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al., Algorithms and algebraic solutions of decay chain differential equations for stable and unstable nuclide fractionation, Computer Physics Communications (2019) 106907, https://doi.org/10.1016/j.cpc.2019.106907.

A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al. / Computer Physics Communications xxx (xxxx) xxx [17] D. Yuan, W. Kernan, J. Appl. Phys. 101 (2007) 094907, http://dx.doi.org/10. 1063/1.2715785. [18] K. Mayer, M. Wallenius, I. Ray, Analyst 130 (2005) 433–441, http://dx.doi. org/10.1039/b412922a. [19] M.J. Kristo, S.J. Tumey, Nucl. Instrum. Methods Phys. Res. B 294 (2013) 656–661, http://dx.doi.org/10.1016/j.nimb.2012.07.047. [20] A. Beiser, Concepts of Modern Physics, sixth ed., McGraw-Hill, New York, 2003, pp. 418–473. [21] S. Hofmann, Nuclear Decay Modes, Institute of Physics, Bristol, 1996, pp. 143–203. [22] X.Z. Li, J. Tian, M.Y. Mei, C.X. Li, Phy. Rev. C. 61 (2000) 24610, http: //dx.doi.org/10.1103/PhysRevC.61.024610. [23] J.K. Tuli, Nuclear Wallet Cards, eighth ed., Brookhaven National Laboratory, Upton, 2011. [24] D.E. Knuth, The Art of Computer Programming – Volume 3: Sorting and Searching, second ed., Addison-Wesley, Reading, 1998. [25] J.M. Hykes, nuclide-data – a Python interface to nuclide data, https: //github.com/jhykes/nuclide-data, 2005 (accessed on 18.09.18).

9

[26] K. Simonov, LibYAML – A C library for parsing and emitting YAML, https: //github.com/yaml/libyaml, 2006 (accessed on 28.07.15). [27] A. Ladshaw, A. Wiechert, C. Tsouris, S. Yiacoumi, Ecosystem Software – A C/C++ library for environmental chemistry and adsorption, https://bitbucket.org/gitecosystem/ecosystem, 2015 (accessed on 15.11.18). [28] B.T. Rearden, M.A. Jessee (Eds.), SCALE Code System, Version 6.2.3, ORNL/TM-2005/39, Oak Ridge National Laboratory, Oak Ridge, 2018. [29] J.M. Hykes, R.M. Ferrer, Proceedings of the 2013 International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering – M&C 2013, American Nuclear Society, Sun Valley, 2013, pp. 1551–1564. [30] P. Maria, Nucl. Sci. Eng. 182 (2016) 297–318, http://dx.doi.org/10.13182/ NSE15-26. [31] D. Gaston, C. Newman, G. Hansen, D. Lebrun-Grandie, Nucl. Eng. Des. 10 (2009) 1768–1778, http://dx.doi.org/10.1016/j.nucengdes.2009.05.021.

Please cite this article as: A. Ladshaw, A.I. Wiechert, Y.-h. Kim et al., Algorithms and algebraic solutions of decay chain differential equations for stable and unstable nuclide fractionation, Computer Physics Communications (2019) 106907, https://doi.org/10.1016/j.cpc.2019.106907.