Lattice parameter of Am, Np bearing MOX fuel: an empirical potential study

An empirical potential study was performed for the americium (Am), neptunium (Np) containing uranium (U) and Plutonium (Pu) mixed oxides (MOX). The configurational space of a complex U1-y-y′-y″PuyAmy′Npy″O2 system was predicted by the rigid lattice Monte Carlo method. Based on the computing time and efficiency performance, the method was found to rapidly converge towards the optimal configuration. From that configuration, the relaxed lattice parameter of Am, Np bearing MOX fuel was investigated and compared with available literature data. As a result, a linear behaviour of the lattice parameter as a function of Am, Np content was observed.


Introduction
Minor actinide [MA (Am, Np), mainly] bearing mixed oxide (MOX) is of interest for use as fuels for fast reactors (FR) and accelerator driven hybrid systems (ADS) [1]. The effect of MAs on the MOX properties, such as thermal conductivity, melting point and lattice parameter, is of paramount importance to the safe, reliable and sustainable operation. Among these material properties, lattice parameters are essential to obtain the theoretical density, which is used to ensure the quality of nuclear fuel during fabrication [2].
Atomic scale simulation techniques, both density functional theory (DFT) and empirical potentials (EP), are promising to evaluate complex oxide systems' properties. Although DFT-based models have been successfully applied to nuclear fuel, the fist principles description of actinide f-electrons is a highly challenging and debated issue [24][25][26].In contrast, empirical potentials are derived from fitting of the individual material's available properties, which can be used for the evaluation of mixed systems. Therefore, this study aims at investigating the lattice parameter of U, Pu, Am and Np mixed oxides as a function of atomic concentrations via EP methods.
The paper is organised based on the following steps: First, the lattice parameter of the simple oxides AnO 2 (An=U, Pu, Am or Np) was calculated as a pre-validation. Second, binary mixture oxides were modelled for various concentrations of actinides in the host UO 2 matrix. The lattice configuration was predicted via a systematic optimization method called rigid lattice method [27]. From that configuration, relaxed system energies and lattice parameters were investigated as a function of composition. Then, the assessment was extended to complex U 1-y-y′ Pu y Am y′ O 2 and U 1-y-y′-y″ Pu y Am y′ Np y″ O 2 systems. Finally, the results of this study were compared to the available experimental data and as well as to the standard ionic radii model.

Atomic scale simulations
Atomic scale simulations were performed according to the supercell method, using the large-scale atomic/ molecular massively parallel simulator (LAMMPS) version 2020 [28]. A 4×4×4 conventional (nonprimitive)fm3m cell under periodic boundary conditions was considered as a starting structure. The lattice parameterwas set to the experimental value of the host matrix before relaxation (see table 2). The forces acting on each ion were obtained by the Busing-Ida type [29,30] empirical potential fields by Arima et al [31][32][33]. We will further denote this as the Arima potential set.
The Arima potential set was developed to reproduce actinide oxide system properties, in particular lattice parameter, bulk modulus and thermal conductivity [31][32][33]; it has been extensively tested in the literature and was therefore selected for this study [31][32][33][34]. The consistent potential set covers the U, Pu, Am, Np and O ion descriptions, which could be used to define pairwise interactions in solids. The potential energy ( ) V r , ij between two ions i and j at a distance apart r ij is expressed by the following form: The parameters a , i b , i c i appearing in the interatomic potential function are given in table 1. The first term of equation (1) simulates the electrostatic interaction between point charges with  , 0 the vacuum permittivity. The electric charge on each atom is assumed to be 67.5% (=x) of the formal charge (Z i ). Classical Ewald summation techniques were used to handle long-range electrostatic interactions in three dimensions [35]. The second term with f , o adjustable parameter expresses the short-range repulsion of overlapping electronic shells, and the third term expresses the van der Waals attraction as a simple inverse power function. In this study the pairwise interactions were evaluated only up to a cut-off distance of 1.1 nm.
In the stoichiometric simple oxide AnO 2 (An=U, Pu, Am or Np), a stable An 4+ is generally considered to maintain the system charge neutral [20][21][22]. In the present study, a perfect system was created, meaning that the oxygen to metal ratio will always remain 2 and the actinide atoms are always sitting in the unit cell of the fluoritelike cation sites. To describe the actinide charge state in a stoichiometric complex mixed oxide system, a common value of (4+) was considered, as confirmed by previous research conducted by Kato and Konichi [23]. However, combinations of other oxidation states e.g. (3+, 5+, 6+), which would also correspond to overall neutrality, were not applied for this study.
The lattice paramteres were investigated via energy minimization techniques at 0 K. The atomic coordinates and lattice vectors were relaxed until the zero strain with the conjugate gradient algorithm implemented in LAMMPS [28]. There are two advantages for the calculation of system properties at 0 K for EP studies. The first one is to check whether the empirical potential is correctly implemented as a pre-validation. Second, a fundamental understanding of the system behaviour can be gained by ignoring the atomic vibrations. To addres the relationship between lattice parameter of Am, Np bearing MOX fuel as a function of atomic fractions (at%), a static approach at 0 K was followed in this study.
However, for the experimental validation 300 K were used to simulate room temperature conditions, while neglecting standard atmospheric pressure of 0.1 MPa. The equation of motion was integrated with velocityverlet algorithm with the unit time step (Δt) of 2 femto-second. The simulations were carried out under NPT (constant number of ions, temperature and pressure) conditions which were controlled by the canonical dissipative Nosé-Hoover thermostat and barostat with relaxation time of 0.1 ps and 0.5 ps, respectively. The time for the system to reach equilibrium was 10 ps and the measured lattice parameter were averaged over the next 10 ps.

Exploration of configuration space
The configurational space for the U 1-y-y′-y″ Pu y Am y′ Np y″ O 2 system was explored in this empirical potential study, also considered as the lattice configuration of a solid solution. The exploration of such an arrangement is a Table 1. Parameters for Busing-Ida type interatomic potential function, as parametrized by Arima et al [31][32][33].
N n N n m m with m=(y′×256), the number of Am atoms distributed in the simulation box, etc. In order to optimize the calculation time, the rigid lattice method was tested [27].
In the rigid lattice approach, the non-relaxed and the relaxed configuration maps were demostrated to be qualitatively similar. To address the energy gain issue due to the permutation of the U, Pu, Am or Np ions the non-relaxed configurations were used. Initially, a host UO 2 lattice was created according to the potential predictions of the atomic coordinates (cf table 3). The Pu, Am or Np ions were randomly distributed as substitutes in the uranium sublattice. The optimal configuration was found then, by coupling the simulation system with a Monte Carlo-Metropolis algorithm, which dictates the swap probabilities to be always in the direction of energy loss [28].
The method was validated for the case of U 0.5 Pu 0.5 O 2 ; ahomogeneous mixture was observed, as illustrated by figure 1; no clustering or separate UO 2 -PuO 2 phases were predicted, which is in line with recent experiments [21]. The systematic configuration space exploration confirms that the minimum energy configuration is indeed reached. The convergence to the energy minimum was observed after about 20,000 montecarlo steps. Based on this computing time gain and efficiency performance, the method is expected to facilitate the search for optimum configurations when more complex oxides such as U 1-y-y′-y″ Pu y Am y′ Np y″ O 2 are at play.

Results and discussions
The lattice parameters for the solid solution of U 1-y-y′-y″ Pu y Am y′ Np y″ O 2 system in the range of y60 at% and y′ (= y″)5 at% were investigated by this empirical potential study. In a first step, the single oxide simulations were performed as a pre-validation to check whether the selected potentials were used properly. Then, the assesment was extended to the binary and more complex oxides. The resulting lattice parameters from this study and the ionic radii model were eventually compared with the available experimental data.
The standard ionic radii model for the lattice parameter (a 0 ) of the f.c.c. cell is expressed as; In equation (2) the body diagonal of the unit-cell is 3 a 0 and one-quarter of the body diagonal is the r c and r a distance. Where r c and r a are the ionic radius of cation and anion, respectively. In order to express the lattice parameter of the U 1-y-y′-y″ Pu y Am y′ Np y″ O 2 system based on ionic radii, equation (3) was derived from equation (2).  It is well known that the ionic radius varies significantly with valence and coordination number (CN) of the ions. In case of a stoichiometric f.c.c. fluorite structure, CNs are 8 for cations and 4 for anions. In this study, ionic radii were taken from the Shannon [38] and Kato and Konashi [23], which are shown in table 2. Although Shannon's set is most widely used, it was assembled in a rather generic way to reproduce the average interatomic distances in solids [38], which was revised and adjusted by Kato and Konashi for a specific AnO 2 (An=U, Pu, Am or Np) system [23]. Table 3 shows the lattice parameters for the simple oxides which were predicted by the consistent empirical potential set of Arima et al [31][32][33]. For UO 2 , the calculated lattice parameter of 5.454 Å agrees with the critical assessment of the interatomic potential for nuclear fuel published by Bertolus et al [39]. This value is in good agreement with 5.455 Å, the 0 K extrapolation of experimental data [40]. At 300 K, both calculations and experiments prove to be in good agreement for all the actinide oxides, as shown in table 3.

Simple oxides
The lattice parameters of AnO 2 (An=U, Pu, Am or Np) were calculated with the standard radii model. Equation (2) was used together with the ionic radii provided by Shannon [38] and Kato and Konashi [23] (cf  table 3). The ionic radii provided by Shannon [38] generally overestimates the actinide oxide system lattice parameter, while Kato and Konashi [23] reproduces good results in the range of±0.001 Å. Therefore, Kato and Konashi data were selected for further investigation [23].

Binary oxide system
Binary oxide U 1-y Pu y O 2 , U 1-y Am y O 2 , U 1-y Np y O 2 lattice configurations were predicted by the rigid lattice method [27]. From that configuration, system lattice parameters were calculated at 0 K. Figure 2(a) shows that the linear relation (also known as Vegard's law) between lattice parameter and content y, for the binary oxide system was observed. To increase the reliability of the current research, the excess energy (or mixing enthalpy at 0 K) was investigated additionally according to the following expression: Where + E A B is the energy of A and B binary solid solution, E A and E B are the standard lattice energies of component A and B, y A and y B are the mole fractions of component A and B, respectively. The whole range of calculated DE excess over the composition is presented in figure 2(b) The calculated excess energies for U 1-y Pu y O 2 , U 1-y Am y O 2 and U 1-y Np y O 2 fit well with a polynomial shape that corresponds to the homogeneous cubic adjustment of the mixed oxide system as a function of composition. Furthermore, the observed trend in DE mix is considered to follow the same trends as the standard ionic radii as tabulated in table 2. In a fixed host matrix  Table 3. Lattice parameters predicted by the consistent potential set of Arima et al [31][32][33]. Equation (2) was used to calculate the lattice parameter of actinide oxides together with the ionic radii provided by Shannon [38] and Kato & Konashi [23]. The two models' results were compared with available experimental data at 300 K.
Lattice parameters (Å) 3.3. Complex U 1-y-y′-y″ Pu y Am y′ Np y″ O 2 system Finally, the lattice parameters were compared with available experimental data for a more complex U 1-y-y′-y″ Pu y Am y′ Np y″ O 2 system. Minor actinide fractions were adjusted up to a share of 5 at% while Pu in the host UO 2 matrix was kept in the range of 5-60 at%. Figures 3(a) and (b) illustrate the effect of minor actinides on the lattice parameter of MOX fuel matrix for a complex U 1-y-y′ Pu y Am y′ O2 (a) and U 1-y-y′-y″ Pu y Am y′ Np y″ O 2 system (b) as predicted by Arima et al at 0 K [31][32][33]. From both figures, it can be seen that the lattice parameter decreases indeed linearly with an increasing amounts of MA dissolved in the MOX fuel matrix. However, for a fixed Pu concentration, changing the Am (y′) -Np (y″) fraction, with y′ = y″ (cf figure 3(b)), results in a sharper decrease than changing the share of Am (y′) in the MOX fuel. This is in line with the experimental work published by Kato and Konashi [23]. Figure 3(c) shows the lattice parameter comparison of this study's EP results to the experimental data [21,23], and the results of ionic radii model as calculated with equation (3). A good agreement between measured and calculated data suggests that the lattice parameter of U 1-y-y′-y″ Pu y Am y′ Np y″ O 2 obeys Vegard's law.

Conclusions
The lattice parameters of up to 5 at% MA containing MOX fuel with maximum 60 at% Pu were addressed by empirical potential based atomic scale simulations. The lattice configurations were predicted by a systematic optimization method called rigid lattice method [27]. The method could be easily applied to other selected Buckingham variants, used to address the ionic oxide crystal [44]. However, one of the advantages of the Arima potential set [31][32][33] is, that it can berelatively easily extended to other oxidation states appearing in nonstoichiometry, e.g. (3+, 5+, 6+), which will be explored in our future works.
The following results were found: • The configurational space exploration of the arrangement of Pu, Am, Np in the host UO 2 matrix was shown to be highly effective with a rigid lattice approach-where no relaxation is performed in a first stage-in combination with a Monte Carlo-Metropolis particle swapping method.
• The lattice parameter of the binary oxide system U  • The lattice parameter of a complex U 1-y-y′ Pu y Am y′ O 2 b) U 1-y-y′-y″ Pu y Am y' Np y' O 2 system decreases linearly with an increasing amount of MA dissolved in the MOX fuel matrix.