Impurity Diffusion as a Possible Metal Chronometer for Pre-Detonation Nuclear Forensics

The ability to determine the age of seized nuclear material—that is, the time that has passed since it was formed— would provide crucial data to be used in its investigation. This paper reviews the methods and mathematical reasoning behind the use of diffusion theory, as previously applied to analysis of metals in ancient artifacts and other objects, to modern investigations in nuclear science. We here examine the time-dependent processes of diffusion, including grain boundary diffusion and discontinuous precipitation, and we assess the utility of examining the profiles of impurity and alloying element concentrations for use as a tool in pre-detonation nuclear forensics.


I. Introduction
Interdiction of metallic nuclear material raises many questions of provenance, only some of which can be satisfactorily addressed today [1-3]. For example, we can analyze for the date of last chemical purification by examining progeny isotopes. This analysis technique develops the approach to answer another key question: when was the specimen cast or formed? We address this question by examining impurity diffusion in the microstructure of the metal.
An enhanced capability in nuclear forensics-one that determines the material's critical parameters, such as its timeline of processing-would support investigations of seized illicit materials. Earlier casework has demonstrated that information important to the investigation can be garnered not just from the nuclear materials themselves, but also from the containers, shields, and assorted other components associated with a seized illicit article [2]. Notably, the first methodology we envision should be extensible for evaluation of the time since forming of any metallic sample, not just potential nuclear device components.
This forensic method of interpretation is not currently developed in the nuclear materials forensics community. In the case of an interdicted sample, all forensics evidence that can be brought to light helps determine possible sources or pathways for the objects in question. Also, the approaches we examine in this article compliment our other work for uranium metals and alloys, providing a cross-check on the results of each, and thus higher confidence in the accuracy of interpretations [4]. In principle, our proposed work is also applicable to metals broadly, including plutonium or other nuclear fuels, including uranium.

II. Background
In the late 1970s, researchers from the Musee d'art et d'historie in Geneva and the Metropolitan Museum of Art in New York City collaborated on a project to investigate changes in the microstructure of Ag-Cu alloys. Their goal was to explain the manner in which normally malleable, ductile silver artifacts become brittle. It has been well documented that ancient silver artifacts are found to be brittle after excavation, whether or not they appear to be corroded [5].
There are three primary causes for embrittlement in ancient silver: corrosion-induced embrittlement, microstructurally-induced embrittlement, and synergistic embrittlement, all of which are discussed in detail elsewhere [6,7]. Schweizer and Meyer investigated discontinuous precipitation (DP) of copper in ancient silver alloys in order to distinguish between forgeries that had precipitated rapidly at high temperatures and copper that had precipitated very slowly over many centuries [5]. They based their conclusions on an observation that the interlamellar distance of a precipitated cell in silver alloys is a function of the temperature at which the precipitation occurred.
In order to apply this principle to verifying the age of a sample, they generated an interlamellar distance versus temperature curve and sought to estimate the rate at which DP occurred at room temperature, in the hopes of approximating the size of the precipitated cells in an ancient alloy that had been aging for several centuries. Schweizer  . The cauldron consists of twelve plates and a bowl, all of 95-98% silver. Chemical analysis has determined copper to be the main alloying (or impurity) element. The cauldron is relatively large, (50 cm in diameter), and due to its high-quality workmanship, it has been the subject of numerous studies.
Wanhill examined the samples using a field emission gun scanning electron microscope (FEG-SEM) and automated electron backscatter diffraction (EBSD) equipment. Samples from the cauldron contained precipitate widths up to 7 µm, well beyond the predicted maximum of 2.1 to 2.2 µm, thus ruling out the use of precipitate widths for authentication.
Furthermore, evidence indicates that the copper precipitate is not lamellar and presents a mottled appearance, eliminating the notion that interlamellar distances could be used to authenticate a sample [9]. A non-uniform growth of DP colonies from different boundaries in the same microstructure is frequently observed in the experimental studies for the following reasons: the physical orientation between the grain boundary (GB) and the precipitate habit planes, and the dynamic properties of grain and interphase boundaries, are dependent on their structure and can vary widely in a polycrystalline aggregate [10]. These investigations of metal in objects formed centuries ago suggest important applications to studies of other metals and contexts-as in the present study of impurity diffusion as a possible metal chronometer for pre-detonation nuclear forensics.

A. Diffusion Fundamentals
Diffusion in solids occurs on microscopic scales through the motion of discrete entities (atoms, molecules, clusters and lattice vacancies). These motions are the result of a statistical distribution of kinetic and potential energy expected among the atoms [11]. The process is thermally activated; relative movement (flux) of atoms or molecules is a response to forces such as gradients in chemical potential or temperature. The process is spontaneous, and must result in a net decrease in free energy [12]. Linear equations are sufficient for relating fluxes on the macroscopic level; however, since each component of a system may be influenced by the gradient in chemical potential of any other component, diffusion in multicomponent systems becomes increasingly complex.

B. Fick's Laws
In the presence of a concentration gradient, impurity atoms will move in the direction of the gradient, from regions of high concentration to regions of low concentration, creating a flux of diffusing atoms. Fick's first law governs this interaction: the flux F per unit area per unit time is proportional to the diffusion coefficient and the one-dimensional concentration gradient: !" (1) where C is the dopant concentration per unit volume, and X is the distance. The diffusion coefficient D depends on the individual type of atom and is strongly temperature dependent: where D 0 is the diffusion coefficient extrapolated to infinite temperature, in units of m 2 /s and E a is the activation energy in eV. The motion of a diffusant is determined by a combination of a material's crystallography, its intermolecular forces, and its defects. Diffusion mechanism is also dependent on factors such as the atomic radius of the diffusing atom compared to the host lattice atoms, and whether or not a defect-assisted process is required.
If an interstitial atom is small enough, when compared to atoms of the host lattice, it can migrate freely from one interstitial site to another. Figure 1 compares bulk diffusion data for several impurity elements in uranium extrapolated to ambient temperatures [13]. Interstitial diffusers have much smaller activation enthalpies than those for self-diffusion.
Unless the experiment attains a steady state, time is also a variable, and a continuity equation must be solved. The time dependence of diffusion t is governed by Fick's second law: This equation considers the time dependence of the diffusion process. However, when D is independent of concentration C, and not a function of position (x), this equation simplifies to: By applying the chain rule to Fick's Second Law (Equation 4), a system of n components in a volume of fixed reference is expressed a follows: Using matrix multiplication, Equation 5 can be rewritten as: or simply as: Where the concentration gradient is defined in Equation 3, and are column vectors and is a matrix of diffusion coefficients. The matrix must always have real, positive eigenvalues, which constrains the values of the elements within the matrix [14]. Often, due to a lack of experimental data, the off-diagonal terms are neglected, and we are forced to assume (at some risk) that they are zero. If n = 2, then Equations 5-7 reduce to the binary diffusion flux, where D 11 is the binary interdiffusion coefficient. The D matrix is usually not symmetrical; however, it can be related to two symmetric matrices, L and G, as follows: where the matrix is the Onsager matrix of kinetic or phenomenological coefficients, and is the thermodynamic matrix. Both matrices have real, positive eigenvalues, and as such, the matrix must meet those conditions.

C. Time-Dependent Diffusion Properties
As we mentioned, scientists have extensively researched diffusion-driven processes such as discontinuous precipitation. Kossowsky determined that, in its simplest form in binary alloys, the linear growth rate of precipitation is obtained from the slope of a line in a plot of the largest cell diameter versus aging time [15]. Malhotra established a linear relationship between time and the thickness of the precipitates [16].
Diffusion in multi-component systems involves the simultaneous flow of more than two components. In these systems, the flux of each component is dependent upon its own concentration and chemical gradient, as well as that of all other diffusing components.
Cellular precipitation does not always occur in multi-component systems such as steel, where the material contains interstitial and substitutional solutes with widely different diffusivities. More importantly, the growth features typically observed during cellular precipitation in binary substitutional systems, such as constant interlamellar spacing, are not present in multicomponent systems [17]. The composition of alloys can impede or promote DP during annealing. Saucedo-Muñoz examined precipitation in steel alloys and found agreement between experimental results and theory that steels with the highest concentration of interstitial solutes (C and N) had the highest kinetics of precipitation, and volume fraction of precipitation, while steel samples with the highest concentration of Mn allowed only carbides to precipitate, as Mn maintains N in the solid solution [17]. In 1991, Kikuchi et al. published a comprehensive review of cellular precipitation in Cr-Ni austenitic steels [18]. Cellular precipitation of Cr 23 C 6 or Cr 2 N often occurs when nitrogen is alloyed to Cr-Ni or Cr-Mn austenitic stainless steel samples. Adding nitrogen to high-chromium austenitic alloys enhances the cellular precipitation of Cr 23 C 6, and when no carbon is used in nitrogen alloyed Cr-Ni and Cr-Mn steels, the dominant precipitate is Cr 2 N. Typically, DP stops in these systems once cells cover 80-90% of the grains, unlike the constant growth rates seen in binary substitutional alloys.
In 2010, Contreras-Piedras et al. examined microstructural evolution and growth kinetics in a magnesium alloy [19]. They sought evidence of GB diffusion and also wanted to evaluate the effects of temperature on cellular spacing. Figure 2 shows the volume fraction of cellular precipitation vs. aging time. The lowest aging temperature yielded the highest volume fraction. The analysis of the plot was done using the Johnson-Mehl-Avrami-Kolmogorov equation, which governs most discontinuous phase transformations: Where X f is the volume fraction transformed, t is time, and k and n are constants whose values depend on the nucleation and growth rate of the transformation product. An n value of 1 corresponds to a boundary, and values obtained experimentally for X f were 1.1, 0.85 and 0.87 for 100 o C, 200 o C and 300 o C, respectively. Cellular spacing increased with aging time, in accordance with the Turnbull theory, and the activation energy indicated a grain boundary diffusion process.
The non-steady-state growth of precipitates can be qualitatively explained by three factors. The first is the transfer of the faster diffusing element, nitrogen. At 1073 K, the volume diffusion coefficient of nitrogen in Cr-Ni austenitic stainless steel is five orders of magnitude larger than that of chromium [18]. Initially, the boundary diffusion of chromium is the rate-controlling factor. However, as the reaction progresses, the growth rate, G (as defined in Equation 10, where s is the interlamellar spacing and D v is the volume diffusion coefficient) becomes larger than the observed growth rate. The difference between the experimentally observed growth rate and the growth rate governed by volume diffusion decreased by an order of magnitude. This is primarily attributed to the interlamellar spacing's inability to increase rapidly enough to match the decreasing migration rate of the moving cell boundary. This suggests that the volume diffusion of chromium may have increased significance in the later stages of cellular precipitation.
The second factor is the deceleration of the moving cell boundary. Nitrogen supersaturation in the untransformed matrix is the driving force for the precipitation of Cr 2 N in Cr-Ni austenitic steel. Once precipitation begins, nitrogen, the faster-diffusing element, flows from the untransformed matrix to the cell via long-range diffusion. This action decreases the amount of nitrogen in the untransformed matrix, which in turn decreases the migration rate of the moving cell boundary. The final factor is the stoppage of cell boundary movement. Under steady-state conditions, where the growth rate is constant, a diffusion zone exists ahead of a moving cell boundary. Its width is defined by the growth rate divided by the volume diffusivity. As the migration rate decreases, the diffusion zone's width increases and is accompanied by a loss in chemical diving force for cellular precipitation. This causes the migration rate to continually decrease, until migration stops altogether.
In 2004, Santhi-Srinivas and Kutumbarao studied DP of a multi-component (Fe-Cr-Mn-N) system [20]. Of note, for temperatures below 700 o C, they observed no DP. And for temperatures above 700 o C, they found that the volume fraction and precipitate cell size increased sigmoidally with time, in agreement with the results shown in Figure 2. Based on their observations and the work of Kikuchi, they proposed a six-step description of the growth mechanism for DP in multicomponent systems. The mechanism outlined in Figure 3 can be described as follows: 1. In the solution-treated sample, the concentration of nitrogen is uniform and is represented by N s in Figure 3(a). 2. After aging, Cr 2 N lamellae form by absorbing nitrogen from the immediate surrounding matrix, which becomes highly depleted of nitrogen. The concentration of nitrogen in the cell matrix (N cm ) decreases with the formation of lamellae while the concentration of nitrogen in the untransformed matrix N um remains the same as N s , as shown in Figure  3(b). 3. This process creates a large concentration gradient between the untransformed matrix and the cell matrix, which results in the flow of nitrogen to the cell matrix. 4. Long-range diffusion of nitrogen from the untransformed matrix is needed for the growth of Cr 2 N lamellae. See Figure 3(c). 5. As aging continues, the concentration of nitrogen in the untransformed matrix (N um ) decreases and nitrogen concentration in the cell (N c ) increases as more nitrogen diffuses to cause growth. As the gradient between the untransformed matrix and the cell continuously decreases, the diffusion rate of nitrogen, and hence the migration rate of the cell boundary, decreases with reaction time. 6. As aging progresses, N um decreases further. Eventually, when chemical equilibrium is reached, N um = N cm , as shown in Figure 3(d), cell growth stops, leading to an incomplete cellular reaction. This six-step process differs from the method proposed by Kikuchi in three ways: long-range diffusion of nitrogen was assumed to take place from the untransformed matrix to the cell; the observed concentration of chromium in the untransformed matrix did not change (emphasis was placed on the volume diffusion of nitrogen and corresponding decrease in the migration rate of the cell boundary), and finally, boundary migration and cell growth were assumed to stop when the cell matrix and the untransformed matrix reached chemical equilibrium.

D. Effective Binary Diffusion Coefficient (EBDC)
Complex diffusion processes in multicomponent alloys are not often solved analytically, and must be dealt with numerically. In order to simplify these calculations, a key assumption is possible to alleviate the need for extensive calculations. The analysis of multi-component diffusion can be simplified by using the concept of effective binary diffusion coefficient (EBDC), treating an alloy as if there were only two components: a solute and a solvent matrix.
The system in Equation 5 required n-1 flux equations in a system of n components. The diagonal terms illustrate the effects of a component's concentration gradient on its flux, and the offdiagonal terms illustrate the ability of the remaining components to influence the flux. The system of equations in Equation 5 can be simplified to: as if there were only two components. D 1 (EB) represents the effective boundary diffusion coefficient of component 1. D 1 (EB) is different for each component, which is not the case in true binary diffusion. If we are to model the diffusion couple as semi-infinite, such that initial concentrations are preserved at sufficiently large distances from the interface (x=0), and that no concentration gradient exists on either side of the couple, we must also assume that the EBDC of a component is independent of distance within the diffusion zone. The solution to the isothermal diffusion equation is as follows [19,21]: where ΔC 0 is the initial difference between the concentrations of the components on each side of the couple, C i (0) is the lower of the two initial values of C i . Calculation of an EBDC requires knowledge of the self-diffusion coefficients of all the diffusing components and the concentration gradients of all but one component. Extensive, current diffusion data can be found in Neumann's Handbook of Experimental Data [13].
Ganguly determined that the true value of Dt, where D is the diffusion coefficient and t is time, is related to the calculated value (Dt) c from the measured concentration profile by the relationship where is the standard deviation of the distribution of intensities of x-rays from the spatial averaging of the microprobe beam.

Influence of Time and Temperature
In 1996 and in several subsequent works, Ganguly sought to use solutions to the diffusion equation to model the time scale of metamorphism in natural garnet diffusion couples, consisting of Ca, Fe, Mg, and Mn [14,21,22]. A diffusion profile clearly exists in Figure 4a. However, the data points for the Mn-profile are too irregular to be of help in the analysis, and the Ca and Fe profiles were used in the determination of Dt. The best fit for the data in Figure 4b using multicomponent diffusion values is given by Dt=7.5x10 -12 cm 2 ; however, it is not required that the Dt values be the same for each best-fit curve.
Diffusion is a function of time due to its dependence on temperature. Non-isothermal diffusion can be reduced to an isothermal diffusion problem (through the period of effective diffusion) by determining the characteristic temperature T Ch as follows: By increasing the number of isothermal steps, one can improve accuracy in calculations. Similarly, one can determine a characteristic diffusion coefficient based on the characteristic temperature in Equation 14 [21,22]. Using the value of Dt=7.5x10 12 cm 2 , and solving Equation 14 for Dt, yields: where Δt= t'-t 0 ,and D i(EB) T (ch) is the effective binary diffusion coefficient at the characteristic  Figure 4a) and modeled concentration profiles with Dt=7.5x10 12 cm 2 (bottom, Figure 4b). [14].
temperature. In order to calculate the EBDC for a component, one must know self-diffusion coefficients for all the diffusing components, and the concentration gradients of n-1 components. Also, as in the case of the Mn profile in Ganguly's geological example [21], one can eliminate diffusing components with poorly defined profiles that would otherwise distort or complicate the solution. Similarly, to average the EBDC of components within the diffusion zone, one can assume a constant D i(EB) and use values from the center of the diffusion zone (about 440 nm in the provided example), as shown in Figure 4b. In order to determine the EBDC for diffusion in a semi-infinite couple, we see that: In order to populate the matrix, Ganguly used Equation 17 for thermodynamically ideal solutions: where ! * is the self-diffusion coefficient of the component i, n is the dependent component, and !" is the Kronecker delta. Ideal values for DFe(EB) and DCa(EB) were calculated from lowtemperature extrapolations. Due to a lack of experimental data, !" * was treated as an unknown, and varied until the values for DFe(EB) and DCa(EB) matched Equation 16. Figure 5 illustrates convergence of ideal EBDCs for Ca and Fe at a Δt of about 4.7 x 10 7 yrs.
In order to adapt this method for the study of metal alloys, the D matrix should be populated by using the Ziebold-Cooper procedure, an extension of Darken's equation of multi-component diffusion as follows [23], in lieu of Equation 17. where: where i is the number of components, D * m is the self-diffusion coefficient, a i is the activity of the i-th component, and δ m is the Kronecker delta function.

A. Sample Calculation
To illustrate the implications of this solution, this paper examines an alloy with known tracer diffusion coefficients (D  [24,25]: Using the Gibbs free energies obtained by Hillert and Waldenstrom [26], and Equation 18, the matrix of thermodynamic coefficients is: For the alloy in question, assuming Cr=1, Ni=2, and Fe=3 (the solvent), the ternary diffusion coefficients are found through the use of Eqn. 22.
Applying the values in Eqns. 20 and 21 to Eqn. 22 yields: For ideal solutions, (where the activity of each atom species is equal to its atom fraction, the atoms species of ideal solutions must be made up of elements that are chemically very similar, if not identical), !" can be replaced with unity [23].
Using the values in the D matrix in Equation 24, and assuming it to be constant, we used the program PROFILER [11] to generate diffusion profiles. The simulation was for 3.6x10 5 s, and exhibited no evidence of uphill diffusion. We assumed the initial concentration differences to be as shown in Figure 6. Given these values for ternary diffusion, one can use an equation similar to Equation 15 and select concentration data from the center of the diffusion zone (about 2100 nm) to ensure an average.
The values were calculated to be: To illustrate the effects of time on isothermal diffusion, Figure 6 shows the effects of varying the aging time at T ch, the characteristic temperature of diffusion, which is based on the experimental data at 1100 o C.

V. Uncertainties
Small errors when calculating T can lead to large errors in calculating D or Δt. An error of ±15°C will cause a significant increase in uncertainty. The activation energy for impurity diffusion in uranium ranges from 52 kJ/mol for Co to 127 kJ/mol for Au [13]. Generally, the higher an atom's activation energy, the more sensitive it is to temperature changes, as shown by the relationship in Equation 26: Figure 7 illustrates the effects of a variance in temperatures and activation energy on calculated D values. The value for D calculated at 300K is 3.75x10 -29 m 2 /s, while at 315K it is 4.24x10 -28 m 2 /s. The greater the activation energy, the steeper the slope of the line and the more significant the variance in temperatures. Another potential source of error is the low-temperature extrapolation of experimental diffusion data according to Arrhenian relations. The lowest experimental temperatures for which data is provided is 683K, nearly 300K greater than ambient temperatures [13]. Whether or not the diffusion mechanism in uranium changes between 300 and 683K is outside the scope of this study, but it bears consideration in future experiments.

VI. Conclusions
The measured properties described in the preceding sections have proven useful in studying the age of geological specimens. Evidence supports the extension of these types of calculations to alloys, to include uranium-based metals. To develop a more reliable method, these calculations must be applied to several impurities/precipitates, such that the estimation of an object's age is determined by solving a system of equations. Also, in the case of the austenitic steels examined, there was no evidence of DP at temperatures less than 700°C. When materials are exposed to ambient temperatures, the nucleation time for DP should be significantly longer than Santhi-Srinivas study [20], and the kinetics sufficiently slow such that the maximum volume transformation occurs not at 10 3 seconds, but over the course of a period more useful to forensic exploration, such as 10s to 100s of years.
We have here presented a method for determining the age (time since casting or forming) of both binary alloys and multicomponent systems. This method involves the diffusion of impurities and alloying elements and the rate of discontinuous precipitation. The ability to validate this method as a useful chronometer depends upon our ability to accurately measure these properties with electron microscopy or similar high-spatial-resolution methods, at much lower temperatures than previously studied.
directs the UT Institute for Nuclear Security. Before joining UT, he spent 20 years at Lawrence Livermore National Laboratory. His research interests include nuclear security applications at the intersection of science, security, and public policy. Professor Hall is a member of the American Nuclear Society, the American Chemical Society, and the Institute of Nuclear Materials Management. He is a Fellow of the American Association for the Advancement of Science and the American Institute of Chemists. Contact: hhall6@utk.edu