# First-principles calculations of uranium diffusion in uranium dioxide

Boris Dorado, David A Andersson, Christopher R Stanek, Marjorie Bertolus,

Blas P Uberuaga, Guillaume Martin, Michel Freyss, Philippe Garcia

► To cite this version:

Boris Dorado, David A Andersson, Christopher R Stanek, Marjorie Bertolus, Blas P Uberuaga, et al.. First-principles calculations of uranium diffusion in uranium dioxide. Physical Review B, 2012, 86 (3), pp.035110. 10.1103/PhysRevB.86.035110. cea-05504034

HAL Id: cea-05504034

https://cea.hal.science/cea-05504034v1

Submitted on 10 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# First-principles calculations of uranium diffusion in uranium dioxide

Boris Dorado,<sup>1,2,*</sup> David A. Andersson,<sup>1</sup> Christopher R. Stanek,<sup>1</sup> Marjorie Bertolus,<sup>2</sup> Blas P. Uberuaga,<sup>1</sup> Guillaume Martin,<sup>2</sup> Michel Freyss,<sup>2</sup> and Philippe Garcia<sup>2</sup>

<sup>1</sup>Materials Science and Technology Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA

<sup>2</sup>CEA, DEN, DEC, Centre de Cadarache, 13108 Saint-Paul-lez-Durance, France

(Received 26 January 2012; published 9 July 2012)

The present work reports first-principles DFT + U calculations of uranium self-diffusion in uranium dioxide (UO₂), with a focus on comparing calculated activation energies to those determined from experiments. To calculate activation energies, we initially formulate a point defect model for UO₂, that is valid for small deviations from stoichiometry. We investigate five migration mechanisms and calculate the corresponding migration barriers using both the LDA + U and GGA + U approximations. These energy barriers are calculated using the occupation matrix control scheme that allows one to avoid the metastable states that exist in the DFT + U approximation. The lowest migration barrier is obtained for a vacancy mechanism along the (110) direction. This mechanism involves significant contribution from the oxygen sublattice, with several oxygen atoms being displaced from their original position. The (110) vacancy diffusion mechanism is predicted to have lower activation energy than any of the interstitial mechanisms and comparison to experimental data for stoichiometric UO₂ also confirms this mechanism.

DOI: 10.1103/PhysRevB.86.035110

PACS number(s): 71.15.Mb, 71.27.+a

## I. INTRODUCTION

Uranium dioxide UO₂ is the most widely used nuclear fuel worldwide and its atomic transport properties are relevant to practically all engineering aspects of the material. During in-reactor operation, the fission of uranium atoms produces a wide variety of fission products that create point defects as they deposit their energy in the surrounding material. Since many of the fuel properties that govern fuel performance are influenced by point defects, it is of high importance to investigate their thermodynamic properties and migration mechanisms.

Activation energies for uranium self-diffusion have been reported for near stoichiometric and stoichiometric UO₂ single crystals ranging from 4.4 eV to 5.6 eV.<sup>1-4</sup> One of our goals here is to use first-principles calculations to provide some physical meaning to these activation energies. Although uranium self-diffusion coefficients have been determined experimentally,<sup>1-5</sup> the migration mechanisms remain unknown and discrepancies are observed even for purportedly stoichiometric compositions, both in relation to the actual value of the self-diffusion coefficient and the associated activation energy. An essential cause of scatter is the enhanced diffusion at the grain boundaries. Sabioni reports a grain boundary diffusion coefficient five orders of magnitude greater than the volume diffusion coefficient between 1773 and 1973 K in a reducing atmosphere.<sup>6</sup> In addition, it has been shown in previous work relating to oxygen self-diffusion<sup>7,8</sup> that one of the issues in obtaining a reliable set of data lies in monitoring the relevant thermodynamic conditions during the experiment, such as oxygen partial pressure, temperature, and the impurity content of the sample.

First-principles calculations based on the density functional theory (DFT) can be used to analyze the experimental studies mentioned above. It is well known, however, that standard approximations to DFT fail to describe accurately uranium dioxide because of the strong correlations among uranium 5f electrons. In order to entirely capture these correlation effects, one needs to use approximations beyond standard DFT, such

as the self-interaction correction,<sup>9</sup> the hybrid functionals,<sup>10-12</sup> or the DFT + U approximation.<sup>13-15</sup>

The DFT + U method has been widely used to study the behavior (formation and migration) of uranium and oxygen point defects in UO₂.<sup>16-22</sup> Despite this, the description of UO₂ remains challenging. It is now recognized that the DFT + U method creates a number of metastable states that make the search for the electronic ground state difficult and can lead to large errors in calculations of defect energies if no care is taken to control the correlated electronic states. To our knowledge, there are currently three methods to circumvent these difficulties: the occupation matrix control (OMC) scheme,<sup>8,21-26</sup> the U-ramping method,<sup>27</sup> and the quasiannealing method.<sup>28</sup> The OMC scheme was used in the present work. It was first developed by Jomard et al.<sup>23</sup> and Amadon et al.<sup>24</sup> then applied to uranium dioxide.<sup>8,21,22,25</sup> Compared to the U-ramping method or the QA approach, the OMC scheme allows us to more precisely control various valence states for uranium atoms in defective UO₂, such as U<sup>3+</sup> and U<sup>5+</sup> cations. We have already applied this scheme in our previous studies,<sup>21,22,29</sup> particularly in relation to oxygen self-diffusion.<sup>8</sup> The results have shown that this method can provide accurate physical properties and defect migration energies and we here extend this work to defect energies involving the cation sublattice.

Assessing defect formation and migration energy calculations against self-diffusion data requires a rigorous point defect model that establishes the point defect concentrations as a function of temperature and oxygen potential. In Sec. II, we establish such a point defect model and specify its range of validity. Analytical expressions for activation energies are then derived in two limiting cases when electronic defects are predominant over anion disorder and vice versa. Section III describes in detail all the energies that are required for the calculation of activation energies. In Sec. IV, we present our first-principles calculations. The energies corresponding to the mass balance equations required for the point defect model are evaluated using both neutral and charged supercell approaches.
Also, different uranium defect migration mechanisms and associated energy barriers are computed. Finally, in Sec. V, we analyze the point defect model using the formation energies calculated from first principles and we assess the model predictions against the experimental data reported in the literature.

## II. POINT DEFECT MODEL FOR UO$_{2\pm x}$

Experimental methods for direct measurement of formation and migration energies are often complex and thus difficult to implement. Self-diffusion properties, however, are directly dependent on these quantities. In order to correlate self-diffusion data with formation and migration energies, a point defect model is needed that enables one to relate the various defect concentrations to prevailing thermodynamic conditions. In this work, we rely on such a point defect model to compare defect energies calculated from first principles to uranium self-diffusion data, as was recently done in the simpler case of oxygen diffusion.⁸ The model developed is similar to others reported in the past³⁰,³¹ and is based on defect chemistry using the Kröger and Vink formalism.³² We assume that charged defect concentrations are governed by a set of simultaneous equations amongst which so-called mass-action laws are the expression of thermodynamic equilibrium.

### A. Applicability

The model presented in this section is only applicable to a limited region of the UO$_{2\pm x}$ phase diagram, typically $x < 2 \times 10^{-3}$ (this validity range will be demonstrated in Sec. V). This is because in the traditional defect chemistry approach that we adopt here, defects are assumed to be far apart, which means both that the configurational entropy is given by the dilute limit approximation (i.e., low defect concentrations) and that strong defect interactions that are known to occur in oxide systems are neglected. Methods have been developed to account for strong defect interactions, as well as the modifications to the configurational entropy due to defect clustering.³³ These approaches, however, are thought to be necessary for defect concentrations greater than approximately $10^{-2}$ mol fraction. Hence the equations described below are applicable to low deviations from stoichiometry only. In addition, the model only accounts for point defects in the strict sense of the term, i.e., it does not explicitly describe the formation of defect clusters that readily appear in UO$_{2}$ as deviation from stoichiometry increases.³⁴ Finally, it is assumed that all point defects are created thermally (i.e., intrinsic regime), which limits the applicability of the model to high temperatures, in practice, above approximately 1273 K. Below 1273 K, the regime is known to be extrinsic with the hole concentration being controlled by impurities.³⁵

### B. Model formulation

Magnetic susceptibility experiments on U$_{4}$O$_{9}$ have shown that this oxide contains a mixture of U$^{4+}$ and U$^{5+}$ ions.³⁶ In addition to this, first principles calculations²² have shown that two U$^{5+}$ ions are more stable than one U$^{6+}$ ion upon introduction of an oxygen interstitial. Hence we follow here the assumption of many other authors³⁵ that electronic disorder

is controlled by the disproportionation of two U$^{4+}$ (5f$^2$) ions to one U$^{5+}$ (5f$^1$) ion and one U$^{3+}$ (5f$^3$) ion, which may be written as follows using the Kröger-Vink notation:

$$2U_U^X \rightleftharpoons h^+ + e^-, \qquad (1)$$

where U$_{U}^X$ designates a uranium atom on a normal uranium lattice site (U$^{4+}$) and h$^+$ and e$^-$ designate holes (U$^{5+}$) and electrons (U$^{3+}$), respectively. Oxygen disorder on the anion sublattice results from the Frenkel equilibrium:

$$O_O^X + V_i^X \rightleftharpoons O_i'' + V_O, \qquad (2)$$

where O$_{O}^X$ is an oxygen atom on a normal oxygen lattice site and V$_{i}^X$ is a vacant interstitial site. The following equilibrium equation describes the incorporation of oxygen atoms from the gas phase into the solid:

$$\frac{1}{2}O_2 + V_i^X + 2U_U^X \rightleftharpoons O_i'' + 2h^+. \qquad (3)$$

Note that in Eq. (3), the two holes (U$^{5+}$) are not bound to the oxygen interstitial. Finally, the Schottky and uranium Frenkel disorder are described as follows:

$$\emptyset \rightleftharpoons V_U'' + 2V_O, \qquad (4)$$

and

$$U_U^X + V_i^X \rightleftharpoons U_i^{\cdots} + V_U^{\cdots}, \qquad (5)$$

In the dilute limit approximation, the configurational entropy terms for each of the chemical potentials of the different defect species are given by the corresponding site fractions. If $\Delta G_\alpha$ designates the Gibbs free energy of one of the above defect formation equations [electron-hole pair ($\alpha = eh$), oxygen Frenkel pair ($\alpha = FP_O$), oxygen interstitial ($\alpha = O_i$), Schottky defect ($\alpha = S$), and uranium Frenkel pair ($\alpha = FP_U$)], then thermodynamic equilibrium determines relations between activities, i.e., site fractions. Site fractions are proportional to defect concentrations and it can be shown that the equilibrium constants may be expressed as a function of defect concentrations normalized to the uranium site concentration. Equilibrium is therefore expressed as five relationships, each corresponding to the five chemical equilibria (1) to (5):

$$K_{eh} = \frac{[e^\prime][h^\prime]}{[U_U^X]^2} = \exp\left(-\frac{\Delta G_{eh}}{k_B T}\right), \qquad (6)$$

$$K_{FP_O} = \frac{[O_i''][V_O^{\ddot{O}}]}{[O_O^X][V_i^X]} = \exp\left(-\frac{\Delta G_{FP_O}}{k_B T}\right), \qquad (7)$$

$$K_{O_i} = \frac{[O_i''][h^\prime]^2}{\sqrt{P_{O_2}} [U_U^X]^2 [V_i^X]} = \exp\left(-\frac{\Delta G_{O_i}}{k_B T}\right), \qquad (8)$$

$$K_S = [V_U^{\ddots}][V_O^{\ddot{O}}]^2 = 4 \exp\left(-\frac{\Delta G_S}{k_B T}\right), \qquad (9)$$

and

$$K_{FP_U} = \frac{[V_U^{\ddots}][U_i^{\cdots}]}{[U_U^X][V_i^X]} = \exp\left(-\frac{\Delta G_{FP_U}}{k_B T}\right), \qquad (10)$$

where the square brackets represent normalized defect concentrations and $P_{O_2}$ the equilibrium oxygen partial pressure. Note that the factor of 4 on the right hand side of Eq. (9) is due to the fact that defect equilibrium constants are expressed as a function of concentrations and not site fractions. There
are nine unknowns to this model and in addition to the five defect equilibrium relationships, three equations express the constraints imposed by the crystalline structure and an additional one guarantees electroneutrality:

$$[U_U^X] + [V_U^{'''}] = 1, \qquad (11)$$

$$[O_O^X] + [V_O^{''}] = 2, \qquad (12)$$

$$[O_i^{''}] + [V_i^X] + [U_i^{''}] = 1, \qquad (13)$$

$$4 [V_U^{'''}] + 2 [O_i^{''}] + [e^+] = 4 [U_i^{''}] + 2 [V_O^{''}] + [h^+]. \qquad (14)$$

A first-order approximation of deviation from stoichiometry $x$ (for small $x$) may then be expressed from the defect concentrations that the model provides:

$$x = \frac{[O_i^{''}] - [V_O^{''}]}{1 - [V_U^{'''}] + [U_i^{''}]} \qquad (15)$$

### C. Analytical expressions for activation energies

Following Kofstad,³⁷ two limiting cases are discussed and a corresponding analytical expression for the uranium activation energy is provided assuming a vacancy or interstitial diffusion mechanism. All expressions are derived under the assumption of proximity to stoichiometric composition.

We assume $[e^+] \approx [h^+]$ (i.e., intrinsic ionisation predominoes) and $[U_U^X] \approx [V_i^X] \approx [O_O^X]/2 \approx 1$. Note that these assumptions will be validated later in Sec. V. In this case, the system of equations (6) to (14) may be solved analytically and in particular $[V_U^{'''}]$ and $[U_i^{''}]$ are expressed as follows:

$$[V_U^{'''}] = \frac{K_{O_i}^2 K_S}{K_{eh}^2 K_{FP_o}^2} p_{O_2}, \qquad (16)$$

$$[U_i^{''}] = \frac{K_{eh}^2 K_{FP_o}^2 K_{FP_u}}{K_{O_i}^2 K_S} \frac{1}{p_{O_2}}. \qquad (17)$$

Note that the above assumption ($[e^+] \approx [h^+]$) implies that the Fermi level is located in the middle of the band gap. This is important because the position of the Fermi level has an influence on the charge states of defects in $UO_2$.¹⁹,³⁸ Crocombette et al. have shown that for a Fermi level located in the middle of the band gap, the most stable defects are $O_i^{''}$ and $V_U^{'''}$ with no local charge compensation around the defects, which is consistent with our model. As for oxygen vacancies, we used a +2 charge state in our model (i.e., with two unbound $U^{3+}$) even though there are indications that they might be charged +1 when the Fermi level is located near the middle of the band gap.³⁸ However, our calculations indicate that the +2 charge state is favored up to a Fermi level positions that are higher than the midgap position. For sake of consistency, we therefore only account for oxygen vacancies in a +2 charge state.

Now, if diffusion proceeds via simple mechanisms, then the expression for the uranium self-diffusion coefficient is given by³⁹

$$D_U = f_{V_U}[V_U^{'''}]D_{V_U} + f_{U_i}[U_i^{''}]D_{U_i}, \qquad (18)$$

where $f_{V_U}$ and $f_{U_i}$ are the correlation factors for the vacancy and interstitial mechanisms. $D_{V_U}$ and $D_{U_i}$ are the vacancy and

interstitial diffusion coefficients, respectively. The latter may be expressed as

$$D_d = D_d^0 \exp \left( - \frac{E_m^d}{k_B T} \right), \qquad (19)$$

where $E_m^d$ ($E_m^{V_U}$ and $E_m^{U_i}$) are the vacancy and interstitial migration barriers, respectively. Combining Eqs. (16) to (19) yields the corresponding activation energies $E_a^{V_U}$ and $E_a^{U_i}$ that can be written as

$$E_a^{V_U} = 2E_{O_i} + E_S - 2E_{eh} - 2E_{FP_o} + E_{pO_2} + E_m^{V_U}, \qquad (20)$$

$$E_a^{U_i} = -2E_{O_i} - E_S + 2E_{eh} + 2E_{FP_o} + E_{FP_U} - E_{pO_2} + E_m^{U_i}, \qquad (21)$$

where $E_a$ are the energies (or strictly speaking enthalpies) associated with the Gibbs free energies $\Delta G_a$ of Eqs. (1) through (5).

Note that the term $E_{pO_2}$ that appears in Eqs. (20) and (21) has been introduced because $D_U$ depends on the oxygen partial pressure via Eqs. (16) and (17). Therefore assessing this theoretical approach to the experimental data would, in principle, require the knowledge of the oxygen potential at which the experiments were carried out. This information is not available. We know, however, that the experiments were carried out under reducing atmospheres generally obtained by introducing a proportion of hydrogen in the carrier gas, which activates the following equilibrium:

$$H_2O \rightleftharpoons H_2 + \frac{1}{2}O_2. \qquad (22)$$

As a result, the actual oxygen potential under which the experiments are carried out has an Arrhenius dependence. We have carried out measurements of oxygen potentials of Ar/5%H₂ at different temperatures and estimated an activation energy of 4.3 eV for the above equilibrium, which is close to the expected value of 5.1 eV that would be expected if the $H_2O \rightleftharpoons H_2 + \frac{1}{2}O_2$ equilibrium were buffering the oxygen partial pressure.

Note that if the uranium vacancy concentration is still predominant with respect to the uranium interstitial concentration in the substoichiometric region of the phase diagram (which is quite possible if $E_{FP_U}$ is high), then the apparent activation energy (i.e., the activation energy of an Arrhenius representation of $D_U$) should show no change of slope, even in the possible case where the material goes from a slightly hyper- to slightly hypo- stoichiometric composition as the temperature increases.

The other limiting case is when internal disorder predominoes (i.e., $[V_O^{''}] \approx [O_i^{''}]$). In this case, it is straightforward to show that the self-diffusion coefficient is no longer dependent upon equilibrium partial pressure, which is in contradiction with experimental results that show a sharp increase in $D_U$ as the deviation from stoichiometry increases. We would therefore expect this situation not to arise if the theory outlined in this work were consistent.
### III. ENERGIES REQUIRED FOR THE CALCULATION OF ACTIVATION ENERGIES

Given the point defect model described above, the following energies are required in order to calculate activation energies for uranium diffusion in UO$_{2\pm x}$: (1) the formation energies of oxygen interstitials ($E_{O_i}$), with two unbound U$^{5+}$ cations. (2) The formation energies of oxygen and uranium Frenkel pairs ($E_{FP_O}$ and $E_{FP_U}$), Schottky defects ($E_S$), and electron-hole pairs ($E_{eh}$). (3) The migration energy of uranium atoms ($E_m^{lu}$ and $E_m^{VU}$).

#### A. Calculations of defect formation energies and associated charge compensations

As mentioned in Sec. II, Eq. (3) describes the insertion of an oxygen interstitial in UO$_{2}$. In this equation, the two holes (U$^{5+}$) are not bound to the oxygen interstitial. This means that the formation energy should be calculated using a supercell in which there is no local charge compensation around the oxygen interstitial. In the same way, Eqs. (1), (2), (4), and (5) assume (i) the charged defects do not interact and (ii) there is no local charge compensation. The corresponding formation energies should therefore also be computed using defective supercells with no charge compensation around the defects.

In a standard DFT calculation (i.e., using a neutral supercell), however, local charge compensations always occur, i.e., there are always U$^{5+}$ or U$^{3+}$ cations created because of the presence of the defects. For instance, an oxygen interstitial always captures two electrons from the surrounding uranium atoms and gets a -2 charge, leaving two U$^{5+}$ cations as charge compensation. In order to be entirely consistent with the model described above, however, it is required to remove the local charge compensations in our calculations. This can be done by using charged defective supercells instead of neutral ones. By adding or removing electrons in the supercell, the charge of the interstitial atom (or vacancy) remains unchanged but the charge of the whole defect (i.e., extra atom or atom deficit along with the local charge compensation on uranium atoms) may be controlled through the following reactions: U$^{5+}$ + e$^{-}$ → U$^{4+}$ and U$^{3+}$ - e$^{-}$ → U$^{4+}$. The difference between a neutral and a charged supercell calculation therefore lies in the presence (or absence) of the local charge compensation.

#### B. Oxygen interstitial formation energy

The accurate determination of the formation energy of an interstitial oxygen remains challenging, both from the experimental and theoretical standpoints. In this section, we describe the various attempts made at calculating this value as well as the difficulties encountered, and we explain why the values obtained from first principles should be considered more reliable than the current experimental ones. We perform a careful comparison of our work with all previous DFT based studies and we infer an approximate value for this formation energy.

##### 1. Determination from experimental measurements

The (effective) formation energy of oxygen interstitials can be calculated from experimental measurements of the oxygen potential as a function of temperature at fixed UO$_{2+x}$

composition. The oxygen potential is given by

$$ \Delta G_{O_2} = RT \ln p_{O_2} = \Delta H_{O_2} - T \Delta S_{O_2}, \quad (23) $$

where $\Delta H_{O_2}/2$ corresponds to the (effective) interstitial formation energy in Eq. (3) calculated from DFT. $\Delta G_{O_2}$ has been measured by several investigators$^{40-48}$ and the results have been analyzed in thermodynamic assessments and reviews of the U-O system.$^{46,49-51}$ For example, based on experimental data from Markin and Bones,$^{41}$ Labroche *et al.* calculated $\Delta H_{O_2}/2$ to range from -1.4 to -1.5 eV between UO$_{2.01}$ and UO$_{2.082}$. $^{51}$ Similarly, Aronson and Belle provided estimates between -1.3 ($x$ = 2.013) and -1.6 ($x$ = 2.203) eV, which is also very close to the oxygen potentials predicted by the assessed CALPHAD thermodynamic model for the same concentration regime.$^{50}$ It should be noted that in this composition range, stoichiometric variations are dominated by clustering phenomena (hence the use of the term *effective*). Similarly, the lowering of the formation energy as function of increasing nonstoichiometry ($x$) follows from increased clustering. There are experimental measurements of $\Delta G_{O_2}$ as the stoichiometry approaches stoichiometric UO$_{2}$ and the corresponding $\Delta H_{O_2}$ values rapidly reaches more negative values,$^{52}$ however, this is not related to lowering of the interstitial formation energy but rather indicates the onset of oxygen vacancies and it is thus related to the transition from UO$_{2+x}$ to UO$_{2-x}$ defect chemistry. Consequently, it is difficult to derive accurate estimates of the interstitial formation energy from measurements of the oxygen chemical potentials. The formation energy of -1.3 eV measured at $x$ = 2.013 represents a lower limit for the interstitial formation energy relevant for the present analysis, but it cannot necessarily be used as an accurate estimate of the reaction in Eq. (3). This leads us to put more emphasis on the oxidation energies derived from DFT.

##### 2. Direct determination from GGA + U calculations

The oxygen interstitial formation energy can be determined directly from GGA + U calculations by evaluating the energy required for the following reaction:

$$ \frac{1}{2}O_2 + V_i^X + 2U_i^X \rightleftharpoons O_i'' + 2h\cdot. \quad (24) $$

This is done by calculating the total energies of the UO$_{2}$ perfect and defective supercells, as well as the total energy of the O$_{2}$ molecule. In our previous work, we have estimated the formation energy to be approximately zero using GGA + U and occupation matrix control,$^{8}$ with holes bound to the oxygen interstitial. In this case, the corresponding equation using the Kröger-Vink notation is

$$ \frac{1}{2}O_2 + V_i^X \rightleftharpoons O_i^X, \quad (25) $$

which is not consistent with the model that assumes randomly distributed and unbound U$^{5+}$ ions [see Eq. (3)]. It is therefore necessary to add to our value the binding energy of two holes with the oxygen interstitial, which was evaluated to be 0.4 eV using charged supercell calculations similar to those applied in Ref. 29. Our oxygen interstitial formation energy with unbound U$^{5+}$ is therefore 0.4 eV, which is surprisingly high given that UO$_{2}$ is known to oxidize easily.
Contrary to what Oxford and Chaka stated,⁵³ the DFT + U approximation is not likely the cause for such a high formation energy. It is rather because of the GGA-PBE functional that fails to accurately describe the O atom and the O₂ molecule, resulting in a calculated O₂ dissociation energy that is off by 30%, 20% of which is due to the PBE functional.⁵⁴ If we correct our value to exactly reproduce the experimental dissociation energy, the formation energy becomes -0.5 eV. Note that the magnitude of our correction is in good agreement with the corrections derived in studies of other metal oxides⁵³,⁵⁵⁻⁵⁷ and molecules.⁵⁴ Based on assessment of a large number of oxidation reactions for transition metal oxides, Wang *et al*.⁵⁶ proposed an empirical correction of -0.7 eV that is meant to capture errors for the O₂ molecule as well as provide a more accurate description of the filling of O 2p orbitals in oxides. If we use Wang's correction, our corrected formation energy becomes -0.3 eV.

### 3. Indirect determination using U₄O₉ compounds

Another way of assessing bounds for the interstitial formation energy is to compare the formation energy of U₄O₉ compounds with the formation energy of single interstitials. By calibrating the calculated data for U₄O₉ to experimental data for the heat of formation, which is a rather well-known quantity for this ordered compound, we are able to estimate a correction that can then be applied to the formation energy of single interstitials. This takes care of the reference energy of the O₂ molecule as well as other uncertainties related to filling the O 2p orbitals of oxygen interstitials.

One of the difficulties of this approach is that the experimental structure of U₄O₉ is very complex and beyond the reach of current DFT calculations. However, recent advances in alternative structure models for this compound based on the ordering patterns proposed by Andersson *et al*.⁵⁸ enable us to determine the U₄O₉ formation energy from DFT with reasonable accuracy.⁵⁹ Using the LDA + U approximation, we calculated the formation energy of U₄O₉ to be -1.4 eV. By comparing our calculated value with the experimental value (-1.8 eV)⁵⁰, we derive a correction of -0.4 eV. Assuming that this correction is constant (measured per interstitial oxygen ion) for all UO₂+x compositions, which is reasonable if the correction is related to the O₂ molecule and filling of the O 2p orbitals, we can apply the same value to the calculated interstitial formation energy. Using LDA + U, the corrected oxygen interstitial formation energy with unbound U⁵⁺ is found to be -0.6 eV after applying the previously derived correction term.

### 4. Assessment of O₁ formation energy

Table I presents a summary of all the DFT based values for the oxygen interstitial formation energy that were discussed in this section. We see from Table I that all DFT based values converged to an oxygen interstitial formation energy between -0.3 and -0.6 eV, i.e., approximately -0.5 eV, which is significantly higher (less negative) than the value derived from oxygen potential measurements away from stoichiometry. Note that for consistency, we did not apply the GGA + U corrections to the LDA + U values, and vice versa. Although a few uncertainties remain, this value probably represents the

<table><caption>TABLE I. Formation energy (in eV) of the oxygen interstitial in UO₂ with unbound U⁵⁺, determined from DFT based calculations. Three different corrections are applied, derived from (i) direct GGA + U calculations, (ii) LDA + U calculations using U₄O₉ compounds, and (iii) Wang's study.⁵⁶</caption><thead><tr><th></th><th>Direct<br>GGA + U</th><th>Wang <i>et al.</i><br>Ref. 56</th><th>U₄O₉<br>LDA + U</th></tr></thead><tbody><tr><td>GGA + U</td><td>-0.5</td><td>-0.3</td><td></td></tr><tr><td>LDA + U<sup>59</sup></td><td></td><td></td><td>-0.6</td></tr></tbody></table>

most accurate estimate of the interstitial formation energy to date. For completeness, we will perform our analysis using four different formation energies ranging from 0 to -1.5 eV (see Table V), thus covering the range of possible values obtained from theory and experiments.

### C. Frenkel pair, Schottky defect, and electron-hole pair formation energies

DFT values of Frenkel pair (uranium or oxygen) and Schottky defect formation energies are more accurate than the oxygen interstitial one because they do not involve any reference state. From the model described in Sec. II, it is seen that the formation energies required are those of isolated Frenkel pairs and Schottky defects, i.e., those with the corresponding defects created in two separate supercells in order to remove the interactions between them. Associated formation energies have been calculated many times before but only those from Ref. 29 are relevant here: in the latter work, Andersson *et al.* indeed used charged supercell calculations in order to remove the creation of U⁵⁺ and U³⁺ cations, which is consistent with our model. We will therefore use the formation energies from Ref. 29 for the oxygen Frenkel pair (3.3 eV), the uranium Frenkel pair (11.2 eV), and the Schottky defect (6.0 eV). It should be noted that these energies were obtained using the LDA functional. We will demonstrate in Sec. IV B, however, that LDA and GGA calculations yield similar results. As for the formation energy of an electron-hole pair [given in Eq. (1)], it has been calculated using charged supercell calculations and estimated to be approximately 1.7 eV. The corresponding equation is

$$2U^{4+} \rightleftharpoons U^{3+} + U^{5+}, \qquad (26)$$

where U⁴⁺ refers to a perfect 96-atom neutral supercell and U³⁺ (respectively, U⁵⁺) refers to a perfect 96-atom charged supercell in which an electron was added (respectively, removed). Note that the electron-hole formation energy that we calculated is consistent with the value of the band gap in UO₂, which is around 2 eV.

Table II displays a summary of all formation energies (E<sub>F</sub>, in eV) needed for the model. The only remaining unknowns in our model are therefore the migration energies for uranium in UO₂. Section IV is dedicated to their determination using both the LDA + U and GGA + U approximations in order to assess the differences between them.
<table><caption>TABLE II. Summary of all formation energies ($E_F$, in eV) needed by the point defect model.</caption><thead><tr><th>Defect</th><th><i>E</i><sub>F</sub> (eV)</th></tr></thead><tbody><tr><td>Oxygen interstitial</td><td>-1.5 to 0</td></tr><tr><td>Oxygen Frenkel pair</td><td>3.3</td></tr><tr><td>Uranium Frenkel pair</td><td>11.2</td></tr><tr><td>Schottky defect</td><td>6.0</td></tr><tr><td>Electron-hole pair</td><td>1.7</td></tr></tbody></table>

## IV. THEORETICAL DETERMINATION OF MIGRATION ENERGIES

### A. Computational details

We only consider the fluorite phase, which is stable above 30.8 K, since it is the experimentally relevant phase. It should be stressed that during the calculation of migration barriers, it is very difficult to maintain the perfect fluorite structure due to the breaking of symmetries that turn some atoms around the defect into the distorted Jahn-Teller phase of UO₂, stable below 30.8 K. However, we have shown in a previous work⁸ that this effect does not significantly affect the calculated migration barriers. All calculations were carried out using the VASP package.⁶⁰⁻⁶² The Kohn-Sham wave functions were calculated within the projector augmented-wave (PAW) formalism⁶³ with the compensation charge $\hat{n}$ truncated up to $L = 6$ (see Ref. 26 for implications of the compensation charge truncation). We used the Liechtenstein approach of DFT + *U* with a double-counting correction in the fully localized limit.¹⁴ For all calculations, the *U* and *J* parameters of the DFT + *U* approximation were set to 4.50 and 0.54 eV respectively, as determined by Yamazaki and Kotani,⁶⁴,⁶⁵ based on the analysis of x-ray photoemission spectra. Migration energies were calculated using the nudged elastic band (NEB) method,⁶⁶ using a 96-atom supercell with a collinear **1k** AFM order as an approximation of the real paramagnetic order. Note that the paramagnetic order is approximated because it cannot be modeled yet with the DFT + *U* approximation due to the prohibitive computational cost of the calculations. In UO₂, the influence on the migration energies of the **1k** AFM order approximation, and of the magnetic ordering in general, is currently unknown. In order to give a rough estimate, we have performed calculations in ferromagnetic UO₂ and we found that the migration energy for the most probable diffusion mechanism was approximately 0.5 eV higher than that in AFM UO₂. The difference is significant but we expect the **1k** AFM order to be closer to the actual paramagnetic state than a ferromagnetic state, because the **1k** AFM order keeps a total magnetic moment of zero (consistent with paramagnetism) though displaying a local ordering of the magnetic moments, instead of a completely random one. A 500 eV cutoff energy was used with a 2 × 2 × 2 Monkhorst-Pack **k**-point mesh. This ensured the convergence of the cell parameters and of the total energy to less than 10⁻³ Å and 10 meV/atom, respectively. For the calculation of fractional occupancies, we used a gaussian smearing with a smearing width of 0.1 eV. Finally, spin-orbit coupling (SOC) is neglected in all calculations. Our previous studies of oxygen point defect formation²² and migration⁸ have shown, however, that the DFT + *U* approximation can

<figure><img src="image_2.png" /><figcaption>FIG. 1. (Color online) Schematics of the interstitial mechanism in the 〈110〉 direction. Light grey and red spheres represent uranium and oxygen atoms, respectively. The dark blue sphere represents the initial interstitial uranium atom.</figcaption></figure>

still quantitatively describe transport phenomena in UO₂, even with the neglect of SOC. We can therefore expect this trend to be the same for uranium self-diffusion.

### B. Migration energies: comparison of LDA and GGA functionals

Five migration mechanisms were considered for uranium self-diffusion.

(1) Direct interstitial mechanism in 〈110〉 direction: a uranium atom at an octahedral interstitial site moves to the nearest octahedral interstitial site, along the 〈110〉 direction (see Fig. 1).

(2) Interstitialcy mechanism: a uranium atom at an octahedral interstitial site first kicks a uranium atom out of a lattice site, which in turn moves to the nearest octahedral interstitial site. There are two interstitialcy mechanisms depending on the direction taken by the kicked atom (see Fig. 2).

(3) Vacancy mechanism in the 〈100〉 direction: A uranium atom moves in the 〈100〉 direction to the nearest uranium vacancy (see Fig. 3). Note that this mechanism involves the uranium ion traveling through an interstitial site.

(4) Vacancy mechanism in the 〈110〉 direction: a uranium atom moves in the 〈110〉 direction to the nearest uranium vacancy (see Fig. 4).

<figure><img src="image_3.png" /><figcaption>FIG. 2. (Color online) Schematics of the two interstitialcy mechanisms, depending on the direction taken by the kicked atom: collinear (top) and noncollinear (bottom). Light grey and red spheres represent uranium and oxygen atoms, respectively. The dark blue sphere represents the initial interstitial uranium atom.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3. (Color online) Schematics of the vacancy mechanism in the (100) direction. Light grey and red spheres represent uranium and oxygen atoms, respectively. The small yellow spheres represent uranium vacancies.</figcaption></figure>

(5) Vacancy mechanism in the (110) direction with contribution of oxygen sublattice: a uranium atom moves in the (110) direction to the nearest uranium vacancy. Compared to the previous mechanism, this mechanism involves a significant contribution from displacements of the oxygen sublattice around the migrating U atom (see Fig. 5).²⁹

Table III reports the migration energies calculated for the mechanisms described above, both in the LDA + U and GGA + U approximations with the equilibrium volumes of the corresponding defective supercells (≈5.45 and ≈5.53 Å for LDA and GGA, respectively). Note that we only considered the noncollinear interstitialy mechanism within both LDA + U and GGA + U. The collinear barrier was calculated within LDA + U to be 3.7 eV. Even though this barrier is lower than the noncollinear mechanism, the conclusions of our work are not affected given the extremely large activation energies of all interstitial mechanisms, irrespective of the mechanism considered (see Sec. V).

We see from Table III that for both approximations the vacancy mechanism in the (110) direction is found to have the lowest migration barrier when assisted by the oxygen sublattice. The calculated migration barriers are 4.8 and 3.6 eV in LDA + U and GGA + U, respectively. The former value is 1 eV higher than the experimental value reported by Matzke,⁴ which is as low as ≈2.4 eV. While it is of the same order of magnitude, we would expect the DFT results to be closer to the experimental value if the correct mechanism had been considered. The experimental migration barrier reported by Matzke for uranium self-diffusion might therefore correspond to diffusion of defect clusters rather than single uranium vacancies.²⁹

As mentioned in a previous work,²⁹ LDA + U always yields higher migration energies than GGA + U. This could be due to the fact that the LDA equilibrium volume is smaller

<figure><img src="image_5.png" /><figcaption>FIG. 4. (Color online) Schematics of the vacancy mechanism in the (110) direction. Light grey and red spheres represent uranium and oxygen atoms, respectively. The small yellow spheres represent uranium vacancies.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5. (Color online) Transition state for the oxygen-assisted vacancy mechanism (left), compared to the regular vacancy mechanism along the (110) direction (right). Light grey and red spheres represent uranium and oxygen atoms, respectively. The dark blue sphere represents the uranium atom migrating to the nearest-neighbor uranium vacancy along the (110) direction. In the assisted mechanism, several oxygen atoms have moved significantly away from their fluorite position.</figcaption></figure>

than that calculated by GGA, which in turn has an influence on the migration energies. To assess the influence of the cell volume, we performed LDA calculations but retained the GGA equilibrium volume. The results are reported in Table IV and are denoted as “LDA + U*.” We see from Table IV that the LDA + U* results are closer to the GGA + U values, indicating that the calculated migration barriers are almost independent of the approximation used for the exchange correlation functional. Therefore, if the GGA + U results reported in Table III are lower than the LDA + U results, it is mainly because the equilibrium volumes differ, even though additional calculations are required to firmly establish this connection. We have also compared point defect energies between LDA + U and GGA + U and, unlike uranium migration, good agreement is found at the respective equilibrium volumes. Apparently, the higher equilibrium volume for GGA + U primarily influences migration barriers of large species such as uranium ions. This emphasizes that we can use the LDA point defect formation energies calculated in Ref. 29 as input data for the model presented in Sec. II.

## V. DISCUSSION

The model presented in Sec. II has been applied with the formation energies defined and calculated in Sec. III. Figure 6 represents the variations of all point defect concentrations as a function of equilibrium oxygen partial pressure calculated at 1873 K, except [Uᵢ⁻] for which concentrations are too low. Deviation from stoichiometry x [as calculated from Eq. (15)] is also represented.

<table><caption>TABLE III. Migration energies (in eV) for the various uranium migration mechanisms considered in the fluorite phase of UO₂ given by the DFT + U approximation.</caption><thead><tr><th rowspan="2">Mechanism</th><th colspan="2">Migration energy (eV)</th></tr><tr><th>LDA + U</th><th>GGA + U</th></tr></thead><tbody><tr><td>Direct interstitial</td><td>8.8</td><td>7.9</td></tr><tr><td>Interstitialy noncollinear</td><td>4.7</td><td>4.1</td></tr><tr><td>Vacancy (100)</td><td>7.6</td><td>7.2</td></tr><tr><td>Vacancy (110)</td><td>6.1</td><td>5.5</td></tr><tr><td>Vacancy (110)(O displacement)</td><td>4.8</td><td>3.6</td></tr></tbody></table>
<table><caption>TABLE IV. Migration energies (in eV) for uranium self-diffusion in the fluorite phase of UO₂ given by the LDA + *U* approximation with the GGA + *U* equilibrium volume.</caption><thead><tr><th rowspan="2">Mechanism</th><th colspan="2">Migration energy (eV)</th></tr><tr><th>LDA + U*</th><th>GGA + U</th></tr></thead><tbody><tr><td>Direct interstitial</td><td>8.1</td><td>7.9</td></tr><tr><td>Interstitialcy noncollinear</td><td>4.3</td><td>4.1</td></tr><tr><td>Vacancy (100)</td><td>7.3</td><td>7.2</td></tr><tr><td>Vacancy (110)</td><td>5.8</td><td>5.5</td></tr><tr><td>Vacancy (110)(O displacement)</td><td>3.8</td><td>3.6</td></tr></tbody></table>

It can be seen from Fig. 6 that defect formation energies lead to uranium defect concentrations that are far below both the majority oxygen defect and electronic defect concentrations over the entire range of oxygen potential calculated. This is a marked improvement over previous point defect models³⁸,⁶⁷,⁶⁸ which, in combination with the formation energies calculated in those studies, have a tendency to predict more stable uranium vacancies in the hyperstoichiometric region of the phase diagram. It should be noted that we used the oxygen interstitial formation energy of -0.5 eV, but the diagram is not significantly affected if this formation energy ranges from -1.5 to 0 eV. The figure also indicates that at compositions near to stoichiometry (i.e., where the model is thought to apply theoretically), the electroneutrality equation that applies is $[e'] \approx [h']$. This is also encouraging since, as pointed out in Sec. II, the relationship $[V_{O}'] \approx [O_{i}']$ is inconsistent with the oxygen potential dependence of the uranium self-diffusion coefficient. In fact, $[e'] \approx [h']$ applies over a range of oxygen partial pressures, the lower and upper bounds of which may be defined as follows:

$$[h'] \approx 2 \times [V_{O}'], \qquad (27)$$

and

$$[e'] \approx 2 \times [O_{i}'], \qquad (28)$$

<figure><img src="image_7.png" /><figcaption>FIG. 6. (Color online) Calculated changes as a function of oxygen potential of various defect concentrations and deviation from stoichiometry *x*.</figcaption></figure>

The corresponding oxygen partial pressures are given by

$$p_{O_2}^{\min} = \frac{K_{eh}^2}{2K_{O_i}^2} \qquad (29)$$

and

$$p_{O_2}^{\max} = \frac{4K_{eh}K_{FPo}^2}{K_{O_i}^2}. \qquad (30)$$

A numerical application of Eqs. (29) and (30) yields values of oxygen potentials of $2 \times 10^{-22}$ and $5 \times 10^{-15}$, respectively, and an upper bound value of deviation from stoichiometry of approximately ±0.002.

We now turn to the calculation of activation energies using Eqs. (20) and (21). Table V reports the activation energies calculated using both the LDA + *U* and LDA + *U*\* migration barriers. Since the GGA + *U* and LDA + *U*\* barriers are very close, we did not create a separate entry for the GGA + *U*. As mentioned in Sec. III B, we use a range of different values for the oxygen interstitial formation energy ($E_{O_i} = -1.5, -1.0, -0.5$, and $-0.05$ eV) given the uncertainties associated with this value. Having said that, based on our earlier analysis we believe $-0.5$ eV to be the most representative. Finally, the values in Table V correspond to activation energies calculated with the experimentally measured activation energy for oxygen partial pressure (4.3 eV), while the values in brackets correspond to those calculated with the value of 5.1 eV that we would expect if the $H_2O \rightleftharpoons H_2 + \frac{1}{2}O_2$ equilibrium were buffering the oxygen partial pressure.

We can draw several conclusions from the results of Table V. First, it is seen that the LDA + *U* activation energies are approximately 1 eV higher than the LDA + *U*\* values. This is significant and highlights the importance of volume effects in the migration process. Second, we notice that our model does not favor low values for the oxygen interstitial formation energy. For $E_{O_i} < -1.00$ eV, the resulting activation energies are lower than all experimental values, especially for LDA + *U*\* results. Finally, interstitial mechanisms may be completely discarded since they yield activation energies for uranium migration which are all higher than 14 eV.

From this study, the oxygen-assisted uranium vacancy mechanism emerges as the most likely mechanism for uranium diffusion in near stoichiometric UO₂, with an apparent activation energy derived from DFT that is close to that reported by Sabioni under a reducing atmosphere (4.4 eV).² In addition, the theoretical approach may be used to rationalize the two studies relating to uranium diffusion in near stoichiometric UO₂ reported by Matzke and Sabioni.²,⁴ Indeed, we see from our study that the term “stoichiometric” covers deviations from stoichiometry of ±0.002. Although these values are small, they correspond to equilibrium oxygen partial pressure variations of over six orders of magnitude. Figure 6 shows that a change in *x* from $10^{-3}$ to 0 (not measurable using standard thermogravimetric techniques) covers a three-order-of-magnitude change in the oxygen partial pressure. Changes in the carrier gas used in different experiments could lead to such variations in the oxygen potential resulting in small undetectable changes in the oxygen content of the material. Hence assuming vacancy assisted uranium diffusion and further assuming that the uranium vacancy concentration
<table><caption>TABLE V. Calculated activation energies for uranium diffusion in UO$_{2}$$_{±}$$_{x}$ (in eV), depending on the uncertainties on E$_{Oi}$ and the mechanism involved for uranium migration. The activation energies in brackets correspond to those obtained with E$_{PO2}$ = 5.1 eV while the others correspond to E$_{PO2}$ = 4.3 eV.</caption><thead><tr><th></th><th colspan="4" style="text-align: center;">E<sub>Oi</sub> (eV)</th></tr><tr><th></th><th style="text-align: center;">-1.5</th><th style="text-align: center;">-1.0</th><th style="text-align: center;">-0.5</th><th style="text-align: center;">0.0</th></tr></thead><tbody><tr><td colspan="5" style="text-align: left;"><strong>LDA + U</strong></td></tr><tr><td style="text-align: left;">Direct interstitial</td><td style="text-align: center;">22.7(21.9)</td><td style="text-align: center;">21.7(20.9)</td><td style="text-align: center;">20.7(19.9)</td><td style="text-align: center;">19.7(18.9)</td></tr><tr><td style="text-align: left;">Interstitialy noncollinear</td><td style="text-align: center;">18.6(17.8)</td><td style="text-align: center;">17.6(16.8)</td><td style="text-align: center;">16.6(15.8)</td><td style="text-align: center;">15.6(14.8)</td></tr><tr><td style="text-align: left;">Vacancy (100)</td><td style="text-align: center;">4.9(5.7)</td><td style="text-align: center;">5.9(6.7)</td><td style="text-align: center;">6.9(7.7)</td><td style="text-align: center;">7.9(8.7)</td></tr><tr><td style="text-align: left;">Vacancy (110)</td><td style="text-align: center;">3.4(4.2)</td><td style="text-align: center;">4.4(5.2)</td><td style="text-align: center;">5.4(6.2)</td><td style="text-align: center;">6.4(7.2)</td></tr><tr><td style="text-align: left;">Vacancy (110)(O displacement)</td><td style="text-align: center;">2.1(2.9)</td><td style="text-align: center;">3.1(3.9)</td><td style="text-align: center;">4.1(4.9)</td><td style="text-align: center;">5.1(5.9)</td></tr><tr><td colspan="5" style="text-align: left;"><strong>LDA + U*</strong></td></tr><tr><td style="text-align: left;">Direct interstitial</td><td style="text-align: center;">22.0(21.2)</td><td style="text-align: center;">21.0(20.2)</td><td style="text-align: center;">20.0(19.2)</td><td style="text-align: center;">19.0(18.2)</td></tr><tr><td style="text-align: left;">Interstitialy noncollinear</td><td style="text-align: center;">18.2(17.4)</td><td style="text-align: center;">17.2(16.4)</td><td style="text-align: center;">16.2(15.4)</td><td style="text-align: center;">15.2(14.4)</td></tr><tr><td style="text-align: left;">Vacancy (100)</td><td style="text-align: center;">4.6(5.4)</td><td style="text-align: center;">5.6(6.4)</td><td style="text-align: center;">6.6(7.4)</td><td style="text-align: center;">7.6(8.4)</td></tr><tr><td style="text-align: left;">Vacancy (110)</td><td style="text-align: center;">3.1(3.9)</td><td style="text-align: center;">4.1(4.9)</td><td style="text-align: center;">5.1(5.9)</td><td style="text-align: center;">6.1(6.9)</td></tr><tr><td style="text-align: left;">Vacancy (110)(O displacement)</td><td style="text-align: center;">1.1(1.9)</td><td style="text-align: center;">2.1(2.9)</td><td style="text-align: center;">3.1(3.9)</td><td style="text-align: center;">4.1(4.9)</td></tr></tbody></table>

is proportional to the equilibrium oxygen partial pressure would explain discrepancies observed in both studies. These uncertainties warrant new uranium self-diffusion studies in which oxygen partial pressure and temperature are monitored carefully.

## VI. CONCLUSION

We report DFT + U calculations of uranium self-diffusion in UO$_{2}$ and carefully compare these calculated values to existing experimental data. To make this comparison possible, a point-defect model based on mass balance equations involving the most probable defects was formulated and the corresponding formation energies were calculated. These calculations show that at stoichimetric composition, electronic defects constitute the predominant defect population, followed by anion and finally cation disorder (mainly in the form of uranium vacancies).

In addition, migration barriers for different mechanisms were calculated and an analytical expression is derived for the uranium self-diffusion activation energy assuming a vacancy and an interstitial mechanism. By comparing the theoretical results to existing data, it is shown that the oxygen-mediated uranium vacancy mechanism is the most probable diffusion mechanism in UO$_{2}$. The lowest migration barrier is obtained for the movement of a vacancy along the (110) direction

involving a substantial concerted distortion of the anion sublattice.

Finally, it is shown that a careful comparison of a theoretical approach to experimental self-diffusion data, although promising, suffers to some extent from a lack of data. These results encourage us to pursue this effort and warrant a thorough investigation of cation self-diffusion under controlled oxygen partial pressures.

## ACKNOWLEDGMENTS

Jean-Paul Crocombette is gratefully acknowledged for discussions and suggestions to improve the original manuscript. Work at CEA/DEN was funded by the MATAV Nuclear Ceramics Basic Research Program. This research is partly supported by the European Commission through the FP7 F-BRIDGE project (Contract No. 211690) and was performed using French HPC resources available at CCRT and CINES. Work at Los Alamos National Laboratory was sponsored by the Nuclear Energy Advanced Modeling and Simulation (NEAMS) program, under the Fuels Integrated Performance and Safety Code (IPSC) project (work package numbers FTLA11MS0603 and MS-12LA060207). Los Alamos National Laboratory, an affirmative action/equal opportunity employer, is operated by Los Alamos National Security, LLC, for the National Nuclear Security Administration of the US Department of Energy under contract DE-AC52-06NA25396.

<footnote>*Corresponding author: boris.dorado@cea.fr</footnote>

<footnote n="1">¹D. K. Reimann and T. S. Lundy, *J. Amer. Ceram. Soc.* **52**, 511 (1969).</footnote>

<footnote n="2">²A. C. S. Sabioni, W. B. Ferraz, and F. Millot, *J. Nucl. Mater.* **257**, 180 (1998).</footnote>

<footnote n="3">³H. Matzke, *J. Nucl. Mater.* **30**, 26 (1969).</footnote>

<footnote n="4">⁴H. Matzke, *J. Chem. Soc. Faraday Trans. II* **83**, 1121 (1987).</footnote>

<footnote n="5">⁵J. F. Marin and P. Contamin, *J. Nucl. Mater.* **30**, 16 (1969).</footnote>

<footnote n="6">⁶A. C. S. Sabioni, W. B. Ferraz, and F. Millot, *J. Nucl. Mater.* **278**, 364 (2000).</footnote>

<footnote n="7">⁷P. Garcia, M. Fraczkiewicz, C. Davoisne, B. Pasquet, G. Carlot, D. Siméone, G. Baldinozzi, and C. Petot, *J. Nucl. Mater.* **400**, 112 (2010).</footnote>

<footnote n="8">⁸B. Dorado, P. Garcia, G. Carlot, C. Davoisne, M. Fraczkiewicz, B. Pasquet, M. Freyss, C. Valot, G. Baldinozzi, D. Siméone, and M. Bertolus, *Phys. Rev. B* **83**, 035126 (2011).</footnote>
⁹L. Petit, A. Svane, Z. Szotek, and W. Temmerman, *Science* **301**, 498 (2003).

¹⁰C. Adamo and V. Barone, *J. Chem. Phys.* **110**, 6158 (1999).

¹¹J. Heyd, G. E. Scuseria, and M. Ernzerhof, *J. Chem. Phys.* **118**, 8207 (2003).

¹²I. D. Prodan, G. E. Scuseria, and R. L. Martin, *Phys. Rev. B* **76**, 033101 (2007).

¹³V. I. Anisimov, J. Zaanen, and O. K. Andersen, *Phys. Rev. B* **44**, 943 (1991).

¹⁴A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, *Phys. Rev. B* **52**, R5467 (1995).

¹⁵S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, *Phys. Rev. B* **57**, 1505 (1998).

¹⁶M. Iwasawa, Y. Chen, Y. Kaneta, T. Ohnuma, H.-Y. Geng, and M. Kinoshita, *Mater. Trans.* **47**, 2651 (2006).

¹⁷H. Y. Geng, Y. Chen, Y. Kaneta, M. Iwasawa, T. Ohnuma, and M. Kinoshita, *Phys. Rev. B* **77**, 104120 (2008).

¹⁸H. Y. Geng, Y. Chen, Y. Kaneta, and M. Kinoshita, *Phys. Rev. B* **77**, 180101(R) (2008).

¹⁹P. Nerikar, T. Watanabe, J. S. Tulenko, S. R. Phillpot, and S. B. Sinnott, *J. Nucl. Mater.* **384**, 61 (2009).

²⁰B. Dorado, M. Freyss, and G. Martin, *Europ. Phys. J. B* **69**, 203 (2009).

²¹B. Dorado, B. Amadon, M. Freyss, and M. Bertolus, *Phys. Rev. B* **79**, 235125 (2009).

²²B. Dorado, G. Jomard, M. Freyss, and M. Bertolus, *Phys. Rev. B* **82**, 035114 (2010).

²³G. Jomard, B. Amadon, F. Bottin, and M. Torrent, *Phys. Rev. B* **78**, 075125 (2008).

²⁴B. Amadon, F. Jollet, and M. Torrent, *Phys. Rev. B* **77**, 155104 (2008).

²⁵B. Dorado, J. Durinck, P. Garcia, M. Freyss, and M. Bertolus, *J. Nucl. Mater.* **400**, 103 (2010).

²⁶B. Dorado, B. Amadon, G. Jomard, M. Freyss, and M. Bertolus, *Phys. Rev. B* **84**, 096101 (2011).

²⁷B. Meredig, A. Thompson, H. A. Hansen, C. Wolverton, and A. van de Walle, *Phys. Rev. B* **82**, 195128 (2010).

²⁸H. Y. Geng, Y. Chen, Y. Kaneta, M. Kinoshita, and Q. Wu, *Phys. Rev. B* **82**, 094106 (2010).

²⁹D. A. Andersson, B. P. Uberuaga, P. V. Nerikar, C. Unal, and C. R. Stanek, *Phys. Rev. B* **84**, 054105 (2011).

³⁰F. A. Kröger, Z. Phys. Chem. **49**, 178 (1966).

³¹K. Park and D. R. Olander, *J. Nucl. Mater.* **187**, 89 (1992).

³²F. A. Kröger and H. J. Vink, *Solid State Physics*, edited by F. Seitz and D. Turnbull (Academic Press, 1956), Vol. III, pp. 307–435.

³³S. Ling, *Phys. Rev. B* **49**, 864 (1994).

³⁴B. T. M. Willis, *Acta Cryst.* **A34**, 88 (1978).

³⁵P. Ruello, *J. Am. Ceram. Soc.* **88**, 604 (2005).

³⁶K. Gotoo, S. Nomura, and K. Naito, *J. Phys. Chem. Solids* **26**, 1679 (1965).

³⁷P. Kofstad, *Nonstoichiometry, Diffusion, and Electrical Conductivity in Binary Oxides*, edited by W. Interscience (John Wiley & Sons Inc., New York, 1972) pp. 38–40.

³⁸J.-P. Crocombette, D. Torumba, and A. Chartier, *Phys. Rev. B* **83**, 184107 (2011).

³⁹J. Philibert, *Atom Movements, Diffusion and Mass Transport in Solids* (EDP Sciences, 1991), pp. 139–145.

⁴⁰S. Aronson and J. Belle, *J. Chem. Phys.* **29**, 151 (1958).

⁴¹T. L. Markin and R. J. Bones, Report AERE-R4042 HL 62/2187 (UKAEA, 1962).

⁴²K. Hagemark and M. Broli, *J. Inorg. Nucl. Chem.* **28**, 2837 (1966).

⁴³T. L. Markin, V. J. Wheeler, and R. J. Bones, *J. Inorg. Nucl. Chem.* **30**, 807 (1968).

⁴⁴Y. Saito, *J. Nucl. Mater.* **51**, 112 (1974).

⁴⁵K. Naito and N. Kamegashira, *Adv. Nucl. Sci. Technol.* **9**, 99 (1976).

⁴⁶T. B. Lindemer and T. M. Besmann, *J. Nucl. Mater.* **130**, 473 (1985).

⁴⁷J. F. Babelot, R. W. Ohse, and M. Hoch, *J. Nucl. Mater.* **137**, 144 (1986).

⁴⁸A. Nakamura and T. Fujino, *J. Nucl. Mater.* **149**, 80 (1987).

⁴⁹P.-Y. Chevalier, E. Fischer, and B. Cheynet, *J. Nucl. Mater.* **303**, 1 (2002).

⁵⁰C. Gueneau, M. Baichi, D. Labroche, C. Chatillon, and B. Sundman, *J. Nucl. Mater.* **304**, 161 (2002).

⁵¹D. Labroche, O. Dugne, and C. Chatillon, *J. Nucl. Mater.* **312**, 21 (2003).

⁵²H. Matzke, *J. Nucl. Mater.* **208**, 18 (1994).

⁵³G. A. E. Oxford and M. Chaka, *J. Phys. Chem. C* **115**, 16992 (2011).

⁵⁴S. Kurth, J. P. Perdew, and P. Blaha, *Int. J. Quantum Chem.* **75**, 889 (1999).

⁵⁵W. Li, C. Stampfl, and M. Scheffler, *Phys. Rev. B* **65**, 075407 (2002).

⁵⁶L. Wang, T. Maxisch, and G. Ceder, *Phys. Rev. B* **73**, 195107 (2006).

⁵⁷A. Jain, G. Hautier, S. P. Ong, C. J. Moore, C. C. Fischer, K. A. Persson, and G. Ceder, *Phys. Rev. B* **84**, 045115 (2011).

⁵⁸D. A. Andersson, J. Lezama, B. P. Uberuaga, C. Deo, and S. D. Conradson, *Phys. Rev. B* **79**, 024110 (2009).

⁵⁹D. A. Andersson *et al.* (unpublished).

⁶⁰G. Kresse and J. Hafner, *Phys. Rev. B* **47**, 558 (1993).

⁶¹G. Kresse and J. Furthmüller, *Comput. Mater. Sci.* **6**, 15 (1996).

⁶²G. Kresse and J. Furthmüller, *Phys. Rev. B* **54**, 11169 (1996).

⁶³P. E. Blöchl, *Phys. Rev. B* **50**, 17953 (1994).

⁶⁴T. Yamazaki and A. Kotani, *J. Phys. Soc. Jpn.* **60**, 49 (1991).

⁶⁵A. Kotani and T. Yamazaki, *Prog. Theor. Phys.* **108** (Suppl.), 117 (1992).

⁶⁶G. Henkelman, B. P. Uberuaga, and H. Jonsson, *J. Chem. Phys.* **113**, 9901 (2000).

⁶⁷J. P. Crocombette, F. Jollet, T. N. Le, and T. Petit, *Phys. Rev. B* **64**, 104107 (2001).

⁶⁸M. Freyss, T. Petit, and J. P. Crocombette, *J. Nucl. Mater.* **347**, 44 (2005).