Hector Balboa

Barthe

teaching and research institutions in France or abroad, or from public or private research centers.

Simulation of thermomechanical properties of U-PuO2

nuclear fuel under irradiation

few microsecond. In order to investigate processes, such as, cation difusion and rare-event annihilation of defects after cascades, another computational tool is required. Atomistic or object kinetic Monte Carlo (KMC) techniques can run for longer timescales than MD. However, for KMC to work accurately, all of the possible inter-state transition mechanisms and their associated rates need to be known a priori. For this reason, the adaptive kinetic Monte Carlo (AKMC) is chosen to overcome these limitations. This method determines the available transition states during simulation. In this way, it takes the system into unforeseen states via complex mechanisms. The power and range of this tool proved to be eicient to discover cation Frenkel pair recombination over a longer periods of time than MD.

Résumé en français

LŠobjectif de cette thèse a été dŠutiliser une approche numérique pour étudier lŠimpact des dommages dŠirradiation sur la microstructure du combustible nucléaire composé dŠun mélange dŠoxyde de plutonium et dŠuranium (MOX). Cette approche, réalisée à lŠéchelle atomique, repose sur lŠutilisation de potentiels empiriques développés dans la littérature pour lŠoxyde mixte (U,Pu)O 2 dans lŠapproximation des interactions dŠions rigides.

Une première étape a été une analyse critique des propriétés structurales, thermodynamiques et mécaniques prédites par 5 potentiels de la littérature. Les calculs de dynamique moléculaire ont été réalisés à lŠaide du code LAMMPS, sur lŠensemble de la gamme de composition de UO 2 à PuO 2 et à des températures comprises entre 300 K et le point de fusion. Nous montrons que les potentiels les plus satisfaisants sont ceux développés par Cooper et Potashnikov. Ces deux potentiels reproduisent correctement la stabilité thermodynamique des phases ainsi que lŠévolution en température des paramètres de maille et de la chaleur spéciĄque. Cependant le comportement mécanique des oxydes est diférent selon le potentiel choisi. Tout dŠabord, on note que les constantes élastiques et le facteur dŠanisotropie obtenus par le potentiel Cooper sont plus Ądèles aux recommandations de la littérature. Ensuite, nous montons que la propagation dŠune fracture induit une transformation de phase en pointe de Ąssure avec le potentiel de Cooper alors quŠune Ąssuration fragile est observée avec le potentiel Potashnikov.

Une seconde étape a été lŠétude des dommages dŠirradiation dans les oxydes mixtes en se limitant à lŠutilisation des deux potentiels de Cooper et de Potashnikov. Des calculs de dynamique moléculaire de cascades de déplacement ont été réalisés à diférentes énergies, température et compositions en plutonium. Ensuite, aĄn dŠévaluer les dommages dŠirradiation sur des temps plus longs que ceux accessibles en dynamique moléculaire, la méthode dŠaccumulation de défauts a été utilisée. Pour les deux potentiels, lŠévolution des dommages primaires avec lŠaugmentation de la dose suit globalement les mêmes étapes que celles trouvées précédemment dans UO 2 . Tout dŠabord, les défauts ponctuels sont créés. Par la suite, ils se regroupent et forment de petites boucles de Frank, qui après se transforment en dislocations parfaites. Cependant, la cinétique de la recombinaison de défauts ponctuels est signiĄcativement plus lente avec le potentiel de Cooper ce qui conduit à la création de petites régions désordonnées pour les cascades dŠénergie élevée. LŠefet de la teneur en plutonium est analysé. Nous montrons en particulier que la densité de dislocations crée diminue lorsque la teneur en plutonium augmente.

Bien que la dynamique moléculaire a été décrite comme un microscope moléculaire en raison de sa capacité à décrire avec précision des systèmes atomiques, elle présente un inconvénient majeur, celui lié aux temps de lŠordre de la femto-seconde nécessaires pour résoudre les vibrations atomiques. Cela limite le temps total de simulation approximativement quelques nanosecondes. AĄn de simuler les processus, tels que la difusion de cations, un autre outil de calcul est nécessaire. Les techniques de Monte Carlo (KMC) atomiques peuvent simuler de plus longtemps que la dynamique moléculaire . Cependant, pour que KMC fonctionne avec précision, il est nécessaire de connaître a priori tous les mécanismes de transition entre les états possibles. Pour cette raison, la methode de Monte Carlo cinétique adaptative (aKMC) est choisie pour surmonter ces limitations. Cette méthode détermine les états de transition disponibles pendant la simulation. De cette manière, elle entraîne le système dans des états imprévus via des mécanismes complexes. La puissance de cet outil sŠest révélée suisamment grande pour étudier la recombinaison de cations sur de plus longues périodes de temps que la DM.

iii

Avant-propos Introduction

The world is more than ever questioning the use and development of the diferent types of energy, especially in the context of global warming, which represents a growing danger for us, the human kind. From all the alternatives, nuclear power seems to be one one of the most reliable energy to achieve our goal of a world free of human CO 2 emissions [START_REF] Shellenberger | Had they bet on nuclear, not renewables, germany & california would already have 100% clean power[END_REF]. However, nuclear power has always been under scrutiny, due to the impact it may have on the environment. This fear comes from accidents such as Chernobyl or Fukushima, or its military use linked to atomic bombs. MOX nuclear fuel cycle has as the advantage of the no proliferation, compared to the uranium cycle. Both cycle are still under debate for choosing the most adequate to the future worldwide need of electricity. MOX fuel cycle would be chosen if the uranium prices increases or/and if the important criteria is the waste disposal. Indeed, it produces less intermediate and high level waste. Last but not least important, MOX is planned to be the fuel for the new generation of reactors, such as ASTRID.

The MOX fuel is submitted to many micro-structural changes while being inside the reactor. One example is the restructuring of the solid matrix. It is due to the large temperatures and temperature gradients. Voids migrate towards higher temperatures and gather in the centre. Conversely, this causes a movement of the solid towards the outer periphery and leads to the creation of a void in the centre. Comparing MOX and UO 2 fuels, we Ąnd that in MOX (mixed oxides) there are higher rates of Ąssion gas release, due to higher linear power levels at the end of the irradiation and higher centre-line temperature (caused by slightly lower thermal conductivity). In general, inside the MOX fuel many complex phenomena take place, leaving more questions than answers. SpeciĄcally, the behaviour at atomic scale is not well understood. Moreover, there is a lack of information related to the entire Pu content domain and temperatures.

An important issue for the future of nuclear power is to ensure safeness and efectiveness 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 eforts gathered in the following reviews [START_REF] Popov | Thermophysical properties of mox and uo2 fuels including the effects of irradiation[END_REF] and [START_REF] Esnii+ | State of the art with a literature review of mox properties[END_REF] , experiments with MOX are diicult 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 uo2 and puo2 by equilibrium molecular dynamics simulations from 300 to 2000k[END_REF], (Arima, Yamasaki, Inagaki, v & Idemitsu, 2006), [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 u0.75pu0.25o2x[END_REF], [START_REF] Cooper | Modelling the thermal conductivity of uxth1xo2 and uxpu1xo2[END_REF] and [START_REF] Li | Molecular dynamics study of thermal conductivities of u0.7xpu0.3amxo2[END_REF]. However, the reliability of the results depends exclusively on the choice of the set of potentials. The potential parameters are usually Ątted 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 have not been Ątted on 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 uo2. part i: static calculations[END_REF], [START_REF] Govers | Comparison of interatomic potentials for uo2: part ii: molecular dynamics simulations[END_REF], [START_REF] Potashnikov | High-precision molecular dynamics simulation of uo2puo2: pair potentials comparison in uo2[END_REF] and [START_REF] Chernatynskiy | Critical assessment of uo2 classical potentials for thermal conductivity calculations[END_REF] but, to our knowledge, not yet for MOX. Therefore, in Chapter 3, we assess available rigid ion model empirical potentials for MOX on the structural, thermodynamics, and mechanical properties. The assessment (Chapter 3) 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. We will present Ąrst the the results related to the thermodynamical properties: 1) lattice parameter, 2) the thermal expansion, 3) the speciĄc heat capacity, 4) a bibliography research about the thermal conductivity, melting point and difusion, 5) and phase stability. Secondly, we present the mechanical properties: 1) the elastic constants, 2) the anisotropy factor, 3) the stress-strain curves under uniaxial deformation, 4) the brittle-to-ductile transition and 5) crack propagation.

It is known that MOX thermomechanical properties vary with respect to irradiation [START_REF] Popov | Thermophysical properties of mox and uo2 fuels including the effects of irradiation[END_REF], [START_REF] Carbajo | A review of the thermophysical properties of mox and uo2 fuels[END_REF] and [START_REF] Esnii+ | State of the art with a literature review of mox properties[END_REF], yet, detailed description of this phenomena is rather sparse due to again the great diiculties of carrying out experiments in such materials (high toxicity and radiation levels). MD has proved to be a useful tool to simulate ballistic damage in UO 2 by means of displacement cascades studies [START_REF] Van Brutzel | Displacement cascade initiated with the realistic energy of the recoil nucleus in uo2 matrix by molecular dynamics simulation[END_REF], (Devanathan, Yu, & Weber, 2009), [START_REF] Martin | Efect of the cascade energy on defect production in uranium dioxide[END_REF] and [START_REF] Crocombette | Molecular dynamics simulations of high energy cascade in ordered alloys: defect production and subcascade division[END_REF]. Moreover, MD simulations using the Frenkel pair accumulation (FPA) technique (Chartier, Meis, Crocombette, Weber, & Corrales, 2005, 2) and [START_REF] Crocombette | Atomistic simulation of amorphization thermokinetics in lanthanum pyrozirconate[END_REF] allows us to go further in picturing the dose efect and have been tested for diferent oxides (Chartier, Catillon, & Crocombette, 2009, 15) and UO 2 [START_REF] Chartier | Early stages of irradiation induced dislocations in urania[END_REF]. Therefore, we will present in Chapter 4 MD ballistic radiation damage characterization of (U,Pu)O 2 solid solution, for several Pu contents, temperatures, and radiation doses. This assessment is shown in Chapter 4 and is divided as follow: the Ąrst part provides some point defect formation energies in order to assess the empirical potentials. The second part examines the results on the lifetimes of Frenkel pairs, which are of importance for the two next following sections. The third part presents the primary damage state obtained via displacement cascades. The fourth part discusses the efect of irradiation dose on the evolution of the microstructure using the Frenkel pair accumulation method. The Ąfth part shows the elastic moduli as a function of the dose. Finally, we assess the simulated XRD (X-ray difraction) patterns as a function of dose.

Generally, the MD approach that solves the motion of atoms using classic mechanics has proven to be successful for a large number of cases. However, we recall that the physical time MD deals with is not suicient to study many important phenomena such as cation difusion and the reaction of some chemical substances, this makes it necessary to Ąnd a way to solve the problem of systems evolving slowly over time. As a solution, we propose to use the aKMC (Adaptive Kinetic Monte-Carlo). Chapter 5 will present the basics of the aKMC method. An example of its application will be illustrated by results obtained over a large period of time for cation Frenkel pair recombination.

Publications and License Information

The author declares that only the personal contributions are used in this thesis. The present PhD work has lead to the publication on the following journal papers:

• H. Balboa, L. Van The present PhD thesis is written with T E X.

Licensees may copy, distribute, display and perform the work and make derivative works and remixes based on it only if they give the author the credits (attribution). The following BibTeX code can be used to cite the current document: @PhdThesis{Balboa:2018, author = {Balboa-Lopez, Hector E.}, title = {{S}imulation of {T}hermo-{M}echanical {P}roperties of (U, Pu)O2 {N}uclear {F}uel under {I}rradiation {U}sing {A}tomistic {S}imulations}, school = {Université Paris-Saclay}, year = {2018}, month = dec, } Interested readers can freely use or adapt the document structure, the title page, etc., to their own needs. 1.6 Microstructure of a fast reactor pelletized MOX fuel pellet [START_REF] Welland | Matter transport in fast reactor fuels. Institute for Transuranium Elements[END_REF]. MOX fuel pellets are often divided into regions deĄned by their restructured state. Going outward from the centreline we Ąnd 1) Central void, 2) Columnar grain region, 3) Equiaxed grain growth and 4) As-fabricated microstructure. 6

List of Figures

1.7 Microprobe mappings of MIMAS MOX. (CEA, 2009). The colour level corresponds to wt % of PuO Edge dislocation, c) Self interstitial atom, d) Vacancy, e) Precipitate of impurity atoms, f) Vacancy type dislocation loop, g) Interstitial type dislocation loop, h) Substitutional impurity atom [START_REF] Foll | Defects in crystals[END_REF]. . . . . . . . . . . . . . . . . 12 1.13 Schematic representation of stalking faults [START_REF] Foll | Defects in crystals[END_REF]. . . . . . . . . . . 13

1.14 Schematic representation of Shockley VS Frank dislocations [START_REF] Radwan | Dislocations in fcc metals[END_REF]. 14

1.15 Schematic representation of a unfaulting process involving a Frank, perfect and Shockley dislocations [START_REF] Foll | Defects in crystals[END_REF] [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with ovito the open visualization tool[END_REF].

2.4 Division of the space by Voronoi cells . . . . . . . . . . . . . . . . . . . . .

2.5

Atomic volume vs count of atoms histogram from the Voronoi cell analysis. .

2.6

Example of how the Burgers vector is created [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with ovito the open visualization tool[END_REF]. . . . . . .

2.7

Construction of dislocation lines. [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with ovito the open visualization tool[END_REF]. . . . . . . . . . . . .

2.8

Graphic representation of the efect of strain in the crystalline structure and its efect on the difraction peaks. a) No strain b) uniform strain c) non-uniform strain is applied. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.1

Evolution with temperature of the lattice parameter for the four potentials studied for pure UO 2 . The results are compared with FinkŠs recommendation [START_REF] Fink | Thermophysical properties of uranium dioxide[END_REF] and experimental work [START_REF] Yamashita | Thermal expansions of npo2 and some other actinide dioxides[END_REF] 3.2 Evolution with temperature of the lattice parameter for the four potentials studied for pure PuO 2 . The results are compared with [START_REF] Yamashita | Thermal expansions of npo2 and some other actinide dioxides[END_REF] The system is then made to go down the minimum energy path (gray lines) on either side of the saddle points which are indicated with *. Here, all four saddle points have a minimum energy path starting at the initial state minimum A, but this does not have to be the case. The rate of each process is then calculated using Equation (5.13). A process is chosen at random using the kinetic Monte Carlo algorithm. In this case (see section Section 5.6), process 1 gets chosen. The system is moved to the Ąnal state of this process, to minimum B. (c) Dimer searches are run from the new minimum, again four distinct saddle points are found. (d) Minimum energy paths are traced out, and the process repeated [START_REF] Henkelman | Long time scale kinetic monte carlo simulations without lattice approximation and predeĄned event table[END_REF] 5.4 Reaction mechanisms with barriers within mk B T of the lowest saddle point energy are considered relevant [START_REF] Xu | [END_REF]. . . . . . . . . . .

5.

5 Force evaluations vs system size for three diferent Al system sizes [START_REF] Xu | [END_REF] Along this chapter we will go through general knowledge about the MOX fuel. Firstly, we will explain what it is used for in Section 1.1. Secondly, we will compare the closed fuel cycle with the open fuel cycle, in order to understand the importance of the MOX fuel in the future of the nuclear power. After, we will describe the behaviour of the fuel inside the reactor and the impact of irradiation. Then, a brief description of the Ćuorite structure and types of defects is included.

Generation of Electricity with Nuclear Power

Nuclear power has played an important role in the production of electricity since its beginnings. It supplies around 30% of the global demand. In Europe, this is especially important, for example, in France it covers more than 75% of their demand. Other countries in Europe, such as, England are still betting on this energy by constructing new reactors.

The nuclear fuel cycle has always been under scrutiny, due to the impact it may have on the environment. This fear comes from accidents such as Chernobyl or Fukushima. Radioactive emissions released in such accidents, can travel around the globe, and therefore cause a vast environmental problem. Moreover, the lack of governmental decisions about the radioactive waste has created a bad reputation for nuclear, worldwide. However, discounting these few major accidents, nuclear energy is important producing low-carbon electricity. This is rather important considering the growing environmental problems related to the global warming. Nuclear power, with its almost non-existent carbon emission is starting again to play an important role in saving the world as we know it.

The generation of nuclear power uses mainly as fuel U 235 . Spontaneous Ąssion of the uranium atoms takes place inside the fuel. This gives energy and releases energetic neutrons. Some of the released neutrons cause further Ąssion processes, releasing more energy. The remaining neutrons are absorbed by control mechanisms or other isotopes, such as, U 238 . Further decays leads to the production of new isotopes, such as, Pu 239 . In certain cases, more Ąssile materials are produced. The energy released from each Ąssion process is extremely energetic (around 200 × 10 6 eV). This great amount of energy, in a small volume, makes this type of fuel one of the most energetic. Part of this energy is transformed in heat. There is a heat exchange between the reactor and the cooling system. The heat is then used to transform water into steam, in order to move turbines and produce electricity.

Why using MOX?

Due to the large number of factors entering the speciĄcations of nuclear fuel, several nuclear fuels have been proposed. The nuclear fuel cycle represents the progression of nuclear fuel from creation to disposal. All of them share the same processes, from mining to fuel being irradiated in the reactor. However, the key factor that diferentiate them starts after the nuclear reactor operation. Some of them consider the option of that the irradiated fuel is reprocessed and the nuclear material is recycled. Choosing the cycle depends on the type of reactor as well (IAEA, 2009).

Mainly, two types of cycles have major acceptance around the globe. The Ąrst one, called ŞThe open fuel cycleŤ, starts by irradiating the nuclear fuel in the reactor. Secondly, after irradiation, the fuel is maintained in ponds in the same nuclear plant for years up to decades until it is sent away (for instance in Spain at ŞEl CabrilŤ ) or it remains in the same ponds awaiting for decisions. This nuclear spent fuel is meant to be conditioned or just packed in special containers and put it into a Ąnal repository (IAEA, 2009). It is important to note, that this fuel cycle strategy is the most popular around the globe. However, some countries that were using this option now are looking forward to re-utilise this valuable material (uranium and plutonium). Another disadvantage is that no Ąnal repository for spent fuel has yet been constructed. The Ąrst will be the Finnish repository at Onkalo. Figure 1.1 shows a schematic representation of this cycle.

The second cycle one is called Şclosed fuel cycleŤ. After the nuclear fuel was irradiated in a reactor and kept cooling during a suicient period of years, the spent fuel is reprocessed to extract the remaining uranium and plutonium (Ąssile and Ąssionable materials), separating the Ąssion products and other actinides. After the reprocessing process, from the original spent fuel 95% of the mass is reprocessed uranium, 4% is nitric acid solution containing Ąssion products and transuranic wastes (TRU), the remaining 1% corresponds to plutonium. It is convenient to remark that a small portion of the radioactivity is released into the atmosphere or into liquid wastes from the reprocessing plant [START_REF] Bunn | The economics of reprocessing versus direct disposal of spent nuclear fuel[END_REF]. The reprocessed uranium and plutonium are then reused in the reactors. This type of cycle has been adopted in many countries, mainly in light water reactors (LWR) in the form of mixed oxide (MOX) fuel. MOX fuel is being used around the world, for instance, 21 reactors are currently using MOX in France, which generate 10% of the total French electricity. Eleven Japanese generating companies were loading MOX fuel in their nuclear generating plants. In the USA, ORANO is working to build the Mixed Fuel Fabrication Facility which will manufacture MOX fuel using ex-military plutonium [START_REF] Areva | Fabrication of recycled mox fuel[END_REF]. Figure 1.2, shows how this cycle can be represented. We now compare both; the open fuel cycle and the closed fuel cycle, it is argued that the former has more advantages related to proliferation, since no material for bombs is separated, and actual costs are favourable for this type of cycle. The latter has advantages concerning resource utilization, due to the fact that recycled actinides reducing the need for enriched uranium. Moreover, another important advantage is that the closed fuel cycle can reduce the requirements of long-term disposal. For example, if actinides are burned in a breeder reactor or even some long-lived Ąssion products, this will impact directly on the half-life and radiotoxicity of the waste that would be disposed [START_REF] Mit | The future of nuclear power[END_REF].

In 2003, the fuel cycle characteristics of the worldwide deployment of nuclear power have been summarized in Figure 1.3. The total capacity of the plants using the once-through enrichment oxide was 325 GWe of electricity. Concerning the plants using MOX in reactors, the total capacity was only about 27 GWe. As we can see, the majority of the plants opted for the once-through cycle. In the case of the closed fuel cycle, the plants using this method opted just for one recycle of the fuel before the spent nuclear fuel is put into disposal [START_REF] Mit | The future of nuclear power[END_REF]. If we take the case where in the year 2050 the deployment under the global growth scenario is achieved by two options. The Ąrst one consists of using the once-through cycle in LWRs or the second one, which consists in using plutonium recycle, where all the spent UO 2 is reprocessed but none of the spent MOX. If new resources of uranium are found or the ones that exist are suicient to maintain reasonable prices, then the closed fuel option could not be economically attractive compared to the open fuel cycle [START_REF] Mit | The future of nuclear power[END_REF]. However, the plutonium recycle option has an advantage concerning Ąnal waste disposal, because it produces less intermediate and high level wastes, although it also produces greater transuranic waste. Contrasting both cycles, both produce almost the same amount of Ąssion products which represent most of the activity in the nuclear spent fuel. In conclusion, by mid-century, if the prices of uranium ore stay reasonable low, the option of recycling plutonium would not be economically attractive. However, if the prices were to become high, the second option could become the most accepted one [START_REF] Mit | The future of nuclear power[END_REF].

MOX inside the reactor

The nuclear fuel pellets for Pressurized Water Reactors (PWRs) are small cylinders of around 5 mm of diameter [START_REF] Nea | State-of-the-art report on innovative fuels for advanced nuclear systems[END_REF]. For MOX, raw powders of UO 2 and PuO 2 are weighted and then mixed to meet the Ąnal fuel composition of 7.38% (for this case) of Pu with a Oxygen to Cation O/C ratio of 1.998% and an average grain size of 10 µm. The powder is then milled and granulated through pre-compaction, crushed and compacted again. Finally, the pellets are sintered. More details of this process can be found in [START_REF] Kaeri | Ceramography analysis of mox fuel after an irradiation test[END_REF]. The plutonium distribution of these pellets are shown in Figure 1.4. As we can observe, the plutonium distribution is quite homogeneous. The maximum diameter of the Pu rich spot is 13 µm with a Pu content of 18%. Under irradiation this inhomogeneity plays an important role. The MOX fuel is submitted to many micro-structural changes while being inside the reactor. One example is the restructuring of the solid matrix. It is due to the large temperatures and temperature gradients. Voids migrate towards higher temperatures and gather in the centre. Conversely, this causes a movement of the solid towards the outer periphery and leads to the creation of a void in the centre. An illustration of the Ąnal pellet state is shown in Figure 1.6, which represents a micrography of a MOX pellet. This phenomena makes the peak in the temperature decrease, since the fuel is now closer to the heat sink. Moreover, it afects the heat transport characteristics of the pellet. In MOX, this results in the migration of oxygen atoms through the pellet and conditions the demixing of U and Pu [START_REF] Welland | Matter transport in fast reactor fuels. Institute for Transuranium Elements[END_REF]. Figure 1.6 Ű Microstructure of a fast reactor pelletized MOX fuel pellet [START_REF] Welland | Matter transport in fast reactor fuels. Institute for Transuranium Elements[END_REF]. MOX fuel pellets are often divided into regions deĄned by their restructured state. Going outward from the centreline we Ąnd 1) Central void, 2) Columnar grain region, 3) Equiaxed grain growth and 4) As-fabricated microstructure.

There is a diference between MOX and UO 2 pellets. In MOX, there is higher rates of Ąssion gas release, due to higher linear power rating and higher centreline temperature (caused by slightly lower thermal conductivity). Furthermore, the heterogeneity of the Ąssile material would play a role as well, to increase the rate of Ąssion gas release. This is shown in Figure 1.7 for MOX fuel fabricated by the MIMAS process. The Pu rich spots are submitted to very high burn-up. This causes restructuring, which leads to further division of grains, precipitation of gas bubbles that go to intra-granular division which make them easily freed when a power transient occurs (CEA, 2009).

In overall, MOX fuel pellets are more complicated to work with, than the UO 2 . Their behaviour inside the reactor has to be carefully addressed, since more processes are involved at the time of the burn-up. Moreover, at the time of producing them, the alpha decays make it diicult to handle and decreases the number of labs capable of sensitize them. Therefore, it is , 2009). The colour level corresponds to wt % of PuO 2 important to understand well what are those micro-structural changes that happen inside the fuel under irradiation and surely the cause as well. This brings the need of new tools to assess the problem. In particular, modelling and computer simulations are expected to contribute to a better understanding of this phenomena.

Physical properties of MOX

The new generation of reactors are expected to work with higher Pu contents. The inĆuence of higher Pu contents are rarely investigated by the private nuclear sector. Furthermore, the efect of stoichiometry is expected to be very important to explain the processes inside the fuel (see Appendix A).

Experiments on MOX are diicult to carry out. Nonetheless, efort has to be made to get more information about this material since it is expected to be used by the new generation of reactors. First of all, a compilation of the knowledge found in the literature about MOX has to be done. This job was accomplished by a group of various European research organization in 1990. They built up a catalogue of properties of MOX, part of which is still used in models of calculation of oxide fuel codes for fast neutron reactors. A new efort was carried out by [START_REF] Esnii+ | State of the art with a literature review of mox properties[END_REF] to create a new catalogue of thermomechanical properties. It aimed to gather MOX fuel data required by the fuel calculation codes as GERMINAL, TRANSURANUS, MACROS, and TRAFIC. SpeciĄcally, input parameters as Temperature, Pu content, O/M ratio, porosity, burnup,. . . relevant for MOX fuels data needed for the reactor prototypes ALFRED (LFR), ALLEGRO (GFR), ASTRID (SFR), and MYRRHA (ADS). Table 1.1 shows the parameters and conditions that this catalogue comprises. Table 1.1 Ű MOX properties catalogue [START_REF] Esnii+ | State of the art with a literature review of mox properties[END_REF].

Although, this catalogue is quite complete and extensive, few gaps have to be completed. This is, for example, expand the Pu content domain and temperatures. Before the ESNII+ attempt, there have been some important works that join simulations and experiments to bring information about MOX. An example is (Devanathan, Van Brutzel, Chartier, & Gueneau, 2010, 10) which gathers information of actinide oxides mainly about crystal structure, thermal expansion and thermodynamic data.

Regarding the mechanical properties, It is well known that UO 2 and MOX mechanical properties are afected by the plutonium (Pu) content, the YoungŠs modulus increases slightly (3 to 10%) with addition of PuO 2 as demonstrated by experiments carried out in the seventies on MOX fuel up to 20% PuO 2 [START_REF] Novion | Mechanical properties of uranium plutonium based ceramics[END_REF], [START_REF] Nutt | Elastic and anelastic response of polycrystalline uo2-puo2[END_REF], [START_REF] Matthews | Mechanical properties and diffusion data for carbide and oxide fuels, ceramics data manual contribution[END_REF]. A recent study performed by Hirooka et al. [START_REF] Hirooka | Sound speeds in and mechanical properties of u,puo2x[END_REF] on MOX fuel containing the full range of Pu content (i.e., from 0 to 100%) conĄrms a 10% increase of the YoungŠs modulus.

New eforts to keep track on this material are being proposed. This is the case, for example, of the INSPYRE project. It stands for Investigation Supporting MOX Fuel Licensing in ESNII Prototype Reactors. The objective is to demonstrate Gen IV fast reactor technologies with closed fuel cycle, to harness European research and industrial capabilities to build advanced fast reactor prototypes and develop supporting infrastructure. Most of this PhD work lies under the task 3 (INSPYRE, 2018) It focuses on improving, developing and applying experimental and modelling methods to reach a better understanding of the underlying mechanisms governing the evolution of mechanical properties. One of the main questions that will be assessed is; "Can we use UO 2 thermomechanical data as MOX data in fuel performance codes?". This PhD research focuses on answering this question by using a computer simulation approach at the atomic scale.

Finally, extensive work has been carried out related to computer simulations (MD and MC). There are a few institution around the world that have undertaken this task. For example, CEA, Imperial College, Los Alamos National Lab, Ural Federal University and Japan atomic energy agency. Their works will be addressed in Chapter 3 and Chapter 4 to compare our results.

Crystallographic structure of MOX

The stoichiometric MOX has a Ćuorite structure. To be more speciĄc, the structure is cubic Fm 3m space group 225. The atomic representation of the structure is shown in Figure 1.8. The ideal crystal adopts a face centred cubic lattice for the cations and inside the anions occupy the tetrahedral sites, which form a simple cubic sub-lattice. Each cation, U or Pu, is then inside a cube where the corners are occupied by O atoms (Figure 1.9) Z U,Pu = 8. In the same manner, the O atoms are in the centre of tetrahedrons with U or Pu in the corners Z O = 4. This structure possesses several symmetries. Thus, it can be described by a unique parameter called "lattice parameter". The principal crystallographic orientation are in Miller index: the ⟨1 0 0⟩, ⟨1 1 0⟩ and ⟨1 1 1⟩ . The MOX Ćuorite structure has large octahedral interstitial holes, in which interstitial ions can be introduced.

Defects in MOX 1.6.1 Point defects

Under irradiation, the crystalline structure may present disorder when atoms are moved from their original crystal position. Defects can be divided in two. The Ąrst type are called intrinsic defects because those are only consisting in atoms from the same material. The extrinsic defects relates to atom type from another material such as Ąssion products. A clear example of this is the doping process. The totality of the work presented here is related to intrinsic defects. The Ąrst one is the Frenkel disorder, in which an atom that is displaced from its initial lattice site onto an interstitial one, leading this way to the creation of a couple vacancy-interstitial. This is shown in Figure 1.10. Both anions and cations can produce this type of defect.

The next type is called Schottky defect and is related to the creation of vacancies but these have to follow stoichiometry. For MOX, it is two Oxygen vacancies and one cation vacancy (Pu or U). Each vacancy species should occupy its own sublattice. A schematic diagram is shown in Figure 1.11. This leads to the third process. The "substitutional" defect is formed when a cation switches site with an anion. The later is not expected in MOX since it involves opposite charge ions and its distorting efect would be higher that just the creation of interstitials and vacancies.

Extended defects

In certain conditions, materials such as MOX, present more complex defects. These defects can be caused, for example, by irradiation after the gathering of point defects Section 1.6. This type of defects can be caused for example, by, irradiation. In general, we can classify defects by their dimension. 0-dimensional defects are point defects (again shown in Section 1.6). 1-dimension defects include perfect dislocation, partial dislocation, dislocation loops, grain boundary. 2-dimension defects are related to stacking faults and grain boundaries. Finally, 3-dimension defects are, for example; precipitates, voids or cavities.

Along this section, it will be brieĆy explained this type of defects. Firstly, Figure 1.12 shows a schematic diagram of point defects and extended defects in a basic lattice. [START_REF] Foll | Defects in crystals[END_REF].

We Ąrst consider dislocations. These defects are characterized by a vector, the Burgers vector, b. It represents the magnitude and direction of the lattice distortion. The movement of dislocations is the main source of plastic deformation in crystals [START_REF] Foll | Defects in crystals[END_REF]. Dislocations with Burgers vectors deĄning a translation of the lattice are called "perfect dislocations". Not all the Burgers vectors deĄne a translation in lattice, as we will see later in this section. In fcc lattices, such as the one belonging to the cation sub-lattice in MOX, the perfect dislocation Burgers vector is normally the shortest translation vector of the lattice given by b = a /2⟨1 1 0⟩. Perfect dislocations leave no change in the crystal structure when the dislocations move.

The fcc lattice can be represented as 3 staking planes formed by cations in a speciĄc way. The Ąrst atomic plane will be noted A, the second B and third C. The atomic plane B that will be stacked after the plane A has to have its atoms right over the depressions of the A (taking each atom as circle that in conjunction with the other will Ąll the whole space). The atomic plane C over B will have to follow the same process as for A to B. This yields an ABCABCABC stacking sequence. Thus, the atomic alignment repeats every third plane. There are cases when this sequence is not followed. For example, A and C planes are neighbours. These type of extended defects are called staking faults. Staking faults are caused by either interstitial or vacancies agglomerating. A schematic representation is shown in Figure 1.13. Intrinsic staking fault is caused by an agglomeration of vacancies. In the same way, extrinsic staking faults are due to interstitials. The Burgers vector of these dislocations is of the type b = a /3⟨1 1 1⟩. Dislocations with this type of Burgers vector are called partial Figure 1.13 Ű Schematic representation of stalking faults [START_REF] Foll | Defects in crystals[END_REF]. dislocations or Frank dislocations. Thus, the agglomeration of point defects in fcc materials may create "Frank loops". Now, if we consider dislocations which do not relate on removing or adding atoms. We have the Shockley partials dislocations. These dislocations are reproduced by cutting on a (1 1 1) plane (in this case between A and B planes). Then, move B to C position. Thus, a staking sequence of the type ABCACABCA will be obtained. The Burgers vector of these type of dislocation is b = a /6⟨1 1 2⟩. Figure 1.14 shows a schematic diagram of these dislocations. [START_REF] Radwan | Dislocations in fcc metals[END_REF].

These three types of dislocations play an important role with respect to the plasticity of the fcc material. They can interact with each other, giving for example, an "unfaulting" process which will be explained next. A Frank partial dislocation formed by a vacancy disc with a Burgers vector of the type b = a /3⟨1 1 1⟩, may split into a perfect dislocation b = a /2⟨1 1 0⟩ and a Shockley partial b = a /6⟨1 1 2⟩. The role of the Shockley partial is to move across the loop in order to remove the staking fault. Thus, a loop bounded by a perfect dislocation is created and is free to move. A simple diagram of these process is shown in Figure 1.15. Figure 1.15 Ű Schematic representation of a unfaulting process involving a Frank, perfect and Shockley dislocations [START_REF] Foll | Defects in crystals[END_REF].

Finally, a cavity is made by joining an important number of vacancies. This is normally found in the centre of displacement cascades as it is going to be shown along this work.

Along this chapter we will describe the background of the principal computational tools that is carried out in this work. Firstly, we will describe the basis of molecular dynamics. Then, we will go through the available interatomic potentials available for MOX. Consequently, a description of simulation techniques is given. Finally, we describe the tools for analysing the atomic structure used in this work.

Molecular Dynamics

Molecular Dynamics is a simulation technique for atomic scales. This technique is described in some reference texts [START_REF] Hoover | Molecular dynamics[END_REF], [START_REF] Goldstein | Classical Mechanics[END_REF]. The theoretical basis concerns the work of great minds of analytical mechanics such as; Euler, Hamilton, Lagrange and of course Newton. MD allows us to build the trajectory of the particles in the phase space of a system composed of a great number of atoms. The simplest form of MD, that is related to structureless particles, involves the NewtonŠs equations solution [START_REF] Rapaport | The art of molecular dynamics simulation[END_REF]. In MD, the atom inputs are the masses m a , position r a and velocities v a see Section 2.1.9.

Simulation box

Particles are in interaction inside a simulation box. The simulation box is usually a parallelepiped form. This box can be described with three vectors x, y and z. The output of a MD simulation are the individual trajectories of each atom along the time described by r a (t) and v a (t) (r a and v a for simplicity). These trajectories and the dimension of the simulation box are required to get thermomechanical properties, such as; temperature, pressure, internal energy, etc. In addition, the MD code needs: 1) the force between atoms, 2) equations that rule the atom trajectories, 3) an integrator that solves numerically these equations with discretization of the temporal domain.

Due to the intrinsic scale of MD, it is often necessary to expand the domain under study. This is done by creating inĄnitely exact replicas of the simulation box. Thus, each atom a in the original simulation box will have inĄnity number of replicas a ′ inside the simulation box replicas. This type of boundary conditions are called periodic. Atoms inside the original simulation box can interact with the simulation box replicas. However, self interaction is prohibited by setting the size of the box twice the distance of the interatomic potential cut-of. Another situation that is allowed to happen is that an atom can cross the original simulation box boundaries. In this case, the atom that have left will be found in one of the neighbour replicas and strictly one atom from the replicas will go inside the original but in the opposite side. Figure 2.1 show a simulation box with its replicas and an atom crossing the limits.

Atomic interactions

Along this work, the atoms are presented as points with a given charge. This representation is usually called rigid ion type. In addition, in this model the interactions between atoms can be taken in terms of potential energy pairs. Rigorously, the total energy of a system of N atoms interacting can be expanded in a many-body expansion:

U(r 1 , r 2 , ...., r N ) = ∑ a U a 1 + ∑ a ∑ b>a U ab 2 + ∑ a ∑ b>a ∑ c>b U abc 3 + .... (2.1)
Each term is important for the next discussion, thus, they are explained below;

1. U 1 is the one-body term. It is related to external force Ąelds or boundary conditions. For example, the wall of a container. 3. U 3 is the three-body term. Here the contribution of a third atom afecting the two body term is taken into account.

This can be easily generalized to any many-body interaction. Thus, for practical reasons we can categorize these interactions in: pair potentials (U 2 ) and many-body potentials (U 3 and higher terms). It is important to highlight that there are examples of many-body potentials in which multi-body efects are included implicitly and can be described through an environment dependence of two-body terms [START_REF] Zhigilei | Introdution to interatomic potentials[END_REF].

In the case of MOX, the potential energy can be divided in two parts: the electrostatic coming from the Coulomb law E ab c and the other part related to intrinsic inter-atomic interactions (excluding of course the electrostatic) noted as φ ab

E a p = E a c + φ a = 1 2 ∑ b (E ab c + φ ab ) (2.2)
We can deĄne the internal energy, such as, the sum of internal energies of the atoms that compose our system under study

U = ∑ a U a = ∑ a (E a c + φ a + E a k ) (2.3)
where E a k is the kinetic energy of atom a. Each component of Equation ( 2.3) will be explained in the following sections. SpeciĄcally, E a c in Section 2.1.6, φ a in Section 2. 1.3 and 2.1.4 and E a k in Section 2.1.9.

Pair potentials

In this subsection we will discuss about the non-coulombic term φ ab . Generally, this term can be divided in two contributions; one corresponding to the attractive van der Waals interaction and the other is of the repulsive kind related to the electronic clouds which forbid two atoms occupy the same volume. The mathematical form of the pair potential is arbitrary. Normally, it is opted to choose a mathematical form that is in concordance with the available physics theory.

There are 3 diferent approaches to obtain a potential function parameters for a particular system:

1. The Ąrst one is to assume a functional form for the potential function. Consecutively, parameters are chosen, in order to reproduce a set of experimental data. These potential are called "empirical". Examples are Lennard-Jones, Morse, etc 2. The second is related to the calculation of the electronic wave function for Ąxed positions. Analytic "semi-empirical" potentials are derived from quantum-mechanical theory.

3. The last is to calculate directly electronic-structure forces using ab initio MD simulations.

The main pair potential forms are the Lennard-JonesŠ potential, BuckinghamŠs potential and MorseŠs potential. The Lennard-JonesŠ potential is directly related to the van der Waals interaction in inert gases and molecular systems. It is often used to model general efects rather than speciĄc properties

φ ab = 4ϵ αβ [ ( σ αβ r ab ) 12 -( σ αβ r ab ) 6 ) (2.4)
where αβ indicates the chemical elements (Pu, U and O for MOX). The Buckingham potential provides a description of the strong repulsion due to the overlapping of closed shell electron clouds.

φ ab = A αβ e -r ab ρ αβ - C αβ (r ab ) 6
(2.5)

The Morse potential (Morse, 1929, 1) is similar to Lennard-Jones but is more suitable for cases when attractive interaction comes from the formation of chemical bonds.

φ ab = D αβ [ e 2B αβ (r ab -r αβ ) -2e -B αβ (r ab -r αβ ) ] (2.6)
It is important to highlight that Equation (2.4), Equation (2.5) and Equation (2.6) go to 0 if r ab increases. This brings the opportunity to establish a cut-of r c (distance from which all the interactions are set to zero) which is important for the MD codes that uses parallelisations. Finally, the force coming from φ ab between two atoms a and b is calculated by

f ab = - 1 r ab
∂φ ab ∂r ab r ab (2.7)

EAM porential

This type of potential includes many-body interactions, which in some case is important for potential energy calculations. EAM stands for embedded atom method (Daw & Baskes, 1984, 12) [START_REF] Finnis | A simple empirical n-body potential for transition metals[END_REF]. These EAM potentials were mostly aimed for metals.

The chemical bonds between pairs of atoms in metals depend signiĄcantly on the nearby atoms caused by quantum efects that describe the inĆuence of the electron gas. The distance between two atoms is no longer suicient to calculate the potential energy. For these potentials the surrounding atomic positions are needed. In general, the EAM potentials have one pair potential contribution that simulates the basic repulsion and attraction of atoms, plus a multi-body term that takes care of the surrounding density of charges.

V E AM (a) = a ∑ b=1 φ(r ab ) + f (ρ a ) (2.8)
here φ is the pair potential, ρ a is the local electron density and f is the embedding function. ρ a depends on the local environment around atom a, and the embedding function f represents how the energy of an atom depends on the local electron density. The electron density is calculated on a pair potential that maps distances between atoms to the corresponding contribution to the local electron density.

The EAM potentials allow a better reproduction of the mechanical properties of metals than pair potentials (Daw & Baskes, 1984, 12). Moreover, it has been applied to actinide oxides [START_REF] Cooper | A manybody potential approach to modelling the thermomechanical properties of actinide oxides[END_REF], [START_REF] Cooper | Modelling the thermal conductivity of uxth1xo2 and uxpu1xo2[END_REF] with Ćuorite structure in order to reproduce the CauchyŠs violation (Axe, 1965, 4A), which cannot be reproduced by only with pair-wise potentials.

ZBL porential

The ZBL potential stands for Zeigler-Biersack-Littmark [START_REF] Ziegler | The stopping and range of ions in matter[END_REF]. It is an universal potential design to reproduce the close repulsion between atoms. [START_REF] Nordlund | Repulsive interatomic potentials calculated using hartree fock and density functional theory methods[END_REF] gives the following description : "It has been constructed by Ątting an analytical function to a large number of repulsive potentials evaluated with the Gombas approximate method for diferent pairs". This potential is important in simulations of energetic ion bombardments. The form of this potential is given as

E ab Z BL = 1 4πϵ 0 Z a Z b e 2 r ab φ( r ab ζ ) ζ = 0.46850 (Z a ) 0.23 + (Z b ) 0.23 φ(x) = 0.18175e -3.19980x + 0.50986e 0.94229x + 0.28022e -0.40290x + 0.02817e -0.20162x (2.9)
where e is the electron charge, ϵ 0 is the electrical permittivity of vacuum, and Z are the nuclear charges. During the displacement cascades, ballistic shocks cause the interatomic distances to be much less than their equilibrium distances. For small interatomic separations, the ZBL potential is used for the distances less than 1.6 Å. The continuity of the cross over between the other potentials and the ZBL potentials is achieved via a Ąfth degree polynomial function chosen such that the potential energy, the forces, and the Ąrst derivative of the forces are continuous at the transition.

Electrostatic force

The electrostatic potential energy for a simulation box with a point charge distribution q b and their corresponding distances to r b with periodic boundary conditions is given by:

E c = ∑ a ∑ b a ∑ l 0 q a q b |r ab -l| (2.10)
here, l is a linear combination of the simulation box given by l = n x x + n y y + n z z and x, y and z the components of the vector between the two atoms under consideration.The easiest way to calculate Equation (2.10) is to propose a cut-of that accomplishes |r ab -l| < r c so any term that does not accomplish this condition in Equation (2.10) will not be taken into account. However, the sum in Equation (2.10) does not converge. Moreover, r c would imply a unreasonable cut-of for the most of the simulation boxes. For this reason, diferent ways of dealing with the electrostatic force have been investigated. Along this work we have used only two. The Ąrst one is a method that is derived from the Ewald sum [START_REF] Ewald | Die berechnung optischer und elektrostatischer gitterpotentiale[END_REF]. We used this method for the thermophysical and mechanical properties calculation. The second is the particle-particle particle-mesh solver (pppm) [START_REF] Hockney | Computer simulation using particles[END_REF]. We used this method for the crack propagation and the irradiation studies. The later, maps atom charge to a 3d mesh. It uses 3d fast Fourier transformation to solve PoissonŠs equation on the mesh, then interpolates electric Ąelds on the mesh points back to the atoms. How to decide between both is merely related to the computational time. The Ąrst scales faster with the number of atoms than the second. However, it is worth to take into account that Ewald is more accurate than pppm. Thus, it depends on the system under study for choosing which one of those methods is going to be used. In this work, we used a r c of 1.2 nm and the accuracy of the pppm was 10 -5 .

MOX potentials

For the mixed oxide compound (U,Pu)O 2 , several interatomic potentials are available in the literature. It exists two main families of potentials. One that considers U and Pu cations as one single entity A, hence they include only three set of parameters (C-C, C-O, and O-O) but depends on the relative percentage of Pu in the MOX. The second one treats explicitly the U and the Pu cations. Therefore, they include six set of parameters (U-U, U-Pu, Pu-Pu, U-O, Pu-O, and O-O) but do not depend on the percentage of Pu. Because we are interested in studying the spatial repartition of both cation sublattices, we will only consider and describe the second type of force Ąeld.

We have found Ąve interatomic potentials that we will coin later on by the name of their Ąrst 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 uo2 and puo2 by equilibrium molecular dynamics simulations from 300 to 2000k[END_REF], Potashnikov [START_REF] Potashnikov | High-precision molecular dynamics simulation of uo2puo2: pair potentials comparison in uo2[END_REF], Tiwary (Tiwary, Walle, & Jeon, 2011, 9), and Cooper [START_REF] Cooper | A manybody potential approach to modelling the thermomechanical properties of actinide oxides[END_REF], [START_REF] Cooper | Modelling the thermal conductivity of uxth1xo2 and uxpu1xo2[END_REF]. These Ąve force Ąelds can be separated according the properties on which they have been Ątted. All potentials have been Ątted 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 Ąrst one followed by Arima and Potashnikov with some improvement at high temperature, up to the melting point. Tiwary potential includes also Ąt on the formation energy of point defects (Frenkel pairs), while Cooper potential focusses on experimental data for single crystal elastic constants.. However, with Tiwary potential it was impossible to run MD simulations, the Ćuorite structure of UO 2 or PuO 2 is not stable, nevertheless it gives good results using static calculations. This is probably due to the fact that the energy landscape of this potential is very rough including a lot of none physical minima. Therefore, we eliminate this potential from our study.

As we can see, the potentials were Ątted on diferent parameters. The results that each potential is able to reproduce depend strongly on what they were Ątted on. All potentials try to simulate the best, the reality. However, as a user, we have to arise some doubts of their capabilities. Generally, at the time of recalculating the properties of which they were Ątted on, all show good results. For instance, Cooper potential was Ątted using elastic constants. One would then, give more conĄdence that this potential will show better results at the time of calculating mechanical properties. Finally, when there is lack of experimental data to compare with, the conĄdence will rely rather on the overall results that they give. Obtaining diferent phenomena when simulating the same system, under the same conditions, for all the potentials, give us more statistics to consider. Thus, we can give an upper and lower limits where the real scenario is hoped to be between them. The way to explain the reality with these type of simulations is then highly dependent on the trends in the results.

The four remaining potentials are described using a general formula with classical shortrange (Buckingham and Morse) and long-range (van der Waals and Coulomb) interactions (see Eq. 2.11). Table 2.1 presents the forms of the four potentials.

U αβ (r) = A αβ e -r/ρ αβ - C αβ r 6 + D αβ [ e -2γ αβ (r-r 0 ) -2 e -γ αβ (r-r 0 ) ] + q α • q β 4πϵ 0 r (2.11)
where A αβ , C αβ , D αβ , γ αβ , q α and r 0 are parameters whose values are given in Table 2.2.

Recalling Section 2.1.3, the Ąrst term in Eq. 2.11 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.

For the Cooper potential, an EAM (Embedded Atom Model) many-body term is added (Eq. 2.12 and Eq. 2.13), in order to reproduce the CauchyŠs violation observed in actinide oxides with the Ćuorite structure (i.e. C 12 C 44 ) (Axe, 1965, 4A), which cannot be reproduced by only pair-wise potentials.

E EAM a = -G α √ ∑ b a σ β (r ab ) (2.12)
where, σ β is computed as follow: .13) erf stands for the error function and η β is a parameter. Tables 2.3 provides the EAM Cooper potential parameters. The functions of the pair energies as a function of the interatomic distances are displayed in Figure 2.2 for all the potentials. The comparison per type of interaction allows us to better visualize the diferences between each potential. In general, all the potential show similar energy proĄles related to the interaction C -O. Diferences rise at interaction C -C. A steeper energy proĄle is given by Cooper and Arima, in comparison with Potashnikov and Yamada.

σ β (r ab ) = 1 2 ( η β (r ab ) 8 ) { 1 + erf [ 20(r ab -1.5) ] } . ( 2

LAMMPS code

In the following studies we used the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code (LAMMPS, 2018) is a classical molecular dynamics code aimed for materials. It can model ensembles of particles in liquid, solid or gaseous state. Moreover, it can model atomic, polymeric, biological, metallic, granular, and coarse-grained systems using a variety of force Ąelds and boundary conditions. LAMMPS is eicient in the way that it can run on single-processor machine up to big clusters using parallelisation. There is no limitation on the number of atoms in the system under study (only computational power limitation). It is an open source code. Thus, it can be modiĄed for personal research objectives. It was originally developed under a US Department Essentially, LAMMPS solves the NewtonŠs equations of motion for atoms or molecules that interact using many diferent potentials available in its packages (or implemented by the user). It increases its eiciency by using neighbour lists to keep track of nearby particles. Using it in parallel machines, LAMMPS decomposes the simulation domain into small 3d sub-domains which are assigned to each available processor. Then, processors communicate and store "ghost" atom information for particles that border their sub-domain [START_REF] Lammps | [END_REF]. In order to achieve maximum eiciency, 3d simulation boxes with uniform density have to be used. Finally, it is friendly to new users, because , it runs from an command-like input.

Simulation techniques

Along this work, we have carried out structure optimization simulations which can be coined static calculations, but also, trajectory calculations at Ąnite temperatures named as dynamic calculations. The static calculations are based in an optimization of an atomic coordinate function (positions and/or velocities). The molecular dynamic techniques are based on the numeric solutions of motion equations. In general, these equations are derived from the analytical mechanic theory formulations. These are the Lagrangian or Hamiltonian formulations. A brief description of these techniques will be presented next. Firstly, the description of the coordinate optimization. Consecutively, a description that simulates a system under thermodynamic ensembles described by the statistical mechanics theory.

Static simulations

These type of simulations are used to obtain information about a system in a stable state. They do not allow us to calculate dynamic quantities. However, they are important in the sense that they complement the dynamic calculations. For instance, we can calculate states with lowest energies. In these kind of states, the atoms are not moving, thus, we can propose that they have temperatures equal to 0. On the other hand, there is an option that allows us to perform the same type of calculations but with system at a certain temperature. This is by using the harmonic approximation in order to calculate the entropy of a system from the phonon spectrum. In general, this type of simulation is much faster than MD.

The static simulations need three elements. The inputs are the atomic positions r a , the simulation box and the potential. The output is an equation of state that can be, for instance, the internal energy, the enthalpy or the free energy. This function will be minimized by changing the system input. There are several algorithms dedicated to this minimization, such as; the conjugate gradient, the Broyden-Fletcher-Goldfarb-Shanno algorithm, the Rational Function Optimization method, etc . Along this work, we have used only the Ąrst, the Polak-Ribiere conjugate gradient version [START_REF] Polak | Note sur la convergence de méthodes de directions conjuguées[END_REF]. The reason of our choice is that this method is one of the most popular ones for solving smooth unconstrained optimization problems due to its simplicity and low memory requirement [START_REF] Yuan | A modiĄed polakribierepolyak conjugate gradient algorithm for nonsmooth convex programs[END_REF].

Microcanonical ensemble

For MD simulations, the equation of motion are solved within thermodynamical ensemble. The micro-canonical ensemble keeps the Volume, the number of particles and the internal energy constant (NVE) during the simulation. We use the Lagrangian formulation.

L = 1 2 ∑ a m a ( ∂r a ∂t ) 2 -Φ (2.14)
where Φ is the potential energy. The motion equations with this formulations is given by the Euler-Lagrange equation ∂ ∂t

( ∂L ∂ q i ) = ∂L ∂q i (2.15)
q i are the generalized coordinates of the system and q i their derivatives with respect of time. From here we can get the second NewtonŠs equation putting these generalized coordinates as the atomic positions.

m ∂ 2 r a ∂t 2 = - ∂Φ ∂r a
(2.16)

The Hamiltonian that in this case corresponds to the internal energy is kept constant during the system evolution

H = 1 2 ∑ a m a ( ∂r a ∂t ) 2 + Φ (2.17)
In order to solve this equation numerically, the time has to be discretized in variable or Ąxed time steps δt LAMMPS uses the scheme called "Velocity Verlet" (Equation (2.18)) to solve the motion equations.The particularity is that the velocities appear explicitly and are calculated at time t. The following two equations are solved:

r a (t + δt) = r a (t) + δt ∂r a ∂t (t) + δt 2 2 ∂ 2 r a δt 2 (t) (2.18) ∂r a ∂t (t + δt) = ∂r a ∂t (t) + δt 2 [ ∂ 2 r a ∂t 2 (t) + ∂ 2 r a δt 2 (t + δt) ] (2.19)

Canonical ensemble (Temperature control)

Now, we will assess the case where the temperature has to be Ąxed to a certain value. This type of situation is important to simulate various systems in nature. For example, phenomena where the systems are surrounded by a thermal bath so the system keeps its temperature constant by exchanging energy with the environment. With these conditions, including maintaining the number of particles (N) and the size of the box constants (V), we can apply the canonical ensemble (NVT) to our system.

In general, two types of algorithms are used to maintain our system at a Ąxed temperature. The Ąrst one consists in a correction of the dynamic variables after the integration step (Berendsen). The second one uses a modiĄed Lagrangian (Nosé-Hoover).

2.1.9.3.1 Berendsen's thermostat Taking into account the equipartition theorem, the temperature is related to the kinetic energy of the particles by

T = 1 3N k b ∑ a m a ( ∂r a ∂t ) 2 (2.20)
Generally in MD, the temperature is deĄned as the quantiĄcation of the thermal agitation, which is not exactly the same as the macroscopic temperature deĄned as the average of Equation (2.20). The Berendsen thermostat [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF] consists then in artiĄcially Ąxing the kinetic energy in order to control the temperature. SpeciĄcally, it makes the temperature of the system T tend to a set temperature T * by changing the velocity of each atom by the factor:

λ 2 = 1 + δt τ ( T * T -1 ) (2.21)
the factor τ controls the eiciency between the system and the thermostat. In the case where it is big enough, we Ąnd the system evolving in a micro-canonical ensemble. Contrarily, the temperature Ćuctuations are heavily pushed to disappear. This thermostat is applied to systems that are far away from the desired temperature. Its practicality makes it a powerful tool. However, we have to notice that the particles trajectories do not strictly follow those of the Newtonian equations.

Nosé-Hoover thermostat

In order to evolve our system close to a real physical behaviour, another alternative from the Berendsen thermostat has to be taken. This is the Nóse-Hoover thermostat. It consists in modifying the Newtonian Lagrangian Equation (2.14) by adding a degree of freedom related to an external thermostat via a change of the time as:

dτ = ξdt (2.22)
This new time τ evolves in a way that the temperature can be controlled. The generalized coordinate ξ is associated with a factious mass Q. With this the Lagrangian is

L = ξ 2 2 ∑ a ( ∂r a ∂τ ) 2 -Φ + Q 2 ( ∂ξ ∂τ ) 2 -N f K B T * ln(ξ) (2.23)
with N f is the degrees of freedom. Following the same process as for the Newtonian Lagrangian, the Euler-Lagrange motion equation are

∂ 2 r a ∂t 2 = 1 m a f a - 1 ξ ∂ξ ∂t ∂r a ∂t (2.24) ∂ 2 ξ ∂t 2 = 1 ξ ( ∂ξ ∂t ) + ξ Q ∑ a m a ( ∂r a ∂t ) 2 - n f k B T * ξ Q (2.25)
These equations difer from Equation ( 2.15) by a factor 1/ξ∂ξ/∂t which is considered as a friction factor. Equation ( 2.24) and Equation ( 2.25) are perfectly reversible in time. The diference yields on that this system does not keep the internal energy constant, rather it conserves the Hamiltonian

H = 1 2 ∑ a m a ( ∂r a ∂t ) 2 + Q 2ξ 2
( ∂ξ ∂t

) 2 + Φ + n f k B T * ln(ξ) (2.26)
This type of system that obeys Equation ( 2.26) will Ćuctuate around a temperature T * . The parameter that can be adjusted to check the eiciency of this thermostat is the Ąctitious mass Q. One can apply the same reasoning as for the factor τ of the Berendsen thermostat. In contrast to the Berendsen thermostat, this method is applied to system close to equilibrium and describes a system that is in a real canonical ensemble. However, it converges slowly and more time steps are necessary to reach the desired temperature (typically 1 ps).

Pressure control (Parrinello-Rahman barostat)

Now, we will assess systems where the pressure has to be controlled. This is the case of systems that are in a isothermal-isobaric ensemble (NPT). For doing this, we will use the Parrinello-Rahman barostat (Parrinello & Rahman, 1980, 14) and the Martyna integrator [START_REF] Martyna | Explicit reversible integrators for extended systems dynamics[END_REF] . It is based of the Nose Hoover thermostat in the sense that it proposes new degrees of freedom on the components of the tensor H associated with generalized coordinates P g and the virtual mass W. The equation of motions including the barostat and thermostat are

∂r a ∂t = p a m a + P g W r a (2.27) p a ∂t = f a - P g W p a - ( 1 N f ) tr(P g ) W p a - ∂ξ ∂t p a (2.28) ∂H ∂t = P g H W (2.29) ∂P g ∂t = det(H)(σ ext -σ) + 1 N f ∑ a (p a ) 2 m a Id - p ξ Q P g (2.30) ∂ξ ∂t = p ξ Q (2.31) ∂p ξ ∂t = ∑ a (p a ) 2 m a + 1 W tr(P t g P g ) -(N f + 9)k B T * (2.32)
Here the adjustable parameters are the masses W and Q. It also involves numerous time steps to converge (typically 10 ps).

Atomic structure analysis

To analyse our system, we used mainly the software OVITO, which is a powerful tool to visualize the raw data coming from MD or MC simulations. This helps the scientist to understand the outcome of his results. It is named as Open Visualization Tool (OVITO). I was created by [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with ovito the open visualization tool[END_REF]. It is an open source program. Consequently, anyone can extend its abilities, resulting in a vast range of applications with many examples online.

From its many functions, OVITO has three major tools that were used along this work for post treating data. These are Wigner-Seitz cell analysis, Voronoi cell analysis and Dislocation Extraction Algorithm (DXA). A brief description is presented below.

Wigner-Seitz cell method

The Wigner-Seitz cell method is used to check if defects have been created in a lattice. It consists on comparing two diferent conĄgurations: the reference state and the displaced conĄguration. The reference state corresponds to the conĄguration including no defects. The displaced conĄguration, as its name proposes, is the conĄguration where changes have occurred and where defects might have been created. The only types of defects that this method can Ąnd are vacancies and interstitials.

The procedure of the Wigner-Seitz cell method is decomposed as follow. First, a reference grid is deĄned. To each point of this reference grid corresponds an atom of the reference state. Subsequently, we determine for each atom of the displaced conĄguration which atom of the reference grid is the closest. Each time that a particular reference grid atom is found as closest, its occupancy number is incremented by one. Hence, the occupancy is deĄned as the number of atoms of the displaced conĄguration sitting within the Wigner-Seitz cell centred on an atom of the reference grid.

Having the occupancy of each atom of the reference grid, we can deĄne a vacancy as the reference grid atom that possesses occupancy equal to 0. Conversely, an interstitial is deĄned as reference grid atom with occupancy superior or equal to 2. LetŠs note that if there is no change in the initial structure, all the occupancies will be equal to 1. Figure 2.3 displays a schematic diagram of how OVITO calculates the occupancy. This tool circumvents the explicit creation of the Voronoi cells and hence it is more computing eicient. Moreover, OVITO allows tracking the type of atom. This means that we can difer from Pu to U for the cation sites. Therefore, we can calculate the number of Pu atoms in substitution of U atoms and vice-versa.

Voronoi cell method

In contrast with the Wigner-Seitz cell method, the Voronoi cell method does not need a reference state conĄguration. It creates a spatial division of the displaced simulation box creating explicitly Voronoi cells. An illustration of this space tessellation is given in Figure 2.4. This method provides two new quantities: the volume of the Voronoi cells and the number of faces attached to each Voronoi cell. The number of faces is equal to the number of the Ąrst-nearest neighbours of each atom and deĄnes the Şcoordination numberŤ. This method has been chosen rather than the Wigner-Seitz cell analysis because numerous disoriented Ćuorite subdomains appear during the simulation. Consequently, the reference grid used for the Wigner-Seitz cell analysis does not apply anymore.

The Voronoi cell volumes have an important role in our post-processing work. OVITO allows the user to get a histogram of the Voronoi cell volumes. Instead, it will increase the Voronoi volume of each of its surrounding atoms. Thus, we take into account those atoms with large Voronoi volumes (peak in the right side in Figure 2.5) in order to calculate the number of vacancies. However, this number of count is divided by the number of Ąrst-nearest neighbours in order to determine the exact number of vacancies. For example, in the FCC cation sub-lattice the number of atoms corresponding to the volume of vacancy count is divided by 12. The cut-of for counting the number vacancies (red lines in Figure 2.5) is chosen such that it matches the number of interstitials created in the Ąrst step of the FPA method (Ąrst introduction of FPs).

Dislocation Extraction Algorithm

The Dislocation Extraction Algorithm (DXA) is implemented in OVITO to search dislocation in the defected conĄguration. It is originally created by the creator of OVITO and a major description can be found in [START_REF] Stukowski | Automated identiĄcation and indexing of dislocations in crystal interfaces[END_REF]. It constructs the Burgers circuits in the simulations box. The method consists in creating a sequence of atom-to-atom steps (∆x )in the dislocated crystal that will be named Burgers circuit C.

= - ∑ C ′ ∆x ′ (2.33)
The Burgers vector b is the closure failure of the path after transferring it to the conĄguration. Since it is impossible to construct all the possible Burgers circuits in the conĄguration under study. The DXA uses the Delaunay tessellation which their edges deĄne the set of elementary ∆x from which Burgers circuits will be constructed. The Delaunay tessellation for a set of points P in a plane is a triangulation DT(P) such that no point in P is inside the circumcircle of any triangle in DT(P). The code then maps each edge of the Delaunay tessellation to a corresponding vector in the perfect reference conĄguration with the help of the Common Neighbour Analysis (CNA) described else where [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with ovito the open visualization tool[END_REF]. The CNA searches for atoms that form a perfect crystal lattice. Therefore, the DXA can diferentiate between a perfect crystal lattice and a defected one. This leads to a division of the space in two: "good zone" and "bad zone". SpeciĄcally, the edges adjacent to atoms that do not form a perfect crystal are considered as "bad". The boundary between these two zones will be called interface mesh. The interface mesh encloses all the defects inside the crystal. The Burgers circuits are built on this surface which encloses all the dislocations. The algorithm counts all Burgers circuits in order of increasing length. This leads to follow the remaining of the dislocation line as shown in Figure 2.7. The analysis of the structure in Section 4.4 evolution is carried out on-the-Ćy. We determined the dislocation density with the DXA algorithm and examined the point defects with Voronoi cell analysis.

Simulation of X-ray powder Diffraction

This expeerimental technique is used to non-destructively analyse, detect and determine the composition of any crystalline phases of samples.

It is based on the principle that each crystalline solid has its unique characteristic X-ray powder pattern which may be used as a "Ąngerprint" for its identiĄcation. This technique uses the generation of X-rays form the bombarding a metal target (Cu, Mo usually) with a beam of electrons emitted from a hot Ąlament (often tungsten). The incident beam of ionized electrons from the K-shell (1s) of the target atom and X-rays are emitted as the resultant vacancies are Ąlled by electrons dropping down from the L (2P) or M (3p) levels (zur Loye, 2013).

Theoretically, X-ray difraction is based on the principle that in order to obtain constructive interference, the path diference between the two incident and the scattered waves, which is 2d sin(Θ), has to be a multiple of the wavelength λ. For this case, the Bragg equation gives the relation between interplanar distance d and difraction angle Θ. In general, from the XRD difractometer λ and Θ are known, thus the distance between any planes of atoms present in any crystalline phases can be detected. However, if the sample has diferent crystalline structures within, the analysis could become complex.

In this work, we use this technique in another way. We use the Debyer code (DEBYER, no date) to create X-ray powder difraction patterns. It uses as an input a Ąle with positions of all the atoms. Difraction patterns are calculated using the Debye scattering equation

I(Q) = ∑ i ∑ j f 1 f j sin(Qr i j) Qr i j (2.34)
where 1. Q is the scattering vector, called also momentum transfer vector Q = |Q| = 4π sin(Θ/λ) being λ the wavelength 2. f i is the atomic scattering factor of i-th atom; in general it depends on Q and we should write it as f (Q) , but we donŠt, to keep the notation simple.

It is convenient to recall some properties that can be extracted from XRD graphs. According to BraggŠs equation, peak position (difraction angle) is a function of the distance (d hkl ) between reĆection planes (hkl) when the wavelength λ is Ąxed. Therefore, the peak position correspond to the distance of reĆection plane. Shifting of difraction peaks, hence, can be related to swelling [START_REF] Zhou | The physical meanings of 5 basic parameters for an x ray difraction peak and their application[END_REF]. Another important property is the "Full width at half maximum intensity" (FWHM). This property is sensitive to the variation in microstructure and stressŰstrain accumulation in the material. The more lattice planes with identical orientation are contributing to a difraction peak, the sharper the difraction pattern will be. In order to clarify more this idea, Figure 2.8 shows the efect of strain in the peaks.

Can Pu atoms be randomly distributed?

Along this work, we have assumed that MOX behaves like an ideal solid solution. This means that we can locate both types of cations (Pu and U) in their respective cation places in the Ćuorite conĄguration without any special restriction. Thus, we have located the Pu atoms randomly distributed in our simulation box. This rises many questions, such as; how reliable is this assumption?, if this assumption does not match the reality, how much this will impact our results?. It is known that for high Pu contents and Oxygen-to-metal ratio diferent than 2, miscibility gaps between UO 2 and PuO 2 can be found [START_REF] Truphemus | Etude des equilibres de phases en fonction de la temperature dans le systeme uo2 puo2 pu2o3[END_REF].

MD simulations are restricted due to their intrinsic time scales. Cation difusion is limited. This means that the cation positions are somehow Ąxed along the simulations. Therefore, the impact that cation heterogeneity (chemical disorder), zones with lower or higher Pu contents may impact on the outcome is not known using these techniques. However, with the help of simulations techniques such as Monte Carlo methods, we can assess this problem. The advantage is that MC can allow exchange one cation (either Pu or U) in its respective cation position for another one and then carry out this process in order to search in the conĄguration space. This will allow us to really understand the thermodynamic equilibrium and study the lower energy conĄguration of MOX solid solution.

This study could not be carried out during this PhD research due to lack of time. However, a CEA colleague named E. Bourasseau has performed this study without publishing it yet. We will brieĆy explain his conclusions. He used Cooper and Potashnikov potentials for a 25% Pu content conĄguration at 1000 K. The simulations box consisting of 96 atoms then periodically duplicated up to 2592 atoms was obtained from ab initio methods. He divided his results in two parts: the Ąrst consists in calculating thermodynamical properties, such as; Enthalpy, SpeciĄc heat, thermal expansion coeicient and molar volume. The second consists in carrying out a research on the cation distribution using the cation exchange technique. The results from the Ąrst part suggest that interchanging cations does not inĆict a signiĄcant change (diference of a few percentage unities) in the thermodynamic properties. The results from the second part show a diference caused by allowing the cation exchange. From the radial distribution function, Potashnikov potential shows a small decrease of the Pu-Pu coordination number with an efect that increases with temperature (e.g., 2.95 instead of 3 for U 7 5Pu 2 5O 2 ). For Cooper potential inverse conclusion is found, with a small increase of the Pu-Pu coordination number (e.gr., 3.1 instead of 3 at 300 K and 3.02 instead of 3 at 1000 K). As we can see, at higher temperatures this potential show less Pu "clustering" than at lower temperatures. Now that we aboarded the clustering. He found as well that there is no signiĄcant clustering of Pu atoms. Instead, they show an average distribution of Pu.

Having in mind the results of our colleague, we can say that the results presented herein with the assumption of a perfect MOX ideal solid solution with random distribution of Pu are somehow supported by the MC method. The impact on the thermodynamic properties is not signiĄcant. However, it does not imply that the Pu atoms will not afect other kind of results, such as; mechanical or the behaviour after irradiation. It is worth to highlight that the MC work only assessed one Pu content. Actually, our work found that there is an impact on the mechanical properties as a function of the Pu content. Furthermore, higher diferences are expected to appear when the stoichiometry is diferent from 2. To contribute to a better understanding of MOX, we have used a computational approach. 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 uo2 and puo2 by equilibrium molecular dynamics simulations from 300 to 2000k[END_REF], [START_REF] Arima | Evaluation of thermal conductivity of hypostoichiometric upuo2x solid solution by molecular dynamics simulation at temperatures up to 2000k[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 u0.75pu0.25o2x[END_REF], [START_REF] Cooper | Modelling the thermal conductivity of uxth1xo2 and uxpu1xo2[END_REF] and [START_REF] Li | Molecular dynamics study of thermal conductivities of u0.7xpu0.3amxo2[END_REF]. However, the reliability of the results depends exclusively on the choice of the set of potentials.

Thermomechanical assessment of empirical potentials

The potential parameters are usually Ątted 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 have not been Ątted on need to be assessed to provide a good insight of advantages and disadvantages of each potential and their range of validity. Thus, we carried out an assessment for MOX on the structural, thermodynamics, and mechanical properties.

The thermomechanical properties shown here were chosen to help the efort to create a new catalogue which is aimed to gather MOX fuel data required by the fuel calculation codes as GERMINAL, TRANSURANUS, MACROS, and TRAFIC [START_REF] Inspyre | Investigation supporting mox fuel licensing in esnii prototype reactors[END_REF]. The properties will be compared to experimental values and computer simulation accordingly.

Along this work, we will compare systematically via molecular dynamics simulations structural properties (lattice parameter, thermal expansion), thermodynamic properties (heat capacity, enthalpy), and mechanical properties (elastic constants, toughness) for four empirical potentials (see Section 2.1.7) over temperatures ranging from 300 K to the melting point (approximately 3000 K) and for compounds from pure UO 2 to pure PuO 2 . This assessment should help us select which interatomic potential to use for studying primary damage in MOX.

The whole study was carried out using the code LAMMPS. Most of the assessment requires relatively small system size: about 7 × 7 × 7 Ćuorite-type unit cells involving 4116 atoms. This size is large enough to avoid self-interaction, because the cut-of distance for the interatomic potentials is equal to 12 Å. Moreover, it ensures enough statistics for the diferent conĄgurations. However, for the crack propagation, the size of the box is increased to 439 × 128 × 7 Ćuorite-type unit cells. In all the cases we are using the periodic boundary conditions to avoid surface efects. The coulomb interactions are calculated with full Ewald summation procedure with exception of the crack propagation test, where we use pppm (see Section 2.1.6). The initial relaxation is done under NPT thermodynamic ensemble using the Parrinello-RahmanŠs algorithm, in order to set pressure and temperature as desired. For the thermodynamic property assessment (Section 3.1), the systems are relaxed 100 ps and all the thermodynamic properties are averaged over the last 20 ps.

This research was limited to study of stoichiometric compounds. Therefore, were can consider the (U 1-y Pu y )O 2 compound as a continuous solid solution. The plutonium atoms are then distributed randomly in the unit cell on the a-Wyckof sites (see related discussion in Section 2.4). The plutonium concentration ranges from 0 to 100% with mainly 15% steps.

Thermodynamical properties

Lattice parameter

The Ąrst structural property is the evolution of the lattice parameter with the temperature. All the interatomic potentials studied should Ąt more or less the experimental results since they were Ątted on this property. However, Yamada and Arima Ątted their potential only up to 2100 K, whereas Potashnikov and Cooper Ątted their potential with values up to 2900 K.

The changes in the calculated lattice parameters for pure UO 2 as a function of temperature for all potentials are shown in Figure 3.1 and compared with some experimental data [START_REF] Yamashita | Thermal expansions of npo2 and some other actinide dioxides[END_REF] and FinkŠs recommendations [START_REF] Fink | Thermophysical properties of uranium dioxide[END_REF]. It is worth to mention that experimental data are really sensitive to the O/U ratio (ESNII+, 2015).

Figure 3.1 Ű Evolution with temperature of the lattice parameter for the four potentials studied for pure UO 2 . The results are compared with FinkŠs recommendation [START_REF] Fink | Thermophysical properties of uranium dioxide[END_REF] and experimental work [START_REF] Yamashita | Thermal expansions of npo2 and some other actinide dioxides[END_REF] As expected, all the resulting lattice parameters calculated follow perfectly the experimental data and recommendations up to 2100 K. At this temperature, a clear bifurcation appears for Yamada potential, which largely underestimates the lattice parameters. Arima potential also underestimates slightly the lattice parameter at high temperatures, while Cooper and Potashnikov potentials follow the recommendation up to the melting point, although Cooper potential seems closer to the recommendation.

The same comparison is made for PuO 2 in Figure 3.2. However, experimental works are scarcely available at high temperatures, thus it is diicult to estimate their reliability. Nevertheless, Cooper, Potashnikov, and Arima potentials Ąt relatively well the experimental data up to 1300 K. For temperatures higher than 2100 K, Arima potential deviates slightly from Cooper and Potashnikov, underestimating consistently the 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 uo2puo2: pair potentials comparison in uo2[END_REF]. With Yamada potential, a phase transition from Ćuorite 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 Ćuorite contrary to experiments [START_REF] Gardner | The plutonium-oxygen phase diagram[END_REF]. Therefore, we eliminate this potential from the following study.

The evolution of the lattice parameters as a function of the temperature for each potential with plutonium content of 25 and 70% respectively are shown in Appendix B. The general behaviour is similar to pure compounds. Hence, we can Ąt 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 Figure 3.2 Ű Evolution with temperature of the lattice parameter for the four potentials studied for pure PuO 2 . The results are compared with [START_REF] Yamashita | Thermal expansions of npo2 and some other actinide dioxides[END_REF].

Table 3.1. The percentage of diference 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. Nonetheless, it is worth to point out that the discrepancy with Yamada potential increases with the increasing amount of plutonium content. Table 3.1 Ű Parameters for the third degree polynomial that Ąts the evolution of the lattice parameters as the function of temperature and plutonium content.

A [Å.K -3 ] B [Å.K -2 ] C [Å.K -1 ] D [Å] m [Å] Arima
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

In order to analyse closely the diferences among the potentials, we display the evolution of the lattice parameters as a function of the plutonium content at 300 K in Figure 3.3. Since the ionic radius of uranium is higher than for plutonium, the lattice parameter of UO 2 is higher than PuO 2 . Between both pure compounds the evolution of the stoichiometric solid solution is expected to follow the VegardŠs law [START_REF] Popov | Thermophysical properties of mox and uo2 fuels including the effects of irradiation[END_REF] and [START_REF] Esnii+ | State of the art with a literature review of mox properties[END_REF]. We recall that VegardŠs law states that the lattice parameters of a solid solution vary linearly with concentration at constant temperature. The lattice parameter is controlled by the relative size of the atoms or species exchanged. It can be expressed with the following equation:

a T = a T UO 2 (1 -y) + a T PuO 2 (y) (3.1)
Where, y, is the mole fraction of Pu and, a T , denotes the lattice parameters at a constant temperature.

The three remaining potentials follow well VegardŠs law at 300 K. However, Potashnikov potential slightly underestimates the lattice parameters for compounds enriched in uranium The VegardŠs law for stoichiometric MOX compounds has been veriĄed up to 2000 K [START_REF] Popov | Thermophysical properties of mox and uo2 fuels including the effects of irradiation[END_REF] and [START_REF] Esnii+ | State of the art with a literature review of mox properties[END_REF]. We investigated this behaviour with all three remaining potentials up to melting points. The percentage of deviation from VegardŠs law as function of temperature for compounds with plutonium content of 25 and 70% are presented in Figure 3.4 and Figure 3.5, respectively.

Both compounds behave similarly for all potentials. From 300 to 2000 K the VegardŠs law is veriĄed within ±0.02%. From 2000 K to melting point, we observe for Potashnikov and Cooper potentials a signiĄcant deviation 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).

Nevertheless, the behaviour for temperatures lower than 2000 K is coherent with both experimental value and recommendation.

Thermal expansion coefficient

A more sensitive quantity to evaluate the structural property with temperature is the linear thermal expansion coeicient (LTEC). This coeicient is calculated as the Ąrst derivative of the lattice parameter with respect of temperature with the following expression: Where, a 0 , is the lattice parameter at 300 K and, ( ∂a ∂T ) P , is the Ąrst derivative at constant pressure calculated by numerical diferentiation from the lattice parameter evolution shown in Section 3.1.1. Figure 3.6, displays the LTEC with Cooper potential for diferent plutonium contents. These values are compared to UchidaŠs experimental data for pure PuO 2 [START_REF] Uchida | Thermal expansion of puo2[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 diference with the plutonium content for this range of temperature. After 2000 K, peaks appear which are dependent of the plutonium content. For Cooper potential they depend on the plutonium content as previously mentioned [START_REF] Cooper | Modelling the thermal conductivity of uxth1xo2 and uxpu1xo2[END_REF]. These peaks, referred to as λ-peak, were studied by Bredig et al. for most of the Ćuorite-like structures [START_REF] Dworkin | Difuse transition and melting in Ćuorite and antiĆuorite type of compounds. heat content of potassium sulĄde from 298 to 1260 degree k[END_REF] and can be attributed to a difusing phase transition that turns compounds into a superionic conductor. Several reviews of experimental works [START_REF] Ronchi | Analysis of recent measurements of the heat capacity of uranium dioxide[END_REF] and (Ralph, 1987, 7) and theoretical works [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] Potashnikov | High-precision molecular dynamics simulation of uo2puo2: pair potentials comparison in uo2[END_REF] and [START_REF] Cooper | Modelling the thermal conductivity of uxth1xo2 and uxpu1xo2[END_REF] show that this transition occurs around 0.8 T m , where T m is the melting temperature. It is often associated with premelting of the oxygen sublattice.

α(a) = 1 a 0 ( ∂a ∂T ) P (3.2)
For the Cooper potential, as already stated by the authors [START_REF] Cooper | Modelling the thermal conductivity of uxth1xo2 and uxpu1xo2[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 conĄrm the existence of a superionic transition. However, according to the theory of 0.8 T m , it should occur around 2400 K since the melting temperature in PuO 2 was recently measured at 3040 K [START_REF] Bohler | Recent advances in the study of the uo2 puo2 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 | Modelling the thermal conductivity of uxth1xo2 and uxpu1xo2[END_REF]. Figure 3.8 displays the results of LTEC for 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 uo2puo2: pair potentials comparison in uo2[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.8 T m .

Overall, all potentials give a good estimation of the LTEC up to 2000 K. Diferences 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 diicult to choose the best potential, which reproduces the best the real material. However, the behaviour of the superionic transition temperature can provide more indication. The Cooper potential, 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 signiĄcant change with the increase of the plutonium content. Finally, Arima potential reproduces no λ-peaks in the range of the temperatures studied. More noticeable changes can be extracted with the analysis of the speciĄc heat. The speciĄc heat capacity coeicients at constant pressure can be calculated directly by numerical diferentiation from the enthalpy increment function with the following relationship:

Enthalpy and specific heat

C p = 1 n ( ∂H ∂T ) P (3.3)
where n is the number of moles. The evolution of C P as a function of temperature for diferent plutonium contents for the three potentials are shown in Figure 3.10, Figure 3.11 and Figure 3.12 . The temperature range investigated herein is above the Debye temperatures calculated with the diferent 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 uo2puo2: pair potentials comparison in uo2[END_REF]. On the other hand, Arima potential seems not to include this transition up to 3000 K.

Thermal conductivity

We did not carry out simulations of thermal conductivity in this work. The reason is that this property has been already highly assessed in the literature. This property is very important for the performance of the fuel inside the reactor. [START_REF] Esnii+ | State of the art with a literature review of mox properties[END_REF] has compiled experimental data about the thermal conductivity until 2015. Early recommendations (until the year 1990) indicate say that MOX can be taken as UO 2 in the range between 500 K < T < melting point, O/M ≤ 0, Pu content 0 < y < 0.30 and porosity 0 ≤ p ≤ 0.4. Later in-pile measurements of the temperature history at the pellet centre during irradiation with a thermocouple by the Halden Reactor Project [START_REF] Esnii+ | State of the art with a literature review of mox properties[END_REF] give a formula to obtain the thermal conductivity for MOX. This formula is obtained by multiplying the phonon term of the formula for UO 2 by a factor 0.92, corresponding to the ratio of the conductivities of fresh UO 2 and LWR MOX. [START_REF] Nichenko | Molecular dynamics study of the mixed oxide fuel thermal conductivity[END_REF] gives a critical overview of the experimental and modelling results until 2013. They say that it is generally accepted that the fresh LWR and FBR (fast breeder reactor) MOX thermal conductivity is lower than that for UO 2 . However, in their overview, it is shown that experimental results exist only within the plutonium concentration range from 0 to 30% and are contradictory. Moreover, the carried out MD simulations for temperatures between 300 to 1700 K using the Arima potential. They obtained a thermal conductivity not constant between 3 to 30% Pu content. The minimum value is found at 45% of the Pu content with a 14% decrease compared to the one of UO 2 .

Ab initio simulations [START_REF] Wang | Thermal conductivity of uo2 and puo2 from Ąrst-principles[END_REF] show that phonon group velocities for PuO 2 are larger that those of UO 2 along the three crystallographic directions. Thus, PuO 2 thermal conductivity is larger than that of UO 2 .

Finally, [START_REF] Cooper | Modelling the thermal conductivity of uxth1xo2 and uxpu1xo2[END_REF] MD work at temperatures between 300 to 2000 K show a small reduction in thermal conductivity in MOX compared to the end members (UO 2 and PuO 2 ). With the minimum at 25% Pu content and a higher thermal conductivity for PuO 2 than UO 2 . In conclusion, Arima and Cooper potentials show good results since the thermal conductivity decreases for MOX compared to UO 2 and PuO 2 .

Melting point

Concerning the melting points for the chosen empirical potentials, most of the potentialsŠ authors calculated this property. Therefore, we did not compute them. However, the method with which they were obtained varies from each other. [START_REF] Potashnikov | High-precision molecular dynamics simulation of uo2puo2: pair potentials comparison in uo2[END_REF] conducted a series of simulations at temperatures where melting may occur then the melting points were detected by sharp changes in density and enthalpy, obtaining this way temperatures around 4000 K. In contrast, Cooper determined the melting points using a moving interface method similar to that used by [START_REF] Govers | Comparison of interatomic potentials for uo2: part ii: molecular dynamics simulations[END_REF]. He reported a melting temperature for UO 2 around 3000 to 3100 K [START_REF] Cooper | A manybody potential approach to modelling the thermomechanical properties of actinide oxides[END_REF] and for PUO 2 around 2800 K [START_REF] Cooper | Modelling the thermal conductivity of uxth1xo2 and uxpu1xo2[END_REF]. The diference between the melting points from both authors arises from their method. Potashnikov states that the cause of too high melting temperatures is that crystals which are MD-simulated under PBC melt in a superheated state (spinodal condition) due to the lack of surface [START_REF] Potashnikov | High-precision molecular dynamics simulation of uo2puo2: pair potentials comparison in uo2[END_REF]. Cooper potential gives melting points closer to the experimental values.

Diffusion in MOX

We did not carry out any difusion calculation. However, we present some literature review about it here. The rates of difusion of diferent ionic constituents are diferent. The difusion mechanism of the constituents are based on the diferent ionic defects. The cationic defects are heavier in terms of their mass than the anionic defects. Therefore, anionic defects are more likely to difuse at lower temperature. Anionic difusion of MOX has been experimentally investigated mostly by the Thermogravimetry experiment. Bayoglu have studied the oxygen chemical difusion extensively with 20% of Plutonium content [START_REF] Bayoglu | Difusion chimique de lŠoxygene dans u, pu o2 x[END_REF], [START_REF] Bayoglu | Difusion chimique de lŠoxygene dans u,pu o2 x[END_REF], [START_REF] Bayoglu | Difusion de lŠoxygene sous gradient de potentiel chimique dans u ce o2 x[END_REF]. These results have shown that the difusion of oxygen varies with the variation of temperature and oxygen to metal ratio. Bayoglu has also proposed several relations between oxygen chemical difusion and temperature for systems with several oxygen to metal ratio. These relations are developed using the Arrhenius deĄnition framework of Difusion coeicient. These relations shows that the activation energy of oxygen chemical difusion is diferent with diferent ratios of oxygen to metal. The variation of activation energy between a hypo-stoichiometric system and a hyper-stoichiometric system is larger than the variation of activation energy between two hypo-stoichiometric systems or, between two hyper-stoichiometric systems. On the other hand, Kato has studied the oxygen chemical difusion of MOX with two diferent Plutonium content-20% Plutonium MOX and 30% Plutonium MOX [START_REF] Kato | Analysis of oxygen potential of u0.7pu0.3o2x and u0.8pu0.2o2x based on point defect chemistry[END_REF]. These results also shows that the oxygen chemical difusion varies with the variation of plutonium content in MOX. Recently, Vauchy et al. determined the oxygen self-difusion coeicient in MOX with a high enrichment of 45% plutonium for two diferent temperatures [START_REF] Vauchy | Oxygen selfdifusion in polycrystalline uraniumŰplutonium mixed oxide u0.55pu0.45o2[END_REF]. Moore and Kato, separately, developed two diferent models for the calculation of oxygen difusion in MOX [START_REF] Moore | Oxygen difusion model of the mixed (u,pu)o2 ± x: assessment and application[END_REF], [START_REF] Kato | Oxygen potentials, oxygen difusion coeicients and defect equilibria of nonstoichiometric u,pu o2 x[END_REF]. MooreŠs sublattice model deĄnes the difusion mechanism based on the migration of interstitial oxygen ion and vacancy and the migration of regular oxygen and vacancy. Whereas, KatoŠs model deĄnes the difusion mechanism as the migration of electrons and holes in the system.

Although, cations have much higher mass than the anions, the cationic difusion also takes place at the higher temperature in MOX. Riemer and Scherf evaluated the MOX with 15% enrichment at the temperature range of 1573 K to 1773 K [START_REF] Riemer | Plutonium difusion in hyperstoichiometric mixed uranium plutonium dioxides[END_REF]. These results are showing the temperature dependence of plutonium difusion in MOX. The rate of difusion increases with the increment of temperature. Matzke and Lambert have studied the self-difusion of plutonium in various MOX system with various enrichment (15%, 18% and 20%) [START_REF] Matzke | The evaporation behavior and metal self-difusion processes in u,pu c and u,pu o2[END_REF], [START_REF] Lambert | The diffusion of plutonium in uranium plutonium mixed oxide single crystals at varying oxygen to metal ratios[END_REF], [START_REF] Matzke | Difusion processes and surface efects in non-stoichiometric nuclear fuel oxides uo2+x, and u, pu o2 x[END_REF]. In these results, it is very evident that the rate of difusion of plutonium depends on the plutonium enrichment as well as oxygen-to-metal ratio.

Anion difusion has been calculated with MD for the three potentials by the authors of each potential (for Cooper potential [START_REF] Cooper | Modelling the thermal conductivity of uxth1xo2 and uxpu1xo2[END_REF]) with exception of Arima potential which was calculated by Potashnikov in [START_REF] Potashnikov | High-precision molecular dynamics simulation of uo2puo2: pair potentials comparison in uo2[END_REF]. For Cooper potential, it can be seen that the super-ionic transition occurs at a higher temperature in UO 2 compared to PuO 2 . The addition of uranium to PuO 2 increases the superionic transition temperature for all the MOX compositions studied. Potashnikov potential shows as well this transition at a temperature around 2600 K for UO 2 . Finally, Arima potential shows this transition at temperatures above 3000 K.

Phase stability (Energy-volume plots)

It is important for the potential assessment to check the relative phase stability of these phases to ensure that Ćuorite structure is the most stable phase for (U,Pu)O 2 compounds. For this reason, we have undertaken energy-volume studies shown in Figure 3.14 with Ąve diferent structures for both pure UO 2 and PuO 2 . The process consists on each step deforming the simulation box by 1% and carry out a coordinate minimization using the Conjugate gradient.

The structures considered in this calculations are: Ćuorite (Fm 3m), cotunnite (Pnma), rutile (P4 2 /mnm), scrutinyite (Pbcn), and marcasite (Pnnm). This study follows previous works on pure UO 2 [START_REF] Chernatynskiy | Critical assessment of uo2 classical potentials for thermal conductivity calculations[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 conĄrm that for both compounds the Ćuorite structure corresponds to the ground state at zero pressure for all three empirical potentials studied. As previously demonstrated with DFT calculations [START_REF] Desai | Stress induced phase transformation in nanocrystalline uo2[END_REF][START_REF] Fossati | Simulation of uranium dioxide polymorphs and their phase transitions[END_REF], we also Ąnd 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 diicult to distinguish which of the two phases is the most favourable. Using the common tangent method and given that P = -∂U ∂V , we calculate the transition pressures from Ćuorite structure to the other phases. These transition pressures are presented in Table 3.2 for pure UO 2 . For all the potentials and pure UO 2 , the lowest transition pressure under tensile loading is found for the transition from Ćuorite to rutile structure, which is coherent with results of the crack propagation.

From Table 3.

Mechanical properties

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 Ątted on those properties. Furthermore, the Cooper potentials has been Ątted 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 Ąnite 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 calculated as initial slopes of stress-strain curves obtained using appropriate components of stress and strain tensors. More speciĄcally, the YoungŠs modulus is deduced by applying uniaxial tension and compression strains individually and calculated as σ i /ε i , where σ i are the components of the stress tensor. Since Ćuorite structure is cubic, we obtain only three none zero independent elastic stifness coeicients: C 11 , C 12 , and C 44 .

Since we are assessing the mechanical behaviour in the (U 1-y Pu y )O 2 solid solution, it is important to check the inĆuence of the Pu content on the elastic properties. Figure 3.15 shows the evolution for the three potentials of the YoungŠs modulus and as a function of the Pu content. For the YoungŠs modulus, apart from the diference of scale, the three potentials show an increase of about 5% with an increase of the Pu content. Recent sound wave measurements carried out by [START_REF] Hirooka | Sound speeds in and mechanical properties of u,puo2x[END_REF] on monocrystal MOX show also an increase of the YoungŠs modulus (about 12%) when the plutonium content increases. Generally, in the litterature (U,Pu)O 2 elastic constants is taken identical to those of UO 2 . However, a moderate increase (about 3%) with the addition of plutonium content is indicated by some authors [START_REF] Novion | Mechanical properties of uranium plutonium based ceramics[END_REF][START_REF] Nutt | Elastic and anelastic response of polycrystalline uo2-puo2[END_REF].

The results that we have obtained for the evolution of the elastic properties as a function of the temperature with MD simulations can be compared to the recommendation based on the review of (ESNII+, 2015), 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 stifness tensor coeicients using the Voigt-Reuss-Hill approximation for randomly oriented polycrystals [START_REF] Hill | The lastic behavior of a crystalline aggregate[END_REF], [START_REF] Mehl | Intermetallic compounds principle and practice i principles[END_REF] and [START_REF] Schreiber | Elastic constants and their measurement[END_REF] . The equations to transform the elastic stifness coeicients found in monocrystal into shear 3.16(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 in agreement with the recommendations in ESNII [START_REF] Esnii+ | State of the art with a literature review of mox properties[END_REF].

Comparing the diferent 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 Ąts well the recommendations around 2400 K up to melting point. The Cooper potential Ąts 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. Therefore, it is expected that it reproduces well the Bulk modulus compared to the monocrystal experimental values of [START_REF] Hutchings | High temperature studies of uo2 and tho2 using neutron scattering techniques[END_REF]. In comparison, Potashnikov potential underestimates the values of the bulk modulus by approximatively 40%. c) is compared with experimental date of [START_REF] Hutchings | High temperature studies of uo2 and tho2 using neutron scattering techniques[END_REF].

However, it is interesting to mention that the absolute values of the YoungŠs modulus without the polycristal correction for Potashnikov potential are in the order of the YoungŠs modulus recommendation of (ESNII+, 2015) for polycrystalline MOX.

Anisotropy factor (Zener's)

Another important criterion is the anisotropy factor that reĆects the diference of bonding with the diferent direction of the crystal. In cubic system the anisotropy can be quantify by the ZenerŠs anisotropy factor deĄned as follow:

Z = (2C 44 ) (C 11 -C 12 ) . (3.5)
When Z is equal to 1, the system is perfectly isotropic, conversely a deviation from 1 shows that the system is anisotropic. Figure 3.17 and Figure 3.18 display 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 , respectively. Figure 3.17 Ű Evolution of the ZenerŠs anisotropy factor as a function of temperature for pure UO 2 for the three interatomic potentials.

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 the crystal is isotropic at lower temperatures and becomes slightly anisotropic around 2400 K where the BredigŠs transition takes place for this potential. The Z value calculated from the experimental elastic constants [START_REF] Hutchings | High temperature studies of uo2 and tho2 using neutron scattering techniques[END_REF] increases from 0.52 to 0.65 in the same range of temperature. It is worth noting that other experimental values calculated at 300 K are closer to values obtained with Cooper potential and range between 0.44 to 0.47 [START_REF] Wachtman | Elastic constants of single crystal uo2 at 25 c[END_REF]) (Berman, Tully, Belle, & Goldberg, no date). Moreover, Z < 1 indicates that the YoungŠs modulus values will be maximum along the ⟨1 0 0⟩ direction and minimum along the ⟨1 1 1⟩ direction. Furthermore, for the three potentials the anisotropic factor at 1600 K (temperature studied for the irradiated case) is still in the order of the one at 300 K. The behaviour of the anisotropy factor with respect of the Pu content at 300 K is shown in Figure 3.19 for the three potentials. For Potashnikov, Z = 1 through the entire range of Pu content, showing no change in the isotropic behaviour. Conversely, for Cooper potential the Z factor decreases quasi linearly as the Pu content increases, yielding to an even more pronounced anisotropy for PuO 2 (i.e., E 111 is 45% of E 100 for UO 2 and E 111 decreases to 35% of E 100 for PuO 2 ). This behaviour is also found with DFT+U calculations (Dorado & Garcia, 2013, 19) but with Z factor decreasing from 0.7 to 0.6 between UO 2 to U 75 Pu 25 O 2 . Arima potential shows no efect due to the Pu content. In all cases the elastic behaviour is quasi-linear with respect to the Pu content. Therefore one intermediate composition between pure urania and plutonia should be suicient to follow the following behaviour with radiation dose.

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 with 7 × 7 × 7 Ćuorite-type unit cells to evaluate the stresses as a function of strain for diferent plutonium contents and at diferent temperatures. The stress-strain 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.

We Ąrst check the stress-strain curves loaded with the three main crystallographic directions of the Ćuorite 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 Figure 3.20. 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 deĄned as the ultimate tensile strength (σ UTS ) at which the system begins to crack. For all the potentials the stifest 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] 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 Figure 3.21.

The ultimate tensile strength decreases almost linearly as the temperature increases for all potentials. Also, there is a slight increase due to the plutonium content. This is in disagreement with recommendations found in the European Commission state of art of MOX fuel report [START_REF] Esnii+ | State of the art with a literature review of mox properties[END_REF] in which it states 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 an increase of around 5% is recommended.

Brittle-to-ductile transition

Another important thermomechanical property is the brittle-to-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 uo2 minus 20 wt puo2 nuclear fuel[END_REF]). This transition is deĄned 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 deĄned 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.

We can deĄne 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 Figure 3.22. For all the potentials, we Ąnd 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 | Brittle fracture of oxide nuclear fuel[END_REF][START_REF] Roberts | Nature of brittle to ductile transition in uo2 minus 20 wt puo2 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.

Figure 3.22 Ű 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 [START_REF] Roberts | Nature of brittle to ductile transition in uo2 minus 20 wt puo2 nuclear fuel[END_REF]).

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.

Crack propagation

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 inĄnite strip containing a semi inĄnite crack[END_REF]. The geometry is illustrated in Figure 3.23. 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 Figure 3.23 Ű Representation of the thin strip geometry of the system used for cracking simulation.

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 <111> 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 Ągures of these snapshots are shown in Figure 3.24. Clear diferences appear between the diferent 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 Figure 3.24) 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 Ąeld necessary for crack propagation and then the volume afected with the secondary phase. But, this phase transition ahead of crack tip has already been observed with MD simulations for pure UO 2 with diferent interatomic potentials [START_REF] Fossati | Simulation of uranium dioxide polymorphs and their phase transitions[END_REF][START_REF] Zhang | Crack tip plasticity in single crystal uo2: atomistic simulations[END_REF]. The stress-strain curves calculated during the crack propagation are displayed in Figure 3.25. The strain at which crack starts to propagate is diferent for the three potentials. Crack initiates Ąrst with the Potashnikov potential at around 4% followed by Arima at 5% and Ąnally 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 3.2. For the Arima and Cooper potentials, σ UTS is greater than their relative transition pressure from Ćuorite 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 Ćuorite 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 GriithŠs theory and experimental results could be described by this behaviour [START_REF] Zhang | Crack tip plasticity in single crystal uo2: atomistic simulations[END_REF].

Figure 3.25 Ű 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 Ćuorite structure.

Conclusions

In this section, 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 the point defects spacial distribution under irradiation, both cations need to be explicitly modelled. Therefore, we found in the literature Ąve interatomic potentials fulĄlling these requirements coined by the name for their Ąrst 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 unstable with MD simulation. Namely, the Ćuorite 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 λ-peak appear, whereas no clear λ-peak 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 stifness constants are best reproduced with the Cooper potential, which has been Ątted 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 stifness 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 Ąnd 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, conĄrming 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 diferent 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 and thermodynamics. It also presents expected brittle behaviour during crack propagation but renders isotropic elastic behaviour. 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. As a consequence, we decide to use Cooper and Potashnikov potentials for our radiation assessment in Chapter 4. Safety analyses are required to ensure that MOX fuel can be burned safely in reactors. They require to take into account the thermomechanical properties of MOX as a function of irradiation, dose and Pu content. It is known that MOX thermomechanical properties vary with respect to irradiation [START_REF] Popov | Thermophysical properties of mox and uo2 fuels including the effects of irradiation[END_REF][START_REF] Esnii+ | State of the art with a literature review of mox properties[END_REF], yet, detailed description of this phenomena is rather sparse due to the great diiculties of carrying out experiments in such materials (high toxicity and radiation levels). In order to circumvent the problems, atomistic simulations can bring valuable information and increase the scope of knowledge in the understanding of radiation-induced defects. This section is decomposed as follow. Section 4.1 provides some point defect formation energies in order to assess the empirical potentials. Section 4.2 examines the results on the lifetimes of Frenkel pairs, which are of importance for the two next following sections. Section 4.3 presents the primary damage state obtained via displacement cascades. Section 4.4 discusses the efect of irradiation dose on the evolution of the microstructure using the Frenkel pair accumulation method. Section 4.5 assesses the evolution of mechanical properties such as Zener ratio and YoungŠs modulus as a function of dose. Finally, Section 4.6 presents simulated XRD patterns as a function of dose.

Damage evolution under irradiation.

Defect formation energies

Ballistic efects of irradiation result primarily in the production of intrinsic point defects, which may further anneal or cluster. It is then important to assess their formation energies. In our MD simulations the number of atoms remains constant. Therefore, we will analyse only Frenkel pairs and Schottky defects. We recall that in MOX, since its structure can be represented as the interconnection between a primitive-cubic anion sublattice embedded in a face-centred cubic cation lattice, two types of FPs can be found: cation and anion. They will be coined hereafter respectively C-FP and O-FP. The second type is the Schottky defect, which implies two oxygen vacancies and one cation vacancy. This defect has been found by positron annihilation spectroscopy to be the most predominant defect detected in the α irradiated UO 2 [START_REF] Wiktor | Coupled experimental and dft u investigation of positron lifetimes in uo2[END_REF].

Authors of the interatomic potentials used herein have already calculated part of these defect formation energies but using diferent approaches and mainly for urania and plutonia. Hence, in order to assess both potentials consistently for diferent Pu contents, we calculate for both potentials formation energies of C-FP, O-FP, and SD defects for 0, 25, 50, and 100% Pu contents.

The formation energy of each defect is calculated by minimization of the energy at 0 K. This method uses an optimization of the ion coordinates of the defected supercell using a conjugate gradient algorithm, more speciĄcally the version of Polak-Ribiere available within the LAMMPS code (referred in Section 2.1.9.1). The defect formation energy is simply calculated by computing the diference between the energy per atom of the defectless system (E bulk ) and the energy per atom of the system containing the defect E defect as follow:

E fp = E defect -E bulk (4.1) E sch = E N-3 defect - N -3 N E N bulk (4.2)
The FP is created by displacing one atom from its lattice position and placing it in a corresponding interstitial site in the supercell. The SD is created by simply removing one cation and two oxygen atoms. Thus, the resulting supercell will be a neutrally charged crystal with periodically copied defects. For each family of defect two types are initiated according to the separation distance between the vacancy and the interstitial and the three vacancies for the FP and the SD respectively. The Ąrst type corresponds to the bounded defect where the point defect (vacancy or interstitial) are connected, namely Ąrst nearest neighbour. Three diferent bounded Schottky defects are possible with the crystallographic orientation of the cubic structure. They are noted SD⟨1 0 0⟩, SD⟨1 1 0⟩ and SD⟨1 1 1⟩. Concerning the bounded FP, it is found that the Ąrst nearest neighbour recombines instantly during the relaxation. For this reason, bounded FP energies are calculated for the second nearest neighbour sites. Conversely, the second type corresponds to the distance between point defects at the "inĄnite". To calculate these values the size of the supercell is gradually increased from n = 4 to n = 18 units cells. For each supercell size, we separate the point defects with the maximal distance within the simulation box. Subsequently, the defect formation energy is calculated as a function of their separation distance. As previously shown by Potashnikov et al. [START_REF] Potashnikov | High-precision molecular dynamics simulation of uo2puo2: pair potentials comparison in uo2[END_REF], the defect formation energies decreases as the inverse of the separation distance. Therefore, one can extrapolate the values of the defect formation energy at an inĄnite separation. These values are noted as FP∞ and SD∞.

In the case of MOX, numerous conĄgurations for the cation arrangements are possible. Instead of studying each conĄguration systematically, we rather assess seven conĄgurations with Pu atoms randomly distributed. The defect formation energy is then averaged over the seven conĄgurations. This provides an indication of the efect of the local surrounding on the defect formation energies. All the defect formation energies with their standard deviation when applicable are reported in Table 4.1. Table 4.1 Ű Formation energies of Frenkel pairs and Schottky defects.

Overall, the defect formation energies are consistent for both potentials and our results show similar trends than previous calculations with the same potentials [START_REF] Potashnikov | High-precision molecular dynamics simulation of uo2puo2: pair potentials comparison in uo2[END_REF] [START_REF] Cooper | Modelling the thermal conductivity of uxth1xo2 and uxpu1xo2[END_REF]. The C-FP formation energies are higher than the SD formation energies which is slightly higher than the O-FP formation energies. Ab initio studies of these formation energies in UO 2 are numerous in the literature. However, there are no studies linked to PuO 2 or MOX. They suggest for UO 2 that Schottky energy is in the range of 5.6 to 10.6 eV [START_REF] Freyss | Point defects in uranium dioxide ab initio pseudopotential approach in the generalized gradient approximation[END_REF] (Geng et al., 2008, 10), the separated oxygen FP ranges from 2.6 to 5.77 eV (Yu, Devanathan, & Weber, 2009) (Dorado, Jomard, Freyss, & Bertolus, 2010, 3) and the uranium FP is 10.6 to 17.2 eV (Geng et al., 2008, 10) (Crocombette, Jollet, Nga, & Petit, 2001, 10). Potashnikov shows formation energy values in agreement with these ab initio values, whereas, Cooper shows values slightly higher. Qualitatively, from these results, both potentials show no drawback for posterior studies since both are in agreement with DFT values.

For the SD defect, the SD⟨1 1 1⟩ conĄguration which possesses the longest distance between the oxygen vacancies shows the lowest energy. This is in agreement with previous MD simulations in UO 2 [START_REF] Cooper | A manybody potential approach to modelling the thermomechanical properties of actinide oxides[END_REF] but in contradiction with ab initio calculations with smaller supercells for which the SD⟨1 1 0⟩ is the most favourable [START_REF] Vathonne | Dft u investigation of charged point defects and clusters in uo 2[END_REF]. Cooper suggests that this diference is caused by an artefact of the relatively small supercell used in the ab initio calculations.

Due to the large Ćuctuation in the formation energies for MOX, the Pu content inĆuence must be handled with care. However, results obtained for pure urania and plutonia show for both potentials that C-FP formation energies increase from urania to plutonia whereas O-FP and SD formation energies decrease. Furthermore, the amplitudes of the standard deviations calculated with Cooper potential are signiĄcantly higher than those calculated with Potashnikov potential with the same MOX conĄgurations. This suggests that Cooper potential is more sensitive to the local conĄguration surrounding the defect. But a more detailed study on the local environments is needed to conĄrm this hypothesis.

Frenkel pair recombination

During displacement cascades numerous FPs with small distances between the interstitial and the vacancy are created. In UO 2 , most of them do not survive and recombine within a few picoseconds [START_REF] Van Brutzel | Displacement cascade initiated with the realistic energy of the recoil nucleus in uo2 matrix by molecular dynamics simulation[END_REF] [START_REF] Martin | Efect of the cascade energy on defect production in uranium dioxide[END_REF]. Assessing the lifetime of such defects is important to interpret the results of displacement cascades and set the parameters for the point defect accumulation method. In this section, we investigate the lifetimes of close C-FPs for both potentials at three temperatures T = 30, 300 and 1600 K, and for Ąve Pu contents: 0, 25, 50, 75, 100%. Previous studies have suggested that the anion sublattice recovers accordingly to the cation sublattice to form perfect Ćuorite structure [START_REF] Crocombette | Atomistic simulation of amorphization thermokinetics in lanthanum pyrozirconate[END_REF]. Also, we have just veriĄed that for the Ąrst and the second neighbour the O-FP recombine spontaneously.

To calculate the recovery time of the defect, the supercell is Ąrst relaxed for 5 ps in the NPT ensemble in order to Ąnd the temperature dependent lattice parameter. Subsequently, a C-FP is created and the system is relaxed under NPT ensemble during 100 ps. If annealing of the C-FP occurs, the time of recombination is recorded. The time of 100 ps for the relaxation is chosen according to the previous studies of FPs recovery in UO 2 , for which most of the lifetimes are found under this time limit (Van Brutzel, Chartier, & Crocombette, 2008, 2) (Devynck, Iannuzzi, & Krack, 2012, 18). It is worth to introduce here that longer periods of time are achieved using aKMC as shown in Chapter 5.

The FPs are classiĄed according to the distance between the interstitial and the vacancy as Ąrst-, second-, third-and fourth-nearest neighbours as shown in Figure 4.1. Due to the Ćuorite symmetry, there are: 6 Ąrst, 8 second, 24 third, and 48 fourth possible nearest neighbours for the C-FPs. For simplicity, we limit the number of third-and fourth-nearest neighbours studied to 8 (all shown in Figure 4.1). For each MOX case, seven diferent conĄgurations with Pu randomly distributed are investigated. The time of recovery is then averaged over 42 for the Ąrst neighbour (7 × 6), and 56 (7 × 8) for the second, third and fourth neighbours. If for a speciĄc microstructure there is no recombination, the case is simply excluded from the average. Therefore, the resulting lifetime corresponds to the lowest limit of the recombination time. Nevertheless, it gives an indication of whether recombination may occur or not. Table 4.2 compiles all the results of the averaged lifetimes with their standard deviation. The symbol, ≪, corresponds to spontaneous recombination i.e., inferior to 1 ps whereas the symbol, >, signiĄes that no recombination takes place within the 100 ps of the simulation for all the cases studied.

1. First-nearest neighbours always recombine instantaneously (recombination occurs in less than 1 ps) for all cases.

2. In contrast, second-nearest neighbours do not recombine for all cases. This is due to the energy barrier imposed by the oxygen laying between the interstitial and the vacancy, as already seen in previous studies in Ćuorite-like structure using rigid ion interatomic potentials [START_REF] Pannier | Molecular dynamics study of frenkel pair recombinations in Ćuorite type compounds[END_REF]. However, investigation with another method (aKMC see Chapter 5) that allows more relaxation time recombination occurs via complex process. 3. For the third-and fourth-nearest neighbours the recombination process is rather complex and involves two or three successive replacements as previously shown in UO 2 (Van Brutzel et al., 2008, 2) (Devynck et al., 2012, 18). Furthermore, there is a clear dependence on the temperature. At higher temperatures, the process of recombination is more probable, indicating a thermal activated process.

4. The recombination process is more probable and occurs at lower temperatures for Potashnikov potential than for Cooper potential.

5. Concerning the Pu content, for Cooper potential it seems that the recombination processes are facilitated with the increase of the Pu content. For Potashnikov potential this tendency is less pronounced especially for the third-nearest neighbours.

6. First and second neighbours for the oxygen FPs present recombination processes that last less than 1 ps.

7. Recombination processes seem to be more efective for Potashnikov potential than for Cooper based on both; the recombination times and temperatures at which they occur.

Overall, recombination of C-FPs falls into two regimes: Ąrst, a spontaneous regime where recombination events occur in less than 1 ps and second, a thermally activated regime. This is important for the parameter of our next study using accumulation of FPs. Therefore, the relaxation time between integration of FPs will be set to 2 ps, which gives enough time for spontaneous recombination.

Primary damage state

This type of simulation helps us to understand the aging of materials. After being irradiated in the reactor, the fuel spends many years in storage. Self-irradiation processes happen inside the fuel. α decay will create during few picoseconds numerous defects. In contrast, β and γ decays will release their energy mainly through electronic efects [START_REF] Robinson | Basic physics of radiation damage production[END_REF]. The damage production can in most cases be divided into two categories: the primary damage that is formed immediately (within a few picoseconds) after the ion/neutron/electron impact by atomic collision processes far from thermodynamic equilibrium, and the long-time scale (nanoseconds to years) damage evolution caused by thermally activated processes. The primary damage state is studied via simulations of displacement cascades. The Ąrst step of the displacement cascade method consists of relaxing the supercell under constant temperature (300 K) and zero pressure during 25 ps. The collision sequence starts by choosing a cation (Pu or U), which is coined hereafter PKA for Primary Knock-on Atom to mimic a recoil nucleus. Subsequently, kinetic energy is given to the PKA. The system is then relaxed under the microcanonical ensemble for around 50 ps. Three PKA energies are studied: 5, 10 and 75 keV. The highest energy correspond to a realistic energy of the recoil nucleus during an α-decay [START_REF] Van Brutzel | Displacement cascade initiated with the realistic energy of the recoil nucleus in uo2 matrix by molecular dynamics simulation[END_REF]. The size of the supercell is chosen so that it contains all the cascade body to avoid cascade self-interaction due to the use of the periodic boundary conditions. For instance, it is 38 × 38 × 38 nm for 75 keV PKA the maximum size reached. Moreover with this box size, at the end of the cascade just a small increase of temperature occurs, less than 150 K. Therefore, no signiĄcant change is expected for the recombination or difusion process of point defects. Hence, no rescaling temperature layer is implemented to release this excess of thermal energy. Due to computational time, only three Pu contents are assessed: 0, 50 and 100%. For 5 and 10 keV PKA energies, 15 diferent cascades are studied for each Pu content. Initial PKA orientation and localization vary between each cascade as well as the micro-structure in order to increase statistics. Pu atoms are randomly distributed in the supercell. Due to the intensive computational time, only 5 cascades with PKA energy of 75 keV for each Pu content are studied.

The end of the cascade is reached when the temperature of the simulation box is constant. At that point, we analyse the number of permanently displaced atoms i.e., atoms that has been moved by more than 0.1 nm from their initial crystallographic position and the number of defects, namely the number of FPs. This is achieved by the Wigner-Seitz cell analysis (Section 2.2.1). With the same software, further analysis is completed via dislocation analysis (DXA) (Section 2.2.3) with the cation sublattice to identify disordered regions i.e., regions where the cation structure is not in fcc structure.

A classical cascade behaviour is observed. Within the Ąrst picosecond, the PKA sufers several collisions releasing its kinetic energy in the form of a large amount of atomic displacements in all directions. These atomic displacements induce a temporary local disorder of the crystal in the core of the cascade. Numerous FPs are created during this stage. Subsequently, the system evolves releasing energy in the form of thermal vibrations, causing a drastic increase of the local temperature. This last stage is usually named the thermal spike. Consequently, numerous recombination events take place, which mainly lead to the reconstruction of the Ćuorite lattice leaving only a few remaining point defects. The error bars correspond to the standard deviation from the average over the diferent cascades. This number is compared to the theoretical number of displaced atoms estimated by the Kinshin-Pease linear law [START_REF] Kinchin | The displacement of atoms in solids by radiation[END_REF]:

N KP disp. = E PKA 2E d (4.3)
herein, E PKA , is the initial kinetic energy of the projectile and E d is the threshold displacement energy in the material. The value of E d is taken as classically equal to 40 eV and 20 eV for the cations and the anions respectively [START_REF] Soullard | High voltage electron microscope observations of uo2[END_REF].

Both potentials show similar results. The number of displaced atoms is higher than the theoretical estimation, except for cations at 5 keV PKA, indicating that the E d for the empirical potentials are probably lower than the ones used in Equation (4.3). Furthermore, the ratio between the number of displaced anions and the number of displaced cations is approximatively equal to 4.3 and 5.5 for Cooper and Posthnikov potentials respectively. Consequently, the anion sublattice is more afected than the cation sublattice by the cascade event. We also Ąnd no inĆuence of the Pu content on those results.

The numbers of remaining FPs at the end of the cascades are reported in Figure 4.3. They compare with the theoretical linear law proposed by [START_REF] Norgett | A proposed method of calculating displacement dose rates[END_REF] (NRT) who rescale the Kinshin-Pease equation by a factor 0.8 as follow:

N NRT FP = 0.4E PKA E d (4.4)
We can observe a slight increment in the number of FP with the increase of the Pu content, especially for Cooper potential. The number of FPs found at the end of the cascades is signiĄcantly lower than the one estimated with the NRT law. This diference resides from the fact that realistically the NRT law does not include the defect recombination during thermal spikes as demonstrated with very high energy cascades in UO 2 by Crocombette et al. [START_REF] Crocombette | Molecular dynamics simulations of high energy cascade in ordered alloys: defect production and subcascade division[END_REF]. However, the number of FPs is signiĄcantly higher with Cooper potential especially at high energy PKA. This is directly linked to a disordered structure that will be mentioned below in which most of the atoms are detected as defects. The ratio between the number of O-FPs and C-FPs is close to 3.1 and 2.2 for Cooper and Potashnikov potentials respectively.These ratios are close to the stoichiometry and lower than those found for the number of displaced atoms. This suggests that the oxygen defects recombine with a higher rate than the cation defects and follow the cation disorder. It also shows that recombination of defects is more efective with Potashnikov potential in agreement with results presented in Section 4.2. All these results conĄrm that the long-time kinetics of defect recovery with the Cooper potential play a crucial role and this explains the unexpected formation of the disordered region found at high energy PKA. Two explanations can be put forward to explain this behaviour. First, the disordered structure is more energetically favourable. This seems unlikely because the Ćuorite structure has been found by [START_REF] Cooper | A manybody potential approach to modelling the thermomechanical properties of actinide oxides[END_REF] to have the lowest cohesive energy compared to some other crystallographic structures (see Section 3.1.7) and because up to 10 keV PKA reconstruction of the Ćuorite structure is observed. Second, the long time kinetics of defects recovery with Cooper potential is too long compared to the time of MD simulation. Hence, the disordered structure could be trapped in a metastable state.

We also analysed how the defects are spatially distributed. Spatial repartition of point defects at the end of the cascade is determined with a cluster analysis. We deĄne a cluster as connected defects of the same type (vacancy or interstitial). The cut-of distances for connected defect is respectively set to 0.237 nm and 0.386 nm for vacancy and interstitial, which correspond to the distance of Ąrst-nearest neighbours in the Ćuorite structure. The number of members of each cluster deĄnes its size. Figure 4.6 shows the ratio in percentage between the number of defects in cluster of a Ąxed size and the total number of defects as a function of cluster size. In Appendix C, we can Ąnd the results of the clustering test for each PKA energy and Pu content for both potentials, including cation and anion separately. Only the results from 10 keV PKA are presented in this section because it was the most comparative option we have, since, the presence of disordered regions at high energy with the Cooper potential does not permit to carry out the clustering test. For both potentials no signiĄcant diference is found with the diferent Pu contents. Therefore, the results of the three Pu content are brought together to increase the statistics. Is it also important to highlight that these Ągures do not show the contribution of single point defect which contribute to the majority of defect. Hence, the total shown in the Ągures does not reach 100%. Overall, both potentials show the same general behaviour even if it is more predominant for Potashnikov potential. Only small interstitial clusters containing up to 4-5 interstitials are found. Conversely, large vacancy clusters up to 22 vacancies corresponding to nanocavities of 1 nm size are found. These large cavities are even more present with Potashnikov potential. One possible explanation for this behaviour comes from visual inspection of the spatial repartition of defects. It shows that vacancy clusters are mainly created in the core of the cascade whereas interstitial clusters are located preferentially at the periphery. Consequently, all remaining vacancies are close to one an other and are more likely to form clusters.

Dose effect

Displacement cascades provide characterization of the primary damage state, which mainly accounts for radiation on single crystal. Therefore, in order to study the dose efect (accumulation of damage) on the MOX microstructure, we use the defect accumulation method. With this method we introduce periodically and consistently FP defects, while the structure is continuously allowed to relax under constant pressure and temperature. The FPs are randomly introduced with minimum distance between the vacancy and the interstitial greater than 1.4 Å. In Section 4.2, we showed that the oxygen sublattice seems to follow the cation disorder [START_REF] Crocombette | Molecular dynamics studies of radiation induced phase transitions in la2zr2o7 pyrochlore[END_REF], consequently only C-FPs (i.e., U or Pu) are introduced in order to reduce the computational time.

A common way to measure irradiation dose on materials is the use of displacement per atom (dpa). This is also a very simple metric from the atomistic simulations point of view. It is deĄned simply by the number of displacements (on average) of one atom, which exactly corresponds to what is simulated with the FPA method. Since only C-FPs are created, we express the dose in displacement per cation (dpc) rather than the classic dpa.

The time between two subsequent C-FP introductions is set to 2 ps. This time corresponds with the Ąrst regime of recombination discussed in Section 4.2. Thus, we allow the system to have spontaneous recombination events. At each C-FP introduction sequence i.e., every 2 ps, 700 FPs are created for Potashnikov potential and 300 for Cooper potential. Therefore, for each C-FP introduction, the dose is increased by 1.17 × 10 -3 and 2.73 × 10 -3 dpc for Cooper and Potashnikov potentials respectively. These numbers have been chosen to minimize the computational time and to ensure that within the time-lapse between two FP introductions the system is relaxed. Namely, the temperature and the pressure are back to the desired values. Because the defect recombination time is higher for Cooper potential (see Section 4.2 ), less FPs can be introduced within the same time interval. An alternative could have been to increase this time interval, but the computational time to reach equivalent dose would have been out-of-range. To verify that this diference in dose rate does not afect the result, we carried out one simulation with Potashnikov potential introducing 300 FPs every 2 ps.

It is worth to mention that experimental dose rates are at least 12 order of magnitudes higher than in the FPA method (10 -4 dpc/s in comparison to our MD of about 10 9 dpc/s). However, A. Chartier suggests in ( [START_REF] Chartier | Early stages of irradiation induced dislocations in urania[END_REF]) a justiĄcation to our method. He argues that dose rates from 10 -3 dpc/s and higher fall into the, so called, recombination regime where mainly athermal nearby FP recombination yields to the restructuring of the Ćuorite structure.

In that regime, no thermal difusion participates to the defect annealing. This is due to the very low uranium self-difusion coeicients, even in presence of defects, which ranges between 10 -31 to 10 -25 m 2 /s. Rapid calculation shows that the time necessary to recombine a FP with a distance between the vacancy and interstitial of 2.4 nm (typical distance between two displaced atoms for a dose rate 10 -3 dpc/s) is greater than 10 3 s. Therefore, FPA simulations model the same mechanism than in the experimental observations. Moreover, [START_REF] Rest | A model for the inĆuence of microstructure, precipitate pinning and Ąssion gas behavior on irradiation induced recrystallization of nuclear fuels[END_REF] shows that in the recombination regime at steady state the dislocation density, ρ, scales with the dose rate K 0 as follows ρ = cK 1/6 0 . Therefore, only rescaling of 10 -3 needs to be applied to compare our MD results on dislocation density with the experimental data.

Due to the high increase in energy arising from each FP introduction, the NPT relaxation is done with the Berendsen algorithm provided by the LAMMPS code Section 2.1.6 in order to reach the desired temperature within the 2 ps. Two diferent temperatures are investigated 300 and 1600 K as well as three Pu contents 0%, 50%, and 100% for both potentials. These two temperatures were chosen in concordance with our study in Section 4.2 that shows a thermal activated recombination at high temperatures. In all the cases studied herein, the size of the system includes 40 × 40 × 40 Ćuorite unit cells, involving 256000 cations. Pu at 1600 K and 300 K . Only, dislocations and main zones where structure departs from perfect Ćuorite structure are reported in these Ągures for clarity. Overall, both potentials show the same general trend for the evolution of the microstructure with the dose at 1600 K. However, diferences appear in the dose at which it occurs. The evolution of the system follows 3 main stages. First, point defect concentration increases without formation of a dislocation (see snapshot 0.2 dpc with cooper potential). Subsequently, at about 0.5 dpc for Cooper and 0.1 dpc for Potashnikov Frank loops with BurgerŠs vector 1 /3⟨1 1 1⟩ nucleate from interstitial clustering. This point is illustrated by the fact that the concentration of interstitials decreases at the dose at which loops nucleate while concentration of vacancies still increases (see Figure 4.7 and Figure 4.8). As the dose increases, these loops transform into perfect loops with BurgerŠs vector 1 /2⟨1 1 0⟩ following this way the "unfaulting" process shown in Section 1.6.2. However, for Cooper potential the number of Frank loops is considerably smaller, which suggests that most of the these loops are highly unstable. Finally, perfect loops grow and self-organise into dislocation lines (last snapshots). The same behaviour was found in UO 2 with MD simulations [START_REF] Chartier | Early stages of irradiation induced dislocations in urania[END_REF] and observed experimentally [START_REF] Onofri | Evolution of extended defects in polycrystalline au irradiated uo2 using in situ tem: temperature and Ćuence efects[END_REF] (Onofri, Sabathier, Palancher, et al., 2016). Howver, for the experimental observations no Frank loops could be detected because of the limitation of the microscope resolution. It is worth noting that for Cooper potential, contrary to what is observed at high energy cascades, no large high disordered region (> 5 nm) is found. Only, small (1-2 nm) regions composed of defect clusters along with prismatic dislocations are found. This leads us to believe that, like for the Potashnikov potential, radiation-induced microstructure likely involves point defects and dislocations instead of highly disordered area. In that case, it conĄrms our hypothesis that with Cooper potential the kinetics of recovery for the high disordered structure is too long compared to the time-frame of the MD simulation.

Temperature plays an important role in the evolution of the extended defects. With Cooper potential no dislocation formation is observed at 300 K whereas some appear for 1600 K. This diference could come from the fact that defect recombination at 300 K for Cooper potential is highly improbable as shown in Section 4.2 For Potashnikov potential some quantitative diferences occur between 300 and 1600 K. Figure 4.11 shows the evolution of the dislocation densities as a function of dose obtained at 300 and 1600 K using the same initial conĄguration (50% Pu content) with Potashnikov potential. For both temperatures, one can observe easily the same sequence of the diferent types of dislocation appearing as the dose increases. First, Frank loops nucleate, peak and decrease abruptly corresponding to the increase of perfect loop density. This transformation is corroborated by the appearance of Shockley partials (BurgerŠs 1 /6⟨1 1 2⟩) at the dose where Frank loops decrease and perfect loops increase. At the higher doses the density of perfect loops stabilizes and a steady state appears. This behaviour is shifted to lower doses as the temperature increases. Peaks of Frank loop and perfect loop density appear successively at 0.4 and 0.8 dpc for 1600 K whereas they appear at 0.6 and 1.2 dpc for 300 K. Moreover, at 1600 K Frank loops are sill present at high doses with density oscillating around 5 × 10 14 m -2 while they are almost non-existent at 300 K. This suggest that at high temperature Frank loops are still created inside the nano-domains delimited by the unfaulted dislocations. This behaviour changes with temperature has been observed experimentally in UO 2 [START_REF] Onofri | Evolution of extended defects in polycrystalline au irradiated uo2 using in situ tem: temperature and Ćuence efects[END_REF]. In this irradiation study performed at -180 • C and 600 • C, the dose at which the loops transform into lines decreases with increasing temperature and loop growth does not occur for the low temperature. The authors suggest that this ofset in the defect changes is attributed to higher defect mobility under irradiation as the temperature rises. We discuss in the following the inĆuence of the Pu content on the primary radiation damage. Figure 4.7 shows the evolution of the density of dislocations as a function of dose obtained with Cooper and Potashnikov potentials at 1600 K. A clear trend is found with Cooper potential. The dislocation density decreases as the Pu content increases and peaks of density shift to higher dose with increasing Pu content. Therefore, the creation of dislocations is easier in urania than in MOX or plutonia. As a hypothesis, we suggest that this behaviour is caused by the fact that for Cooper potential the recombination processes are facilitated with the increase of the Pu content (see Section 4.2). Thus, more defects in urania are able to create dislocations, such as, Frank dislocations and initiate the unfaulting process that Ąnish in perfect dislocations. In contrast, defects being annealed faster in plutonia will create less Frank dislocations and therefore less perfect dislocations. It is worth noting that for this potential, few Frank loops nucleate implying that they are highly unstable. With Potashnikov potential, the Pu inĆuence on the dislocation density is less signiĄcant, even if above 2 dpc the density of 1 /2⟨1 1 0⟩ dislocations seems slightly higher for the MOX than for plutonia.

Dependency on Pu content can also be found with the evolution of the number of point defects (vacancy and interstitial). Figure 4.8 displays this evolution carried out at 1600 K for both potentials as a function of dose. For both potentials, a similar trend is observed. The number of vacancies and interstitials slightly increases as Pu content increases. Since interstitial defects feed dislocations, this behaviour corroborates previous conclusion stating that less dislocations are formed as Pu content increases. A similar clustering test to the one in Section 4.3 for Ąnding clusters was carried out for this section. Results show that there is no clustering of cation species (U clustering with U or Pu clustering with Pu). As maximum, clusters of around 7 cations of the same specie can be found along the whole process.

Elastic moduli vs dose

Radiation induced damage afects the structure and the thermomechanical properties of all the nuclear core components. In particular, pellet-cladding interaction can increase due to the increasing stresses induced by the fuel swelling onto the cladding [START_REF] Michel | 3.22 modeling of pellet cladding interaction[END_REF] [START_REF] Baurens | 3d thermo chemical mechanical simulation of power ramps with alcyone fuel code[END_REF]. Therefore, to ensure the stability and durability of the nuclear reactor core in these operating conditions, a solid knowledge of the thermomechanical properties under irradiation of the fuel is needed. For this reason, we assess some of the mechanical properties as a function of dose in this section. The elastic mechanical properties are then evaluated as a function of the irradiation dose. Figure 4.12 and Figure 4.13 compare respectively the evolution of the bulk modulus and of the Zener ratio as a function of dose calculated at 1600 K for pure UO 2 and PuO 2 and for a solid solution containing 50% Pu.

For both potentials we observed the same behaviour for the bulk modulus. It decreases rapidly for doses less than 0.1-0.2 dpc, and becomes almost constant for the highest dose levels. The decrease is sensibly the same for all Pu contents indicating no real inĆuence of the Pu content on the elastic moduli after accumulation of damage, even with Cooper potential. This decrease at low doses is strongly correlated with the rapid increase of the point defects for the same doses. This indicates that point defects contribute for the most part in the reduction of the bulk modulus at low doses. This behaviour has already been reported by (Gao & Weber, 2004, 22) for another ceramic material, SiC. Furthermore, the analysis of simulated XRD pattern in UO 2 [START_REF] Chartier | Early stages of irradiation induced dislocations in urania[END_REF] shows that the lattice expansion at low doses is mainly related to point defects accumulation, while the onset of lattice contraction starts when perfect dislocations nucleate. The bulk modulus being a measure of the resistance to compressibility, it is not surprising that a decrease of the bulk modulus promotes swelling.

Contrarily, the evolution of the anisotropy factor really difers between the two potentials. For Cooper potential, the Zener ratio Ąrst exhibits a rapid increase and saturates for higher doses, corresponding to a decrease of the anisotropy. The Pu content also has an impact; the richest Pu content shows higher anisotropy. This behaviour is similar to the one found for bulk modulus and indicates that point defects play a major role in this evolution.

Conversely, for Potashnikov potential, which renders isotropic feature, the Zener ratio drops rapidly to a minimum of Z = 0.87 before increasing around 0.2-0.3 dpc and returns to almost perfect isotropic behaviour at high doses. This short anisotropy increase is related to point defect creation, but also to the creation of Frank loops, which almost reaches its peak density at that dose. After 0.5-0.6 dpc, Frank loops transform into perfect dislocations and at that point Zener ratio almost reaches back to its value for non-irradiated (U,Pu)O 2 .

X-ray powder Diffraction

Simulated X-ray difraction (XRD) have been generated with the code DEBYER as exposed in section Section 4.6, for snapshots corresponding to diferent doses for both potential at 1600 K. Each conĄguration has been relaxed 50 ps under NPT to allow extra annealing of the defects. The atomic positions are then taken as input for the DEBYER code. Subsequently, XRD patterns are generated and analysed.

As an example, Figure 4.14 shows the evolution of the XRD patter for UO 2 as a function of dose obtained with Cooper potential. We observe easily a big shift of the peaks as the dose increases from 0 to 0.3 dpc. This is correlated to the increase of the lattice parameter and hence to the swelling to the system. It is worth noticing that this dose corresponds as well to the dose at which the number of point defects is the highest (see Figure 4.8). Hence, we can correlate the point defects contribution directly to the swelling. This is congruent with studies in UO 2 about swelling by [START_REF] Chartier | Early stages of irradiation induced dislocations in urania[END_REF].

Another interesting feature of the XRD pattern is the analysis of the full width at half maximum intensity (FWHM) of the peaks. It provides indication about the disorder and the non-uniform strain inside materials. However, it is not usually related directly to a particular defect. Since our MD simulations ofer a detail description of the microstructure we can overcome this diiculty.

We analysed the FWHM of the Ąrst peak (1 1 1) because it is the most intense. We can easily notice from the evolution of the FWHM peaks that it correlates with the perfect dislocation densities along the whole range of dose for both potentials. The link between the broadening of the difraction peaks and a dislocation density is well known and the theory dates back to the works of Wilkens [START_REF] Wilkens | The determination of density and distribution of dislocations in deformed single crystals from broadened x ray difraction proĄles[END_REF]. However, our results indicates that the broadening is only due to the perfect dislocation density, implying that the contribution of Frank loops is negligible. Two points can be put forward to explain this observation. First, the magnitude of the Burgers vector of the Frank loop b = a /3⟨1 1 1⟩, is lower than that of the perfect dislocations b = a /2⟨1 1 0⟩. Because the strains Ąelds generated by dislocations is proportional to the magnitude of the Burgers vector, a smaller broadening is expected with smaller Burgers vectors. Second, the strain Ąeld generated by a dislocation loop decreases much faster with the distance (1/r 3 decay) than for a straight dislocation (1/r decay).

Conclusions

In this section, we assessed the primary radiation damage in (U 1-y ,Pu y )O 2 solid solution with molecular dynamics simulations. Two interatomic potentials were studied, as coined by Cooper and Potashnikov. These potentials were chosen according to the conclusions of our previous thermomechanical assessment in Chapter 3. The radiation assessment consisted of six diferent studies: 1) defect formation energies; 2) Frenkel pair recombination; 3) Primary state damage via displacement cascades; and 4) study of radiation dose efect with Frenkel pair accumulation method; 5) elastic moduli vs dose and 6) X-ray powder Difraction. In each For the defect formation energies, both potentials give approximately the same values and corroborate with existing data for urania. Concerning the FP recombination, the recombination of C-FPs falls into two regimes: Ąrst, a spontaneous regime where recombination events occur in less than 1 ps and second, a thermally activated regime. Recombination processes seem to be more efective for Potashnikov potential that for Cooper based on both; the recombination times and temperatures at which they occur.

From the simulation of displacement cascades, we obtained the classical cascade behavior; within the Ąrst picosecond, the PKA sufers several collisions releasing its kinetic energy in the form of a large amount of atomic displacements in all directions. These atomic displacements induce a temporary local disorder of the crystal in the core of the cascade. Numerous FPs are created during this stage. Subsequently, the system evolves releasing energy in the form of thermal vibrations, causing a drastic increase of the local temperature (thermal spike). Consequently, numerous recombination events take place, which mainly leads to the reconstruction of the Ćuorite lattice leaving only a few remaining point defects. Diferences start to appear at 75 keV. For the cooper potential, the recovery stage does not lead to full reconstruction of the Ćuorite structure. Instead, a large disordered region with no crystallographic atomic structure of several nanometres is found at the cascade core. This behaviour is probably linked to the slow kinetics of Frenkel pair recombination. In contrast, Potashnikov potential presents a quasi-total reconstruction of the Ćuorite structure.

The FP accumulation study was carried out to investigate the dose efect. Whatever the Pu content, Potashnikov potential presents a similar evolution in three main stages as the dose increases: accumulation of point defects that cluster and form Frank loops which in turn transform into perfect loop and dislocation lines. This behaviour is diferent for Cooper potential. With this potential at 300 K no dislocation is found. Nevertheless, at 1600 K dislocations form but only few Frank loops are observed suggesting a direct mechanism for the creation of perfect loops or that they are highly unstable. Hence, the temperature plays a major role.This diference between both temperatures 300 and 1600 K could come from the fact that defect recombination at 300 K for Cooper potential is highly improbable as shown in the FP recombination study. Concerning the Pu content, a clear trend is found. The dislocation density decreases and dislocations appear for higher doses as the Pu content increases.

Concerning the elastic moduli vs dose, for both potentials we observed the same behaviour for the bulk modulus. It decreases rapidly for doses less than 0.1-0.2 dpc, and becomes almost constant for the highest dose levels. This decrease at low doses is strongly correlated with the rapid increase of the point defects for the same doses. Moreover, point defects are as well the most important factor for changes in the Zener ratio, since it drops rapidly to a minimum of Z = 0.87 where the point defects are at their maximum.

The X-ray powder difraction study shows that a high number of point defects produce the most notorious shift in the difraction peak. Hence, we can correlate the points defects contribute the most to swelling. In contrast, the evolution of the FWHM is highly correlated with the perfect dislocation densities along the whole range of dose for both potentials.

Overall, this assessment shows that both potentials have advantages and disadvantages when studying radiation-induced damage in (U,Pu)O 2 solid solution. Cooper potential seems to be more sensitive to Pu content than Potashnikov potential. However, with Cooper potential, displacement cascades at high energies lead to very high disordered microstructure, which seems to be contradictory with experimental results that show no amorphisation for UO 2 or MOX even at high doses. Nevertheless, microstructure evolution with radiation dose studied with FPA method shows no amorphisation. We hypothesize that this diference of behaviour could be explained by the kinetics of the defect recovery. Therefore, this potential could be more suitable for fast kinetic radiation phenomena or for methods, like kinetic Monte-Carlo, where time-frame is less constrained.

Our journey to understand the behaviour of systems along time is not a simple task. Atomistic scale systems can show a large gamut of phenomena during their time evolution. In all the cases, the conĄguration of atoms evolves along time caused by either just thermal agitation or another more complex phenomena. How the atoms interact between each other has always been a question for the human kind. We have tried to answer this question since we rationalized the concept of an atom. During this long process, some new human tools have been created such as computers. This new tool has given us the capability to search along a vast quantity of scenarios that a pen and paper have not been able to get us close to. The increasing computational power allowed us to solve numerically mathematical forms that represent the interaction between atoms. The higher the computational power is, the longer the systems being analysed can be. This is the case of MD which is presented in Chapter 2. Generally, the MD approach that solves the motion of atoms using classic mechanics has proven to be enough for a large number of cases. However, one of the strongest properties of MD which is showing the "exact" position and velocity of each atom is as well one of its weak points. Thermal vibration frequencies are in the order of 10 13 Hz, thus, in order to reach laboratory times or to simulate "rare events" using MD, it would take hundreds of human lifespans using the most powerful super computer available nowadays. This leads us to look for another way to deal with the problem of systems evolving slowly a long time.

Transition state theory

This theory focuses on rare events. Atomic systems usually can be found trapped in a local potential energy minimum for long periods. Due to the complexity of those systems (various atomic conĄgurations, temperatures, pressures, etc), we can expect that there are many other local potential energy minimums available. The rate of shifting from one minimum to another is very important to understand the behaviour of systems along time. This is the aim of the Transition State Theory (TST). It has been a decisive theory in the physics, chemistry, engineering and biology domains for a long time. Henry Eyring is one of the progenitors of this theory. This can be found in his works about chemical reactions published around the year 1935 [START_REF] Eyring | The activated complex in chemical reactions[END_REF] and [START_REF] Glasstone | Fracture Mechanics: Fundamentals and Applications[END_REF]. This theory explains how in a chemical reaction, the reactants are in stable state with a low potential energy. After, the reactants form a complex which has higher energy and is unstable. The complex decays rapidly forming products with lower potential energy. An application to solids was developed by (Wert & Zener, 1949, 8) and [START_REF] Vineyard | Frequency factors and isotope efects in solid state rate processes[END_REF]. The next step is to try to get a mathematical expression for the transition rate from one state to another k from basic concepts. First, letŠs note that each atomic conĄguration of the system under study has a corresponding potential energy. Thus, we can deĄne a potential energy surface such as the one deĄned by the value of the potential energy at each atomic conĄguration. This potential energy surface is usually multidimensional. It would be convenient to make a projection of this multidimensional space in order to work just with one scalar which is called reaction coordinate. A reaction coordinate can be for example the distance between two atoms or molecules.Figure 5.1 shows a one dimensional potential as a function of a reaction coordinate with the most important points highlighted.

It will be assumed that the reaction coordinate x in this case is the distance between two particles. We can deduce a mathematical expression for the rate k using the Boltzmann-Maxwell distribution [START_REF] Toda | Statistical Physics I[END_REF]. At equilibrium, the probability distribution as a function of positions x and velocities v is given by

P(x, v) = exp - ( 1 kT ( 1 2 mv 2 + V(x)) ) ∫ ∫ exp - ( 1 kT ( 1 2 mv 2 + V(x))
) dxdv

(5.1)

LetŠs focus our attention to the particles ubicated close to the potential energy maximum E + . This maximum will be called saddle point. More strictically we deĄne a saddle point such as a point on a potential surface at which the force is zero and at which there is one negative curvature. Since any particle with an x that surpasses this energy will fall into the potential energy V b and in the other sense it will fall into the potential energy V a . The latter can be expressed with the next inequalities.

x + -v∆t > x < x + v > 0 (5.2)
The Boltzman distribution will help us to calculate the probability P c (∆t) of crossing the saddle point in a time ∆t.

P c (∆t) = ∫ ∞ 0 dv ∫ x + x + -v∆t dx exp -1 kT ( 1 2 mv 2 +V(x)) ∫ ∫ dvdx exp -1 kT ( 1 2 mv 2 +V(x))
(5.3) using Equation ( 5.2) we have

P c (∆t) = ∫ ∞ 0 dv exp -1 kT ( 1 2 mv 2 ) ∫ x + x + -v∆t dx exp -1 kT V(x + ) +O(x -x + ) ∫ ∫ dvdx exp -1 kT ( 1 2 mv 2 +V(x))
(5.4) solving

P c (∆t) = ∫ ∞ 0 dv exp -1 kT ( 1 2 mv 2 ) v∆t exp -1 kT V(x + ) ∫ ∞ -∞ dv exp -1 kT ( 1 2 mv 2 ) ∫ ∞ -∞ dx exp -1 kT V(x)
(5.5) The probability that the particle is initially in the state a is

Ąnally P c (∆t) = √ kT 2πm ∆te -1 kT V(x + ) ∫ ∞ -∞ dx exp -1 kT V(x) (5
Pa = ∫ Ω dx∆te 1 kT V(x) ∫ ∞ -∞ dx∆te -1 kT V(x)
(5.8) where ∫ Ω is the unidimensional surface integral in the state A. Using Equation (5.7), Equation (5.8) and Equation (5.6) we have that

k(T) = √ kT 2πm e E + -Ea kT ∫ Ω e -1 kT V(x)-E a (5.9)
Here is when the harmonic approximation will be taken in to account by using a Taylor approximation of V(x) around x a . This approximations afects directly the prefactor in Equation (5.13). Its advantages and disadvantages are discussed in the next paragraphs.

V(x) = V(x a ) + 1 2 mw 2 (x -x a ) 2 + O(x -x a ) 3 ≈ Ea + 1 2 mw 2 (x -x a ) 2 (5.10)
w is the angular frequency. The next step is to approximate the integral over Ω to a integral over all the space, which can be done since the resulting Gaussian function will decay at large distances from x a . Then

∫ Ω dxe -1 kT (V(x)-(E + -E a )) ≈ ∫ ∞ -∞ dxe -mw 2 (x-xa ) 2 kT = √ 2kT mw 2 ∫ ∞ -∞ dye y 2 = √ 2πkT mw 2
(5.11)

Using this with Equation ( 5.9) we have a simple expression for the rate k as a function of temperature

k(T) = w 2π e -(E + -Ea ) kT
(5.12)

This process was generalized for N dimensions by [START_REF] Vineyard | Frequency factors and isotope efects in solid state rate processes[END_REF]. He used the assumption that the surface can be considered as a heat bath of harmonic oscillators as we did in Equation (5.10). The system is then described by N normal modes with angular frequencies w a n in the state corresponding to the minimum potential energy E a for all the 3N -6 normal coordinates, and w + n are the angular frequencies of vibrations at the saddle point corresponding to the maximum of potential energy E + , however, the mode corresponding to traversing the reaction coordinate which connects the saddle point with the reactants and products has an imaginary frequency and has to be excluded. Therefore, there are 3N-7 normal modes at E + .

k(T) = 1 2π ∏ 3N-6 n=1 w a n ∏ 3N-7 n=1 w + n e -(E + -Ea ) kT
(5.13)

As we can notice, all the quantities can be evaluated directly from the potential energy surface without dynamical calculations and entropic or thermal efect was included through the harmonic partition function. However, the calculation of the prefactors gives an additional computational cost which we are not strictly necessary forced to take since the relevant phonon modes are not expected to vary much from site to site whereas the exponential term can vary several orders of magnitude. Thus, practically it is often chosen to only use an order of magnitude for the prefactor. This can be done using the Debye temperature Θ D of the substrate as Γ 0 ≈ k B Θ D /h, or with the VdenskyŠs harmonic approximation Γ 0 ≈ 2k B T/h (Clarke & Vvedensky, 1987, 21). Usually, the prefactor ranges between 10 12 to 10 13 s -1 . Recent studies about the importance of the prefactor show than it has more importance that it was believed. This is the case of some difusion mechanisms that would not be accessible if using a constant prefactor. [START_REF] Lazauskas | InĆuence of the prefactor to detect motion in alpha iron during long time scale simulations[END_REF] working with the inĆuence of the prefactor on motion of defects in α -Fe found that the calculated prefactor increases the relative interstitial-vacancy difusion rates by an order of magnitude compared to the constant one. Moreover, the rate table changes, as well. He compared his results with MD and TAD (Temperature Accelerated Dynamics).

Accelerated dynamics

Now that we have a theory that gives us the transition rates, the problem is then located at identifying the states of a system and Ąnding the mechanism from a current state to a new one. This is not a simple task. The system can be composed of a large number of atoms implying a large number of degrees of freedom. Search through this complex potential energy for saddle points will spend most of computational time available. The uncertainty coming from the Ąnal states of low lying saddle points is considerable with these complex potential energy systems. There are some methods available in the literature that will help us to face this problem. These methods were invented to increase the chances of the rare events to occur, increasing in this way the real time of simulations while keeping the same computational time.

One method is to add a repulsive potential energy (bias potential) to the actual potential in order to increase the probability of rare events. This bias potentials will increase the potential energy minima in such a way that it will not modify the potential of unknown transition states. Therefore, the dynamics of this speciĄc system would not be changed. Voter has proposed a formulation of this method (Voter, 1997a) (Voter, 1997b, 20). In order to achieve this, the bias potential must vanish at the transition state. Consequently, under the harmonic TST, a energy rim near the saddle points should be stated as the cut-of for the bias potential. The evolution of the system with the added potential is in concordance with the one without in the sense that the probability of having any particular sequence of states is the same for both. The problem with this method comes at the time of building the bias potential itself. Figure 5.2 shows a bias potential constructed in such a way that it replaces the real potential by a constant equal to an energy called "boost energy" when the real potential energy is below this boost energy. Therefore, this boost energy will help to overcome the energy barriers. It is worth to emphasize that the bias potential has to be lower than the saddle point in order to maintain the dynamics untouched. Finally, the time is calculated as for regular MD but multiplied by an instantaneous boost factor, the inverse Boltzmann factor for the bias potential at that point.

The second method is called Parallel Replica Dynamics (Voter, 1998, 22). The main objective of this method is to replicate a system into several processors in order to explore more extensively the phase space. SpeciĄcally, a system is copied to N number of processors and for each case the momenta is periodically randomized to eliminate correlations between replicas. When an escape pathway is detected, a signal is sent to the other processors to stop the current simulation. After, the system where the transition happened is copied to the available processors and the cycle starts again. The simulation time is advanced taking into account all the accumulated trajectory time summed over all replicas. Therefore, this time corresponds to the total time spent at the basin until an transition of state happens. An advantage using this method is that there is no corruption in the probabilities of the possible escape paths.

The last method is the one known as Temperature Accelerated Dynamics (TAD) [START_REF] Sorensen | Temperature accelerated dynamics for simulation of infrequent events[END_REF]. It consists in increasing the temperature from a temperature of interest T l to Figure 5.2 Ű Bias potential constructed in such a way that it replaces the real potential by a constant equal to an energy called "boost energy" T h in order to drive the system to higher transition probabilities. However, transitions that should not have happened at T l have to be removed while evolving the system. The harmonic TST is used to describe the system evolution. The process consists in letting the system evolve with T equal to T h and when a transition is detected the new state basin is recorded. Consecutively, the corresponding saddle point is calculated using the Nudged Elastic Band method [START_REF] Henkelman | A dimer method for Ąnding saddle points on high dimensional potential surfaces using only Ąrst derivatives[END_REF][START_REF] Henkelman | A climbing image nudged elastic band method for Ąnding saddle points and minimum energy paths[END_REF] [START_REF] Henkelman | Improved tangent estimate in the nudged elastic band method for Ąnding minimum energy paths and saddle points[END_REF]. A list of escape paths and attempted escape times for the system at T h is recorded. Using Equation (5.13) it can be extrapolated the escape time at T l . Finally, the process with shortest time at T l is chosen to be the correct transition and the system clock is advanced by the shortest time at T l .

Other methods to explore longer periods of time, for instance k -Dynamics [START_REF] Lu | Communication: kapp dynamics an exact method for accelerating rare event classical molecular dynamics[END_REF] are equally important. For more information about these methods, a complete and comparative study of this method and the last three as well as the improvement on these techniques can be found elsewhere [START_REF] Uberuaga | Atomistic simulation methods for long time dynamics in materials for nuclear energy systems[END_REF].

The dimer method

The dimer method has the objective of Ąnding the saddle points of a system. The dimer concept comes from the fact that the process is focused on two images or two diferent replicas. Thus, the dimer is the ensemble of these two replicas. Both replicas atom coordinates difer by a small displacement. The purpose is to Ąnd the path in the potential energy surface from the local minima up towards a saddle point. A rotation is done along the way to Ąnd the lowest curvature mode at the point where the dimer is located [START_REF] Henkelman | A dimer method for Ąnding saddle points on high dimensional potential surfaces using only Ąrst derivatives[END_REF]. The lowest energy orientation is located along the lowest curvature mode which allows us to avoid calculating the Hessian matrix. This method was proposed by Voter (Voter, 1997b, 20). The process consists on letting free the dimer to rotate. Then a rotational force is calculated with the diference of the force of both replicas. A minimization is carried out to align the dimer with the lowest curvature mode. Now that dimer will orient itself along the lowest rotational mode, we need to translate it to the saddle point. A saddle point is located at the maximum along the lowest curvature mode direction and a minimum in all other directions. The dimer quantities are calculated from the energy and force acting at the two images. Translate the dimer up in the potential surface can be achieved by deĄning a efective force on the dimer at the middle point between the replicas which is simply the average between both the force of both images. The later force has pulled the dimer towards a minimum. Thus, a modiĄed force is proposed having the force component along the dimer inverted. A minimization of this modiĄed force will move the dimer to a saddle point. It was found that is most eicient to give a tolerance for the rotational force rather than fully converge the dimer orientation. This means that the dimer will rotate a few times until the rotational force reaches the cut-of saving this way some computational time. Figure 5.3 shows an example of the dimer method process for Ąnding saddle points.

One important feature of this method is that when the system size is increased, given this way an increase at the number of degrees of freedom, the number of force evaluations rises slowly. The reason is that it avoids the calculation of the second derivative Hessian matrix.

The dimer searches are started close to a local minima by moving the atom coordinates slightly. The dimension of this change is important for the eiciency of Ąnding the saddle points. There are three modes of displacing the atom coordinates. The Ąrst one is a very generic concept. It displaces all the atoms by a normally distributed random amount. The reason of doing this is that by displacing them randomly, we can hope that the new positions will result to a value of the potential energy closer to a saddle point. This does not assure us that we will Ąnd them all. However, the problem becomes just to sample them. The second method relies on the fact that we can save some computational power if we have certain information a priori of the kinetic events. Knowing beforehand this, we can focus our dimer searches to the atoms that we know are participants of the main kinetic events. This is done by choosing the "important" atoms and displacing atoms in the local region. Bear in mind that those displacement do not constrain the events to the local minimum [START_REF] Xu | [END_REF]. The Ąnal displacement mode is related to a choice of a local region. In other way, the user chooses the area of interest. The later has a meaning when addressing surface difusion and phenomena happening just in a certain part of the system. Identifying which of these three displacement mode is better to engage the eiciency at the time of Ąnding the saddle points is not usually a evident issue. Problems can arise when the system is totally unknown a priori. The user can wrongly bias the search, given this way undesirable transitions. Finally, the magnitude of the initial displacement is important. Small displacements will lead to a poor distribution of saddle points around the local minima. Large ones may bring searches at forbidden energies or others that lead to a saddle points that do not lead to the initial minimum. For this reason, minimisations are carried out for both directions around the saddle to point. If the local minimum is not reached, the saddle point is discarded.

One question of the optimal number of saddle point searches is now discussed. This questions is quite important for a good simulation of the system under study. A successful search of saddle point will result in an accurate rate table that will be used by the Monte Carlo method (see Section 5.6). A stopping criterion is proposed to evaluate when there were The system is then made to go down the minimum energy path (gray lines) on either side of the saddle points which are indicated with *. Here, all four saddle points have a minimum energy path starting at the initial state minimum A, but this does not have to be the case. The rate of each process is then calculated using Equation (5.13). A process is chosen at random using the kinetic Monte Carlo algorithm. In this case (see section Section 5.6), process 1 gets chosen. The system is moved to the Ąnal state of this process, to minimum B. (c) Dimer searches are run from the new minimum, again four distinct saddle points are found. (d) Minimum energy paths are traced out, and the process repeated [START_REF] Henkelman | Long time scale kinetic monte carlo simulations without lattice approximation and predeĄned event table[END_REF] enough saddle points found. It is based on the history of previous searches to evaluate that an important saddle point has been missed. The advantages are that this criterion gives a level of conĄdence to the accuracy and it saves computational time by relating the number of searches with the number of processes that each state has. For instance, states with only a small number of processes require only a small number of searches and vice versa. If it is assumed that all saddle points are found with an equal probability as they are found, we have

C = (1 - 1 N r ) (5.14)
where N r is the number of sequential searches that Ąnd relevant but redundant (non-unique) processes. How to deduce this formula can be found elsewhere [START_REF] Xu | [END_REF]. As an example, if we set C = 95% then N r = 20 so the computer will complete saddle searches until N r = 20 without Ąnding a new, unique saddle point. The assumption that all saddle points have an equal probability to be found is quite strong. However, it can be smoothed by introducing a parameter α which is the relative probability of Ąnding the saddle point that is least likely to be found as compared to the ideal case above. Thus, we have the next correlation:

C = (1 - 1 N r α ) (5.15)
Here α runs from 0 to 1. It is easy to note that if α is equal to 1 we get the same as Equation (5.14). This parameter is added because there is some non-uniform probability distribution of Ąnding the saddle points.

It is important now to set how relevant is a kinetic event occurring in the system under study. We can recall that the number of saddle points increases exponentially with the system size. Many of these kinetic events are in a energy range that is extremely high for the type of phenomena under study. It is expected that an accurate rate table would contain all low energy saddle points with high rates. For this reason, a parameter has to be introduced in order to let us truncate the saddle point energy range of the system. The range of relevant kinetic events can be found within mk B T of the lowest saddle point energy, for an m accordingly to the case. Figure 5.4 shows a schematic diagram of these energy window. [START_REF] Xu | [END_REF].

Proposing an m equal to 20 means that an event with higher energy is e 20 as likely to occur as the lower energy event, if the two have similar prefactor (see Equation (5.13)). From this an important assumption for this work will take place. Given that the exponential grows faster than any other function, along this work it will be set a constant prefactor for all the cases under study even if they can vary many order of magnitude.

Recycling saddle points and super basins

Up to the moment, we have set all the saddle points that are going to be included in our rate table with a certain conĄdence criteria. It is the job of the Monte Carlo method to advance the system to its new state. At this new state, a new rate table has to be calculated taking in account all these new probable events. One could hope that some of these new events are similar to the ones of the last state. This would save us an important amount of computational time at the time of building the new event table. Mainly, this is aimed to processes that are localized to a speciĄc area and the rest of the system remains largely unchanged. Therefore, a distance parameter dR is introduced for deciding when to recycle or not. Atoms that move by more of this parameter in the chosen process are set in the Ąnal conĄguration, whereas, atoms that move less than this parameter are placed in the conĄguration of the saddle points of the recycled process. [START_REF] Xu | Adaptive kinetic monte carlo for Ąrst principles accelerated dynamics[END_REF] carried out a test for Al atoms system with diferent sizes. Figure 5.5 shows their results on the saddle points recycling. Figure 5.5 Ű Force evaluations vs system size for three diferent Al system sizes [START_REF] Xu | [END_REF].

In order to end our discussion about saddle points, a discussion about the issue of states connected by low energy saddle points will be undertaken. This topic is important since the KMC eiciency is directly impacted by the choice of saddle points. For instance, if the rate table has mainly rates coming from low energy (compared with minimum) saddle points, the system evolution through time will be poor. The KMC will spend must of its computation time evolving the system in states connected closely in time which is the contrary of the main purpose of investigating systems in the largest possible time periods. Figure 5.6 shows a couple of states connected by a low energy saddle point compared to those that will lead the system to leave the superbasin. It can be proposed that if the system evolves quickly from A to B, an equilibrium will be soon reached, thus, these two states can be considered as a single superbasin. When the system is found in a new state, the rate of each event is compared to the total rate to leave the state, then a ratio between both is calculate and set as a parameter to Ąnd superbasins. This way, if a state has a dominant process which leads to another state that is under this later parameter, then those are grouped in a superbasin. Consecutively, a new rate table is built with all the processes that lead out the superbasin. Finally, if two low barrier states are visited a posteriori, then a combined rate table is created to go out of the superbasin.

Monte Carlo

Monte Carlo (MC) methods are related with the application of random sampling to problems of applied mathematics.Table 5.1 shows a general possible classiĄcation of diferent numerical methods [START_REF] Lambert | Monte carlo methods and simulation techniques[END_REF]. As we can notice, the Monte Carlo methods are stochastic solutions to exactly formulated problems. [START_REF] Amar | The monte carlo method in science and engineering[END_REF] [START_REF] Gardiner | Stochastic methods[END_REF] have been used for many applications, such as, materials science, nuclear power, economy, etc. Its use started around 1940 and 1950 with the apparition of computers [START_REF] Metropolis | The beginnig of the monte carlo method[END_REF]. The name is linked with the statical properties of random events found at gambling in the casinos of Monte Carlo.

MC has a mathematical basis that will be slightly shown in this work. The MC algorithm starts with the speciĄcation of a sampling procedure to estimate characteristic of a stochastic process. Normally, every MC computation is dedicated to estimate the value of an integral 5.16) here f (x) is a known function of x and p(x) is a known probability density function deĄned over the space M. Now a sampling procedure is proposed in order to get an unbiased estimate θ ˆof θ. Equation ( 5.16) can be interpreted using the deĄnitions of the expectation operator E as θ = E[ f (x)] (5.17)

θ = ∫ M f (x)p(x)dx (
where the expectation E is with respect to p(x). The mean θ in Equation ( 5.17) is estimated by the sampling mean θ N , deĄned by

θ N N -1 N ∑ j f (x j ) (5.18)
Using the notation x j [x j ; j = 1, 2, ..., N] for a random sample of size N drawn from the probability density function p(x). As we can notice the subscript denotes that the estimate is based on a Ąnite sample of size N. Finally, we have that the random variable θ N is one stochastic solution to the exactly formulated integral Equation (5.16). The mathematical basis of MC is supported by two main theorems [START_REF] Handschin | Monte carlo techniques for filtering and prediction of nonlinear stochastic processes[END_REF]):

1.-The strong law of large numbers Let x j denote a random sample of size N drawn from the probability density function p(x) to compute the random variable θ N using Equation (5.18). If the independent variates ( f (x j ); j = 1, 2, ...., N) have a common distribution and if the integral of Equation ( 5.17) exists in the ordinary sense, then θ N converges with probability one to θ as N tends to inĄnity.

2.-The central limit theorem

On the premises that the second order moment θ (2) E[ f 2 (x)], deĄned by 5.19) exits and the sample size N is large, the probability that the event E, deĄned by E: θδ θ N θ + δ, occurs is asymptotically independent of the exact nature of f (x) or p(x). Indeed, the probability depends only on N and the variance of f (x), deĄned by 5.20) We can subtract from the theorems that the probability is highly dependent of the sampling size N. The more sampling we make, the better. Moreover, theorem 2 postulates that the exact nature of f(x) or p(x) is not important.

E[ f 2 (x)] ∫ X f 2 (x)p(x)dx (
var( f (x)) E[( f (x) -θ) 2 ] = ∫ X ( f (x) -θ) 2 p(x)dx (
The next mathematical step that will help us to understand the Monte Carlo type used in this work (Kinetic Monte Carlo) is the concept of Markov Chains. First, we denote as uncorrelated chains of a particular sequence of a system which transits from its initial state x i to its Ąnal state x f whose probability of occurrence is statically uncorrelated

P(x i , x i+1 , ....., x f -1 , x f ) = P 1 (x i )P 1 (x i+1 )....P 1 (x f ) (5.21)
where P 1 (x) is independent probability of occurrence for the state x. The Markov chain is a correlated sequence of random events whose probability densities at a time interval depend upon the previous states of the system. It can be deĄned using the transition probability W(x a → x b ) which successfully gives state x b from x a , thus, the probability of having a sequence of states x is

P(x i , x i+1 , ....., x f -1 , x f ) = P 1 (x i )W(x i → x x+1 )W(x x+1 → x x+2 )...W(x f -1 → x f ) (5.22)
From here we can set the time evolution of the transition probability P(x, t) of the system in state x given by the Markovian master equation [START_REF] Levi | Theory and simulation of crystal growth[END_REF]

(K. & W, 2010) ∂P(x f , t) ∂t = - ∑ x f W(x f → x i )P(x f , t) + ∑ x i W(x i → x f )P(x i , t) (5.23)
P(x, t) is the probability that a system is in state x at time t and W(x i → x f ) is the transition probability per unit time that the system will transit from state x i to state x f . The Ąrst term on the right hand side of the Markovian Equation represents all processes where one moves away from the considered state x i at the moment in time t and hence its probability is decreased, while the second term contains all reverse processes which hence lead to an increase of the probability of Ąnding x i . Here the controlling factor is the transition probability. In many cases the system is in a steady state. For example, it occurs when the system is at thermal equilibrium. in this case the occupation probability is Ąxed according to the Boltzmann distribution

P i ∝ e -E i /k B T (5.24)
where E i is the energy of state i, T is the temperature, and k b is the BoltzmannŠs constant.

The steady state occurs when the time derivative of the Markovian Equation is zero, that is ∂P(x f ,t) ∂t = 0. This implies that the sum of all the transitions into a particular state x i equals the sum of all the transitions out of a particular state x f . Therefore, the steady state can be stated as ∑

x f W(x f → x i )P(x f , t) = ∑ x i W(x i → x f )P(x i , t) (5.25)
To be consistent with the Boltzmann distribution in equilibrium, the detailed-balance criterion is usually imposed on the MC transition probabilities. It can be written as

W(x i → x f )P(x i , t = ∞) = W(x f → x i )P(x f , t = ∞)
(5.26) Equation ( 5.26) implies that the ratio of the transition probabilities for a move x f → x i and the inverse move x i → x f depends only on the energy change. It is worth to highlight that the detailed balance is a necessary but not suicient condition for thermodynamic equilibrium. This type of MC that evolves form state to state is usually called Kinetic Monte Carlo or KMC. It was developed around 1970s for atomistic systems of crystal growth using probabilistic rules to govern deposition, difusion and other phenomena. It was Ąrst adopted by [START_REF] Abraham | Computer Simulation of Vapor Deposition on Two Dimensional Lattices[END_REF] and [START_REF] Gilmer | Simulation of Crystal Growth with Surface Difusion[END_REF]. The KMC is usually applied to simulate the evolution of a physical system through numerical sampling of (Markovian) stochastic processes. In contrast the traditional MC method is applied to systems close to equilibrium, while KMC allow us to advance our system under study along time, giving us this way the opportunity to study phenomena out of equilibrium. The basic principle of KMC is to carry out elementary jumps on a grid, in order to represent real atomic jumps on a crystal. The code Ąrst make a random selection to a process with the use of transition rates for all possible processes. Then, the selected transition is performed and the corresponding change of the state vector. Finally, the time counter and transitions probabilities are updated.

Adaptive Kinetic Monte Carlo

The Adaptive Kinetic Monte Carlo (aKMC) method is a relatively new type of MC. It surpasses its predecessors by calculating its rate or event table every time the system is advanced to a new state [START_REF] Henkelman | Long time scale kinetic monte carlo simulations without lattice approximation and predeĄned event table[END_REF]. In contrast, the KMC has to have an accurate event table a priori. Therefore, only events that were included in the beginning of the simulation can be used to advance the system. Moreover, the KMC uses an on-lattice approximation that prohibits the simulation of materials that have gone through large deformations such as nuclear materials inside the reactor.

Having the ability to calculate on-the-Ćy the event table allows the aKMC to search into unforeseen states that would arise from complicated mechanisms. Thus, it beneĄts of longer periods of time while maintaining an atomic representation of the system in evolution under the force coming from an interatomic potential.

The aKMC has already been tested for nuclear materials. The Nuclear National Laboratory has used the open source DL/ / AK MC code developed within the ScientiĄc Computing group at the Science and Technology Facilities Council (STFC) Daresbury laboratory (Gunn, Allan, & Purton, 2014, 33) to study oxygen difusion in UO 2 . The results showed that aKMC predicts anion difusion coeicient close to the ones calculated with MD [START_REF] Platts | Modelling the evolution of radiation induced defects in nuclear fuel using adactive kinetic monte carlo simulation[END_REF]. Moreover, aKMC has been used to study the evolution of radiation damage near grain boundaries in fcc Ni and a Ni-Cr alloy [START_REF] Tooq | Modelling radiation damage at grain boundaries in fcc nickel and ni based alloy using long time scale dynamics techniques[END_REF].

The newest versions of DL/ / AK MC include both parallelisation fo the transition states searches and saddle point re-cycling. Up to the moment the code has succeeded in demonstrating a good scaling up to many hundreds of CPU cores. Figure 5.7 shows the schematic of the parallelisation and simulation process implemented in DL/ / AK MC.

Long term recombination

As we stated in Section 4.2, during displacement cascades numerous FPs with small distances between the interstitial and the vacancy are created. In the same section we analysed them using MD. In the present section, we have used the DL/ / AK MC code to analyse again the possible FPs that displacements cascades produce. This study will allow us to search for longer periods of time compared to MD. (Gunn, Allan, & Purton, 2014, 33) Almost the same process is followed as in Section 4.2. To calculate the recovery time of the defect, the supercell is Ąrst relaxed for 5 ps in the NPT ensemble in order to Ąnd the temperature dependent lattice parameter using MD. A C-FP is created after this relaxation. This defected structure is used as input for the DL/ / AK MC code. If annealing of the C-FP occurs, the time of recombination is recorded. The FPs are classiĄed according to the distance between the interstitial and the vacancy as second-, third-and fourth-nearest neighbours as shown in Figure 4.1. First neighbours are excluded since they anneal instantly. Due to the Ćuorite symmetry, there are: 8 second, 24 third, and 48 fourth possible nearest neighbours for the C-FPs. For simplicity, we limit the number of third-and fourth-nearest neighbours studied to 8 (all shown in Figure 4.1).

Due to the high computational power needed to carry out MC simulations, we could only analyze 5 MOX, 2 UO 2 and 1 PuO 2 structures. MOX structures have Pu randomly distributed. If for a speciĄc structure there is no recombination, the case is simply excluded from the average. Therefore, the resulting lifetime corresponds to the lowest limit of the recombination time. Nevertheless, it gives an indication of whether recombination may occur or not. Table 5.2 compiles all the results of the averaged lifetimes with their standard deviation obtained with the aKMC. Only Potashnikov potential could be implemented. Eforts to implement Cooper potentials are being made. Only one temperature 1600 K was assessed again due to the high computational power required. We chose this temperature because it is the most important in Section 4.2 since at higher temperatures, the process of recombination is more probable, indicating a thermally activated process.

The supercells consisted of 3 × 3 × 3 Ćuorite unitary cells. The dimer method (see Section 5.4 is chosen to carry out the search of saddle points. A minimum activation energy of 1 eV for a transition to be allowed was chosen in order to limit the search of transitions only for cations . The energy parameter m is chosen at 160 (see Figure 5.4). This gives a range of energy of 1 to 22 eV. This range was chosen in this way to focus only on cation movements. The Ąrst important point that comes to light from Table 5.2 is that aKMC allows to search for much longer periods of time compared with MD. In Section 4.2, we limited our research to 100 ps. Here we reached up to almost milliseconds. Secondly, we could limited our research on only rare events. This is cation movements. As a consequence, we found that 2 nd are able to anneal at this speciĄc conditions for times around 10 µs. The process that gives annealing for the 2 nd neighbour is quite complex. It is not a direct recombination, meaning that the interstitial does not go directly to the vacancy. Instead, two surrounding cations to the vacancy plus the interstitial needed to move in order to have annealing of the 2 nd neighbour. This is not surprising due to the energy barrier imposed by the oxygen laying between the interstitial and the vacancy, as already seen in previous studies in Ćuorite-like structure using rigid ion interatomic potentials [START_REF] Pannier | Molecular dynamics study of frenkel pair recombinations in Ćuorite type compounds[END_REF].

The recombination time for third neighbours is quite long compared to our study in Section 4.2. This is probably due to our energy window. Since we limited the minimum activation energy of 1 eV for a transition to be allowed. This could block recombination paths that give short recombination times found in Section 4.2. Finally, only one recombination event for the fourth neighbours is found. This could mean that recombination at this interstitial distance from the vacancy involves processes with energy below 1 eV. The later is in concordance with fourth-neighbour activation energies for UO 2 that are around 0.6 eV (Van Brutzel et al., 2008, 2).

Conclusions

As a conclusion, we proved that aKMC is a powerful tool to overcome the time of MD simulations. It allows to search for much longer periods of time. Important results due to this longer period of times are achieved, for example, we found that 2 nd are able to anneal at this speciĄc conditions for times around 10 µs. In contrast, we could not Ąnd any annealing for the second neighbour in MD. Moreover, results show that recombination is possible to happen at times around milliseconds. However, the number of recombination events is lower that those that take place before 100 ps. As a consequence, we propose that recombination of C-FPs falls into three regimes: Ąrst, a spontaneous regime where recombination events occur in less than 1 ps and second, a thermally activated regime and Ąnally, rare-event recombination events taking place in longer periods of time (up to milliseconds in this study).

Both techniques working together could give broader insights than each one working separately could give. More eforts have to be made to increase the limits of this research. For example. the implementation of Cooper potential or analyze extended defects with aKMC such as, dislocations or cavities.

Conclusions

Nuclear fuel inside the reactor sufers micro-structural changes. One example is the restructuring of the solid matrix mainly due to high temperatures and large temperature gradients. Voids migrate towards higher temperatures and gather in the centre. Conversely, a movement of the solid matrix takes place towards the outer periphery. There is a diference between MOX and UO 2 pellets. In MOX, there are higher rates of Ąssion gas release, due to higher linear power levels at the end of the irradiation and higher centre-line temperature (caused by slightly lower thermal conductivity). Furthermore, the heterogeneity of the Ąssile material would play a role as well. The Pu rich spots are submitted to very high burn-up. This causes restructuring, which leads to further division of grains, precipitation of gas bubbles that go to intra-granular division which make them easily freed when a power transient occurs. In general, MOX fuel pellets are more complicated to work with compared to UO 2 . Their behaviour inside the reactor has to be carefully addressed, since more processes are involved at the time of the burn-up. Moreover, at the time of producing them, the alpha decay make them a diicult task and decreases the number of capable labs that are able to sensitize them. Therefore, it is important to understand well what are the micro-structural changes which happen inside the fuel under irradiation and surely the cause as well. This brings the need of new tools to assess the problem. For example, this PhD work research focuses on answering this question by using a computer simulation approach at the atomic scale.

One of the purpose of this doctoral research is to contribute to the knowledge concerning MOX and its thermomechanical properties. This will help to build updated catalogues of MOX properties for the new reactor prototypes. Moreover, our work is part of the European project INSPYRE which focuses on the investigation of fast reactor and MOX fuel. Most of this PhD work lies under the work package 3, which focuses on improving, developing and applying experimental and modelling methods to reach a better understanding of the underlying mechanisms governing the evolution of mechanical properties.

Our Ąrst computational tool was MD. The reliability of the results obtained with this tool depends on the choice of the set of potentials. For this reason, we assessed empirical potentials for the (U 1-y ,Pu y )O 2 solid solution. We found in the literature Ąve interatomic potentials coined by the name for their Ąrst 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 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 included ranges of temperatures and compositions as well as properties that have not been studied by the authors.

Our Ąrst study was to compute the evolution of the lattice parameter. This helped us to Ąnd that Tiwary potential is unstable with MD simulations. For this potential the Ćuorite structure collapses after few steps of the simulation. Thus, we removed this potential for the following tests. Furthermore, this test allowed us to eliminate, as well, Yamada potential, since it shows large discrepancies on the Ćuorite structure at high temperatures (>2100 K). The three remaining potentials seem to behave in a similar way for the thermodynamic properties (thermal expansion and heat capacity).

For the mechanical properties, results show that the elastic stifness constants and brittleto-ductile transition are best reproduced with the Cooper potential. The Potashnikov potential gives fairly good agreement with experimental data while Arima potential overestimates largely the elastic stifness constants. However, anisotropy is presented even at low temperature for the Cooper and Arima potentials whereas it appears only at high temperatures for Potashnikov potential. Therefore, the main diference between potentials concerning the mechanical properties is the anisotropy.

Concerning the Pu content, only Cooper potential shows a small sensitivity on the thermodynamics properties. This sensitivity is linked to a greater accuracy to predict melting temperature according the Pu content.. However, all the mechanical properties show little dependence on the plutonium content, conĄrming the assumption that mechanical properties of UO 2 can be, to a large extent, applied to MOX.

Overall, with the structural, thermodynamics, and mechanical properties assessment, the three remaining potentials show advantages and disadvantages. 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. It also presents the expected brittle behaviour during crack propagation but renders isotropic elastic behaviour. 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. As a consequence, we have decided to keep Cooper and Potashnikov potentials for our radiation assessment.

The radiation assessment is performed is four steps. First, the calculation of the formation energy of point defects that assesses the potentials on their capability to reproduce correctly or not the creation of point defects. Second, the analysis of the Frenkel pair recombination time, which is used as set parameter for the last step. Third, simulations of displacement cascades that render the type and number of defect created by one single radiation event. And last, simulations of Frenkel pair accumulation, which is designed to study the dose efect.

From defect formation energies study, both potentials give approximately the same values and corroborate with existing data for urania.

Concerning the FP recombination. There is a clear dependence on the temperature. At high temperatures (1600 K), the process of recombination is more probable, indicating a thermal activated process. Recombination processes seem to be more efective for Potashnikov potential than for Cooper based on both; the recombination times and temperatures at which they occur.

For displacement cascades the same general behaviour is observed for both potentials. Within the Ąrst picosecond, the PKA sufers several collisions releasing its kinetic energy in the form of a large amount of atomic displacements in all directions. These atomic displacements induce a temporary local disorder of the crystal in the core of the cascade. Numerous FPs are created during this stage. Subsequently, the system evolves releasing energy in the form of thermal vibrations, causing a drastic increase of the local temperature (thermal spike). Consequently, numerous recombination events take place, which mainly leads to the reconstruction of the Ćuorite lattice leaving only a few remaining point defects. However, diferences between potentials start to appear at 75 keV. For the Cooper potential, the recovery stage does not lead to full reconstruction of the Ćuorite structure. Instead, a large disordered region with no crystallographic atomic structure of several manometers is found at the cascade core. This behaviour is probably linked to the slow kinetics of Frenkel pair recombination. In contrast, Potashnikov potential presents a quasi-total reconstruction of the Ćuorite structure.

The Frenkel Pair accumulation study was carried out to investigate the dose efect. Potashnikov potential presents the following process as the dose increases: accumulation of point defects that cluster and form Frank loops which in turn transform into perfect loop and dislocation lines. This behaviour is diferent for Cooper potential. At 1600 K, dislocations are formed but only few Frank loops are observed suggesting either a direct mechanism for the creation of perfect loops or that they are highly unstable.

Temperature plays also an important role in the evolution of the extended defects. With Cooper potential no dislocation formation is observed at 300 K whereas some appear at 1600 K. This diference could come from the fact that defect recombination at 300 K for Cooper potential is highly improbable as shown in the FP recombination study. For Potashnikov potential some quantitative diferences occur; the full process is shifted to lower doses as the temperature increases. About the Pu content, the dislocation density decreases as the Pu content increases and peaks of density shift to higher dose with increasing Pu content. Therefore, the creation of dislocations is easier in urania than in MOX or plutonia. As a hypothesis, we suggest that this behaviour is caused by the fact that for Cooper potential the recombination processes are facilitated with the increase of the Pu content. Thus, more defects in urania are able to create dislocations, such as, Frank dislocations and initiate the process that ends in perfect dislocations. In contrast, defects being annealed faster in plutonia will create less Frank dislocations and therefore less perfect dislocations. This is in agreement with the number of defects, which shows that the number of vacancies and interstitials slightly increases as Pu content increases. Finally, we could not Ąnd any inĆuence of the Pu content on the elastic moduli as a function of the dose and the X-ray powder difraction studies. In general, Cooper potential is more sensitive to changes in the Pu content.

The spatial repartition of defects was analysed after displacement cascades occurs and Frenkel pair accumulation simulations. In general, we could not Ąnd an important role of the Pu content on the morphology of defects. This means that we did not Ąnd any trend of Pu atoms clustering with Pu atoms (same can be applied to U). However, we could Ąnd important clustering in point defects at the end of the cascades. Interstitials cluster less than vacancies. Vacancies can cluster until a point they form nanocavities of 1 nm size.

We continued our study on the efect of dose. The elastic moduli as a function of dose and the X-ray powder difraction studies agreed on one point. The points defects were one of the most important factors that played a role when accumulating dose. They caused a sharp decrease in both; the Bulk modulus and the Zener ratio. In contrast, the evolution of the FWHM is highly correlated to the perfect dislocations.

Finally, we used aKMC, which is a powerful tool to overcome the time of MD. It allowed to search for much longer periods of time compared with MD. Important results due to this longer period of times were achieved, for example, we found that 2 nd neighbours were able to anneal after about 10 µs at 1600 K. In contrast, we could not Ąnd any annealing for the second neighbours by MD.

In general, along the whole PhD work, we systematically carried out studies with the Pu content, temperature and diferent potentials as variables of interest. We gave pros and contras of each potential. We tested our results when available with the scarce experiment values in the literature. Therefore, the trends found through our results are believed to be representative of the material under study; MOX.

There is still an important amount work left to do in order to better understand MOX, starting with exploring or creating new potentials that are able to simulate phenomena out of stoichiometry. Secondly, the local conĄguration around defects has to be analysed in order to evaluate the inĆuence of changing U for Pu and vice versa on the formation energy of defects. In the same way, an study of the impact of the diference between cation species could bring a insight of the annealing path in the FP recombination study. Thirdly, the spatial repartition of defects could be used as input for the FP accumulation method so each step of the latter will simulate the Ąnal state of a cascade (including nanocavities). Finally, the aKMC method has a huge potential to expand our studies to longer time periods.

1. 1 1 . 2 5 1. 5

 11255 Schematic presentation of the open fuel cycle (MIT, 2003) . . . . . . . . . 3 Schematic presentation of the closed fuel cycle (MIT, 2003) . . . . . . . . 3 1.3 Fuel cycle characteristics of current plants (MIT, 2003) . . . . . . . . . . . 4 1.4 X-ray Mappings of a newly fabricated MOX pellet (KAERI, 2009). The colour level corresponds to wt % of PuO 2 . . . . . . . . . . . . . . . . . . Micrographs of MOX fuel pellet. (KAERI, 2009) . . . . . . . . . . . . . . 5

  2 . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.8 Representation of the MOX Ćuorite structure. . . . . . . . . . . . . . . . . 9 1.9 Coordinance tetrahedron and cubic in the MOX Ćuorite structure. . . . . . 10 1.10 Diagram showing a showing 2D projection onto (1 0 0) plane of a cation Frenkel defect. The red circles represent cations (Pu or U) and the blue ones, anions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.11 Diagram showing 2D projection onto (1 0 0) plane of a Schottky defect. . . . 11 1.12 Schematic lattice network showing defects. a) Interstitial impurity atom, b)

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Evolution of the lattice parameter as a function of Pu content at 300 K. . . . . 3.4 Percentage of deviation from VegardŠs Law for (U 0.75 Pu 0.25 )0 2 . . . . . . . . 3.5 Percentage of deviation from VegardŠs Law for (U 0.30 Pu 0.70 )0 2 . . . . . . . . 3.6 Linear thermal expansion coeicient obtained with Cooper potential for diferent plutonium contents. . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Linear thermal expansion coeicient obtained with Potashnikov potential for diferent plutonium contents. . . . . . . . . . . . . . . . . . . . . . . . . .3.8 Linear thermal expansion coeicient obtained with Arima potential for diferent plutonium contents. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Evolution of the enthalpy increment as a function of temperature calculated with Potashnikov, Cooper, and Arima potentials. These values are compared with FinkŠs recommendation . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Evolution of the SpeciĄc hear C P as a function of the temperature for diferent plutonium contents for Cooper potential. . . . . . . . . . . . . . . . . . . 3.11 Evolution of the SpeciĄc hear C P as a function of the temperature for diferent plutonium contents for Potashnikov potential. . . . . . . . . . . . . . . . . . xiv 3.12 Evolution of the SpeciĄc hear C P as a function of the temperature for diferent plutonium contents for Arima potential. . . . . . . . . . . . . . . . . . . . . 3.13 Comparative study of MooreŠs model, KatoŠs model and VauchyŠs experimental data for 1273 K (a) and 1073 K (b) respectively. Taken from (Chakraborty, 2017). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Energy-volume plots shown with Ąve diferent structures for pure UO 2 . . . . . 3.15 Evolution of monocrystal YoungŠs modulus as a function of the plutonium content calculated with Arima, Potashnikov, and Cooper potentials and are compared with the experimental values in[START_REF] Hirooka | Sound speeds in and mechanical properties of u,puo2x[END_REF] . . . . .3.16 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 (ESNII+, 2015) and (c) is compared with experimental date of (Hutchings, 1987). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.17 Evolution of the ZenerŠs anisotropy factor as a function of temperature for pure UO 2 for the three interatomic potentials. . . . . . . . . . . . . . . . . 3.18 Evolution of the ZenerŠs anisotropy factor as a function of temperature for pure PuO 2 for the three interatomic potentials. . . . . . . . . . . . . . . . . 3.19 Evolution of the ZenerŠs anisotropy factor as a function of the Pu content at 300 K for the three interatomic potentials. . . . . . . . . . . . . . . . . . . 3.20 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 Ćuorite 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 . . . . . . . . . 3.21 UTS as a function of temperatures for the three potentials. y = 0 corresponds to urania and y = 1 corresponds to Pu . . . . . . . . . . . . . . . . . . . . 3.22 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 (Roberts & Wrona, 1971). . . . . . . . . . . . . . . . . . . . 3.23 Representation of the thin strip geometry of the system used for cracking simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.24 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 Ćuorite structure. The inserts in each Ągure depict a close up of the crack front where phase transition occurs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv 3.25 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 Ćuorite structure. . 4.1 Schematic diagram of the ⟨1 0 0⟩ and ⟨1 1 0⟩ planes showing the possible crystallographic position for the cation interstitials in a Ćuorite-like lattice. The positions are classiĄed as 1 st , 2 nd , 3 r d , and 4 th nearest neighbour according to their distance to the vacancy. . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Evolution of the number of (a) cation and (b) anion displaced atoms at the end of displacements cascades as a function of PKA energy for both potentials (C for Cooper and P for Potashnikov) and the three Pu contents. The black doted line corresponds to the theoretical number of displaced atoms estimated with the Kinshin-Pease law. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Evolution of the number of (a) cation and (b) anion FPs present at the end of displacements cascades as a function of PKA energy for both potentials (C for Cooper and P for Potashnikov) and the three Pu contents. The black doted line corresponds to the theoretical number of FPs estimated with the NRT law. 4.4 Snapshot of the defect analysis at the end of the cascade initiated with 75 keV PKA for Cooper potential for 50% Pu content. The red zone marked the regions where structure departs from perfect Ćuorite structure. The insets represent close-up visualisation of the atoms in the main disordered region. 4.5 Snapshot of the defect analysis at the end of the cascade initiated with 75 keV PKA for Potashnikov potentials for 50% Pu content. The red zone marked the regions where structure departs from perfect Ćuorite structure. The insets represent close-up visualisation of the atoms in the main disordered region. 4.6 Distribution of point defects in cluster (interstitial and vacancy) as a function of their size for Cooper and Potashnikov potentials calculated from 10 keV PKA cascades. The error bars correspond to the standard deviation of the diferent values averaged over the diferent microstructures. . . . . . . . . 4.7 Evolution of dislocation densities as a function of dose (expressed in dpc) obtained with (a) Cooper and (b) Potashnikov potentials at 1600 K for 0%, 50%, and 100% Pu. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Evolution of the number of point defects (vacancy and interstitial) as a function of dose (expressed in dpc) obtained with (a) Cooper and (b) Potashnikov potentials at 1600 K for 0%, 50%, and 100% Pu. . . . . . . . . . . . . . . . 4.9 Snapshots of the evolution of extended defects as a function of dose (expressed in dpc) for Cooper potential at a) 300 K and b) 1600 K for 50% Pu. The red zone marked the regions where structure departs from perfect Ćuorite structure. The cyan lines are 1 /3⟨1

5. 1

 1 One dimensional potential as a function of a reaction coordinate with the most important points highlighted. . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Bias potential constructed in such a way that it replaces the real potential by a constant equal to an energy called "boost energy" . . . . . . . . . . . . . . 5.3 Dimer method in a two-dimensional potential surface. The system is initially in state A. (a) Ten dimer searches are started from random positions around the minimum. They converge on four distinct saddle points (two of the searches practically overlap). (b)

1

 1 Point defects . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.6.2 Extended defects . . . . . . . . . . . . . . . . . . . . . . . . . . 10

Figure 1 . 1 Ű

 11 Figure 1.1 Ű Schematic presentation of the open fuel cycle (MIT, 2003)

Figure 1 . 2 Ű

 12 Figure 1.2 Ű Schematic presentation of the closed fuel cycle (MIT, 2003)

Figure 1 .

 1 Figure 1.3 Ű Fuel cycle characteristics of current plants (MIT, 2003)

Figure 1 . 5 Figure 1 . 5 Ű

 1515 Figure 1.4 Ű X-ray Mappings of a newly fabricated MOX pellet (KAERI, 2009). The colour level corresponds to wt % of PuO 2An irradiation test is undertaken with an average burnup of 50 MWd/kgHM with an irradiation time of 1020 EFPD. The peak in the temperature in the fuel centre was estimated around 1673 K. The results of the irradiation show many Pu rich spots and metallic Ąssion product precipitates. The Pu rich spots were found mainly in middle and peripheral regions, in contrast to the centre where the density was lower. Pu in the spots difuses out to the UO 2 matrix in the central region[START_REF] Kaeri | Ceramography analysis of mox fuel after an irradiation test[END_REF]. Pu spots are surrounded by a ring of small pores called "Halo" as shown in Figure1.5

Figure 1 .

 1 Figure 1.7 Ű Microprobe mappings of MIMAS MOX. (CEA, 2009). The colour level corresponds to wt % of PuO 2

Figure 1 .

 1 Figure 1.8 Ű Representation of the MOX Ćuorite structure.

Figure 1 .

 1 Figure 1.9 Ű Coordinance tetrahedron and cubic in the MOX Ćuorite structure.

Figure 1 .

 1 Figure 1.10 Ű Diagram showing a showing 2D projection onto (1 0 0) plane of a cation Frenkel defect. The red circles represent cations (Pu or U) and the blue ones, anions.

Figure 1 .

 1 Figure 1.11 Ű Diagram showing 2D projection onto (1 0 0) plane of a Schottky defect.

Figure 1 .

 1 Figure 1.12 Ű Schematic lattice network showing defects. a) Interstitial impurity atom, b) Edge dislocation, c) Self interstitial atom, d) Vacancy, e) Precipitate of impurity atoms, f) Vacancy type dislocation loop, g) Interstitial type dislocation loop, h) Substitutional impurity atom[START_REF] Foll | Defects in crystals[END_REF].

Figure 1 .

 1 Figure 1.14 Ű Schematic representation of Shockley VS Frank dislocations (Radwan, 2018).

Figure 2 . 1 Ű

 21 Figure 2.1 Ű Periodic boundary conditions

Figure 2 .

 2 Figure 2.3 Ű Diagram showing how OVITO calculates the occupancy (Stukowski, 2010).

  Figure 2.5 shows an example of an atomic volume histogram as a function of the count of atoms from the Voronoi cell analysis. In this Figure, one can observe that the Voronoi volume associated to interstitials is smaller than the volume related to atoms in a regular lattice position because it has extra neighbours closer to him. Conversely, a vacancy should have no Voronoi volume associated.

Figure 2 .

 2 Figure 2.4 Ű Division of the space by Voronoi cells

  Figure 2.6 shows an example of how the Burgers vector is created.

Figure 2 .

 2 Figure 2.6 Ű Example of how the Burgers vector is created (Stukowski, 2010).

Figure 2 .

 2 Figure 2.7 Ű Construction of dislocation lines. (Stukowski, 2010).

Figure 2 .

 2 Figure 2.8 Ű Graphic representation of the efect of strain in the crystalline structure and its efect on the difraction peaks. a) No strain b) uniform strain c) non-uniform strain is applied.

Figure 3

 3 Figure 3.3 Ű Evolution of the lattice parameter as a function of Pu content at 300 K.

Figure 3 . 4 Ű

 34 Figure 3.4 Ű Percentage of deviation from VegardŠs Law for (U 0.75 Pu 0.25 )0 2 .

Figure 3

 3 Figure 3.6 Ű Linear thermal expansion coeicient obtained with Cooper potential for diferent plutonium contents.

Figure 3

 3 Figure 3.7, displays the results of the LTEC for Potashnikov potential. The values from 300 to 2000 K are higher than the experimental data and the recommendation. The λ-peaks are present but less pronounced than for the Cooper potential. Furthermore, no clear diference appears on the increment with plutonium content.

Figure 3

 3 Figure 3.7 Ű Linear thermal expansion coeicient obtained with Potashnikov potential for diferent plutonium contents.

Figure 3

 3 Figure 3.8 Ű Linear thermal expansion coeicient obtained with Arima potential for diferent plutonium contents.

Figure 3

 3 Figure 3.9, shows the evolution as a function of the temperature up to the melting point of the enthalpy increment calculated with Potashnikov, Cooper, and Arima potentials. The enthalpy values up to 2000 K are almost identical and Ąt perfectly the FinkŠs recommendation. Discrepancy appears at 2400 K, where Arima potential still gives linear feature, while FinkŠs recommendation and both Potashnikov and Cooper potentials deviates from linearity. For these two last potentials discrepancy with FinkŠs recommendation appears around 2600 K, which is the temperature of the λ-peak.

Figure 3

 3 Figure 3.9 Ű Evolution of the enthalpy increment as a function of temperature calculated with Potashnikov, Cooper, and Arima potentials. These values are compared with FinkŠs recommendation

Figure 3 .

 3 Figure 3.11 Ű Evolution of the SpeciĄc hear C P as a function of the temperature for diferent plutonium contents for Potashnikov potential.

Figure 3 .

 3 Figure 3.12 Ű Evolution of the SpeciĄc hear C P as a function of the temperature for diferent plutonium contents for Arima potential.

  A comparative study of MooreŠs model, KatoŠs model and VauchyŠs experimental data of Oxygen self difusion is presented in Figure 3.13.

Figure 3 .

 3 Figure 3.13 Ű Comparative study of MooreŠs model, KatoŠs model and VauchyŠs experimental data for 1273 K (a) and 1073 K (b) respectively. Taken from (Chakraborty, 2017).

  Figure 3.14 Ű Energy-volume plots shown with Ąve diferent structures for pure UO 2 .

Figure 3 .

 3 Figure 3.15 Ű Evolution of monocrystal YoungŠs modulus as a function of the plutonium content calculated with Arima, Potashnikov, and Cooper potentials and are compared with the experimental values in[START_REF] Hirooka | Sound speeds in and mechanical properties of u,puo2x[END_REF]

Figure 3 .

 3 Figure 3.16(a), Figure 3.16(b), and Figure 3.16(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 in agreement with the recommendations in ESNII[START_REF] Esnii+ | State of the art with a literature review of mox properties[END_REF].

Figure 3 .

 3 Figure 3.16 Ű 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 (ESNII+, 2015) and (c) is compared with experimental date of[START_REF] Hutchings | High temperature studies of uo2 and tho2 using neutron scattering techniques[END_REF].

Figure 3 .

 3 Figure 3.18 Ű Evolution of the ZenerŠs anisotropy factor as a function of temperature for pure PuO 2 for the three interatomic potentials.

Figure 3 .

 3 Figure 3.19 Ű Evolution of the ZenerŠs anisotropy factor as a function of the Pu content at 300 K for the three interatomic potentials.

  Figure 3.20 Ű 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 Ćuorite 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 3 .

 3 Figure 3.21 Ű UTS as a function of temperatures for the three potentials. y = 0 corresponds to urania and y = 1 corresponds to Pu

  Figure 3.24 Ű 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 Ćuorite structure. The inserts in each Ągure depict a close up of the crack front where phase transition occurs.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

Figure 4 . 1 Ű

 41 Figure 4.1 Ű Schematic diagram of the ⟨1 0 0⟩ and ⟨1 1 0⟩ planes showing the possible crystallographic position for the cation interstitials in a Ćuorite-like lattice. The positions are classiĄed as 1 st , 2 nd , 3 r d , and 4 th nearest neighbour according to their distance to the vacancy.

Figure 4 .

 4 Figure 4.2 shows the number of displaced cations and anions as a function of PKA energy.The error bars correspond to the standard deviation from the average over the diferent cascades. This number is compared to the theoretical number of displaced atoms estimated by the Kinshin-Pease linear law[START_REF] Kinchin | The displacement of atoms in solids by radiation[END_REF]):

Figure 4 . 2 Ű

 42 Figure 4.2 Ű Evolution of the number of (a) cation and (b) anion displaced atoms at the end of displacements cascades as a function of PKA energy for both potentials (C for Cooper and P for Potashnikov) and the three Pu contents. The black doted line corresponds to the theoretical number of displaced atoms estimated with the Kinshin-Pease law.

Figure 4 . 3 Ű

 43 Figure 4.3 Ű Evolution of the number of (a) cation and (b) anion FPs present at the end of displacements cascades as a function of PKA energy for both potentials (C for Cooper and P for Potashnikov) and the three Pu contents. The black doted line corresponds to the theoretical number of FPs estimated with the NRT law.

Figure 4

 4 Figure 4.4 Ű Snapshot of the defect analysis at the end of the cascade initiated with 75 keV PKA for Cooper potential for 50% Pu content. The red zone marked the regions where structure departs from perfect Ćuorite structure. The insets represent close-up visualisation of the atoms in the main disordered region.

Figure 4 . 5 Ű

 45 Figure 4.5 Ű Snapshot of the defect analysis at the end of the cascade initiated with 75 keV PKA for Potashnikov potentials for 50% Pu content. The red zone marked the regions where structure departs from perfect Ćuorite structure. The insets represent close-up visualisation of the atoms in the main disordered region.

Figure 4

 4 Figure 4.6 Ű Distribution of point defects in cluster (interstitial and vacancy) as a function of their size for Cooper and Potashnikov potentials calculated from 10 keV PKA cascades. The error bars correspond to the standard deviation of the diferent values averaged over the diferent microstructures.

Figure 4 .

 4 Figure 4.9 and Figure 4.10 display snapshots of the damage as a function of dose expressed in dpc obtained with Cooper and Potashnikov potentials respectively for MOX with 50%Pu at 1600 K and 300 K . Only, dislocations and main zones where structure departs from perfect Ćuorite structure are reported in these Ągures for clarity. Overall, both potentials show the same general trend for the evolution of the microstructure with the dose at 1600 K. However, diferences appear in the dose at which it occurs. The evolution of the system follows 3 main stages. First, point defect concentration increases without formation of a dislocation (see snapshot 0.2 dpc with cooper potential). Subsequently, at about 0.5 dpc for Cooper and 0.1 dpc for Potashnikov Frank loops with BurgerŠs vector 1 /3⟨1 1 1⟩ nucleate from

Figure 4

 4 Figure 4.7 Ű Evolution of dislocation densities as a function of dose (expressed in dpc) obtained with (a) Cooper and (b) Potashnikov potentials at 1600 K for 0%, 50%, and 100% Pu.

Figure 4

 4 Figure 4.8 Ű Evolution of the number of point defects (vacancy and interstitial) as a function of dose (expressed in dpc) obtained with (a) Cooper and (b) Potashnikov potentials at 1600 K for 0%, 50%, and 100% Pu.

  Figure 4.9 Ű Snapshots of the evolution of extended defects as a function of dose (expressed in dpc) for Cooper potential at a) 300 K and b) 1600 K for 50% Pu. The red zone marked the regions where structure departs from perfect Ćuorite structure. The cyan lines are 1 /3⟨1 1 1⟩ Frank loops, the green lines are 1 /6⟨1 1 2⟩ Shockley partial dislocations, and the dark blue lines are 1 /2⟨1 1 0⟩ perfect dislocations.

Figure 4 .

 4 Figure 4.11 Ű Evolution of dislocation densities as a function of dose (expressed in dpc) obtained with Potashnikov potential at 300 and 1600 K for 50% Pu. (circle) 1 /3⟨1 1 1⟩ Frank loops, (triangle) 1 /6⟨1 1 2⟩ Shockley partial dislocations, (star) 1 /2⟨1 1 0⟩ perfect dislocations.

Figure 4 .

 4 Figure 4.12 Ű Evolution of the bulk modulus as a function of irradiated dose expressed in displacement per cation (dpc) for (a) Cooper and (b) Potashnikov potentials at 1600 K.

Figure 4 .

 4 Figure 4.13 Ű Evolution of the Zener ratio as a function of irradiated dose expressed in displacement per cation (dpc) for (a) Cooper and (b) Potashnikov potentials at 1600 K.

  Figure 4.15 and Figure 4.16 show simultaneously the FWHM and dislocation densities as a function of dose calculated at 1600 K for Cooper and Potashnikov potentials, respectively.

Figure 4 .

 4 Figure 4.14 Ű Evolution of the XRD pattern for UO 2 for Cooper potential as a function of dose

Figure 4 .

 4 Figure 4.16 Ű Evolution of the FWHM and dislocation densities for Potashnikov potential.

  .6)We want to know what is the relation or rate of particles that will cross the saddle point in comparison with the particles in the state a, thus

Figure 5 .

 5 Figure 5.3 Ű Dimer method in a two-dimensional potential surface. The system is initially in state A. (a) Ten dimer searches are started from random positions around the minimum. They converge on four distinct saddle points (two of the searches practically overlap). (b)The system is then made to go down the minimum energy path (gray lines) on either side of the saddle points which are indicated with *. Here, all four saddle points have a minimum energy path starting at the initial state minimum A, but this does not have to be the case. The rate of each process is then calculated using Equation (5.13). A process is chosen at random using the kinetic Monte Carlo algorithm. In this case (see section Section 5.6), process 1 gets chosen. The system is moved to the Ąnal state of this process, to minimum B. (c) Dimer searches are run from the new minimum, again four distinct saddle points are found. (d) Minimum energy paths are traced out, and the process repeated[START_REF] Henkelman | Long time scale kinetic monte carlo simulations without lattice approximation and predeĄned event table[END_REF]

Figure 5 .

 5 Figure 5.4 Ű Reaction mechanisms with barriers within mk B T of the lowest saddle point energy are considered relevant[START_REF] Xu | [END_REF].

Figure 5 .

 5 Figure 5.6 Ű Superbasin composed by two states (a and b) and the two barriers (∆E 1 , ∆E 2 ) (Xu & Henkelman, 2008)

Figure 5 .

 5 Figure 5.7 Ű Schematic of the parallelisation and simulation process implemented in DL_AK MC(Gunn, Allan, & Purton, 2014, 33)

Figure A. 1 Ű

 1 Figure A.1 Ű Ternary phase diagram of the O-Pu-U system at 1500 K[START_REF] Guéneau | Thermodynamic modelling of advanced oxide and carbide nuclear fuels: description of the uŰpuŰoŰc systems[END_REF]

  Brutzel, A. Chartier, and Y. Le Bouar, "Assessment of empirical potential for MOX nuclear fuels and thermomechanical properties", Journal of Nuclear Materials, 495 (2017), 67 -77.

• H. Balboa, L. Van Brutzel, A. Chartier, and Y. Le Bouar, "Damage characterization of (U,Pu)O 2 under irradiation by molecular dynamics simulations", Journal of Nuclear Materials, accepted in press, doi:https://doi.org/10.1016/j.jnucmat.2018.07.056. • L. Van Brutzel, H. Balboa, and A. Chartier, "Mechanical properties of (U,Pu)O 2 under irradiation studied by Molecular Dynamics simulations", submitted to Journal of Nuclear Materials.

  Periodic boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Details of the pair potential . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Diagram showing how OVITO calculates the occupancy

. . . . . . . . . . . . . . . . . . . . 14 xiii 2.1

  Snapshots of the evolution of extended defects as a function of dose (expressed in dpc) for Potashnikov potential at a) 300 K and b) 1600 K for 50% Pu. The red zone marked the regions where structure departs from perfect Ćuorite structure. The cyan lines are 1 /3⟨1 1 1⟩ Frank loops, the green lines are 1 /6⟨1 1 2⟩ Shockley partial dislocations, and the dark blue lines are1 /2⟨1 1 0⟩ perfect dislocations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0⟩ perfect dislocations. . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Evolution of the bulk modulus as a function of irradiated dose expressed in displacement per cation (dpc) for (a) Cooper and (b) Potashnikov potentials at 1600 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	4.11 Evolution of dislocation densities as a function of dose (expressed in dpc)
	obtained with Potashnikov potential at 300 and 1600 K for 50% Pu. (circle)
	1 /3⟨1 1 1⟩ Frank loops, (triangle) 1 /6⟨1 1 2⟩ Shockley partial dislocations, (star)
	1 /2⟨1 1

1 1⟩ Frank loops, the green lines are 1 /6⟨1 1 2⟩ Shockley partial dislocations, and the dark blue lines are 1 /2⟨1 1 0⟩ perfect dislocations. xvi 4.10 4.13 Evolution of the Zener ratio as a function of irradiated dose expressed in displacement per cation (dpc) for (a) Cooper and (b) Potashnikov potentials at 1600 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Evolution of the XRD pattern for UO 2 for Cooper potential as a function of dose 4.15 Evolution of the FWHM and dislocation densities for Cooper potential. . . 4.16 Evolution of the FWHM and dislocation densities for Potashnikov potential.

  Distribution of point defects in cluster (interstitial and vacancy) as a function of their size for Cooper potentials calculated from 5 keV PKA cascade and 50% Pu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.9 Distribution of point defects in cluster (interstitial and vacancy) as a function of their size for Cooper potentials calculated from 5 keV PKA cascade and 25% Pu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.10 Distribution of point defects in cluster (interstitial and vacancy) as a function of their size for Potashnikov potentials calculated from 5 keV PKA cascade and 25% Pu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	List of Tables
	1.1 MOX properties catalogue (ESNII+, 2015). . . . . . . . . . . . . . . . . .	8
	2.1 Type of physical function for the diferent interatomic potentials studied. . .	22
	xix

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.8 2.2 Parameters for the interatomic potentials. . . . . . . . . . . . . . . . . . . 23 2.3 Parameters for Cooper EAM term. . . . . . . . . . . . . . . . . . . . . . . 23 3.1 Parameters for the third degree polynomial that Ąts the evolution of the lattice parameters as the function of temperature and plutonium content. . . . . . . 40 3.2 Transition pressure at T = 0 K from Ćuorite 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. . . . . . . . . . . . . . . . . . . . . . . 51 4.1 Formation energies of Frenkel pairs and Schottky defects. . . . . . . . . . . . 67 4.2 Recombination times in ps for C-FPs for both Cooper and Potashnikov potentials at diferent Pu contents and temperatures. The symbol, ≪, corresponds to a spontaneous recombination (less than 1 ps) and the symbol, >, signiĄes that no recombination takes place within the 100 ps of the simulation. The number in parenthesis corresponds to the number of recombination that occurs over the averaged cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.1 ClassiĄcation of diferent numerical methods. . . . . . . . . . . . . . . . . 98 5.2 Recombination times in µs for C-FPs for Potashnikov potential at diferent Pu contents for the second-, third-and fourth-nearest neighbours. The number in parenthesis corresponds to the number of recombination that occurs over the averaged cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 xxi 1 MOX Fuel overview

Table 2 .

 2

	Potential	Coulomb Buckingham Morse EAM
	Yamada
	Arima
	Potashnikov
	Cooper

1 Ű Type of physical function for the diferent interatomic potentials studied.

Table 2 .

 2 2 Ű Parameters for the interatomic potentials. Details of the pair potential of Energy CRADA (Cooperative Research and Development Agreement) between two DOE labs and 3 companies. It is distributed by Sandia National Labs.

				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 0.78093 Pu-O 0.564005	--	--	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 2.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 3 .

 3 2 Ű Transition pressure at T = 0 K from Ćuorite 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. PuO 2 UO 2 PuO 2 UO 2 PuO 2 UO 2 PuO 2

					Ćuorite (Fm 3m) to
		Cotunnite	Rutile	Scrutinyite	Marcasite
		Pnma	P4 2 /mnm	Pbcn	Pnnm
	UO 2 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

Table 4 .

 4 2 Ű Recombination times in ps for C-FPs for both Cooper and Potashnikov potentials at diferent Pu contents and temperatures. The symbol, ≪, corresponds to a spontaneous recombination (less than 1 ps) and the symbol, >, signiĄes that no recombination takes place within the 100 ps of the simulation. The number in parenthesis corresponds to the number of recombination that occurs over the averaged cases.

	Pu	Temperature		Cooper			Potashnikov
	Content	(K)	1 st 2 nd	3 r d	4 th	1 st 2 nd	3 r d	4 th
		30	≪ >	>	>	≪ >	>	>
	0% Pu	300	≪ >	>	>	≪ >	>	>
		1600	≪ >	29.7±36.6 (8) 45.6±37.8 (2) ≪ >	2.5±3.6 (6)	10.0±7.0 (2)
		30	≪ >	>	>	≪ >	>	>
	25% Pu	300	≪ >	>	>	≪ >	5.5±8.2 (29) 51.7±68.5 (26)
		1600	≪ > 28.1±12.1 (50) 29.6±22.6 (18) ≪ >	6.5±5.7 (54) 18.5±17.4 (20)
		30	≪ >	>	>	≪ >	>	>
	50% Pu	300	≪ >	>	>	≪ > 11.7±15.6 (26) 19.8±6.0 (9)
		1600	≪ >	21.7±6.2 (55) 28.6±16.1 (28) ≪ >	20.2±9.1 (50) 21.5±18.1 (17)
		30	≪ >	>	>	≪ >	>	>
	75% Pu	300	≪ >	>	>	≪ >	7.3±6.0 (34) 15.6±21.6 (10)
		1600	≪ >	10.3±3.8 (56) 28.5±13.8 (34) ≪ >	8.0±5.4 (54) 29.8±22.0 (27)
		30	≪ >	>	>	≪ >	>	>
	100% Pu	300	≪ >	>	>	≪ >	5.6±4.1 (6)	15.0±0.0 (2)
		1600	≪ >	5.0±6.2 (8)	26.3±37.5 (2) ≪ >	5.6±7.0 (6)	0.7±0.1 (2)

Table 5 .

 5 1 Ű ClassiĄcation of diferent numerical methods.

			PROBLEM FORMULATION
			exact	stochastic
	METHOD OF	exact	Numerical analysis	Probability theory
	SOLUTION	stochastic Monte Carlo methods Simulation methods
	MC algorithms

Table 5 .

 5 2 Ű Recombination times in µs for C-FPs for Potashnikov potential at diferent Pu contents for the second-, third-and fourth-nearest neighbours. The number in parenthesis corresponds to the number of recombination that occurs over the averaged cases.

	Pu content	2 th	3 r d	4 th
	0%	28±29 (4) 40±54 (4)	-
	50%	16±36 (9) 16±54 (12) 0.8 (1)
	100%	17 (1)	29±50	-

Remerciements

A

Phase diagram of MOX

In order to understand better which states the MOX fuel goes through at diferent temperatures and Pu contents, the ternary phase diagram O-Pu-U will be assessed. This is shown in Figure A.1 for a system at 1500 K.

Here the legend "C1 M OX2" corresponds to the Ćuorite structure, PUO161 to PuO 1.61 and the rest is self-explicit. Above 1000 K, MOX forms a solid solution across their entire composition range. In the sub-stoichiometric range, the loss of O atoms is compensated by the reduction of Pu +4 to Pu +3 [START_REF] Moore | Development of a thermo-kinetic diffusion model for uo 2 and (u, pu)o 2 oxide fuels using the dictra code[END_REF]. On the other hand, the over-stoichiometric region consists on the extra oxygen interstitial and the oxidation of U +4 to U +5 at temperatures below 1000 K. The highest stoichiometric ratio is 2.2 for Pu content in the range 60 % > x > 32 % [START_REF] Guéneau | Thermodynamic modelling of advanced oxide and carbide nuclear fuels: description of the uŰpuŰoŰc systems[END_REF]. The deviation in stoichiometry as a function of Pu content at 1500 K is shown in Figure A.2.

At 1500 K, the stoichiometry changes as a function of the plutonium content. The sub-stoichiometry goes from a UO x slightly below 2 all the way up to PuO 1.61 at higher plutonium contents. The hyper-stoichiometry decreases as well, going from a value slightly above 2.2 to less that 2.1.

We now present the stability of the system for a wider temperature range. [START_REF] Markin | The uranium-plutonium-oxygen ternary phase diagram[END_REF]. The red dashed lines correspond to calculated equilibria involving Thermo-Calc Software and the thermodynamic model proposed by [START_REF] Guéneau | Thermodynamic modelling of advanced oxide and carbide nuclear fuels: description of the uŰpuŰoŰc systems[END_REF]. Taken from [START_REF] Vauchy | Efect of cooling rate on achieving thermodynamic equilibrium in uraniumŰplutonium mixed oxides[END_REF]

B

Extended lattice parameter

C

Cluster analysis

This appendix shows distribution of point defects in cluster (interstitial and vacancy) as a function of their size for Cooper and Potashnikov potentials calculated from all PKA energies. Cation (Pu and U) vacancy and interstitial will me named U for simplicity.