Hector Balboa

Laurent Van Brutzel
email: laurent.vanbrutzel@cea.fr

Alain Chartier

Yann Le Bouar

Assessment of empirical potential for MOX nuclear fuels and thermomechanical properties

Keywords: Molecular dynamics, nuclear fuel, mechanical, assessment

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Nowadays, uranium dioxide (UO 2 ) is being used as the standard nuclear fuel in fission nuclear reactors and has been extensively studied since the sixties. In parallel, nuclear fuel containing a mixture of uranium and plutonium oxides (MOX) as principle components provides an alternative, due to the fact that: [START_REF] Popov | Thermophysical Properties of MOX and UO 2 Fuels Including the Effects of Irradiation[END_REF] it allows large quantities of fissile isotopes produced in spent nuclear fuel from light water reactors to be recycled, [START_REF]ESNII+ European commission. Deliverable D7.11 "state of the art with a literature review of MOX properties[END_REF] it can be seen as a more efficiently way of using the uranium dioxide, since the abundant 238 U found in natural uranium is a Pu producer, (3) it can be taken as a solution for the increasing stockpile of Pu around the globe coming from either nuclear weapons and commercial reactors, and (4) it is designated as the most probable fuel for future fast breeder reactors, such as ASTRID (Advanced Sodium Technological Reactor for Industrial Demonstration).

A very important issue for the future of nuclear power is to ensure safeness and effectiveness during processes involving MOX fuel such as fabrication, operation and recycling. Yet, the toxicity of plutonium and high radiation levels make experiments less viable. Nonetheless, beside previous experimental efforts gathered in the following reviews [START_REF] Popov | Thermophysical Properties of MOX and UO 2 Fuels Including the Effects of Irradiation[END_REF][START_REF]ESNII+ European commission. Deliverable D7.11 "state of the art with a literature review of MOX properties[END_REF], experiments with MOX are difficult to perform, especially at high temperatures and under irradiation condition. For these reasons, numerical approaches can be chosen to bring some insight on basic physical phenomena that take place in the fuel matrix. For instance, over the last decade, several atomistic approaches using molecular dynamics (MD) simulations have been carried out to study thermal conductivity properties in (U,Pu)O 2 [START_REF] Yamada | Evaluation of thermal properties of mixed oxide fuel by molecular dynamics[END_REF][START_REF] Kurosaki | Molecular dynamics study of mixed oxide fuel[END_REF][START_REF] Terentyev | Molecular dynamics study of oxygen transport and thermal properties of mixed oxide fuels[END_REF][START_REF] Arima | Evaluation of thermal properties of UO 2 and PuO 2 by equilibrium molecular dynamics simulations from 300 to 2000 K[END_REF][START_REF] Arima | Evaluation of thermal conductivity of hypostoichiometric (U,Pu)O 2-x solid solution by molecular dynamics simulation at temperatures up to 2000 K[END_REF][START_REF] Arima | Equilibrium and nonequilibrium molecular dynamics simulations of heat conduction in uranium oxide and mixed uranium?plutonium oxide[END_REF][START_REF] Nichenko | Molecular Dynamics study of the mixed oxide fuel thermal conductivity[END_REF][START_REF] Ma | Molecular dynamical study of physical properties of (U 0.75 Pu 0.25 )O 2-x[END_REF][START_REF] Cooper | Modelling the thermal conductivity of (U x Th 1-x )O 2 and (U x Pu 1-x )O 2[END_REF][START_REF] Li | Molecular dynamics study of thermal conductivities of (U 0.7-x Pu 0.3 Am x )O 2[END_REF].

However, the reliability of the results depends exclusively on the choice of the set of potentials. The potential parameters are usually fitted to reproduce a few physical properties, typically the lattice parameter, the cohesive energy, and complementary the elastic constants, which comes from experimental values or if not available from ab initio calculations. Therefore, each potential has its domain of validity. Subsequently, others physical properties for which the set of parameters are not been fitted need to be assessed to provide a good insight of advantages and disadvantages of each potential and their range of validity. This type of study has already been carried out in the case of UO 2 [START_REF] Govers | Comparison of interatomic potentials for UO 2 . Part I: Static calculations[END_REF][START_REF] Govers | Comparison of interatomic potentials for UO 2 . Part II: Molecular dynamics simulations[END_REF][START_REF] Potashnikov | High-precision molecular dynamics simulation of UO 2 -PuO 2 : Pair potentials comparison in UO 2[END_REF][START_REF] Chernatynskiy | Critical assessment of UO 2 classical potentials for thermal conductivity calculations[END_REF] but, to our knowledge, not yet for MOX. Therefore, in this study, we assess available rigid ion model empirical potentials for MOX on the structural, thermodynamics, and mechanical properties. The assessment is performed over the full range of plutonium composition, from pure UO 2 to pure PuO 2 and for temperatures ranging from 300 K to melting point.

With MD method, actinide atoms are usually simulated in the approximation of rigid ions and pair interactions. For the mixed oxide compound (U,Pu)O 2 , several interatomic potentials are available in the literature. There exists two main families of potentials. One that considers U and Pu cations as one single entity and hence they include only three set of parameters (A-A, A-O, and O-O) but depends on the relative percentage of Pu in the MOX [START_REF] Basak | A Novel Pseudo-Ion Approach in Classical MD Simulation: A Case Study on (U 0.8 Pu 0.2 )O 2 Mixed Oxide[END_REF]. The second one treats explicitly the U and Pu cations. Therefore, they include six set of parameters (U-U, U-Pu, Pu-Pu, U-O, Pu-O, and O-O) and do not depend on the percentage of Pu. Because we are interested in studying the spatial repartition of both cation and anion sublattices, we will only consider and describe the second type of force field.

Five rigid ion model potentials have been found in the literature and are tested herein. They will be coined by the name of the first author: Yamada [START_REF] Yamada | Evaluation of thermal properties of mixed oxide fuel by molecular dynamics[END_REF], Arima [START_REF] Arima | Evaluation of thermal properties of UO 2 and PuO 2 by equilibrium molecular dynamics simulations from 300 to 2000 K[END_REF], Potashnikov [START_REF] Potashnikov | High-precision molecular dynamics simulation of UO 2 -PuO 2 : Pair potentials comparison in UO 2[END_REF], Tiwary [START_REF] Tiwary | Interatomic potentials for mixed oxide and advanced nuclear fuels[END_REF], and Cooper [START_REF] Cooper | A many-body potential approach to modelling the thermomechanical properties of actinide oxides[END_REF][START_REF] Cooper | Thermophysical properties and oxygen transport in the (U x ,Pu 1-x )O 2 lattice[END_REF]. We will present first the method used, then we will discuss the results obtained for the lattice parameter, the thermal expansion, the specific heat capacity, the elastic constants, the stress-strain curves under uniaxial deformation, and crack propagation.

Computational method

These five force fields can be separated according the properties on which they have been fitted. All potentials have been fitted to reproduce correctly the thermal expansion up to the maximum temperature available by experiments at the time, which is about 2100 K. Historically, Yamada was the first one followed by Arima and Potashnikov with some improvement at high temperature, up to the melting point. Tiwary potential includes also fits on the formation energy of point defects (Frenkel pairs), while the Cooper potential focuses on reproducing accurately the mechanical behaviour and has been fitted on experimental data for single crystal elastic constants.

However, the Tiwary potential, did not provide a stable fluorite structure for either UO 2 nor PuO 2 using finite temperature MD simulations. Nevertheless it gives good results using static calculations (0 K). Therefore, we drop this potential from our assessment.

The four remaining potentials are described with classical short-range (Buckingham and Morse) and long-range (van der Waals and Coulomb) interactions (see Eq. 1). Table 1 presents the forms of the four potentials. 6 +D αβ e -2γ αβ (r-r 0 ) -2 e -γ αβ (r-r 0 ) + q α • q β 4πǫ 0 r (1) where A αβ , C αβ , D αβ , γ αβ , q α and r 0 are parameters whose values are given in Table 2. The first term in Eq. 1 is designed to reproduce the repulsion originating from the Pauli's exclusion principle. The second term is the attractive van der Waals' interaction. The third term (Morse) is designed to describe more accurately covalent bond and the vibrations in molecules. Finally, the last term represents the long-range Coulomb interaction.

U αβ (r) = A αβ e -r/ρ αβ - C αβ r
For the Cooper potential, an EAM (Embedded Atom Model) many-body term is added (Eq. 2 and Eq. 3), in order to reproduce the Cauchy's violation observed in actinide oxides with the fluorite structure (i.e. C 12 C 44 ) [START_REF] Axe | Long-Wave Lattice Dynamics of the Fluorite Structure[END_REF], which cannot be reproduced by only pair-wise potentials.

E EAM i = -G α j i σ β (r i j ) (2)
where, σ β is computed as follow:

σ β (r i j ) = 1 2        η β r 8 i j        1 + erf 20(r i j -1.5) . (3)
erf stands for the error function and η β is a parameter. Tables 3 provides the EAM Cooper potential parameters.

Potential

Coulomb Buckingham Morse EAM Yamada Arima Potashnikov Cooper In this study, we will compare systematically via MD simulations structural properties (lattice parameter, thermal expansion), thermodynamic properties (heat capacity, enthalpy), and mechanical properties (elastic constants, stress-strain curves) for the four empirical potentials remaining over temperatures ranging from 300 K to the melting point (∼3000 K) and for compounds from pure UO 2 to pure PuO 2 .

The whole study was carried out using the code LAMMPS [22]. Most of the assessment requires relatively small system sizes: about 7 × 7 × 7 fluorite structure unit cells involving 4116 atoms. Moreover, it ensures enough statistics for the different configurations. For the crack propagation, the size of the box is increased to 440 × 110 × 7 fluorite structure unit cells. In all the cases we are using periodic boundary conditions to avoid surface effects. The Coulomb interactions are calculated with the full Ewald summation procedure. The relaxation is done under NPT thermodynamic ensemble using the well-known Parrinello-Rahman's algorithm [START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF]. The systems are relaxed 100 ps and all the thermodynamic properties are averaged over the last 20 ps.

Since the force fields are based on rigid ion model, it is not possible to study none stoichiometric compounds. Therefore, we can consider the (U 1-y Pu y )O 2 compound as a continuous solid solution as shown by the assessment of Belin et al. [START_REF] Belin | In situ high temperature X-Ray diffraction study of the phase equilibria in the UO 2 -PuO 2 -Pu 2 O 3 system[END_REF]. The plutonium atoms are then distributed randomly in the unit cell on the a-Wyckoff sites (fluorite-like cation sites). The plutonium concentration ranges from 0 to 100 at.% with mainly 15 at.% steps. However, for clarity most of the figures presented herein will only display compounds containing 0, 25, 75, and 100 at.% of plutonium atoms or in some cases it will only be shown for the most relevant case.

Results and discussion

Lattice parameter

The first structural property to fit is the evolution of the lattice parameter with the temperature. Therefore, all the interatomic potentials studied should fit more or less the experimental results. However, Yamada and Arima fitted their potential only up to 2100 K, whereas Potashnikov and Cooper fitted their potential with values up to 2900 K. It is worth to mention that experimental data are really sensitive to the O/U ratio hence we keep our comparison with the strict stoichiometric compound. Furthermore, experimental works are scarcely available at high temperatures (close to the melting point), thus it is difficult to estimate the statistical dispersion for these temperatures. As expected, the changes in the calculated lattice parameters for pure UO 2 as a function of temperature for all potentials are in good agreement with experimental data [START_REF] Yamashita | Thermal expansions of NpO 2 and some other actinide dioxides[END_REF] and Fink's recommendations [26] up to 2100 K (not shown). We have performed a similar comparison for pure PuO 2 and the result is presented in Fig. 1. Cooper, Potashnikov, and Arima potentials fit relatively well the experimental data up to 1300 K. For temperatures higher than 2100 K, Arima potential underestimates consistently the recommendation lattice parameters. On the other hand, Yamada potential shows a large discrepancy with no continuous evolution for the full range of temperature. This behaviour has already been pointed out by Potashnikov et al. in their potential assessment [START_REF] Potashnikov | High-precision molecular dynamics simulation of UO 2 -PuO 2 : Pair potentials comparison in UO 2[END_REF]. With the Yamada potential, a phase transition from fluorite to rutile-like structure occurs spontaneously for PuO 2 at high temperatures. This demonstrates that for Yamada potential this new phase is more stable than the fluorite structure contrary to experiments [START_REF] Gardner | The plutonium-oxygen phase diagram[END_REF]. Therefore, this potential is not considered for the following study.

The evolution of the lattice parameters as a function of the temperature for each potential, excluding Yamada's due to the anomalies discussed above, with respect of the plutonium content is similar than for the pure compounds. Hence, we can fit MD data with a general expression of the form: a(T, y) = AT 3 + BT 2 + CT + D + my, where y denotes the plutonium content (0<y<1) and T is the temperature (300 K to melting point). The parameters are reported in Table 4. The percentage of difference with Fink's recommendation for pure UO 2 (y=0) is less than 1% for all remaining potentials, which denotes a good structural evolution of the potentials. Table 4: Parameters for the third degree polynomial that fits the evolution of the lattice parameters as the function of temperature and plutonium content. We also check if the lattice parameter evolves as the Vegard's law for stoichiometric MOX compounds up to the melting point. This behaviour has been verified experimentally up to 2000 K [START_REF] Popov | Thermophysical Properties of MOX and UO 2 Fuels Including the Effects of Irradiation[END_REF][START_REF]ESNII+ European commission. Deliverable D7.11 "state of the art with a literature review of MOX properties[END_REF]. The percentage of deviation from Vegard's law as function of temperature for compounds with different plutonium contents behaves similarly for all potentials. From 300 to ∼2500 K the Vegard's law is verified within ±0.02%. From 2500 K to the melting point, the Potashnikov and Cooper potentials significantly deviate from Vegard's law, whereas Arima potential seems to follow Vegard's law for the full temperature range. The discrepancy with Potashnikov and Cooper potentials is usually attributed to a high sublattice oxygen disorder (see next chapter for more details).

A [Å.K -3 ] B [Å.K -2 ] C [Å.K -1 ] D [Å] m [Å]

Thermal expansion coefficients

A more sensitive quantity to evaluate the structural property with temperature is the linear thermal expansion coefficient (LTEC). This coefficient is calculated as the first derivative of the lattice parameter with respect to the temperature using the following expression:

α(T ) = 1 a 0 ∂a ∂T P (4)
where, a 0 denotes, in our cas,e the lattice parameter at 300 K, and (∂a/∂T ) P is the first derivative at constant pressure calculated by numerical differentiation from the lattice parameter evolution. Fig. 2 displays the LTEC with Cooper (C ending), Potashnikov (P ending) and Arima (A ending) potentials for different plutonium contents (0<y<1). These values are compared to Uchida's experimental data for pure PuO 2 [START_REF] Uchida | Thermal expansion of PuO 2[END_REF] and Fink's recommendation for pure UO 2 up to 2000 K. As expected, the LTEC follows the experimental data and the recommendation from 300 to 2000 K. There is also no noticeable difference with the plutonium content for this range of temperature. After 2000 K, peaks appear for Potashnikov and Cooper potentials. For Cooper potential they depend on the plutonium content as previously mentioned [START_REF] Cooper | Thermophysical properties and oxygen transport in the (U x ,Pu 1-x )O 2 lattice[END_REF]. These peaks, referred to as λ-peak, were discovered and characterized by Bredig et al. for most of the fluorite-like structures [START_REF] Dworkin | Diffuse Transition and Melting in Fluorite and Anti-Fluorite Type of Compounds: Heat Content of Potassium Sulfide from 298 to 1260 K[END_REF] and can be attributed to a diffusing phase transition that turns compounds into a superionic conductor. Several reviews of experimental works [START_REF] Ralph | Specific Heat of UO 2 , Tho 2 , PuO 2 , and the Mixed Oxides (Th x U 1-x )O 2 , and (Pu 0.2 U 0.8 )0 1.97 by Enthalpy Data Analysis[END_REF][START_REF] Ronchi | Analysis of recent measurements of the heat capacity of uranium dioxide[END_REF] and theoretical works [START_REF] Potashnikov | High-precision molecular dynamics simulation of UO 2 -PuO 2 : Pair potentials comparison in UO 2[END_REF][START_REF] Cooper | Thermophysical properties and oxygen transport in the (U x ,Pu 1-x )O 2 lattice[END_REF][START_REF] Yakub | Molecular dynamics simulation of premelting and melting phase transitions in stoichiometric uranium dioxide[END_REF][START_REF] Lunev | A classical molecular dynamics study of the correlation between the Bredig transition and thermal conductivity of stoichiometric uranium dioxide[END_REF][START_REF] Korneva | Atomistic Simulation of a Superionic Transition in UO 2[END_REF] show that this transition occurs around 0.8T m , where T m is the melting temperature. It is often associated with premelting of the oxygen sublattice. For the Cooper potential, as already stated by the authors [START_REF] Cooper | Thermophysical properties and oxygen transport in the (U x ,Pu 1-x )O 2 lattice[END_REF], the λ-peak for pure UO 2 occurs around 2600 K in concordance with the experimental value of superionic transition temperature of 2670 K [START_REF] Hiernaut | Premelting Transition in Uranium Dioxide[END_REF]. For the case of pure PuO 2 , there is no experimental data available which could confirm the existence of a superionic transition. However, according to the theory of 0.8T m , it should occur around 2400 K since the melting temperature in PuO 2 was recently measured at 3040 K [START_REF] Böhler | Recent advances in the study of the UO 2 -PuO 2 phase diagram at high temperatures[END_REF]. Therefore, the superionic transition temperature of 2300 K found with the Cooper potential for PuO 2 is maybe slightly underestimated but still in the range. Moreover, λ-peaks decreases with the increase of the plutonium content, which is in agreement with previous studies [START_REF] Cooper | Thermophysical properties and oxygen transport in the (U x ,Pu 1-x )O 2 lattice[END_REF].

In the case of the Potashnikov potential, the values from 300 to 2000 K are slightly higher than the experimental data and the recommendation. The λ-peaks are present but less pronounced than for the Cooper potential. Due to the lack of statistics with our results, no clear difference appears with plutonium content. Finally, for the Arima potential, like for the two others, the LTEC behaves linearly up to 2000 K. At higher temperatures, although there are small oscillations, no clear λ-peak can be noticed for all plutonium content up to 3000 K. However, previous studies on pure UO 2 [START_REF] Potashnikov | High-precision molecular dynamics simulation of UO 2 -PuO 2 : Pair potentials comparison in UO 2[END_REF] show a clear λ-peak with the Arima potential around 3500 K. The melting point with this potential is found to be around 4500 K, which overestimates the experimental data. Nevertheless, it still follows the theory of the 0.8T m .

Overall, all potentials give a good estimation of the LTEC up to 2000 K. Differences between potentials appear at higher temperatures, i.e. in condition for which the experimental data for MOX and PuO 2 are lacking. It is then somehow difficult to choose the best potential, which reproduces most the real material. However, the behaviour according the superionic transition temperature can provide more indication. The Cooper poten-tial, gives a clear distinction between the superionic transition temperatures with the plutonium content and seems to follow best the theoretical λ-peak temperatures while the Potashnikov potential does not show significant change with the increase of the plutonium content. Finally, Arima potential reproduces no λ-peaks in the range of the temperatures studied.

Enthalpy and specific heat

We have calculated the enthalpy evolution as a function of temperature up to the melting point for Arima, Potashnikov, and Cooper potentials (not shown). The enthalpy values up to 2000 K are almost identical and fit perfectly the Fink's recommendation. Discrepancy appears at 2400 K, where Arima potential still provides a linear feature, while Fink's recommendation and both the Potashnikov and Cooper potentials deviate from linearity. For these two last potentials discrepancy with Fink's recommendation appears around 2600 K, which is the temperature of the λ-peak.

More noticeable changes can be extracted with the analysis of the specific heat. The specific heat capacity coefficients at constant pressure can be calculated directly by numerical differentiation from the enthalpy increment function with the following relationship:

C P = 1 n ∂H ∂T P ( 5
)
where n is the number of moles. All the evolutions of C P as a function of temperature for different plutonium contents are shown in Fig. 3. The temperature range investigated herein is above the Debye temperatures calculated with the different potentials, which are between 350 and 480 K. Overall, the behaviour is similar to that of the LTEC. However, λ-peaks appear more clearly than with the LTEC in the case of the Potashnikov potential even though no clear dependence with the plutonium content can be drawn. This demonstrates that this potential can also reproduce the superionic transition, as already been shown in the case of pure UO 2 [START_REF] Potashnikov | High-precision molecular dynamics simulation of UO 2 -PuO 2 : Pair potentials comparison in UO 2[END_REF]. On the other hand, Arima potential seems not to include this transition up to 3000 K.

Elastic constants

The assessment of the elastic constants as a function of temperature and plutonium content is important because the Potashnikov and Arima potentials were not fitted on those properties. Furthermore, the Cooper potentials has been fitted only with experimental data at room temperature. It is then an important checkpoint on the reliability of these potentials for further studies of the mechanical properties under irradiation. In addition, evolution of elastic properties with temperature in MOX is rather scarce in the literature, thus, atomic simulation could bring some insight in the subject.

The elastic constants, at finite temperature, are computed in a stepwise fashion. First, small deformations (0.1%) are imposed on the box for all the six Voigt components of the strain tensor ǫ i . This step is followed by an equilibration at constant temperature and constant volume for 10 ps. Elastic moduli were then The same colour code than for the LTEC is used. calculated as initial slopes of stress-strain curves obtained using appropriate components of stress and strain tensors. More specifically, the Young's modulus is deduced by applying tension and compression uniaxial strains individually and calculated as σ i /ǫ i , where σ i are the components of the stress tensor.

The results that we have obtained with MD simulations will be compared to the recommendation based on the review of Martin et al. [START_REF] Martin | The elastic constants of polycrystalline U0 2 and (U, Pu) mixed oxides: a review and recommendations[END_REF], who assesses experimental data from mainly polycrystalline materials. In this review, the authors concluded that the elastic constants of stoichiometric (U,Pu)O 2 should be taken identical to those of UO 2 . Also, recommendation shows only the evolution of Young's and shear moduli. In consequence, we will compare computed Young's (noted E) and shear (noted G) moduli calculated from elastic stiffness tensor coefficients using the Voigt-Reuss-Hill approximation for randomly oriented polycrystals [START_REF] Hill | The Elastic Behaviour of a Crystalline Aggregate[END_REF][START_REF] Mehl | Intermetallic Compounds: Principle and Practice[END_REF][START_REF] Schreiber | Elastic Constants and Their Measurement[END_REF]. The equations to transform the elastic stiffness coefficients found in monocrystal into shear and Young's polycrystalline moduli for a cubic system are presented below:

G polycrystal = (G Reuss + G Voigt )/2 with G Reuss = 5(C 11 -C 12 )C 44 [4C 44 + 3(C 11 -C 12 )] and G Voigt = (C 11 -C 12 + 3C 44 )/5 E polycrystal = -9BG polycrystal (3B + G polycrystal ) (6)
where B = (C 11 + 2C 12 )/3 is the bulk modulus.

Fig. 4(a), 4(b), and 4(c) show the evolution as a function of temperature and plutonium content of Young's, shear, and bulk moduli, respectively. All moduli for all compositions decrease with increasing temperature. Between 300 and 1600 K, the decrease seems to be linear for the three potentials, whereas at higher temperatures, they decrease more rapidly as per the recommendation. Moreover, the increment in the plutonium con-tent causes a moderate increase, of about 5%, for the Arima and Postashnikov potentials and a moderate decrease, of about 5%, for the Cooper potential for both the shear and Young's moduli. Even if the recommendation states that (U,Pu)O 2 elastic constants should be taken identical to those of UO 2 , a moderate increase (about 3%) with the addition of plutonium content is indicated by some authors [START_REF] De Novion | Mechanical Properties of Uranium Plutonium Based Ceramics[END_REF][START_REF] Nutt | Elastic and Anelastic Response of Polycrystalline UO 2 -PuO 2[END_REF].

Comparing the different potentials, it appears clearly that the Arima potential gives the highest elastic constant overestimating largely the recommendation for the Young's modulus. The Potashnikov potential underestimates the Young's and shear moduli at low temperatures but fits well the recommendations around 2400 K up to melting point. The Cooper potential fits really well the recommendations over the full range of temperature. We recall here that the Cooper potential has been designed to reproduce the elastic constant of experimental data for single crystal. Furthermore, the correction for polycrystalline materials turns to be a good estimate.

Another important criterion is the anisotropy factor that reflects the difference of bonding with the different directions of the crystal. In cubic system the anisotropy can be quantified by the Zener's anisotropy factor defined as follow:

Z = (2C 44 ) (C 11 -C 12 ) . ( 7
)
When Z is equal to 1, the system is perfectly isotropic, conversely a deviation from 1 shows that the system is anisotropic. We calculate the evolution for the three potentials of the Zener's factor as a function of the temperature for pure UO 2 and pure PuO 2 . The crystals produced with the Arima and Cooper potentials are anisotropic for the full range of temperature (Z ∼ 0.5), whereas with the Potashnikov potential is isotropic at lower temperatures and becomes anisotropic around 2400 K where the Bredig's transition takes place for this potential. Therefore, one can expect a strong influence of the strain direction on the mechanical properties.

Stress-strain curves

In order to assess the empirical potentials on the cracking behaviour of (U,Pu)O 2 solid solution, we carried out simulations to evaluate the stresses as a function of the strain for different plutonium contents and at different temperatures. The stressstrain curves are calculated by imposing an uniaxial deformation on the box with a constant strain rate (10 8 /s) until the system cracks in two. During the simulation the stress component corresponding to the direction of deformation is recorded while relaxing the other components of the stress tensor under NσT. In such a way, it is possible to investigate the anisotropy of the system.

We first check the stress-strain curves for the three main crystallographic directions of the fluorite crystal (i.e. 1 0 0 , 1 1 0 , and 1 1 1 ). The results for the three potentials at 300 K for pure UO 2 and pure PuO 2 are displayed in the Fig. 5. For all cases, the shapes of the stress-strain curves are consistent with a classical brittle fracture. Indeed, we observe a linear increase of the stress corresponding to the elastic deformation followed by an abrupt decrease indicating the loss of the crystal structure. The highest stress point is then defined as the ultimate tensile strength (σ UTS ) at which the system begins to crack. For all the potentials the stiffest direction appears clearly to be the 1 0 0 direction. The ultimate tensile strengths in the 1 1 0 and 1 1 1 directions are very close to each other in the case of the Arima and Potashnikov potentials. However, with the Cooper potential, σ UTS is slightly lower in the 1 1 1 direction. This in agreement with theoretical [START_REF] Zhang | Crack tip plasticity in single crystal UO 2 : Atomistic simulations[END_REF][START_REF] Fossati | Simulation of uranium dioxide polymorphs and their phase transitions[END_REF] and experimental results [START_REF] Rapperport | Deformation modes of single crystal uranium dioxide from 700 • C to 1900 • C[END_REF][START_REF] Robins | Uranium Oxide Cleavage[END_REF] that show that crack propagates mainly along the {1 1 1} planes. Therefore, for the rest of the investigation, we will only focus on the 1 1 1 direction.

The impact of the temperature on σ UTS obtained from the stress-strain curves ranging from 300 K to the melting point for pure UO 2 and pure PuO 2 for the Arima, Potashnikov, and Cooper potentials is shown in Fig. 6.

The ultimate tensile strength decreases almost linearly as the temperature increases for all potentials. Also, there is no remarkable difference concerning the plutonium content. This is in agreement with recommendations found in the European Commission state of art of MOX fuel report [START_REF]ESNII+ European commission. Deliverable D7.11 "state of the art with a literature review of MOX properties[END_REF] in which it stats that the yield and ultimate stress for MOX are on a precautionary basis the same as UO 2 , whatever the Pu content, and up to 1500 K. Our MD study suggests that this behaviour could be extrapolated to higher temperatures.

Another important thermo-mechanical property is the brittleto-ductile transition, which occurs at high temperature around half the melting point in both UO 2 and MOX fuel [START_REF] Roberts | Nature of brittle-to-ductile transition in UO 2 -20 wt% PuO 2 nuclear fuel[END_REF][START_REF] Roberts | Brittle fracture of oxide nuclear fuel[END_REF]. This transition is defined as the temperature (T C ) where yield stress and ultimate strength deviate from one another. Namely, below T C fuel behaviour is brittle whereas above T C fuel behaviour exhibits some plastic features. This temperature can be obtained as well from the stress-strain curves. In our simulations, yield stress is defined as the boundary in between elastic, reversible strain, and plastic, irreversible deformation using the following procedure. Regularly during the uniaxial deformation under tensile mode, the system is unloaded along the same direction until zero stress is reached. Then, the point where the resulting strain of the unloaded system is not equal to zero is taken as the yield stress. Hence, we can define a ratio between the yield stress and σ UTS , which equals 1 at low temperature where the yield stress is combined with σ UTS . The results of this ratio as a function of the temperature is reported in Fig. 7. For all the potentials, we find the same behaviour. As expected, from low temperatures up to around 1500 K for the Copper potential and around 1900 K for the Potashnikov and Arima potentials, yield stress and σ UTS are the same; the ratio is equal to 1. Above these temperatures, the ratio decreases rapidly until a plateau is reached. The transition between these two regimes is the brittle-to-ductile temperature. First, all potentials are able to reproduce this brittle-to-ductile behaviour. Cooper potential gives T C very close to the experimental value (∼1673 K [START_REF] Roberts | Nature of brittle-to-ductile transition in UO 2 -20 wt% PuO 2 nuclear fuel[END_REF][START_REF] Roberts | Brittle fracture of oxide nuclear fuel[END_REF]) while the others overestimate it by about 300 K. After the transition, the ratio stabilizes at 0.85, 0.80, and 0.75 for the Potashnikov, Arima, and Cooper potentials, respectively. This indicates that with the Cooper potential the material is inclined to experience plastic deformation for lower temperatures than with the other potentials.

However, conclusion about the reliability of the potentials based on this property needs to be taken with care. Indeed, the experimental values for the brittle-to-ductile transition are obtained from polycrystal samples and under strain rates far lower than the one of MD simulations. The physical phenomena underneath this transition are rather complex involving thermal creep which can not be captured by MD simulations. Therefore, the accuracy between MD and experiment could be a coincidence. Nonetheless, it still provides an indication on the behaviour at high temperatures.

To complete the assessment on the mechanical properties, we investigate crack propagation behaviour. The crack propagation simulations require large systems in order to include an initial crack and enough material to analyse its propagation. We use the thin strip geometry, where a constant strain is applied perpendicularly to the initial crack. The advantage of this geometry is that the energy release rate does not depend on the crack length, and can be found analytically by considering the energetics of an advancing crack. This is applicable if the system length (x-direction) to height (y-direction) ratio is at least 4 (i.e., L x ≥ 4L y ) [START_REF] Buehler | Hyperelasticity governs dynamic fracture at a critical length scale[END_REF][START_REF] Knauss | Stresses in an Infinite Strip Containing a Semi-Infinite Crack[END_REF]. Therefore, we use a system including 4 × 10 6 atoms with initial box size roughly equal to 240 × 60 × 4 nm in the x, y, and z directions respectively. The initial notch is designed as an ellipse to ensure maximal stress concentration at the crack tip and equals 40×10 nm in the x and y directions respectively. Such simulations are computationally expensive due to the size of the system. Therefore, only one simulation per interatomic potential is carried out. We choose a system with 25 at% of plutonium at 300 K and the load is applied along the 1 1 1 crystal orientation, which is considered the weakest. Loading is applied with constant strain-rate (10 8 /s) until complete cracking of the system occurs. During the simulation, the stress tensor is recorded.

The mechanism of crack propagation can be determined by analysing the snapshots during the simulations. All the figures of these snapshots are shown in Fig. 8. Clear differences appear between the different potentials. For the Potashnikov potential, the crack propagates classically with cleavage-like behaviour, the crack opens straight with steady velocity. This behaviour is expected for a pure brittle material. Conversely, for the Arima and Cooper potentials the crack propagates at the interface or within a secondary phase that forms ahead of the crack tip. One can also observe small cavities forming which grow and coalesce with the main crack. This secondary phase (marked in grey colour in Fig. 8) is of rutile-like structure. In the case of Cooper potential, it covers relatively a large zone before crack actually opens up. This type of propagation denotes an unusual plastic-like behaviour at low temperature. However, the high strain-rate used here could cause to overestimate the stress field necessary for crack propagation and then the volume affected with the secondary phase. But, this phase transition ahead of crack tip has already been observed with MD simulations for pure UO 2 with different interatomic potentials [START_REF] Zhang | Crack tip plasticity in single crystal UO 2 : Atomistic simulations[END_REF][START_REF] Fossati | Simulation of uranium dioxide polymorphs and their phase transitions[END_REF].

Therefore, it is important for the potential assessment to check the relative stability of these phases to ensure that fluorite structure is the most stable phase for (U,Pu)O 2 compounds. For this reason, we have undertaken energy-volume studies (not shown) with five different structures for both pure UO 2 and pure PuO 2 . These structures are of type: fluorite (Fm 3m), cotunnite (Pnma), rutile (P4 2 /mnm), scrutinyite (Pbcn), and marcasite (Pnnm). This study follows previous works in pure UO 2 [START_REF] Chernatynskiy | Critical assessment of UO 2 classical potentials for thermal conductivity calculations[END_REF][START_REF] Zhang | Crack tip plasticity in single crystal UO 2 : Atomistic simulations[END_REF][START_REF] Fossati | Simulation of uranium dioxide polymorphs and their phase transitions[END_REF] carried out with other empirical potentials and DFT calculations. We confirm that for both compounds the fluorite structure corresponds to the ground state for all three empirical potentials studied. As previously demonstrated with DFT calculations [START_REF] Zhang | Crack tip plasticity in single crystal UO 2 : Atomistic simulations[END_REF][START_REF] Fossati | Simulation of uranium dioxide polymorphs and their phase transitions[END_REF][START_REF] Desai | Stress-induced phase transformation in nanocrystalline UO 2[END_REF], we also find that cotunnite structure is the most stable at high isotropic pressure (compression) and that either rutile or scrutinyite structures are the most stable under negative isotropic pressure (tensile load in all the three directions). These two phases are almost energetically degenerate, hence it is impossible to distinguish which phase is the most favorable. Using the common tangent method, we calculate the transition pressures from fluorite structure to the other phases. These transition pressures are presented in Table 5 for both pure UO 2 and pure PuO 2 . For all the potentials and pure UO 2 , the lowest transition pressure under tensile loading is found for the transition from fluorite to rutile structure, which is coherent with results of the crack propagation. The stress-strain curves calculated during the crack propagation are displayed in Fig. 9. The strain at which crack starts to propagate is different for the three potentials. Crack initiates first with the Potashnikov potential at around 4% followed by Arima at 5% and finally Cooper at 6%. However, the corresponding σ UTS is the highest for Arima potential at 10 GPa followed by Cooper at 8 GPa, and Potashnikov at 5.5 GPa. These σ UTS are much lower than for the bulk case due to the presence of the initial crack, which concentrates the stresses at the crack tip. These values can also be related to the lower transition pressures in Table 5. For the Arima and Cooper potentials, σ UTS is greater than their relative transition pressure from fluorite to rutile structure. Therefore, the transition may occur at the crack tip where the stresses are the highest, explaining the plastic-like behaviour observed. However, it is noteworthy that these secondary phases are less stable than the fluorite structure and disappear behind the crack front as the crack advances. This could explain why this secondary phase cannot be directly observed experimentally. However, discrepancies between Griffith's theory and experimental results could be described by this behaviour [START_REF] Zhang | Crack tip plasticity in single crystal UO 2 : Atomistic simulations[END_REF]. Therefore, it is difficult to draw a conclusion about the accuracy of the empirical potentials for crack propagation.

Conclusions

In this paper, we assess empirical potentials for the (U 1-y ,Pu y )O 2 solid solution. To date, only empirical potentials using rigid ion model are available. Since we are interested in studying for our future study on the mechanical behaviour under irradiation the point defects distribution, both cations need to be explicitly modeled. Therefore, we found in the literature five interatomic potentials fulfilling these requirements coined by the name for their first author: Yamada, Arima, Potashnikov, Tiwary, and Cooper. In our assessment, the structural (lattice parameter, relative phase stability) and thermodynamics (thermal expansion, Heat capacity) properties are systematically calculated for the full range of temperature from 300 K to melting point, and for the full range of plutonium content from pure UO 2 to pure PuO 2 . We also investigate the potentials through their mechanical properties (elastic and crack propagation). Thus, this assessment includes ranges of temperatures and compositions as well as properties that have not been studied by the authors.

Tiwary potential turns quickly to be instable with MD simulation. Namely, the fluorite structure collapses after a few steps of simulation. Therefore, we eliminate this potential from our study. Yamada potential shows large discrepancies on the lattice parameter at high temperatures (> 2100 K) with Fink's recommendation. Therefore, we also eliminate this potential from the rest of the assessment. For the three potentials remaining (i.e. Arima, Potashnikov, and Cooper), thermal expansion and heat capacity show good agreement up to 2000 K. Nevertheless, at higher temperature with the Potashnikov and Cooper potentials so called λ-peaks appear, whereas no clear λ-peaks appears with the Arima potential up to 3000 K. These peaks are usually associated with a superionic transition corresponding to the premelting of the oxygen atoms sublattice and have been observed experimentally around 0.8 of the melting point in UO 2 . The assessment shows as well that temperatures and the behaviour with plutonium content for the superionic transition seem better reproduced with the Cooper potential.

The mechanical elastic properties are also investigated as a function of temperature and plutonium content. The results show clearly that the elastic stiffness constants are best reproduced with the Cooper potential, which has been fitted on the elastic constant of experimental data for single crystal. The Potashnikov potential gives fairly good agreement with experimental data while Arima potential overestimates largely the elastic stiffness constants. However, anisotropy is present even at low temperature for the Cooper and Arima potentials whereas it appears only at high temperatures for Potashnikov potential. Hence, Cooper potential appears to be the best potential to reproduce the mechanical elastic properties.

Analysis of stress-strain curves obtained with uniaxial loading shows that the 1 1 1 crystallographic direction gives lowest ultimate tensile strength, in agreement with experimental observations. We also find a brittle-to-ductile transition for the three potentials with transition temperature in good agreement with experimental values for Cooper potential and relatively close for the two others. Furthermore, all these mechanical properties show little dependence on the plutonium content, confirming the assumption that mechanical properties of UO 2 can be, to a large extent, applied to MOX.

Finally, the behaviour during crack propagation simulations is very different between the three potentials. For the Cooper and Arima potentials the crack propagates through secondary phase of rutile-like structure that appears ahead of the crack tip leading to an unexpected plastic-like behaviour. Conversely, for Potashnikov potential crack propagates by cleavage, which is typical of a brittle-like behaviour. However, it is strenuous to conclude which potential reproduces best the reality since no direct experimental observation is available.

Overall, with the structural, thermodynamics, and mechanical properties assessment the Cooper interatomic potential reproduces the best the Fink's recommendation, yet it renders an unexpected plastic-like behaviour during crack propagation. The Potashnikov potential gives fairly good agreement for structural, thermodynamics, and elasticity properties. It also presents expected brittle behaviour during crack propagation. Finally, the Arima potential gives good results for structural and thermodynamics properties under 2100 K, but presents discrepancies at high temperatures and gives average results for the mechanical properties.

Figure 1 :

 1 Figure 1: Evolution of the lattice parameter as a function of the temperature for the four potentials studied for pure PuO 2 . The results are compared with Yamashita's experimental work (Yama 97 [25]).

Arima 1 .

 1 03 × 10 -12 3.89 × 10 -9 4.99 × 10 -5 5.450 -0.074 Potashnikov 1.80 × 10 -12 1.82 × 10 -9 6.12 × 10 -5 5.441 -0.065 Cooper 1.90 × 10 -12 2.82 × 10 -9 5.16 × 10 -5 5.453 -0.072

Figure 2 :

 2 Figure 2: Evolution of the linear thermal expansion coefficient obtained for different plutonium contents as a function of temperature (the value of y indicates the plutonium content, (U 1-y Pu y )O 2 ). Arima potential (A ending) is referred with blueish colours and star symbols, Potashnikov potential (P ending) is referred with greenish colours and triangle symbols, and Cooper potential (C ending) is referred with reddish colours and square symbols.

Figure 3 :

 3 Figure 3: Evolution of the C P as a function of the temperature for different plutonium contents (the value of y indicates the plutonium content, (U 1-y Pu y )O 2 ). The same colour code than for the LTEC is used.

Figure 4 :

 4 Figure 4: Evolution of Young's (a) shear (b) and bulk (c) polycrystalline moduli as a function temperature and plutonium content calculated with Arima, Potashnikov, and Cooper potentials (the value of y indicates the plutonium content, (U 1-y Pu y )O 2 ). (a) and (b) are compared with the recommendations in [26] and (c) is compared with experimental date of Hutchings et al. [48].

Figure 5 :

 5 Figure 5: Stress-strain curves for the three interatomic potentials ((a) Arima, (b) Potashnikov, and (c) Cooper) at 300 K strained in the main three crystallographic directions of the fluorite crystal (i.e. 1 0 0 , 1 1 0 , and 1 1 1 ). For each potential, it is shown the results for pure UO 2 and pure PuO 2 .

Figure 6 :

 6 Figure 6: Evolution of the ultimate tensile strength as a function of the temperature and plutonium content (the value of y indicates the plutonium content, (U 1-y Pu y )O 2 ).

Figure 7 :

 7 Figure 7: Ratio between the Yield stress and the ultimate tensile strength (σ Yield /σ UTS ) as a function of temperature. These values are compared with experimental data found in [51].

Figure 8 :

 8 Figure 8: Snapshots during crack propagation in (U,Pu)O 2 system with 25 at.% of plutonium at 300 K for (a) Arima, (b) Potashnikov, and (c) Cooper potentials. The load is applied along the 1 1 1 crystallographic direction in the fluorite structure. The inserts in each figure depict a close up of the crack front where phase transition occurs.

Figure 9 :

 9 Figure 9: Stress-strain curves for the three interatomic potentials during crack propagation in (U,Pu)O 2 system with 25 at.% of plutonium at 300 K. The load is applied along the 1 1 1 crystallographic direction in the fluorite structure.

Table 1 :

 1 Type of physical function for the different interatomic potentials studied.

Table 2 :

 2 Parameters for the interatomic potentials.

				Yamada	Arima	Potashnikov	Cooper
			U	2.4	2.7	2.74492	2.2208
	q [e]		Pu	2.4	2.7	2.74492	2.2208
			O	-1.2	-1.35	-1.37246	-1.1104
	A [eV]	U-U U-Pu 1752.102 7.83068×10 13 442.161 2.48128×10 13 Pu-Pu 32606.8 2.80460×10 14 U-O 1018.46 55892.6	0 --873.107	18600 18600 18600 448.779
			Pu-O	5329.83	57425.2	871.79	527.516
			O-O	2345.9	978.718	50211.7	830.283
			U-U	0.32	0.072	-	0.2747
			U-Pu	0.24	0.0685	-	0.2691
	ρ [Å]	Pu-Pu U-O	0.16 0.32	0.065 0.202	-0.3592	0.2637 0.387758
			Pu-O	0.24	0.1985	0.3561	0.379344
			O-O	0.32	0.332	0.18446	0.352856
	C [eV.Å 6 ]		O-O	4.146	17.3544	74.7961	3.884372
				Yamada Arima Potashnikov Cooper
	D [eV]		U-O Pu-O 0.564005 0.78093	--	--	0.6608 0.70185
	γ [Å -1 ]	U-O Pu-O	1.25 1.56	--	--	2.05815 1.98008
	r 0 [Å]		U-O Pu-O	2.369 2.339	--	--	2.38051 2.34591

Table 3 :

 3 Parameters for Cooper EAM term.

		G α [eV.Å 1.5 ]	η β [Å 5 ]
	U	1.806	3450.995
	Pu	2.168	3980.058
	O	0.69	106.856

Table 5 :

 5 Transition pressure at T = 0 K from fluorite structure to secondary phases in GPa for the three interatomic potentials for pure UO 2 and pure PuO 2 . The negative sign denotes tensile loading.

					fluorite (Fm 3m) to
		Cotunnite	Rutile	Scrutinyite	Marcasite
		Pnma UO 2 PuO 2 UO 2 PuO 2 UO 2 PuO 2 UO 2 PuO 2 Pbcn Pnnm P4 2 /mnm
	Arima	44.3 106.9 -8.3	-9.3	-9.6	-7.4	-	-
	Potashnikov 16.0 19.6 -10.6 -11.8 -12.1 -13.3	-	-13.5
	Cooper	3.0	4.0	-6.8 -10.3 -8.0	-8.2 -16.0 -15.8

Acknowledgements

This work was granted access to the HPC resources of [TGCC] under the allocation 2016-mtt7073 made by GENCI. This research contributes to the joint programme on nuclear materials (JPNM) of the European energy research alliance (EERA).