Method of False Rates for Radioactive Decay Chain Calculation

Many engineering applications involving radioactive materials requires the time history of the radioactivity of nuclides within a decay chain. The system of differential equations with initial conditions or initial value problem describing radioactive decay of a parent and daughter nuclides was posited by Ernest Rutherford who was awarded a Noble Prize in 1910 for this work. Harry Bateman (1910) provided an analytic solution to the radioactive decay chain problem with the constraint that initial inventory of all daughter elements is zero. Required data for the decay chain calculation consists of the parent and all daughters’ radioactive half-life. The half-life for essentially all radionuclides has been established and is available from multiple sources. Solutions other than Bateman’s (non-zero initial conditions) can be computed analytically but become unwieldy for longer decay chains. For this reason, many applications use a numerical solution. However, a numerical solution can require constraints on the time step size. The proposed method of false rates provides a unique algorithm for the decay chain activities. The method treats the decay chain with arbitrary initial conditions and the calculation is analytic or exact. The method is unexpectedly simplistic. An example decay chain calculation compares the solutions by the method of false rates with a numerical method. The comparison is a verification of the method of false rates calculation. The method of false rates is easily coded as a stand-alone application or as a sub-module of a more general code such as a contaminant transport model.


Introduction
Radioactive materials have applications in many disciplines. A short list includes a) Nuclear medicine b) Reactor fuel for power generation c) Inspection of manufacture metal materials and integrity of welds d) Radiocarbon dating e) Production of nuclear weapon materials f) Environmental remediation of nuclear waste Computer models for radioactive processes are found in the private sector [1][2][3]. In most cases the private sector models are proprietary and require the purchase of a license. Generally, the radioactive decay chain calculation of the private sector codes is only a minor portion of the possible algorithmic solutions of numerous mathematical problems. For the private sector models the source code is not provided and specific numerical algorithms are not disclosed. The government laboratories [4][5][6] provide numerous radioactive decay calculators. Information on the use of these decay calculations is provided, but again the numerical algorithm and source code are not disclosed. Further, most radioactive decay chain calculators are specific to a stated application. The decay chain calculator proposed introduces no model constraints and is unusually simplistic.
The mathematical model proposed by Ernest Rutherford for the decay chain mass is a system of first order linear homogeneous constant coefficient differential equations , 0 where nuclide mass, nuclide decay rate, K=decay chain length.
In Eq. 2 the first term in the derivative expression represents the rate of mass loss to daughter decay, while the second term represents the rate of mass gained from ingrowth of the preceding element in the decay chain. The decay chain ends if the last element in the chain is stable (nonradioactive) or if the application regards further daughters as unimportant to the investigation.
In 1910 Harry Bateman [7] applied the method of Laplace Transform to provide a solution to Rutherford's model. Bateman's solution requires zero initial mass for all daughters in Eq. 2. Bateman's solution for the kth element in the chain is An algebraic solution of Bateman's equations requiring an eigenvalue/eigenvector calculation and matrix inversion is demonstrated in [8]. More recent work has focused on error estimations and a power law approximation for Bateman's solution [9,10]. Much of the analysis for decay chain calculations focus on the Bateman's equations. The method of false rates is not limited by the initial concentration constraint imposed by Bateman's solution.
For nuclear waste management regulatory limits are expressed in terms of nuclide activity, ) = , which represents a mass rate of change. In certain applications, such as the remediation of nuclear waste sites, the initial mass of daughter elements can be greater than zero. The system of differential equations for activity with non-negative initial conditions is These results and a thorough treatment of the physics of radioactivity are provided in [11].

Half-Life and Decay Rate
The half-life of a radionuclide element, denoted /, , is determined by solving Eq. 3 for the time at which the radionuclide activity is equal half the initial activity. The nuclide activity decays according to

Solution for the half-life is
If the half-life is known the decay rate is

Decay of Parent and Single Daughter
Consider the system of differential equations for nuclide activity with non-negative initial conditions, Eqs. 3 and 4. For a parent and single daughter the system is The equation for the parent is separable with solution Substitute ) into Eq. 6 yields The solution for the daughter activity is This analytic approach for the parent and first daughter activity is used to construct a solution for all remaining daughters in the decay chain.

Method of False Rates for Multiple Daughters
Suppose the decay chain has two or more daughters. Assume it is of interest to calculate the nuclides decay time history. Introduce a time discretization, < < , < ⋯ < 0 < 0 , where is initial time and 0 is end-time. If n=1, the decay occurs over a single time step. Assume the solution is advanced to . Then the solution at 5 over time step Δ = 5 − for the parent and first daughter activity follows from Eqs. 7 and 8 ) , 5 = ) Now focus on the solution for the second daughter activity. The parent and first daughter solution at time 5 are computed from Eqs. 9 and 10. Assume that the first daughter's decay masquerades as a parent decay with false rate 7 , . Then ) , 5 ) , !"-7 , Δ .
Solution of Eq. 11 for the false decay rate yields 7 , 8 9: ; With the false decay rate representing the first daughter's decay as a parent, the second daughter decay is modeled by Eq. 10 with nuclide index incremented by one For a third daughter, the second daughter's decay is assumed to decay as a false parent, which provides a false decay rate for the second daughter. This false decay rate is used with nuclide index incremented in Eq. 14 to calculate the third daughter's activity, ) C 5 . This iterative calculation is continued to the end of the decay chain. Subsequently, this method is referenced as the method of false rates.
For many applications the decay history is required. For example, if the decay calculation is included in a contaminant transport model and the contaminant transport is a discrete model, say a finite difference model, the transport model takes relatively small time steps compared to the total run time in order to control truncation error. Consequently, the decay calculation is performed within each time step of the transport calculation.

Example of Decay Chain Calculation
Consider the decay chain involving the following radioactive isotopes: The half-life and initial condition for each decay-chain element is provided in Table 1. Note Ra-228 is radioactive but assume for this example the daughters of Ra-228 are not of interest. The decay chain and half-life data are obtained from [12]. For comparison a numerical solution is calculated with the 4 th order Runge-Kutta (RK) method [13]. The decay chain models (MFR and RK) are coded in Excel ® . Decay is modeled over the time interval [0, 200 yr] with 200-time steps or Δ 1 TF .
Assume the unit for activity is Curie [Ci]. Since the initial activity of U-236 is zero, the false rate for U-236 on the initial time step is calculated as the following variant of Eq. 12 The comparison of the two methods (MFR and RK) at 200 yr is shown in Table 2. The results for Cm-244 show a difference in the 8 th significant digit, while other nuclides agree to 8-significant digits. Figure 1 is the time history of the nuclide activity as calculated by the method of false rate (MFR) and Runge-Kutta (RK). The comparison shows very good agreement.  Consider the qualitative behavior of results in Figure 1. In the 200 yr run time the Cm-244 decays by approximately 11 half-lives, which results in more than three orders of magnitude decay. The Pu-240 decays over 0.3 half-lives. The Pu-240 decay together with the ingrowth from Cm-244 shows a slight decay response (see Table 1). The half-lives of U-236 and Th-232 are several orders of magnitude greater than 200 yrs. The response of U-236 with initial activity zero is due to ingrowth of Pu-240. From Table 1 results, Th-232 show almost no decay. Finally, the decay of the Ra-228 demonstrate an interesting behavior in which the activity increases and is asymptotic to Th-232 activity. For the Ra-228 calculation W XY,,R 0.12/TF and 7 C 7 Z ,>, 2.5\ 11/TF , where 7 C is approximately 10 orders of magnitude less than W . The result for Ra-228 activity from Eq. 13 with incremented nuclide index is The estimate 7 C ≪ W provides the approximation The 2 nd term in the sum decays to zero with rate W . This demonstrates that the activity for Ra-228 is asymptotic to the Th-232 activity. This behavior is called secular equilibrium and occurs when the decay rate of the daughter (Ra-228) is much greater than the decay rate of the daughter's preceding element in the decay chain (Th-232).

Other Decay Modes
The former analysis and example considered a decay chain for which each daughter nuclide has a single predecessor and a single successor. Radioactive decay can also occur in the following modes: A single nuclide can decay to two distinct nuclides (bifurcate).
Two nuclides can decay to the same nuclide (confluent). The decay chain shown in Figure 2 demonstrates both the above cases. In Figure 2 the nuclides are indexed for reference.
The decay of Cm-243 bifurcates to both Am-243 and Pu-239 with stoichiometry coefficients 0.0024 and 0.9976, respectively. The bifurcation from Cm-243 to Am-243 and Pu-239 requires the solutions of the rate equations for Am-243 and Pu-239 include the corresponding stoichiometric coefficient as a scaler multiplier of the Cm-243 decay rate.
The decay from Am-243 to Np-239 requires a false rate for decay of Am-243. The decay of Cm-243 (1) and Np-239 (3) to Pu-239 (4) is confluent. The analysis for the first daughter equation, Eq. 8, does not include this case. The differential equation for activity of Pu-239 (4) analogous to Eq. 4 is If Np-239 (3) is treated as a false parent with false rate, then the solution of Eq. 14 yields a result analogous to Eq. 8 As indicated in Eq. 15, the calculation for Pu-239 (4) activity requires a Np-239 (3) false rate. The calculation order would necessarily agree with the order of indexing.

Conclusion
For decay chain calculations the advantages of the method of false rates is its simplicity and the method of false rates provide an exact solution for any time step. After calculation of the activity for the parent and first daughter, the algorithm requires evaluation of the false decay rate for the first daughter and the application of daughter activity, Eq. 13, for the activity of the 2 nd daughter. This step is repeated for successive daughters. The method takes advantage of the structure of the decay chain system of differential equations and is easily implemented.
For the other decay modes, the impact of the mode on the rate equations and the corresponding solutions are discussed. With the analysis of the treatment of the 1 st daughter solution for the confluence case, the method of false rates is applicable to both the bifurcation and confluence modes.