## 1. Introduction

Nuclear forensics for denuclearization is considered essential for nuclear non-proliferation efforts to prevent nuclear security incidents, which may include illegal trading, theft, loss, and terrorism using nuclear materials. To manage nuclear security threats, reliable technology for estimating the time of nuclear activity is required for nuclear forensics. To verify nuclear non-proliferation, sample analysis and assessment techniques for nuclear materials and byproducts are necessary.

The International Atomic Energy Agency (IAEA) conducts inspections of facilities capable of developing nuclear weapons to prevent and respond to nuclear security incidents in member countries through the Nuclear Nonproliferation Treaty. As part of such activities, nuclear forensic technologies are being developed by the IAEA through international cooperation under the leadership of the US and EU [1, 2]. Advanced nuclear forensics and radio-chronometry technologies include various nuclear tracking technologies, such as quantitative analysis of related substances, isotope analysis, constituent element analysis, and shape analysis using mass spectrometers (US). Research and development have been conducted in the EU on the radio-chronometry of nuclear materials using isotopes of thorium, and the UK possesses technology to predict the production year accurately through radiation spectroscopy and isotope analysis based on mass spectrometry for samples containing plutonium.

Uranium radio-chronometry using nuclear forensic techniques can estimate the time of nuclear activities, such as enrichment and reprocessing of nuclear materials. Such estimation can be performed by identifying the proportions of radionuclide signatures, such as uranium and its daughter nuclides, and analyzing their characteristics. If the characteristics of radionuclides in the uranium found in nuclear materials and byproducts with their corresponding decay chains are accurately evaluated for each radioactive equilibrium ratio, the time of specific nuclear activity can be identified based on each radionuclide ratio. Each type of radionuclide signature and the ratio of each radionuclide in the uranium decay chain exhibit unique characteristics. Moreover, the time of production or enrichment can be estimated from the measured ratio of the radioactive equilibrium relationship between uranium and other radionuclides in the uranium decay chain.

Based on radioactive equilibrium relationships, a nuclear forensic methodology for estimating the time of nuclear activity was investigated. To this end, the radioactive equilibrium relationship, radio-chronometry, and a mathematical model of nuclear materials were analyzed, a related novel algorithm was implemented, and a comparative evaluation was performed. Thus, a radio-chronometry algorithm for nuclear materials or activities was developed as a basic technology for nuclear forensics to strengthen the capability of controlling nuclear materials and bolster global nuclear nonproliferation efforts.

## 2. Methodology and Result

### 2.1 Radioactive Equilibrium Relationships of Nuclear Materials and Radio-Chronometry

Uranium and thorium are representative natural radionuclides with half-lives exceeding 100 million years. They undergo decay chains and become stable through continuous decay reactions; during this process, radioactive equilibrium appears if the half-life of the parent nuclide is longer than that of the daughter nuclide and the decay time is sufficient.

Radio-chronometry for nuclear materials uses the radioactive equilibrium between nuclides in a decay chain. The Bateman equation includes simultaneous differential equations for the branching decay ratio, decay constant, and initial condition of each nuclide, which are the main factors of the decay chain and radioactive equilibrium, and it represents the ratio of nuclides over time when various types of nuclide decay simultaneously. The radioactivity of the parent and daughter nuclides in the radionuclide decay process can be estimated using the Bateman equation.

### 2.2 Mathematical Models for Radio-Chronometry

#### 2.2.1 Review of Radioactivity Calculation Models for Radioactive Decay Chain

Radio-chronometry for nuclear materials in nuclear forensics consists of solving simultaneous differential equations describing the radioactive isotopes that exist in nuclear materials.

Nuclear materials that have a decay chain, such as actinide nuclides in natural radionuclides and spent nuclear fuels, can be modeled using the Bateman equation, which is a linear first-order differential equation that can be expressed as follows [3, 4]:

If the radioactive decay chain is “1^{st} → 2^{nd} → 3^{rd} → … → *i ^{th}* → … →

*n*,,”

^{th}

where, *λ _{i}* is the decay constant of

*i*nuclide

^{th}*b _{i}* is the branching ratio of

*i*− 1

^{th}nuclide decay into

*i*

^{th}If *N*_{1} (0) ≠ 0 and *N _{i}* (0) = 0, and

*i*> 1, the radioactivity of the

*n*nuclide can be expressed as follows [4]:

^{th}

where, ${a}_{i}={\displaystyle \prod _{\begin{array}{c}j=1\\ j\ne i\end{array}}^{n}\frac{{\lambda}_{j}}{{\lambda}_{j}-{\lambda}_{i}},\left(i=1,n\right)}$

Here, B is the branching ratio term, $B={\displaystyle \prod _{k=1}^{n-1}{b}_{k}}$

*b _{k}* is the branching ratio of

*i*− 1

^{th}nuclide decay into

*i*

^{th}According to Cetnar et al., when some nuclides appear in the decay chain more than once, the analytical solution becomes more complicated. In this case, a general analytical solution can be expressed as follows [4]:

where, ${a}_{i}={\displaystyle \prod _{\begin{array}{c}j=1\\ j\ne i\end{array}}^{n}{\left(\frac{{\lambda}_{j}}{{\lambda}_{j}-{\lambda}_{i}}\right)}^{{m}_{j}}}$

Here, B has the same definition as that for Eq. (2).

Moral et al. obtained the numerical analysis solution of the Bateman equation that includes branching decay by using a matrix exponential function as follows [5, 6]:

where (in the case of a homogeneous system)

The solution of Eq. (4) can be expressed as

where, $\left[\Lambda \right]=Diag\left[{e}^{{\Lambda}_{1}t},{e}^{{\Lambda}_{2}t},\cdots \cdots ,{e}^{{\Lambda}_{n}t}\right]$

If one defines

If branching decay is included,

where

The main models to calculate the radioactivity of the nuclides in a radionuclide decay chain are listed in Table 1.

#### 2.2.2 Current Status of Radio-Chronometry Calculation Models

Regarding radio-chronometry, the Lawrence Livermore National Laboratory (LLNL) defined “model age” and “sample age” [10]. The sample age of nuclear materials is defined as the difference between the analysis (reference) date and the production date (or the time of nuclear activity). This can be applied to a closed system with no loss or increase in the radioactive isotopes in the decay chain. The model age can be calculated when there are no daughter nuclides of interest in the decay chain during the production of nuclear materials. The amount of daughter nuclide in the investigated nuclear materials is determined by the amount of parent nuclide and the decay rate. The maximum value of the generated time of the target nuclear material was estimated using the model and sample ages. The Nuclear Forensics International Technical Working Group (ITWG) also assumed a pure parent nuclide and closed system, similar to the guidelines of the LLNL. Based on these assumptions, Eq. (7) can be deduced by summarizing the solution of the Bateman equation for time t:

The radio-chronometry of the materials of interest, such as uranium and plutonium, using the equations proposed in the guidelines of the LLNL and ITWG, were used in monitoring research on the uranium production of nuclear weapons and the estimation of the production time of standard sources for nuclear forensics [13-16]. However, in Eq. (7), the parent nuclide at the time of production (or refinement) must be in a pure state without any daughter nuclide or nuclide of interest, which is only valid in a closed system. Moreover, calculations between nuclides beyond one generation, such as parent and granddaughter nuclides, are impossible. Therefore, it is impossible to obtain analytically the solution for the elapsed time t in nuclear forensics, in which the radioactivity ratio of environmental samples or daughter nuclides in an open system is important. This ratio is generally applied in reality. To overcome the calculation limitations of more than one generation and closed systems, radio-chronometry should be performed by repeated calculation of the Bateman equation or inverse calculation using a radioactivity ratio table or graph [11, 17].

### 2.3 Decay Chain Radioactivity Calculation Method

In this study, to estimate the production time (model age) of nuclear materials, the radioactivity calculation method for radionuclides in a decay chain was implemented first. Since Bateman suggested a solution for calculating radioactivity in a radioactive decay chain, various approaches for calculating the accurate radioactivity of radionuclides in complex decay chains have been proposed for uranium enrichment and nuclear forensics, focusing on nuclear materials and spent nuclear fuels. The calculation of radioactivity for radionuclides in the decay chains of nuclear materials and spent nuclear fuels based on the Bateman equation includes the aforementioned general solution method, transmutation trajectory analysis method, matrix exponential method, and Monte Carlo method.

However, because the conventional radioactivity calculation method is optimized for use, some modifications are required for radio-chronometry. Therefore, in this study, the methods for efficient calculation of radioactivity implemented using a symbolic–numerical solution and matrix exponential function approaches, which have been investigated recently, were reviewed.

#### 2.3.1 Method Using Symbolic–numerical Solution

In this study, to calculate radioactivity more accurately, which includes analyzing the radioactivity of radionuclides in complex decay chains, the algebraic numerical solution was obtained through Eq. (1) using MATLAB software (R2021a Update 5).

First, the uranium decay chain was selected to obtain symbolic solutions for calculating the radioactivity of each nuclide in the uranium decay chain using MATLAB. To solve Eq. (1) efficiently, similar solutions for each nuclide were obtained using the MATLAB built-in “dsolve” function, as shown in Table 2. The “dsolve” function, which is a built-in function of such numerical analysis programming tools as MATLAB and Mathematica, is used to obtain a symbolic solution of ordinary differential equations.

To calculate the symbolic solution numerically, the radioactivity curves for each nuclide were obtained using MATLAB. The radioactivity calculation results for each nuclide developed for the uranium decay chain are summarized in Fig. 1 and Table 3.

The above results reveal that this approach obtained a more accurate value than the other methods, but with a longer computation time. Thus, it is not appropriate for applications because of the long computation time necessary for analyzing composite nuclear materials with complex decay chains and similar spent nuclear fuels.

#### 2.3.2 Method Using Matrix Exponential Function

This study used the matrix exponential function, which is known to calculate the radioactivity of radionuclides in complex decay chains efficiently. Eqs. (5) and (6), obtained based on the above-mentioned matrix exponential function for the Bateman equation, were employed to implement the radioactivity calculation method using MATLAB. The results of the radioactivity calculation for each nuclide applied to the uranium decay chain using the methods proposed by Moral and Pacheco [5] are summarized in Fig. 2 and Table 4.

The method proposed by Moral and Pacheco fails to calculate accurately the radioactivity of the daughter nuclide that undergoes branching decay, as indicated in yellow in Fig. 2. To improve its accuracy, a correction method, such as excluding short-lived nuclides [5], is necessary.

In contrast, the radioactivity calculation method based on the matrix exponential function proposed by Levy can be used to calculate the radioactivity even when branching decay is included. The method proposed by Levy obtains [F], where the lower triangular matrix [AM] of the decay chain differential equation has the relationship [AM] [F] = [F] [AM] with the lower triangular matrix [F] of the decay chain general solution. Matrix [F], defined in Eq. (6), can be obtained using the Schur–Parlett algorithm and is implemented using the MATLAB built-in “funm” function [6, 18]. The radioactivity calculation results for each nuclide applied to the uranium decay chain using this method are shown in Fig. 3 and Table 5.

### 2.4 Radio-Chronometry Algorithm and Influencing Factors

A numerical calculation algorithm for the generated time corresponding to elapsed time t was developed using three algorithms that were implemented to calculate the decay chain radioactivity. To calculate the radioactivity ratio for elapsed time t, the Monte Carlo method was used. It was implemented using the MATLAB built-in “vpasolve” function [19], which is used to solve polynomial or nonpolynomial equations. Because the equation that calculates the radioactivity ratio for each nuclide of interest in this study is in a nonpolynomial form, a random initial value within the search range was chosen using the MATLAB built-in random number generator, and a solution that meets the boundary condition was obtained.

Fig. 4 shows a graph of the function for the radioactivity ratio in a decay chain for 10 years after the generation of uranium nuclear materials. The point marked as “Age dating” indicates the position at time t when the radioactivity ratio of the parent and daughter nuclides is 0.75. The algorithm for the algebraic numerical solution and the Levy method obtained the same results: 0.1319672 years for ^{234}Th and 0.1324292 years for ^{234m}Pa. However, the Moral method led to a different result for ^{234m}Pa, although the same result of 0.1319672 was obtained for ^{234}Th. This could be because the branching decay calculation must be corrected, as mentioned above.

The factors influencing the age dating of the nuclear materials consist of the decay constant of each nuclide, radioactivity measurement uncertainty of the target nuclides of nuclear forensics (parent nuclide, daughter nuclide), and sampling uncertainty of the samples for nuclear forensics, which serve as the main input factors of Eq. (7) [20, 21].

### 2.5 Comparison of Radio-Chronometry Algorithms

#### 2.5.1 Comparison of Radioactivity Calculation Algorithms

A comparative evaluation was performed to verify the accuracy of the proposed algorithm for decay chain nuclides in radio-chronometry. The results of the Chebyshev rational approximation method (CRAM) solver and WISE Uranium decay chain calculator used for the radioactivity calculation of decay chain nuclides, the symbolic–numerical solution algorithm, and the matrix exponential function algorithm based on the Moral and Levy methods are summarized in Table 6. In this study, to compare the implemented algorithms with commonly used algorithms, such as the CRAM solver and WISE Uranium decay chain calculator under equivalent conditions, the radioactivity was calculated for ^{241}Am and its daughter nuclides.

All three algorithms implemented in this study for the radioactivity calculation of nuclides in a radioactive decay chain showed the same values, thereby validating the Moral method as well if there is no branching. The Moral method also showed an approximate value compared with the CRAM and WISE Uranium decay chain calculator algorithms. The computation time was 2.34 s for the moral method, 2.42 s for the Levy method, and 173.05 s for the algebraic numerical solution. Thus, the method using the matrix exponential function yielded the best results. In particular, the algorithm using the algebraic numerical solution consumed dozens of minutes during radioactivity calculation after the seventh generation when branching decay was included.

This requires additional research, such as optimization. Regarding the Moral and Levy methods, the former was better in terms of computation time; however, it requires additional correction when branching decay is included. Therefore, the Levy method was found to be the most efficient in this study.

#### 2.5.2 Comparison of Radio-Chronometry Calculation

To compare the radio-chronometry algorithms, the calculation results obtained using the implemented radiochronometry algorithms with the Levy method based on reference data [24], which measured the nuclide ratios of ^{234}U and ^{230}Th, are shown in Fig. 5. In addition, based on Eq. (7), a summary of the calculation results obtained by LLNL, Los Alamos National Laboratory (LANL, US), Japan Atomic Energy Agency (JAEA), and Commissariat a l’Energie Atomique et aux Energies Alternatives (CEA, France) are presented in Table 7.

The radioactivity ratio was first obtained using the ratios of ^{234}U and ^{230}Th, which were analyzed by each agency and then applied to each calculation. The calculation period was assumed to be 100,000 years. The comparison results showed that the results of the model ages, which were the radio-chronometry results of each agency and the implemented radio-chronometry results using the Levy method in this study, were approximate. In particular, the result of LLNL demonstrated an accuracy for the degree of difference owing to rounding to the third decimal place of the calculation result. Regarding the results of CEA, which only presented the first decimal place, rounding to the second decimal place made it more approximate. Therefore, it can be estimated that most of the differences in Table 7 contribute to the accuracy of the input used for the calculation, as well as the decimal processing of the output value.

## 3. Conclusion

In this study, a novel radio-chronometry algorithm for age dating of nuclear materials was implemented based on the radioactive equilibrium relationship with mathematical modeling as a basic technology in nuclear forensics to strengthen the capability of controlling nuclear materials for denuclearization. Algorithms using a symbolic–numerical solution and matrix exponential function were implemented based on the Bateman equation to calculate the radioactivity of each radionuclide in the decay chain of nuclear materials, such as uranium or spent nuclear fuel.

The proposed algorithms were confirmed to obtain the same level of calculation accuracy as commonly used radioactivity calculation algorithms, such as the CRAM solver and WISE Uranium. However, the symbolic–numerical solution algorithm suffered from a large computation time during the increment of generations of radionuclide. In contrast, the Moral method failed to calculate the radioactivity of the daughter nuclide accurately when branching decay occurred.

To verify the accuracy of the implemented radio-chronometry algorithm, the calculation results were compared with those of the following four agencies: LLNL, LANL, JAEA, and CEA. The results verified that the implemented algorithm could perform chronometry calculations with the same level of performance. However, to improve the utilization of the implemented radio-chronometry algorithm, follow-up studies, such as comparisons with calculation results reflecting complex decay chains and nonhomogeneous conditions, as well as the optimization of calculation algorithms, are required.

The proposed radio-chronometry algorithm can be used as a research and development tool for identifying radionuclide signatures and developing related technologies in nuclear forensics, as well as a tool for comparatively verifying the calculation results of operators for independent verification of spent nuclear fuel management.