domain: chim phys

text:
C Takoukam-Takoundjou

E Bourasseau

M J D Rushton

V Lachet

Optimization of a new interatomic potential to investigate the thermodynamic properties of hypo-stoichiometric mixed oxide fuel U1-yPuyO2-x

The behaviour of stoichiometric U1-yPuyO2 compounds used as nuclear fuel is relatively well understood. Conversely, the effects of stoichiometry deviation on fuel performance and fuel stability are intricate and poorly studied. In order to investigate what affect these have on the thermophysical properties of hypo-stoichiometric U1-yPuyO2-x mixed oxide fuel, new interaction parameters based on the many-body CRG (Cooper-Rushton-Grimes) potential formalism were optimized. The new potential has been fitted to match experimental lattice parameters of U0.70Pu0.30O1.99 (O/M = 1.99) and U0.70Pu0.30O1.97 (O/M = 1.97), where M represents the total amount of metallic cations, through a rigorous procedure combining classical molecular dynamic and classical molecular Monte Carlo simulation methods. This new potential provides an excellent description of the U1-yPuyO2-x system. Concerning lattice parameter, although fitted on only one Pu content (30%) and two stoichiometries (1.99 and 1.97), our potential allows good predictions compared to available experimental results as well as to available recommendations in wide ranges of O/M ratio, Pu content and temperature. For the U0.70Pu0.30O2-x hypo-stoichiometric system (30% Pu content and O/M ratio ranging from 1.94 to 2.00), some direct properties (lattice parameter and enthalpy) and some derivative properties (linear thermal expansion coefficient and specific heat) were systematically investigated from room temperature up to the expected melting temperatures and a good agreement with experiments is found. Moreover, our potential shows good transferability to the plutonium sesquioxide Pu2O3 system.

Introduction

Uranium-plutonium mixed oxides (MOX) with a plutonium content below 9 wt% are currently used as nuclear fuels in Pressurized Water Reactors in France. Within this range, U1-yPuyO2 compounds have been extensively studied, both experimentally and theoretically. The temperature dependence of the physical properties of MOX below this plutonium content are well known [1]- [3]. Homogeneous MOX fuels with higher amounts of plutonium (Pu content more than 20%) are also considered as reference fuels for IV th generation Sodium cooled Fast neutron Reactor (SFR). However, the properties of this homogeneous MOX fuel are poorly understood especially for Pu contents between 25-45% where a miscibility gap was experimentally observed [3]- [5]. Therefore, further research is essential to better understand phase equilibria as a function of temperature and MOX composition since they govern both the manufacturing process and the behaviour under irradiation. By understanding and predicting these phase equilibria over an important range of conditions, in particular a large temperature domain, a deep knowledge of the fuel behaviour can be achieved.

The evolution of the nuclear fuel under irradiation or in storage conditions is affected by many phenomena that cannot always be isolated experimentally. Computer simulations are a way to avoid such complications, as separated processes can be simulated. However, the models and theories upon which they are based could affect the accuracy of simulations. Atomistic simulation techniques have been used for many years to predict the elementary mechanisms that govern the observable behaviour of mixed oxides using both empirical interatomic potentials and ab-initio calculations [6]- [9]. The use of empirical potentials implies that the interatomic forces are computed in the framework of classical mechanics whilst for ab-initio methods these are calculated using quantum physics. The complex chemistry associated with these mixed oxide systems has meant that, ideally, the interatomic forces governing their properties would be described using quantum mechanical methods. However, due to their significant computational cost, the system sizes that can be considered in these calculations are limited to several tens of atoms. For this reason, methods employing classical potentials are often used to allow the calculation of properties on larger length and time scales.

Interatomic potentials are most often used in two complementary simulation methods: molecular dynamics (MD) and molecular Monte Carlo (MC). During a MD simulation, the time evolution of atomic trajectories is predicted by solving Newton's equations of motion, allowing information related to the temporal evolution of the properties of the system to be obtained. By comparison, the MC method which is based on the Metropolis-Hastings algorithm [10], consists of performing random atom moves for a system of atoms, to build a statistical ensemble of configurations from which averages can be made. The MC method seems to be particularly relevant for mixed oxide fuel investigation since it can enable the exchange of uranium and plutonium atoms during a simulation, which allows exploring the space of the cationic configurations in order to achieve the real thermodynamic equilibrium. Such exchanges, whilst possible in MD, are very unlikely to occur within the very limited simulation time accessible through most MD simulations.

Interatomic potentials commonly used to model interactions involved in nuclear oxide fuels fall into the general class of Born-Meyer-Huggins (BMH) potentials used to describe ionic crystals. This class of potentials typically consists in three contributions: Coulomb interactions, core repulsions, and van der Waals dispersion [11], [12]. In the literature, several authors have proposed potentials based on the Buckingham form [13], [14] whilst some others improve it with the Morse form [15], [16] which is applied only to cation-anion pairs. Comparing several such pair potentials, Govers et al. [17], [18] have shown that while a specific potential may well reproduce a subset of thermomechanical properties, none are able to reproduce an extensive set of properties. Moreover, these simple pairwise interatomic potentials failed to reproduce the Cauchy violation observed in actinide oxides with the fluorite structure (i.e. C12≠C44 where Cij are elastic tensor elements in Voigt notation).

Recently Cooper, Rushton and Grimes (CRG) developed a different potential model, which combines Buckingham and Morse forms added to the Embedded Atom Method (EAM) [19] to take into account many-body effects. Thus, the CRG model successfully captures the Cauchy violation, which is necessary to reproduce accurately the elastic constants and consequently, the mechanical behaviour of fluorite oxides. This robust potential provides an excellent description of the interatomic interactions involved in stoichiometric MOX (where the oxygen to metal ratio is equal to 2) and in other fluorite oxides [20]- [22]. It is capable of accurately predicting a large number of the thermophysical properties of stoichiometric MOX and fluorite oxides systems from room temperature up to the melting point.

Since MOX are modelled as ionic compounds, the atomic charge representation is particularly important. In simulations using BMH potentials, charged species are typically represented as constant point charges, as atoms in constant oxidation states. Investigation of stoichiometric U1-yPuyO2 MOX considers the oxidation state +4 for the uranium and plutonium cations. This oxidation state of the cations is not established for the case of hypo-stoichiometric MOX and different cation oxidation states are used in the literature. For instance, Arima et al. [11] have investigated hypo-stoichiometric U1-yPuyO2-x mixed oxide using a partially ionic BMH potential, maintaining the charge neutrality through the appropriate charge balance for the number of Pu 3+ and Pu 4+ ions 1 . In their study, Arima et al. [11] represented each metal oxidation state with different interatomic potential parameters, based on the analysis of different phases of plutonium oxides. Similar approaches have also been applied including specific M 3+ and M 5+ ions in the system in order to maintain charge neutrality to simulate other offstoichiometric nuclear fuels [23]- [25].

Up to now, the Arima et al. [11] potential is the only one available in the literature to study hypostoichiometric U1-yPuyO2-x MOX fuel. This potential was first optimized at room temperature by adjusting the Pu 3+ cation potential parameters against the experimental Pu2O3 cubic bixbyite structure (𝐼𝑎3 ̅ ) at room temperature. However, according to the Pu-O phase diagram, this structure is not stable at room temperature [3], [26]. In order to obtain a more reliable potential above room temperature, Arima et al. assumed that UO2 and PuO2 with the fluorite type structure (𝐹𝑚3 ̅ 𝑚) can easily form solid solution with M2O3 (where M represents rare earth: Er, Gd, Nd or Y). The effect of temperature was then introduced into the potential by optimizing interaction parameters on thermal expansion and compressibility data based on those of M2O3 since these physical properties are unknown for the cubic Pu2O3 phase.

The CRG potential developed by Cooper et al. gradually became a reference potential for the simulation of actinide oxides and their mixed oxides in the last years. Rushton et al. [27] have extended the CRG potential to include interactions concerning several trivalent rare earth species (Yb, Er, Ho, Dy, Gd, Sm, Nd, and La) typically used as dopants within CeO2. Their potential model allows a good description of the thermal expansion and elastic properties of ceria with an excellent agreement with experimental data. The CRG potential has also been applied to hyper-stoichiometric UO2+x, through the addition of U 5+ interaction terms [28]. Combining classical MD and DFT force matching approach, potentials for gas-oxygen and gas-actinide interactions that are consistent with the many-body CRG potential model set have also been developed in order to investigate the behaviour of fission gases in mixed oxides [29]. Nevertheless, up to now, no potential parameters are available for hypo-stoichiometric systems in the CRG framework. A recent study by Maxwell et al. [30] on hypo/hyper-stoichiometric systems used the CRG potential, together with an approach inspired by the charge optimization-based potentials from Li et al. [31] in which the atomic charge of metal ions around point defects is rescaled in order to maintain charge neutrality. However, this method does not allow the displacement of defects; rather, this method likely pins the defects to the points in the lattice wherein the charges are perturbed.

In the present work, we combined classical MD and MC simulation methods to optimize additional parameters in the CRG potential framework to include Pu 3+ interactions. Thus, we build a rigorous empirical fitting procedure to derive the interactions between the trivalent plutonium and the other ions in the system, keeping unchanged the other interactions already included in the original CRG potential.

The structure of this paper is as follows: a description of the methods employed in this work and of the fitting procedure is presented in sections 2.1 to 2.4. The optimization results for the Pu 3+ interactions compatible with the original CRG potential are presented in section 3.1. Then, the validation of our new potential is performed and its transferability to Pu2O3 system is discussed in section 3.2. In sections 3.3 and 3.4 we investigated the effect of the O/M ratio (where M represents the total amount of metallic cations, i.e. U 4+ , Pu 4+ and Pu 3+ ) on the thermophysical properties of hypo-stoichiometric U0.70Pu0.30O2-x MOX fuel within a wide temperature range.

Method

Potential model

The empirical potential model used to describe interactions between atoms in this work is the CRG potential model [20], which provides an excellent description of both UO2 and PuO2 oxides as well as uranium plutonium mixed oxide (MOX). The CRG potential has been detailed in previous works [20], [22], [32]; thus only the essential features will be recalled here. This model combines a pair potential description with the many-body Embedded Atom Model (EAM) [19] from Daw and Baskes. The interaction potential energy 𝐸 𝑖 between an atom 𝑖 with respect to all other atoms is given by:

𝑬 𝒊 = 𝟏 𝟐 ∑ 𝑼 𝜶𝜷 (𝒓 𝒊𝒋 ) 𝒋 -𝑮 𝜶 (∑ 𝝈 𝜷 (𝒓 𝒊𝒋 ) 𝒋 ) 𝟏 𝟐 Eq. 2.1
This equation has two distinct components. 𝑈 𝛼𝛽 (𝑟 𝑖𝑗 ) represents the pairwise interaction between two atoms 𝑖 and 𝑗 separated by a distance 𝑟 𝑖𝑗 . 𝛼 and 𝛽 are used to label the species of atom 𝑖 and 𝑗 respectively. The pairwise interaction term 𝑈 𝛼𝛽 (𝑟 𝑖𝑗 ) is described using both long-range (Coulomb, 𝑈 𝐶 (𝑟 𝑖𝑗 )) and short-range (Buckingham, 𝑈 𝐵 (𝑟 𝑖𝑗 ) and Morse, 𝑈 𝑀 (𝑟 𝑖𝑗 )) contributions given by the equations (Eq. 2.2 to Eq. 2.4) below. The Coulomb and Buckingham contributions added together give the well-established "rigid ion potential" 𝑈 𝑅𝐼𝑃 (𝑟 𝑖𝑗 ). The Ewald summations were used as long-range electrostatic correction. Concerning short-range interactions, a cut off equal to 11 Ȧ has been used as proposed in the original CRG potential. In this equation, a short-range cut-off as an error function is applied at 1.5 Å that reduces the EAM component gradually to prevent unrealistic forces occurring at short separations; ensuring that there is no discontinuity in the interatomic energy, which would arise from an abrupt cut-off as pointed out by Cooper et al. [32].

Simulation methods

Molecular Monte Carlo calculations

Monte Carlo (MC) simulations with cation exchanges were carried out within the isobaric-isothermal (NPT) ensemble: both the atomic coordinates and the cell dimensions were allowed to vary during the simulation. The cation exchange moves allow random exchange of uranium and plutonium atoms. This move, specific to the MC method, was applied only between two cations of different types in order to account for the effect of the cationic distribution [22], [33]. The temperature effect is included by allowing random translations of randomly selected atoms, and the pressure effect was sampled using random volume changes of the simulation box. The attempt probabilities attributed to these three MC moves were 79 %, 1 %, and 20 % for translations, volume changes and exchanges respectively. A move attempt was accepted or rejected following the standard Metropolis-Hastings scheme based on the potential energy difference between the old and the new configurations [10]. The maximum variations in the atomic displacements and in the volume changes were adjusted automatically during the simulation to maintain an acceptance/rejection ratio of approximately 0.5. Simulations were performed with the MC GIBBS package [START_REF] Ungerer | Applications of Molecular Simulation in the Oil and Gas Industry: Monte Carlo Methods[END_REF], [START_REF] Gibbs | [END_REF] from room temperature to the expected melting temperature; O/M ratios varying from 2.00 to 1.94 were investigated. Initial configurations for a given stoichiometry were built using a particular method described in the first step of our fitting process (see section 2.4) and the values of each thermodynamic property were averaged over 450 million configurations after a 50 million step equilibration run. Periodic boundary conditions were used to avoid limited size effects.

Molecular dynamics calculations

Molecular dynamics calculations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF] whereby Newton's laws of motion are numerically integrated to predict atom positions and velocities as a function of time using the forces obtained from the potential model described in section 2.1. This simulation technique was mainly used in stages ii) and iii) of the fitting procedure (see section 2.4). Concerning the fitting process, the MD calculations were carried out under standard constant pressure-temperature (NPT) ensemble at zero external pressure from 300 K to 1900 K using Berendsen thermostat and barostat with relaxation times of 0.1 ps and 0.5 ps respectively. The fitting process implies mainly the calculation of lattice parameters, MD simulations were carried out on supercells of U1-yPuyO2-x hypo-stoichiometries MOX with a fixed timestep of 2 fs (these supercells were built by performing a 6×6×6 replication of a 12 atom cubic cell of the MO2 fluorite structure). These large supercells were chosen in order to avoid system size effects. We used the same initial configurations for MD and MC calculations. Each MD simulation was carried out for 15 ps at a given temperature before being averaged over 5 ps, thereby ensuring enough time for the system to reach equilibrium. Again, periodic boundary conditions were employed throughout, allowing a continuous bulk crystal to be simulated. Moreover, we used the PPPM (Particle-Particle, Particle-Mesh) method available in the LAMMPS package to speed up the Ewald summations.

Choice of the optimization database

The optimization of the original CRG potential was based on experimental elastic constants at room temperature as well as lattice parameters at different temperatures of the pure oxide systems. In the literature, we did not find any experimental measurement showing the effect of stoichiometry on the elastic constants of MOX. We did however find several papers giving experimental lattice parameter data and associated recommendations for hypo-stoichiometric MOX at room temperature. The most reliable references in the literature are from Schmitz et al. [38], Philipponneau [START_REF] Philipponneau | Paramètre cristallin et homogénéité dans le combustible MOX[END_REF], Harding et al. [START_REF] Harding | thermophysical and thermochemical properties of fast reactor materials[END_REF], Benedict et al. [START_REF] Benedict | [END_REF], Omichi et al. [42] and Kato et al. [43]. Belin et al. [44] have measured lattice parameters of hypo-stoichiometric MOX from 298 K to 1750 K using in situ high temperature X-ray diffraction (XRD) for Pu contents of 0.14, 0.24, 0.35, 0.46, 0.54 and 0.62 but the O/M ratio of their samples during the measurements are unknown. These authors highlight that no method is available to determine the O/M ratio during in situ high temperature-XRD experiments. Moreover, back at room temperature, the mass of samples retrieved is also too low to allow for O/M determination using traditional gravimetric methods. To our knowledge, it is very difficult to control the variations of the O/M ratio during experiments at high temperature. The effect of stoichiometry on thermal expansion of MOX has been investigated using thermal dilatometry measurements in a recent study from Kato et al. [43] for Pu contents of 0.30 and 0.48, in a temperature range from 1400 K to 1950 K and for three O/M ratios (1.99, 1.98 and 1.97). Kato et al. [43] have controlled oxygen partial pressure by adjusting the hydrogen partial pressure to moisture partial pressure ratio in a Ar/H2/H2O gas mixture to hold a constant oxygen to metal ratio in the U1-yPuyO2-x during the measurement. The O/M ratio in their study was stated to be stable during the measurement, and uncertainty was estimated to be ± 0.004 for the variation of the O/M ratio during measurements. Based on these measurements, Kato et al. have derived thermal expansion correlations, which take into account the effect of O/M ratio, Pu content and temperature on thermal expansion of hypo-stoichiometric U1-yPuyO2-x.

From the thermal expansion correlations of Kato et al. [43], we determined the lattice parameters as a function of temperature for three O/M ratios (1.99, 1.98 and 1.97) used as reference data, the lattice parameters at room temperature being provided by the same authors in a previous work [45]. These lattice parameters were used in our fitting procedure and are summarized in the first part of Table 1. Note that only plutonium content of 30% and stoichiometries of 1.99 and 1.97 were used within our fitting process.

Fitting procedure

Parameters to be optimized

In the present work, the CRG potential is extended to include interactions involving Pu 3+ . The parameters needed to take into account these interactions, and as a consequence to investigate hypostoichiometric U1-yPuyO2-x MOX fuel, are not available in the original version of the CRG potential. Therefore, we build a fitting procedure to derive the interactions between the trivalent plutonium and the other ions in the system. The goal of the fitting process is to minimize an objective (or merit) function measuring the discrepancies between the reference data presented above and the calculations made by the new candidate potential.

We chose to keep the existing interactions from the original CRG model unchanged (interactions concerning Pu 4+ , U 4+ and O 2-ions). This is particularly important because it maintains the consistency between the new potential and the original CRG potential. Thus, only the interactions involving Pu 3+ ions were optimized i.e. Pu 3+ -Pu 3+ , Pu 3+ -Pu 4+ , Pu 3+ -U 4+ and Pu 3+ -O 2-interactions. As starting parameters for Pu 3+ interactions, we adopted the parameters given for Pu 4+ interactions.

Moreover, we choose to exclude some Pu 3+ parameters from the optimisation. In previous works [20], [28] concerning the extension of the CRG potential to other elements, parameters A and C of the Buckingham contribution of the cation-cation interactions were fixed at 18600 eV and 0 eV.Å 6 respectively for all cations. So we kept these values for the corresponding Pu 3+ interactions. The C parameter of Pu 3+ -O 2-interaction is also fixed to zero. Likewise, the covalent Morse term of cationcation interactions was set to zero in these previous works, so we also set it to zero for Pu 3+ . Concerning the EAM interaction parameters [Pu 3+ ] and G[Pu 3+ ] we chose to kept them identical to those of Pu 4+ (3980.058 Å 5 and 2.168 eV.Å 1.5 respectively). In addition, in order to maintain electro neutrality of the system, the effective charge of Pu 3+ was set to 𝑞 𝑃𝑢 3+ = -3 2 𝑞 𝑂 2-= 1.6656 with 𝑞 𝑂 2-= -1.1104 being the effective charge of O 2-ions. This approach considerably reduces the number of parameters to adjust and therefore simplifies the problem and its numerical cost. The remaining parameters to be adjusted are the three cation-cation ρ interaction parameters (Pu 3+ -Pu 3+ , Pu 3+ -Pu 4+ , Pu 3+ -U 4+ ) and the five components A, ρ, D, ɣ and r0 of the cation-anion interaction (Pu 3+ -O 2-).

Fitting procedure of the parameters

Figure 1 is an illustration of the fitting procedure we have used. The goal of this procedure is to minimize a merit function and to determine the parameters of the potential which give the best agreement in regard with the reference data set presented above. Our fitting procedure can be separated in three main steps: i)

First we built a supercell containing 2592 atoms (864 cations and 1728 anions) by duplicating the Special Quasi-random Structure (SQS) [46] supercells of 96 atom in the three spatial directions. This supercell correspond to the 666 replication of a 12 atom cubic system of the MO2 fluorite structure. The arrangements of U and Pu cations from von Pezold et al. [47] were used for the corresponding Pu content. Then we relaxed the 666 U1-yPuyO2 bulk supercells during a 60 million MC step equilibration run at 300 K within the NPT ensemble with the Monte Carlo GIBBS package [START_REF] Ungerer | Applications of Molecular Simulation in the Oil and Gas Industry: Monte Carlo Methods[END_REF], [START_REF] Gibbs | [END_REF]. This run is performed using the cation exchange move in order to ensure a correct distribution of U and Pu cations in our initial configurations [22]. This Monte Carlo relaxation appears to be relevant a posteriori since it allowed us to have a good improvement of the quality of our optimization results. Then, we created oxygen vacancies in the relaxed bulk supercells to obtain U1-yPuyO2-x systems. Oxygen vacancies are randomly located but must have at least two plutonium atoms among their first neighbours. For each created vacancy, two Pu 4+ first neighbours were indeed turned into Pu 3+ in order to accommodate the charge variation caused by the oxygen vacancies. This procedure has been chosen based on DFT predictions [48] showing that, when an oxygen vacancy is created in U1-yPuyO2, the minimum energy is found for configurations where two atoms of plutonium first neighbours of the vacancy are reduced to Pu 3+ while all the uranium atoms remain U 4+ . Figure 2 presents a U1-yPuyO2-

x structure with a 30 % Pu content and a O/M ratio of 1.97 obtained from the protocol described above.

ii)

Once the structures of U1-yPuyO2-x had been created, we adjusted the Pu 3+ interaction parameters to reproduce the lattice parameter of U1-yPuyO2-x at room temperature using a global minimizer based on the Differential Evolution Algorithm (DEA) of Storn and Price [49]. The lattice parameters of the various targeted systems were then calculated using MD LAMMPS calculations at room temperature. This step is completed when the value of the merit function reaches 10 -4 . The potential parameters adjusted during this step constitute the initial parameter set for the next step of our fitting.

iii)

The structures and parameters obtained during the previous step were then used in a third step. This step consists in MD calculations of the thermal expansion of our systems using again the LAMMPS package. The calculations are carried out for a plutonium content of 30%, stoichiometries of 1.99 and 1.97, and temperatures of 300 K, 600 K, 900 K, 1400 K and 1900 K. This last step uses the Nelder-Mead (or Simplex) local minimisation algorithm [50] to find the parameter set which reproduces the thermal expansion. This step is completed when the value of the merit function reaches 10 -3 when the experimental data of Kato et al. are considered [43].

Note that we tried to simultaneously adjust all the interaction parameters without success. We finally started to fit r0 and ρ together, keeping the other parameters fixed. Then we fitted A and D parameters only, and finally the ɣ parameter alone. The final parameters obtained during this procedure are presented in section 3.1.

Figure 1: Illustration of the fitting procedure used for the optimization of the Pu 3+ interaction parameters in the CRG potential. The final parameters are summarized in Table 2.

Figure 2: Simulation box obtained after MC relaxation and the removal of some oxygen atoms in order to create oxygen vacancies (yellow) with the Pu 4+ ions (light grey), Pu 3+ (green). The remaining U 4+ and O 2-ions are dark grey and red respectively.

Results

Optimisation results

The optimized Pu 3+ parameters are summarized in Table 2. Despite the fact that the fitting procedure is empirical, the evolution of fitted Pu 3+ parameters seems to respond to a physical logic. This is for instance the case for the values of the parameters ρ and r0 of the cationanion interactions. These two parameters represent the equilibrium distance of the Buckingham and the Morse potential respectively. The values of these parameters for the original CRG potential are ρ[Pu 4+ -O 2-] = 0.3793 Å and r0[Pu 4+ -O 2-] = 2.346 Å while those optimized in this work for the Pu 3+ -O 2-interaction are ρ[Pu 3+ -O 2-] = 0.384 Å and r0[Pu 3+ -O 2-] = 2.408 Å. These two parameters increase from Pu 4+ to Pu 3+ . This is consistent with the fact that the reduction of Pu 4+ to Pu 3+ corresponds to the capture of an electron, which leads to a larger ionic radius. Experimentally, the ionic radius of Pu 4+ and Pu 3+ are 0.96 Å and 1.00 Å respectively [45], [51]. This is a satisfying feature of our optimization. Table 1 presents the experimental data of lattice parameters of U0.70Pu0.30O2-x from Kato et al. [43] used in the fitting process along with the computed values using our new potential. A perfect match is found between experimental and calculated data. Although the fitting procedure did not include data for a stoichiometry of 1.98, the new model is able to predict properties with a high level of accuracy at this stoichiometry for all the temperatures considered. This constitutes a first validation of our potential. Figure 3 presents the differences between the reference and the calculated data with the new potential. These differences are very small (always under 0.05 %), they increase with increasing temperature and they reach their maximum of around 0.045 % at a temperature of 1900 K and a stoichiometry of 1.97.

Interaction

Potential validation on U1-yPuyO2-x and transferability to Pu2O3

To validate the fitted potential on data not included in the fitting process, we investigated the evolution of the lattice parameter of U1-yPuyO2-x as a function of the deviation from stoichiometry (x) for various plutonium contents (y). Figure 4 shows the evolution at room temperature of U1-yPuyO2-x lattice parameter as a function of the O/M ratio for three different plutonium contents y=0.125, y=0.30 and y=0.50. These calculations were carried out using Monte Carlo simulations with the newly developed potential. A comparison is made between our calculations (diamonds) and the measurements from Schmitz et al. [38] (lines) which constitute a reference in the literature, as well as the recommendation proposed by Philipponneau [START_REF] Philipponneau | Paramètre cristallin et homogénéité dans le combustible MOX[END_REF] (dotted lines) largely used in fuel performance codes. The calculation results obtained by molecular dynamics from Arima et al. [11] are also presented (circles). In this figure, the black line represents the reduction limit based on the measurements of Schmitz, i.e. the upper limit above which no further reduction of Pu 4+ is possible because all the Pu 4+ in the system have been reduced to Pu 3+ . The evolution at room temperature of the U1-yPuyO2-x lattice parameter as a function of the O/M ratio is linear for all plutonium contents. Although we fitted our potential only on the thermal expansion of U0.70Pu0.30O1.99 and U0.70Pu0.30O1.97, we found a very good agreement between the predictions of our potential and the experimental results from Schmitz et al. [38] as well as the recommendations from Philipponneau [START_REF] Philipponneau | Paramètre cristallin et homogénéité dans le combustible MOX[END_REF]. Conversely, it can be seen that the results from Arima et al. [11] do not provide as good a match as our results particularly at high Pu content. [52] for details) are too small to be seen (about 5×10 -6 Å).

Table 3 presents the average slope (da/d(O/M)) of the evolution curves of the lattice parameter (a) as a function of the stoichiometry O/M. Our MC results are compared to the MD results from Arima et al. [11] obtained with their empirical potential as well as experimental measurements and recommendations taken in the literature. The slope provided by the potential developed in this work is in better agreement with experimental measurements and recommendations than the one provided by the potential of Arima et al. A key property of an interatomic potential is its ability to accurately predict material properties that were not in the training dataset. This is generally referred to as transferability. To extend the validation of our potential, we checked its transferability to another system involving Pu 3+ interactions: Pu2O3. All the plutonium atoms in this system are trivalent. Actinide oxides crystallized in three different structures (hexagonal, monoclinic and cubic). As far as Pu2O3 plutonium sesquioxides are concerned, only two phases have been synthesized: cubic (α-Pu2O3) and hexagonal (β-Pu2O3). Due to the lack of experimental data on the hexagonal β-Pu2O3 phase, our investigation focused only on the cubic α-Pu2O3 phase.

Table 4 presents the lattice parameter (a), the volume (V), the density (ρ), the specific heat (Cp) and the bulk modulus (B0) of the α-Pu2O3 phase computed with our potential using Monte Carlo simulations at room temperature. We applied the calculation procedure described in our previous work [22]. The results we have obtain with our new potential have been compared to literature values obtained using molecular dynamics (MD) and molecular statics (MS) with a BMH type potential fitted on α-Pu2O3 structure. These aimed to reproduce the structure's experimental cell parameter, density and space group at room temperature. Experimental measurements and ab-initio calculations are also compared to our results in Table 4. The values of lattice parameter, volume and density in this table are all consistent whatever the method. Our results are in good agreement with the experimental measurements and previous theoretical works. This demonstrates the transferability and validity of our new potential. Concerning specific heat and bulk modulus of the α-Pu2O3 system, they are slightly overestimated. Calculated specific heat lies within 4% of the experimental data while bulk modulus lies within 13% of the abinitio data. These results remain satisfactory since the experimental errors on these properties are generally high. [START_REF] Haire | Advance in plutonium chemistry 1967-2000[END_REF]. b ref [START_REF] Prodan | [END_REF]. c ref [56]. d ref [57]. e ref [58].

Thermophysical properties of the U0.70Pu0.30O2-x system

According to the U-Pu-O phase diagram, the hypo-stoichiometric region of U1-yPuyO2-x mixed oxide is known to exhibit a large biphasic domain at room temperature [3], [26] depending on the Pu content. This biphasic domain leads to a miscibility gap, what can have direct consequences on the fuel properties. However, the studies show that this miscibility gap exists for stoichiometry under 1.97 at 300 K and its size decreases with increasing temperature. At temperatures higher than ~1000 K, the hypo-stoichiometric U1-yPuyO2-x mixed oxide exhibits only one phase independently of Pu content and stoichiometry.

In order to verify the capacity of our potential to predict properties of hypo-stoichiometric U1-yPuyO2-x mixed oxide outside and inside the miscibility gap, we have carried out studies for stoichiometry ranging from 2 to 1.94 from room temperature up to the expected melting temperatures.

Direct properties

Lattice parameter

In Figure 5 the lattice parameter of the hypo-stoichiometric U0.70Pu0.30O2-x is plotted as a function of temperature for various O/M ratios alongside the correlation derived from experimental data by Kato et al. [43]. These results show that the lattice parameter increases with temperature for all stoichiometries up to the expected melting temperature. From room temperature up to 2600 K, the evolution of the lattice parameter as a function of stoichiometry is almost regular and shows an increase with a decrease of the O/M ratio. The behaviour observed in this region is similar to the one obtained previously by molecular dynamics with the potential from Arima [11], [25]. Moreover, our results obtained with our new potential fit perfectly the Kato's experimental correlation. According to Watanabe et al. [24] the substitution of 4+ ions by larger 3+ ions [45], [51] associated with the presence of vacancies reduces the resulting Coulomb interaction strength and tends to cause an increase in the lattice parameter.

Figure 5: Lattice parameter of U0.70Pu0.30O2-x as a function of temperature for different stoichiometries. Solid lines are correlations derived from experimental data by Kato et al. [43] and data points are MC simulations using our new potential. The error bars calculated from the block averaging method (see ref. [START_REF] Grossfield | Chapter 2 Quantifying Uncertainty and Sampling Quality in Biomolecular Simulations[END_REF] for details) are too small to be seen (about 5×10 -6 Å).

Above 2600 K, the effect of the O/M ratio on the lattice parameter gradually decreases with increasing temperature. At high temperatures, above 2900 K, this effect becomes almost negligible, the lattice parameters of hypo-stoichiometric U0.70Pu0.30O2-x being identical to those of stoichiometric U0.70Pu0.30O2 mixed oxide. It is commonly observed that some crystalline materials begin to acquire significant ionic mobility above a specific onset temperature [59], [START_REF] Golunski | Why use platinum in catalytic converters?[END_REF] termed the "Tammann temperature". This specific onset temperature has been identified in superionic materials. Concerning UO2, MD studies [61]- [63] showed significant oxygen disorder at temperatures above 2000 K in agreement with neutron scattering experiments [64]- [67]. These temperatures also acknowledge a change of the material dynamics from that of a "simple" crystal to some kind of highly defective crystalline structure. In this condition, the oxygen sublattice is disordered and it is very difficult to identify oxygen vacancies. The hypo-stoichiometric U0.70Pu0.30O2-x oxygen sublattice should be therefore relatively similar to the corresponding stoichiometric U0.70Pu0.30O2 sublattice at high temperature (above 2600 K). This particular behaviour at high temperature probably reduces the effect of the O/M ratio on the thermodynamic properties and leads to the same lattice parameters for stoichiometric and hypostoichiometric mixed oxides under these conditions.

Enthalpy

Enthalpy is defined as H = U + PV = 3 2 NkT + U conf + PV, with N=(3-x) where x describes the stoichiometry. U, P, and V are respectively the internal energy, the pressure and the volume of the system. The energy U can also be written as the sum of the configurational energy of the system Uconf (the energy coming from the potential) and a kinetic contribution equal to 3 2 NkT . Enthalpy is plotted in Figure 6.a as a function of temperature for various O/M ratios. We can see that the increase of the O/M ratio leads to a decrease of the enthalpy, whatever the temperature up to the expected melting point. The computed enthalpy increment i.e. H(T) -H(300 K) of hypo-stoichiometric U70Pu30O2-x is displayed in Figure 6.b for various O/M ratios in comparison with experimental measurements from Lorenzelli et al. [68], Orgard et al. [START_REF] Ogard | High temperature heat content and heat capacity of uranium dioxide and uranium dioxide-plutonium dioxide solid solution[END_REF], and Gibby et al. [START_REF] Gibby | [END_REF] for MOX with three different plutonium contents (25 %, 20 %, and 19.8 %). Figure 6.b shows the evolution of the U1-yPuyO2-x enthalpy increment as a function of the temperature: it increases linearly below 2400 K while for temperatures above 2400 K, we note a "bump". This behaviour is observed for all investigated stoichiometries. It could be attributed to the phase transition due to the anionic sublattice instability as previously pointed out for the stoichiometric MOX [22], [71]. Our calculation results are consistent with the experimental measurements up to 1800 K. Above 1800 K, our potential tends to underestimate the enthalpy increment. This underestimation has also been observed for stoichiometric MOX in previous works using original the CRG potential [START_REF] Balboa | [END_REF], [22], [72]. [69] and Gibby et al. [START_REF] Gibby | [END_REF]. Diamonds are predictions from our new potential for 30 % Pu content and plus stand for experimental measurements for three Pu contents (20; 19.8 and 25 %). The error bars calculated from the block averaging method (see ref. [START_REF] Grossfield | Chapter 2 Quantifying Uncertainty and Sampling Quality in Biomolecular Simulations[END_REF] for details) are too small to be seen (about 7×10 -2 J.mol -1 ). a) b)

Derivative properties

Linear thermal expansion coefficient

The linear thermal expansion coefficient (α) is defined as the first derivative of the lattice parameter with respect to the temperature, as shown in the following equation: 3) [22], [73].

𝜶 = 𝟏 〈𝑽〉𝒌𝑻 𝟐 (〈𝑽𝑯 ̂〉 -〈𝑽〉〈𝑯 ̂〉) Eq. 3.2
where 𝑉 is the volume of the system, 𝑘 is Bolztmann constant, 𝐻 ̂ is the configurational enthalpy, 𝐻 ̂= 𝑈 𝑐𝑜𝑛𝑓 + 𝑃𝑉 where 𝑈 𝑐𝑜𝑛𝑓 is the potential energy. We then deduce the linear thermal expansion coefficient dividing the volume thermal expansion coefficient by a factor of 3.

To our knowledge, experimental data on the linear thermal expansion coefficient for hypostoichiometric mixed oxide in the literature are scarce. Moreover, the available data from Leblanc et al. [START_REF] Leblanc | EURATOM/USA Rep[END_REF], Lorenzelli et al. [68], Roth et al. [START_REF] Roth | [END_REF] and Gibby [START_REF] Gibby | [END_REF] are scattered and old. Such data exist only for plutonium contents between 19.1% and 25% and for temperatures between 293 K and 1273 K. Nonetheless, these data have been used to establish Martin's recommendations for all Pu contents up to 1800 K [76].

Table 5 summarises the averaged linear thermal expansion coefficient (αL) of hypo-stoichiometric U1-yPuyO2-x mixed oxide in the temperature range from 300 K to 1270 K. Our calculated values are presented with the previous ones from MD as well as experimental measurements and recommendations. For the whole range of stoichiometry investigated, our results show no significant effect of the O/M ratio on the averaged linear thermal expansion coefficient (α) of U1-yPuyO2-x mixed oxide. By comparison others show a large effect of the O/M ratio [76], [68]. Nevertheless, recent thermal dilatometry measurements of linear thermal expansion coefficient in an oxygen partial pressurecontrolled atmosphere from Kato et al. [43] are consistent with our results. Kato measured at 1573 K, the linear thermal expansion coefficient equal to 3.79, 3.84, 3.90 for O/M ratio equal to 1.99, 1.98 and 1.97 respectively. These measurements deviate by less than 3 % from our results. In the latter experiments, the oxygen partial pressure was controlled by adjusting the hydrogen partial pressure to moisture partial pressure ratio in an Ar/H2/H2O gas mixture. This leads to hold a constant oxygen to metal ratio in the U1-yPuyO2-x during the measurement, which was not necessarily the case in older experiments. They pointed out that the published data on which Martin's recommendations were based had been obtained from measurements carried out in Ar/H2 gas mixture, which could lead to a change of the O/M ratio during the measurement. Therefore, they considered that the O/M change during the measurement is the cause of the larger O/M dependency of Martin's recommendations compared to their data. Conversely, the O/M ratio in the study from Kato et al. [43] was stated to be stable during the measurement with an uncertainty estimated to be around ± 0.004. We can add that the data used within Martin's recommendations were obtained using XRD measurements. According to Belin et al. [44] no direct method is available to determine the O/M ratio during in situ HT-XRD experiments.

Moreover, the mass of samples cooled back to room temperature is also too low to allow for O/M determination using traditional gravimetry methods. Consequently, we trust our calculations results even if they are only consistent with Kato's experimental results.

Table 5: Averaged linear thermal expansion coefficient (αL in 10 -6 .K -1 ) of hypo-stoichiometric U1-yPuyO2-x mixed oxide in the temperature range from 300 K to 1270 K. Our calculated values by Monte Carlo simulation are presented along with previous results from molecular dynamics (MD) as well as experimental measurements and recommendations. An estimation of error bars calculated from the block averaging method (see ref.

[52] for details) is presented in Figure 7. Two main regions appear in this graph: the first region concerns temperatures ranging from 300 K to 2300 K, where 𝛼L values show a regular increase with temperature for all O/M ratios. Martin's recommendations also suggest an increase in linear thermal expansion coefficient with O/M ratio in this region. Concerning our results, it is difficult to highlight the effect of the O/M ratio due to large uncertainties, however, as mentioned above, our results are consistent with recent measurements from Kato et al. [43] indicating no significant effect of O/M ratio up to 1923 K. Moreover, according to Ma et al. [25] the impact of oxygen vacancy concentration on the linear thermal expansion coefficient of U0.75Pu0.25O2-x is negligible when the temperature falls below 1300 K. Ma et al. [25] highlight that the presence of oxygen vacancies probably have a larger effect on the linear thermal expansion coefficient in the high temperature range.

O/

Concerning the second region visible in the Figure 7 (temperatures from 2300 K up to the expected melting temperature), our results show that the linear thermal expansion coefficient of the hypostoichiometric mixed oxide increases significantly with temperature up to a maximum before decreasing. This defines a peak. The location of this peak seems to vary for each O/M ratio, but in our results the effect of the O/M ratio is difficult to detect taking into consideration the error bars. Nevertheless, the global behaviour observed for the linear thermal expansion coefficient calculated with our potential is similar to the previous molecular dynamics simulations results on stoichiometric mixed oxide using the original CRG potential [START_REF] Balboa | [END_REF], [21], [22]. The authors highlight that the sudden rise followed by a decrease observed in the high temperature region may be related to the Bredig transition, which occurs at the temperature Tλ (see the section devoted to the specific heat for more details). Unfortunately, no experimental data exist on hypo-stoichiometric mixed oxides in the high temperature domain and high oxygen deficiency conditions.

Specific heat

Compared to the enthalpy increment, more information can be extracted from analysing specific heat. As already pointed by Jorgensen et al. [77], the specific heat can also be separated into two contributions residual (𝑪 𝒑 𝒓𝒆𝒔 ) and ideal (𝑪 𝒑 𝒊𝒅 ) parts. The ideal part can be considered as constant above the Debye temperature and obtained by considering the compound as an ideal gas. The residual part is computed by analysing the fluctuations of the energy within a Monte Carlo simulation using the following formula. More details on this computational method are available in ref. [22], [73]. The specific heat of U1-yPuyO2-x mixed oxide is shown in Figure 8 as a function of temperature for different O/M ratios alongside experimental data from Gibby et al. [START_REF] Gibby | [END_REF], Affortit et al. [START_REF] Affortit | Mesure de la chaleur spécifique des métaux jusqu'à leur température de fusion[END_REF] and ANL [79] coming from the IAEA Therpro database [START_REF]Therpro database from IAEA[END_REF]. As expected, the calculated specific heat shows the same trend and the same behaviour as previously observed for the linear thermal expansion coefficient, with two main regions.

The first region concerns temperatures below 2300 K where specific heat regularly increases with the temperature, and without any noticeable effect of the O/M ratio. In this region, Kato et al. [START_REF] Kato | [END_REF] found that the specific heat at constant volume (CV) slightly decreases with an increasing deviation x in PuO2x due to a decrease of the number of atoms per molecule, and the change of Cp with O/M ratio from 2.00 to 1.92 was small and within 2 J/mol/K. They highlight that this effect of the O/M ratio on Cp was negligible in the O/M range of 2.00-1.92. These observations are consistent with our results on specific heat of U1-yPuyO2-x.

The second region goes from 2300 K up to the expected melting temperature where specific heat shows the characteristic peak attributed to the superionic transition of the fluorite-type structure. This transition is largely discussed in the literature (see for examples the work from Cooper et al. [20] and the one from Zhang et al. [START_REF] Bathellier | Heat capacity calculation of mixed-oxide (MOX) fuel using classical molecular dynamics[END_REF]). Although the previous studies on stoichiometric U1-yPuyO2 mixed oxide do not show any effect of plutonium content in the temperature range from room temperature to around 2300 K, its effect is noticeable at high temperature (around 2600 K) showing an evolution of the peak maximum with the plutonium content [START_REF] Bathellier | Heat capacity calculation of mixed-oxide (MOX) fuel using classical molecular dynamics[END_REF]. In the present study, the effect of stoichiometry is hardly noticeable in the second region (high temperature region) although we should expect peak maxima depending on the O/M ratio. Following what has been done by Bathellier et al. [START_REF] Bathellier | Heat capacity calculation of mixed-oxide (MOX) fuel using classical molecular dynamics[END_REF] and Cooper et al. [20] to highlight the effect of the Pu content on heat capacity of stoichiometric MOX, we plan to perform MD calculations through numerical derivatives of the enthalpy increment, in order to investigate the effect of the O/M ratio in this region. Nevertheless, the superionic transition peak temperature (Tλ) of fluorite-type materials structure is related to the melting point (Tm) by the equation (𝑇 𝑚 = 𝑇 𝜆 0.85

) [START_REF] Zhang | [END_REF]. As the melting temperature decreases with the deviation from stoichiometry [84], this helps explain the dependence of the peaks of the Bredig transition in function of stoichiometry. However, as already observed for stoichiometric oxide and mixed oxide systems, all peaks are around 2400 K -2600 K.

Conclusion

In this work, we built an original procedure combining classical molecular dynamics and Monte Carlo simulation methods in order to optimize an extension of the many-body CRG potential including Pu 3+ interactions. This new potential has been fitted to match the most recent dilatometry measurements of the lattice parameters of U0.70Pu0.30O1.99 and U0.70Pu0.30O1.97 from 300 K to 1900 K [43], [45]. This new potential has proven to be very effective for accurate description of many thermophysical properties of hypo-stoichiometric mixed oxide systems. Concerning the lattice parameter, although we fitted our potential only on one Pu content and two stoichiometries, we found a very good agreement between the predictions based on our potential and the available experimental results as well as the recommendations in wide ranges of O/M ratios, Pu contents and temperatures. Furthermore, the potential shows good transferability to the plutonium sesquioxide α-Pu2O3. This new validated potential has then been used to investigate properties of hypo-stoichiometric U0.70Pu0.30O2-x mixed oxide fuels from room temperature up to the expected melting temperatures. Both direct properties, namely lattice parameter and enthalpy, and derivative properties, namely linear thermal expansion coefficient and the specific heat, have been considered. Our results show a significant effect of the degree of hypo-stoichiometry on the direct properties, with a gradual increase of the lattice parameter and the enthalpy as the O/M ratio decreases, and a good agreement with experiments is found. Conversely, as far as derivative properties are concerned, no significant effect of the O/M ratio can be brought out. However, around the Bredig transition temperature, the obtained values are too scattered to allow reliable conclusions concerning a possible effect of the hypo-stoichiometry degree.

As outlooks of this work will be the calculation of these thermophysical properties for other Pu contents, in order to have a global vision of the coupled effects of Pu content and stoichiometry on the properties of MOX. Moreover, the empirical interatomic potential developed in this work will allow us to investigate the existence of a supposed miscibility gap in the hypo-stoichiometric domain [3]- [5]. The procedure developed in this study will also be used to optimize other interactions parameters such as Am 3+ in order to investigate the effect of americium content in mixed oxide fuel.

Figure 3 :

 3 Figure 3: Differences between reference data and calculated data with the new potential. The reference data corresponding to the stoichiometry of 1.98 were not used in the fitting procedure.

Figure 4 :

 4 Figure 4: Evolution at room temperature of the U1-yPuyO2-x lattice parameter as a function of stoichiometry (O/M ratio) for different plutonium contents y=0.125 in red; y=0.30 in green and y=0.50 in blue. Predictions using the new Pu 3+ potential (diamonds) are compared to the experimental measurement from Schmitz et al. [38] (lines), molecular dynamics data from Arima et al. [11] for y=0.10, 0.20 and 0.30 (circles) as well as Philipponneau's recommendation [39] for y=0.125, 0.30 and 0.50 (dotted lines). The error bars calculated from the block averaging method (see ref. [52] for details) are too small to be seen (about 5×10 -6 Å).

Figure 6 :

 6 Figure 6: Enthalpy a) and enthalpy increment b) of U1-yPuyO2-x mixed oxide as a function of temperature for different O/M ratios alongside experimental measurements from Lorenzelli et al. [68], Orgard et al.[69] and Gibby et al.[START_REF] Gibby | [END_REF]. Diamonds are predictions from our new potential for 30 % Pu content and plus stand for experimental measurements for three Pu contents (20; 19.8 and 25 %). The error bars calculated from the block averaging method (see ref.[START_REF] Grossfield | Chapter 2 Quantifying Uncertainty and Sampling Quality in Biomolecular Simulations[END_REF] for details) are too small to be seen (about 7×10 -2 J.mol -1 ).

1 where.

 1 𝑎 stands for the lattice parameter. Note that the linear thermal expansion coefficient is equal to one third of the volumetric thermal expansion which is defined as 𝛼 = In our Monte Carlo simulations, we first compute the thermal expansion coefficient by analysing volume fluctuations, as shown in equation (Eq.3.

Figure 7 :

 7 Figure 7: Linear thermal expansion coefficient (α) of U0.70Pu0.30O2-x mixed oxide as a function of temperature for different O/M ratios compared to experimental Martin's recommendations. Data points are predictions from the new potential and dashed lines are experimental Martin's recommendations for hypo-stoichiometric mixed oxide. Solid line represents Martin's recommendation for UO2.

  〈𝑼𝑯 ̂〉 -〈𝑼〉〈𝑯 ̂〉) + 𝑷 𝒌𝑻 𝟐 (〈𝑽𝑯 ̂〉 -〈𝑽〉〈𝑯 ̂〉) -𝑵𝒌 Eq. 3.3

Figure 8 :

 8 Figure 8: Specific heat evolution of U1-yPuyO2-x mixed oxide as a function of temperature for different O/M ratios alongside experimental measurements from Gibby et al. [70], Affortit et al. [78] and ANL [79] taken in the IAEA Therpro database [80]. Diamonds are predictions from our new potential for 30 % Pu content and lines stand for experimental measurements for two Pu contents (20 and 25 %).

  𝑼 𝜶𝜷 (𝒓 𝒊𝒋 ) = 𝑼 𝑹𝑰𝑷 (𝒓 𝒊𝒋 ) + 𝑼 𝑴 (𝒓 𝒊𝒋 ) 𝑞 𝛽 , 𝐴 𝛼𝛽 , 𝜌 𝛼𝛽 , 𝐶 𝛼𝛽 , 𝐷 𝛼𝛽 , 𝛾 𝛼𝛽 and 𝑟 𝛼𝛽 0 are adjustable parameters.

								Eq. 2.2
	𝑼 𝑹𝑰𝑷 (𝒓 𝒊𝒋 ) = 𝑼 𝑪 (𝒓 𝒊𝒋 ) + 𝑼 𝑩 (𝒓 𝒊𝒋 ) =	𝒒 𝜶 𝒒 𝜷 𝟒𝝅𝜺 𝟎 𝒓 𝒊𝒋	+𝑨 𝜶𝜷 𝒆𝐱𝐩 (-	𝒓 𝒊𝒋 𝝆 𝜶𝜷	) -	𝑪 𝜶𝜷 𝒓 𝒊𝒋 𝟔	Eq. 2.3
	𝑼 𝑴 (𝒓 𝒊𝒋 ) = 𝑫 𝜶𝜷 [	𝐞𝐱𝐩 (-𝟐𝜸 𝜶𝜷 (𝒓 𝒊𝒋 -𝒓 𝜶𝜷 𝟎 )) -𝟐𝐞𝐱𝐩 (-𝜸 𝜶𝜷 (𝒓 𝒊𝒋 -𝒓 𝜶𝜷 𝟎 )) ]		Eq. 2.4
	where 𝑞 𝛼 , The second component of equation Eq. 2.1 represents the many-body EAM interaction term added by
	Cooper et al. to improve the potentials based on pair interaction models. This term allows to successfully
	capture the Cauchy-violation (i.e. C12≠C44 where Cij are elastic tensor elements in Voigt notation), which
	is necessary for accurate reproduction of the elastic constants and the mechanical behaviour of fluorite
	type structure. The parameter 𝜎 𝛽 in this equation is given by the equation Eq.2.6., in which the function
	erf [32] stands for the error function and 𝜂 𝛽 and 𝐺 𝛼 the adjustable parameters for the many-body EAM
	interaction term.
	𝝈 𝜷 (𝒓 𝒊𝒋 ) =	𝟏 𝟐 (	𝜼 𝜷 𝒓 𝒊𝒋 𝟖 ) {𝟏 + 𝒆𝒓𝒇[𝟐𝟎(𝒓 𝒊𝒋 -𝟏. 𝟓)]}		Eq. 2.5

Table 1 :

 1 Experimental lattice parameters (Å) of U0.70Pu0.30O2-x hypo-stoichiometric from Kato et al.[43] used in the fitting process and lattice parameters (Å) predicted with the new Pu 3+ potential. Note that the potential has not been fitted on the data corresponding to the stoichiometry of 1.98. The agreements between the predictions of new Pu 3+ potential and the experimental values are given as a percentage in the brackets.

	Stoichiometry	1.99	1.98	1.97
	T (K)		Experiment
	300	5.450	5.454	5.457
	600	5.467	5.470	5.474
	900	5.484	5.488	5.492
	1400	5.516	5.521	5.526
	1900	5.554	5.560	5.565
		This study (using our new potential)
	300	5.450 (0.000 %)	5.453 (0.006 %)	5.456 (0.020 %)
	600	5.467 (0.008 %)	5.470 (0.001 %)	5.473 (0.019 %)
	900	5.485 (0.012 %)	5.489 (0.009 %)	5.492 (0.014 %)
	1400	5.518 (0.020 %)	5.522 (0.011 %)	5.525 (0.019 %)
	1900	5.557 (0.029 %)	5.560 (0.006 %)	5.563 (0.046 %)

Table 2 :

 2 Pu 3+ interaction parameters obtained in this work as an extension of the original CRG potential.

Table 3 :

 3 Slopes

		Method Slope da/d(O/M)	Refs.
	Experiment	XRD	2.75	Benedict et al. [41]
		XRD	3.45	Schmitz et al. [38]
		XRD	3.3	Omichi et al. [42]
		XRD	3.2	Duriez et al. [53]
	Recommended	-	3.20	Philipponneau [39]
		-	3.15	Harding et al. [40]
	Calculations	MD	4.35	Arima et al. [11]
		MC	3.15	This study

da/d(O/M) of the evolution curve of the U1-yPuyO2-x lattice parameter (a) as a function of the stoichiometry O/M. Our Monte Carlo (MC) results are presented alongside MD results and experimental X-ray diffraction measurements and recommendations.

Table 4 :

 4 Lattice

parameter (a), volume (V), density (ρ), specific heat (Cp) and bulk modulus (B0) of α-Pu2O3 calculated with our new potential using Monte Carlo simulations at room temperature. These results are compared to results from molecular dynamics (MD), molecular statics (MS) as well as experimental measurements and ab-initio calculations.

α-

Pu2O3 Exp a This study Ab-initio b,c MD d MS d,e

	a (Å)	11.04	10.92	11.17	11.05	11.04
	V (Å 3 )	1345.6	1303.14	1393.7	1349.2	1345.6
	ρ (g/cm 3 )	10.4	10.73	10.04	10.39	10.39
	Cp (J/mol/K) 120.35	126.06	-	117.6	105.3
	B0 (GPa)	-	145.73	128	-	143
	a ref

  Figure7shows the evolution of the linear thermal expansion coefficient of U0.70Pu0.30O2-x mixed oxide as a function of temperature for different O/M ratios alongside Martin's experimental data[76]. The change in linear thermal expansion coefficient as a function of temperature shows a similar trend, no matter what the O/M ratio is.

					M ratio
			1.99	1.98	1.97	1.96	1.94	Method	Refs.
		20	3.85	-	-	-	4.36	Exp.	Lorenzelli et al. [68]
	% Pu	0 25	3.60 -	3.74 3.76	3.87 -	4.01 -	4.29 3.89	Recom. MD	D.G. Martin [76] J. Ma et al. [25]
		30	3.73	3.67	3.66	3.68	3.63	MC	This study

For simplification, the following notation will be adopted. Any ion X of total charge [±C] will be noted X ±C instead of X[±C] although its total charge is not always equal to its effective charge. For instance, trivalent plutonium Pu[+III] = Pu 3+ .

Acknowledgements

We thank M.W.D. Cooper for his practical advice. The authors wish to express their gratitude to M. Bertolus, G. Jomard, M. Freyss, D. Batellier, J. Lechelle, R. Dianzinga, N. Chauvin, M. Lainet and P. Martin for fruitful discussions. We used in this work the GIBBS Monte Carlo code co-developed by IFPEN and CNRS. High-performance resources from Grand Equipement National de Calcul Intensif (GENCI) [Centre Commun de recherche et Technologie (CCRT)] were used for this investigation. This research is part of the Investigations Supporting MOX Fuel licensing in ESNII Protocol Reactors (INSPYRE) project, which has received funding from Euratom research and training program 2014-2018 under grant 754329.