Is excess faster than deficient? A molecular-dynamics study of oxygen-interstitial and oxygen-vacancy diffusion in CeO2

The diffusivity of oxygen interstitials (Di) and of oxygen vacancies (Dv) in fluorite-structured CeO2 was studied by means of classical molecular dynamic simulation techniques. Simulations were performed on cells that were either oxygen abundant or oxygen deficient at temperatures 1500 ≤ T / K ≤ 2000 for defect site fractions 0.18% ≤ ni/v ≤ 9.1%. In general, we found that at a given temperature T and defect site fraction ni/v the vacancy diffusivity Dv was higher than the interstitial diffusivity Di. Isothermal values of Di and Dv were constant at low defect site fractions (ni/v < 0.91%), but the behaviour diverged at higher ni/v: whereas Dv decreased at higher nv, Di increased at higher ni. The analysis also yielded, as a function of ni/v, activation enthalpies (ΔHmig) and entropies (ΔSmig) of vacancy migration and of interstitial migration. A constant value of Δ H mig , v ≈ 0.6 eV was found for low nv, with increases in Δ H mig , v observed for nv > 0.91%. For low ni a constant value of Δ H mig , i ≈ 1.4 eV was found, with a surprising decrease in Δ H mig , i for ni > 0.91%. The effect of dopants on the behaviour of the defect diffusivities was also studied. Doping with Gd3+ had a detrimental effect on vacancy diffusion, with a slight decrease in Dv and an increase in Δ H mig , v being observed. Donor doping with Nb5+, in contrast, was beneficial, resulting in higher Di and a decrease in Δ H mig , i . We suggest that the migration mechanism of oxygen interstitials in CeO2, non-collinear interstitialcy, is responsible for the lower defect diffusivity and higher migration barrier.

At the most basic level the process of ion migration in an ionic solid consists of a charged particle moving within a lattice of other charged particles. On this basis there appears to be no reason why oxygen interstitials should on principle be either far less mobile or far more mobile than oxygen vacancies [14]. One class of materials, for which both cases have been studied, is the family of AO 2 fluorite-structured oxides, such as ZrO 2 , HfO 2 , CeO 2 , ThO 2 , UO 2 and PuO 2 . Most of the published computational studies on these materials, however, only report activation barriers for oxygen-vacancy migration (from static calculations either based on densityfunctional-theory calculations [23][24][25][26][27][28][29] or based on empirical pair-potentials calculations [30][31][32][33][34][35]). The activation barrier for vacancy migration is predicted to be in the range of 0.4eV to 0.6 eV. For the few studies, in Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. which oxygen-interstitial migration has been examined, higher values are reported, in the range of 0.8 eV to 1.1 eV [25,[36][37][38][39][40][41] (with the exception of HfO 2 , for which a barrier of 1.8 eV is reported [42]).
There have been attempts in the literature to compare interstitial and vacancy transport in oxides [1,11,12,14,16,43], but they deal with the subject neither systematically nor in depth. In this study, taking the cubic (high-symmetry, isotropic) oxide CeO 2 as our model AO 2 system, we examine systematically the question of which point-defects (oxygen interstitials or oxygen vacancies) exhibit the higher diffusivity. To this end, we employ molecular dynamics (MD) simulation techniques, as they give us, from data obtained as a function of temperature, not only the activation enthalpies (ΔH mig ) and entropies (ΔS mig ) of migration, but also absolute values for the diffusion coefficients. In particular, we compare the defect diffusivities D i v as a function of defect site fraction n i/v : since D i v should be independent of n i/v for dilute solutions, deviations from this behaviour allows us to study defect-defect self-interactions and to determine the defect concentration up to which the dilute solution treatment is still applicable, as far as point-defect migration is concerned. Lastly, we look at the influence of dopants on the oxygen-interstitial and oxygen-vacancy diffusivity. To this end we introduce Nb 5+ or Ta 5+ as donor dopants; or we introduce Gd 3+ as an acceptor dopant; and we compare the results with those obtained for the dopant-free cells.

Computational section
The system is described by short-range interactions with a Buckingham potential and long-range coulombic interactions. Full ionic charges are assumed. The set of pair-potential parameters derived by Sayle et al [44] for CeO 2 have already been used to simulate structures and properties of extended defects [44,49], activation barriers of oxygen-vacancy migration in pristine [30] and in strained [32] CeO 2 and non-linear oxygen-vacancy mobility in CeO 2 [35]. The cation-anion parameters used in this study were all derived for Catlow's oxygenoxygen interaction [45] and should therefore be compatible with one another. The inter-cationic interactions, with the exception of coulombic interactions, were assumed to be zero. The set of potential parameters used in our simulations is given in table 1.
In this study, simulations of bulk CeO 2 employed a cell of dimensions´á a a 14 14 14 (where a is the lattice parameter of cubic CeO 2 ), containing 32928 ions, i.e. Ce 10976 O 21952 . We performed two different sets of simulations. In the first set, we randomly introduced vacancies or interstitials throughout the cell and compensated the charge with a uniform background charge. We refer to these simulation cells as bgcc (background charge-compensated); they allow us to study the diffusion behaviour of the two oxygen defects, free from the interactions with compensating defects [29,[50][51][52]. The second set of simulations were performed on cells containing dopant ions to charge compensate the respective oxygen point-defects; donor dopants (Nb Ce • or Ta Ce • ) for the cells containing oxygen interstitials (  O i ) and acceptor dopants ( ¢ Gd Ce ) for the cells with oxygen vacancies (V O •• ). The dopants were randomly distributed over the possible cation sites replacing regularĆe Ce . Between 20 to 1000 oxygen interstitials or oxygen vacancies were introduced into the simulation cell: this corresponds to defect site fractions between 0.18% to 9.11% for the interstitials and 0.09% to 4.56% for the vacancies.
We utilized the LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [53,54] code for our MD simulations. In all simulations a NpT ensemble was used to describe the system. In this ensemble the temperature T, the pressure p and the particle number N, of the system are fixed, while its energy E and volume V are allowed to vary. A Nosé-Hoover thermostat (barostat), as implemented in LAMMPS [54] controlled the system's temperature (pressure). For each calculation, the simulation cell was first brought into an equilibrium state. Equilibrium was deemed to have been reached when U and V were observed to only fluctuate around constant values. The time required to attain equilibrium varied between 150 ps and 500 ps depending on the number of oxygen defects. As a rule of thumb: the more time was needed, the smaller the number of migrating particles. After equilibration, the simulation was run for 1.4 ns to produce transport data. The software package OVITO [55] was used to visualize the simulation results.

Oxygen-interstitial and oxygen-vacancy in CeO 2
The same three-step procedure was used to determine the defect diffusivity, regardless of whether the results were obtained from oxygen-deficient or oxygen-excess simulation cells. First, we computed from the appropriate MD simulation the mean-squared-displacement of the oxygen ions, á ñ r O 2 , as a function of time t. Such data are shown in figure 1 for a defect site fraction of 0.91%, for the two different defect species at various temperatures. In both cases, á ñ r O 2 increases linearly with time, indicating that diffusion took place. Second, we applied the standard Einstein relation, to these data to obtain the tracer diffusion coefficients of oxygen, D O *. It is worth noting that even without calculating exact values of D O * one recognises that, at T=2000 K for example, á ñ r O 2 is much higher for oxygendeficient than for oxygen-excess cells. This is a qualitative indication that vacancies are more mobile than interstitials. Third, knowing the defect site fraction n i/v and the tracer correlation factor f i v * for the two oxygen defects, we extracted the respective defect diffusivity, f * depends on the migration mechanism and the symmetry of the sublattice on which diffusion takes place. For the fluorite structure, migration of anion interstitials can either occur by an interstitial mechanism or by a noncolinear interstitialcy mechanism. All our MD simulations confirmed that the trajectories of the anions were consistent with the second possibility. For the interstitials, we thus have 58,59]. Defect site fractions of oxygen interstitials and oxygen vacancies ranged, as noted previously, from ca.0.1% to 1.8%; the obtained defect-diffusivity data are shown in figure 2.
Comparing D v with D i , for a given temperature T and defect site fraction n i/v , we see that the vacancy diffusivity is always higher than interstitial diffusivity. For low n i/v both isothermal values of D v and D i are constant, albeit with different values. With increasing n i/v the defect diffusivities start to deviate from the values obtained for low n i/v . D v decreases with increasing n v , whereas D i increases with increasing n i . Since no other defects are present in these bgcc simulation cells, these deviations can be attributed to defect-defect selfinteractions. For both defect species the deviations from the regions of constant diffusivity values begin at roughly , oxygen point-defects should not be treated as a dilute solution, as far as diffusion is concerned, but rather as a concentrated solution.
Assuming, for simplicity, that the defects execute a random walk in three dimensions, one can express the defect diffusivity in terms of the jump distance d ; i v the number of possible lattice sites to which the defect can jump, Z ; i v the attempt frequency ν 0 ; and the activation entropy ΔS mig and activation enthalpy ΔH mig of migration [13,60,61].
where a(T) is the lattice parameter taken from the MD simulation at temperature T. For vacancy migration, Z v =6 and d v =a(T)/2. Assuming a common ν 0 value (of 10 13 s −1 ) for both interstitialcy and vacancy migration, we find, by applying equation (3) to the data in

Oxygen-interstitial diffusion in donor-doped CeO 2
In figure 4 we compare the results for our Nb Ce • -doped cell with the undoped cell and two differently doped cells. The first case, figure 4(a), compares the values for D i of the bgcc cell (with n i =0.45%), with values from a cell containing the same site fraction of oxygen interstitials, but also containing Nb Ce • , such that the charge of the interstitials is exactly compensated by the charge of the dopant. There is a significant difference between the two cases, with D i in the Nb Ce • -doped cell being higher and exhibiting a smaller activation enthalpy of interstitial migration, ΔH mig,i =(1.22±0.04) eV, compared with ΔH mig,i =(1.42±0.03) eV for the undoped cell. These results clearly show that Nb Ce • has a beneficial effect on the diffusion of oxygen interstitials. In the second case ( figure 4(b)) both cells contained a site fraction of 0.9% Nb Ce • ions, but in the one case the dopants were charge-compensated only by oxygen interstitials and in the second case they were compensated mostly by ¢ Ce Ce with only a few  O i . The reason for this comparison is that, according to defect chemical modelling [62], Nb Ce • as a donor dopant is compensated primarily by ¢ Ce Ce for a wide range of conditions. Here, both sets of MD simulations indicate similar interstitial diffusion coefficients. The data for the second set display more scatter, and this is due to the low number of oxygen interstitials. The activation enthalpy of interstitial migration for this case is ΔH mig =(1.08±0.12) eV, with the large error coming from the poor statistics, but it is comparable to the data for the cell containing only  O i to compensate ) . As the ¢ Ce Ce ions do not seem to influence the diffusion of oxygen interstitials, we will not consider ¢ Ce Ce in the following simulations and therefore only focus on donor dopants and oxygen interstitials.
In the third case (figure 4(c)) the results for the two different donor dopants (Nb Ce • and Ta Ce • ) are compared. Both simulation sets show very close agreement for the absolute values for D i , and the activation enthalpies of interstitial migration are similar, too: ΔH mig =(1.23±0.04) eV for the Nb-doped cell and ΔH mig =(1.22±0.03) eV for the Ta-doped one. Evidently, it does not matter which of these two pentavalent dopants is used, as they both yield essentially the same results, at least according to these simulations.
Given the positive effect of Nb on the interstitial diffusivity, we performed MD simulations of oxygeninterstitial diffusion as a function of Nb content, with the donors' charge being compensated exactly by the interstitials' charge. For each concentration three simulations were conducted to enhance the quality of the statistics. The results are plotted in figure 5(a). There is a clear trend of D i increasing and ΔH mig,i decreasing with increasing Nb concentration. Defect-defect interactions between oxygen interstitials can be excluded as the sole reason for this behaviour, as the deviation in the bgcc cells was smaller (see figure 2(a)) and the activation enthalpy of interstitial migration was unaffected for the dilute solution case, opposed to the behaviour shown in figure 5(a). Figure 5(b) shows the activation enthalpy of interstitial diffusion obtained from the simulations plotted against the mean Nb-Nb spacing. With a decrease in the mean Nb-Nb spacing the migration enthalpy

Oxygen-vacancy diffusion in acceptor-doped CeO 2
For the sake of completeness, we also examine the interactions between dopant cations and mobile defects in cells with oxygen vacancies. It can be seen in figure 7 that the presence of the acceptor-dopant ¢ Gd Ce leads to a decrease in the oxygen-vacancy diffusivity and to an increase in the activation enthalpy of vacancy migration, In figure 8(a) values of D v for various oxygen-vacancy site fractions are plotted against inverse temperature. In comparison with results for oxygen-vacancy diffusion in the bgcc simulation cell (see figure 2(b)), the presence of Gd 3+ leads to an overall decrease in diffusivity by a factor of roughly 2, with no further modification beyond the already observed one due to oxygen vacancy-vacancy interactions. The ¢ Gd Ce site fraction does not seem to have an effect on the activation enthalpy of vacancy migration, beyond the increase of DH mig,v in comparison to the obtained values for the undoped cell. For all ¢ Gd Ce site fractions looked at, the obtained values for DH mig,v seem to fall around (0.82±0.05) eV, in contrast to activation enthalpies of migration reported by Faber et al (green circles), which show a slight decrease with increase in ¢ Gd Ce . A possible reason for the difference between the reported values by Faber et al and our results is that Gd was randomly distributed in our simulation cells, whereas in the samples of Faber et al, Gd was mobile at the sintering temperature, and thus able to form clusters and nanophases, as predicted by Žguns et al [63].
In figure 9 the oxygen trajectories for the undoped and ¢ Gd Ce -doped cell, with the same n v are shown. For the undoped cell (shown in figure 9(a)) the oxygen ions have moved homogeneously through the cell, as indicated  by the homogeneous distribution of the trajectories. In figure 9(b) the oxygen trajectories for the ¢ Gd Ce -doped cell are shown. The trajectory distribution is less homogeneous than for the undoped cell, indicating an interaction between the dopant and the vacancies. In the proximity of the dopant ions, the trajectory distribution is more dense, indicating more oxygen jumps. The effect is similar to the observation for the Nb Ce • -doped cell, albeit not as pronounced, with the fundamental difference that in the acceptor-doped cells the overall oxygen-vacancy diffusivity is reduced in comparison to the undoped cells. The oxygen-vacancies seem to prefer to reside in the proximity of ¢ Gd Ce , similar to the behaviour of the oxygen-interstitials in the donor-doped case. But in contrast to the behaviour of oxygen-interstitials, here the interaction between the dopant and the  oxygen vacancy leads to a higher migration barrier and through this to a decrease in the oxygen vacancy diffusivity. The oxygen vacancies are less mobile in the proximity of the dopant ions.

Literature discussion
In figure 10 we compare our computational results with experimental and computational results for dilute solutions of oxygen defects and dopants in CeO 2 and other AO 2 systems. We take our results for doped cells, since the experimental systems are doped at similar levels. In table 2 we give a comprehensive list of experimentally obtained and theoretically predicted activation enthalpies of migration for both dilute and concentrated systems.
In general, we see in figure 10 that the diffusion coefficients for oxygen vacancies are higher than those for oxygen interstitials for all fluorite-structured oxides. The temperature dependence of the defect diffusivity is stronger for oxygen interstitials than for oxygen vacancies, that is, the activation enthalpy of migration is higher for oxygen interstitials than for oxygen vacancies.
Comparing our predicted D i values with reported values for oxygen-interstitial diffusion in CeO 2 and two other AO 2 oxides (UO 2 , PuO 2 ), we find all of them to be of comparable order of magnitude. In particular, we find a good agreement between the experimentally found values for Ce 0 Ta: 0.75±0.08 exp. The activation enthalpy of interstitial migration DH mig,i in AO 2 oxides lies in the range of ΔH mig,i =(0.9 to to 1.3) eV [25,39,[65][66][67]  for CaF 2 . Similar to the oxygen migration in AO 2 the activation enthalpy of vacancy migration is lower than the activation enthalpy of interstitial migration.
The reason for such behaviour, we propose, is rooted in the crystal structure of AX 2 fluorite-structured materials. The vacancy mechanism does not require the migrating ion to come into direct contact with other anions, whereas the interstitialcy mechanism specifically requires this. Comparing the activation enthalpies of vacancy migration for AO 2 and ¢ A F 2 materials, we recognise that the fluorine diffusion in ¢ A F 2 has a lower migration enthalpy than oxygen diffusion in AO 2 , and the same holds true for the diffusion of interstitials. The reason for this could be that the migration of single charged F − is less hindered than the migration of the double charged -O 2 . Genreith-Schriever et al [28] examined using computational methods the migration behaviour of various anions in CeO 2 and found that a higher negative charge on the migrating ion leads to a higher migration barrier. This would support our explanation for the observed difference between the F − and the -O 2 migration.

Conclusions
(i) At a given temperature and defect site fraction, oxygen vacancies were found to exhibit a higher diffusivity than oxygen interstitials in Ce , the activation enthalpy of vacancy migration was substantially lower than that of interstitial migration, D < D H H mig,v mig,i . (iii) Defect-defect self-interactions are different in nature for vacancies and interstitials. Going beyond the dilute limit leads to D i increasing but D v decreasing. This is accompanied by DH mig,i decreasing, but DH mig,v increasing.
(iv) The fundamental difference between oxygen-vacancy and oxygen-interstitial migration in fluortitestructured materials is attributed to the crystal structure.