Structural dynamics of Schottky and Frenkel defects in CeO2: a density-functional theory study

Cerium dioxide CeO2 (ceria) is an important material in catalysis and energy applications. The intrinsic Frenkel and Schottky defects can impact a wide range of material properties including the oxygen storage capacity, the redox cycle, and the ionic and thermal transport. Here, we study the impact of Frenkel and Schottky defects on the structural dynamics and thermal properties of ceria using density functional theory. The phonon contributions to the free energy are found to reduce the defect formation free energies at elevated temperature. The phonon dispersions of defective CeO2 show significant broadening of the main branches compared to stoichiometric ceria. Phonon modes associated with the defects are identifiable in the infrared spectra through characteristic shoulders on the main features of the stoichiometric fluorite structure. Finally, the presence of Frenkel and Schottky defects are also found to reduce the thermal conductivity by up to 88% compared to stoichiometric CeO2.


Introduction
Cerium dioxide (CeO 2 or ceria) is a lanthanide oxide that has been studied for decades due to its use in a wide range of industrial applications. CeO 2 is employed as an oxygen storage material in catalysis [1,2], a nanozyme in biomedical applications [3], a solid electrolyte in solid oxide fuel cells [4][5][6][7], and as a surrogate material for actinide oxide nuclear fuels [8][9][10][11].
CeO 2 is often used as a simulant material for PuO 2 , as both Ce and Pu can switch between 3+ and 4+ oxidation states and share the cubic CaF 2 fluorite structure (Fm3m), which CeO 2 retains until melting at 2753 K [12,13]. Experimentally, the toxicity and radioactivity of PuO 2 make using CeO 2 as a simulant material attractive. Even for computational studies PuO 2 requires computationally-demanding spin-orbit coupling and non-collinear magnetism to describe its longitudinal 3k antiferromagnetic ground-state [14], making non-magnetic [15] CeO 2 an ideal proxy.
The performance of nuclear fuels is closely linked to their thermal conductivity. The difference in thermal conductivity between CeO 2 , PuO 2 and UO 2 is a relatively small 1-2 W m −1 K −1 across a wide temperature range of ∼300-2000 K [9,10,16]. The thermal conductivities of mixed oxide (MOX) fuels such as U 0.92 Pu 0.08 O 2 (Pu-MOX) and U 0.92 Ce 0.08 O 2 (Ce-MOX) are also very similar above 800 K and differ by ∼1 W m −1 K −1 at 400 K [10].
The presence of defects arising from the extreme operating conditions in nuclear reactors and other high-temperature applications impact the thermophysical properties, including the thermal conductivity. Classical molecular dynamics simulations found that both anion and cation vacancies significantly reduce the thermal conductivity of UO 2 and PuO 2 by 16.5%-25.2%, even at low concentrations of ∼0.1%, with cation vacancies having a larger impact as they provide a greater disruption to the structural dynamics [17]. In hyper-stoichiometric structures, oxygen ions occupy the fluorite octahedral interstitial sites, resulting in additional phonon scattering and a reduction in thermal conductivity of UO 2+x , with laser-flash measurements at 873 K identifying a reduction of ∼1 W m −1 K −1 when the oxygen content of UO 2+x increases from x = 2.002 to x = 2.033 [18].
At 900 K and p(O 2 ) = 10 −10 atm, the concentrations of intrinsic defects such as Frenkel and Schottky defects are lower than those of point defects such as oxygen vacancies by around 20 orders of magnitude [19]. However, it is feasible that large, non-equilibrium local concentrations could build up due to the energy release associated with nuclear events in working reactors, and hence the impact of complex defects on the thermal conductivity is of great interest. Whilst the defect formation energies of Frenkel and Schottky defects have been calculated, their impact on the thermal conductivity has yet to be reported. In the similar fluorite-structured ThO 2 , Frenkel and Schottky defects were shown to broaden the phonon dispersion and reduce the thermal conductivity by up to 90% [20]. In insulators such as ceria, the low frequency acoustic phonons (<12 THz), which have the largest group velocities and longest lifetimes, dictate the magnitude of the heat transport [8,21]. Assessing the impact of the defects on this part of the phonon dispersion is therefore important, but data is available only for individual point defects in CeO 2 [22], and this study did not perform band unfolding nor simulate other routine phonon spectra to facilitate better comparison to experimental data.
In this work, we use density-functional theory (DFT) to establish the impact of the intrinsic Frenkel and Schottky defects in CeO 2 on the structural dynamics. We compute reference phonon spectra, including identifying spectral signatures in the infrared (IR) spectra that would allow the defects to be identified and quantified, we quantify the effect of the phonon free energy on the defect formation energies, and we also estimate the impact of the defects on the thermal conductivity.

Computational methodology
A 96-atom supercell of stoichiometric CeO 2 was built from a 2 × 2 × 2 expansion of the 12-atom conventional cell, from which anion Frenkel, cation Frenkel and Schottky defect models were generated with defect concentrations of 3.1%.
DFT simulations were performed using the Perdew-Burke-Ernzerhof generalised-gradient approximation (GGA) exchange-correlation functional revised for solids (PBEsol) [23] within the Vienna Ab initio Simulation Package (VASP) code [24]. Spin-polarized calculations were performed, as are standard in calculations on ceria. However, we do not see any spin polarization in our calculations, as the defective systems we study are charge-neutral and comprised of closed-shell O 2− and Ce 4+ ions. A Hubbard U parameter of U eff = 5 eV was applied to the Ce 4 f states to account for self-interaction error [15] as suggested in previous literature [25]. However, the defects we examine do not produce Ce 3+ ions and we do not have issues with localization. A plane-wave basis set with a cut-off energy of 500 eV, which gives total energies converged to 4 meV per CeO 2 formula unit compared to a 20% larger cutoff, was used with projector augmented wave pseudopotentials [26,27] with frozen cores of [Kr] and [He] on the Ce and O atoms, respectively, to model the valence electronic structure. A Γ-centred Monkhorst-Pack k-point grid [28] with 4 × 4 × 4 subdivisions was used to sample the Brillouin zone (BZ) of the stoichiometric and defective systems. The electronic and ionic tolerances were set to 10 −9 eV and 10 −3 eV Å −1 respectively. The stoichiometric supercell was minimised under constant pressure, while the defective structures were minimised with the volume and cell shape fixed to those of the stoichiometric model. In general, the intrinsic defect concentrations are expected to be low and fixing the volume to that of the host lattice is appropriate. Furthermore, the band-unfolding algorithm works better if the reference and distorted cells have similar lattice vectors. Relaxing the structures at constant pressure resulted in a small increase in the unit-cell volume of <0.5% for the majority of the structures and ∼2% for the cation Frenkel defect model. The phonon frequencies usually decrease with volume, leading to lower free energies and lattice thermal conductivities. However, across all the defects we would only expect the volume change in the cation Frenkel model to have a material impact on the properties we calculate, and since this is predicted to be the least stable defect this is not too concerning.
The Phonopy code [29] was utilised to calculate 2nd-order harmonic interatomic force constants (IFCs) using the finite-displacement method, with VASP employed as the force calculator. The phonon dispersion of stoichiometric CeO 2 was determined using a transformation matrix to the primitive cell. The band unfolding methodology [30,31] in Phonopy was used to obtain the phonon dispersions of the defective models projected onto the primitive stoichiometric cell for comparison. This procedure was previously used successfully for ThO 2 [20]. The density functional perturbation theory (DFPT) [32] routine in the VASP code was employed to calculate the Born effective charges and high-frequency dielectric constants required to apply a non-analytical correction to the phonon frequencies. The atom-projected phonon density of states (PDOS) curves for the stoichiometric supercell and defective models were calculated on 10 × 10 × 10 Γ-centred q-point grids. To obtain IR spectra of the models the Phonopy-Spectroscopy code [33], which implements the methodology described in [34], was used.
The Phono3py code [35] was employed to calculate the 3rd-order IFCs, which were used with the single-mode relaxation-time approximation (RTA) method to calculate the lattice thermal conductivity of stoichiometric CeO 2 . Convergence testing of the phonon group velocities, phonon lifetimes and thermal conductivity identified a 32 × 32 × 32 sampling mesh to be sufficient for computing the modal properties for the primitive cell. The thermal conductivity of the defective models was estimated using a q-point mesh with 10 × 10 × 10 subdivisions and one of two approximate methods [36]. The two-phonon weighted joint density of states (w-JDOS) functions were also generated for the stoichiometric and defective models using the same 10 × 10 × 10 mesh. A larger q-point mesh was not feasible due to memory restrictions, which may result in small errors in the calculated properties.
The Frenkel defect is an interstitial-vacancy pair, where an atom migrates from a lattice site to an interstitial site. Formation of the anion (O) and cation (Ce) Frenkel defects (figures 1(B) and (C)) are defined by equations (1) and (2) respectively: Two configurations of anion Frenkel defects were optimised with the O ′ ′ i and V •• O pair aligned along the ⟨111⟩ and ⟨133⟩ directions, which we denote OF⟨111⟩ and OF⟨133⟩ respectively. The distance between the interstitial-vacancy pairs in the OF⟨111⟩ and OF⟨133⟩ models are 7.01 and 5.80 Å. The OF⟨133⟩ structure was found to be more stable than the OF⟨111⟩, with formation energies per point defect of 1.96 eV and 2.01 eV respectively, so we discuss the OF⟨133⟩ in this work. The cation Frenkel defect was created with the Ce •••• i and V ′ ′ ′ ′ Ce pair aligned along the ⟨111⟩ direction with an interstitial-vacancy pair distance of 4.36 Å and is referred to as CeF⟨111⟩. Based on previous literature, the interstitial-vacancy pairs in both types of defect were placed as 2 nd nearest neighbours to prevent spontaneous recombination during optimisation [20,42].
It should be noted that in the literature some authors refer to the anion-Frenkel defect in equation (1) as an 'anti-Frenkel' , but the equations they present describe the formation of the equivalent oxygen interstitial-vacancy pair [43][44][45].
To form the Schottky defect, a complete formula unit is removed from the structure leaving one V ′ ′ ′ ′ Ce and two V •• O as described in equation (3): Three configurations of the Schottky defect were considered with the two V •• O aligned along the ⟨111⟩, ⟨110⟩ and ⟨100⟩ directions with respect to the central V ′ ′ ′ ′ Ce , which we denote S⟨111⟩, S⟨110⟩ and S⟨100⟩ respectively ( figure 1(D)). Based on previous literature, the vacancies were positioned as first nearest-neighbours [20]. The distance between the V ′ ′ ′ ′ Ce and V •• O is 2.45 Å in each Schottky defect and the distance between the 2 V •• O are 4.90, 4.00 and 2.83 Å in the S⟨111⟩, S⟨110⟩ and S⟨100⟩ models, respectively. The Helmholtz free energies of formation ∆A f for the Frenkel and Schottky defects were determined according to equations (4) and (5) respectively: where A Stoich , A Frenkel , A Schottky and A CeO2 are respectively the free energies of the stoichiometric supercell, the supercell models containing Frenkel and Schottky defects and a single CeO 2 formula unit. The prefactors of  (4) and (5) using the phonon contribution to free energy A vib (T) as described in equation (6): where k B is the Boltzmann constant, Z vib (T) is the vibrational partition function, the sums are taken over phonon modes qj with wavevectors q, band indices j and frequencies ω qj , and N is the number of q included in the summation. The A vib (T) were obtained using the ω qj at wavevectors q across the full BZ. The internal energy U vib (T) and entropy S vib (T) are calculated according to equations (7) and (8) respectively: Summing the lattice internal energy U 0 , equated to the DFT total energy, and the A vib (T) gives the total Helmholtz free energy of the system A (T) as in equation (9): Table 1 compares our calculated defect formation energies with other modelling studies using DFT [19,43,44,46] and potential models (PMs) [42,47]. Variations in the literature values for the formation energies of Frenkel and Schottky defects could be due to how they are calculated, either by summing the energies of the individual point defects [43,44,47], or by considering a model with defects in a cluster, which also accounts for the binding energies between point defects [19,42,46]. For the Schottky defect, theoretical values for the defect formation energies calculated using DFT range from 1.19-1.64 eV per point defect using a cluster model, whereas summing the energies of individual point defects predicts formation energies in the range of 1.20-2.29 eV per point defect [19,43,44,46]. We note that formation energies can be reported as either energies per point defect [43,47] or energies per disorder [19,42,44,46]. We adopt the former Table 1. Calculated athermal (0 K) defect formation energies (eV per point defect) for the anion Frenkel, cation Frenkel and Schottky defects examined in this work compared to other theoretical literature. Defect formation energies marked by asterisks were originally reported as energy per disorder and converted to energy per point defect by dividing by two (Frenkel) or three (Schottky). Energies marked by daggers were calculated by summing the energies of individual point defects. The size of the supercell and method used to calculate the energies are given in the top row (PM-potential model). OF∞, CeF∞ and S∞ denote formation energies for the Frenkel and Schottky defects in the limit of infinite dilution.
Considering the energies per point defect yields a stability ordering of Schottky > anion Frenkel > cation Frenkel, with energies per point defect of 1.44-1.64, 1.96 and 5.51 eV, respectively. The substantially larger formation energy of the cation Frenkel suggests it is likely to be present at insignificant concentrations, and can be attributed to the larger disruption to structure that this defect causes. There is disagreement within the theoretical literature as to whether the Schottky defect is the most stable, as found in the present work [19,42,44], or the anion Frenkel is the most stable [43,47]. When considering energies per disorder, we find a different stability order of anion Frenkel > Schottky > cation Frenkel, with energies per disorder of 3.92, 4.33-4.92 and 11.02 eV, respectively. This is mainly due to the fact that a Schottky defect creates more point defects than a Frenkel defect, and the change in stability ordering is also reported in [42]. However, no matter how one defines the energy, i.e. either per point defect or per disorder, there are still discrepancies in the literature. As an example, Zacherle et al [43] found the anion Frenkel defect to be consistently the most stable, whereas Keating et al [19] and Nakayama and Martin [44] found the Schottky defect to be most stable. Differences in methodology apparently do not account for this discrepancy, as Zacherle et al [43] and Keating et al [19] both used PBE + U with U = 5 eV and accounted for the chemical potential of the species, yet obtained different outcomes.
For the actinide oxides ThO 2 , UO 2 , NpO 2 , PuO 2 , AmO 2 and CmO 2 , the order of stability in terms of energies per point defect follows the trend Schottky > anion Frenkel > cation Frenkel [42]. This is consistent with the present study.
The ∆A f of the anion Frenkel, cation Frenkel and Schottky defects are plotted as a function of temperature from 0-2750 K in figure 2. As temperature increases, a reduction in ∆A f is observed in all cases, suggesting that defect formation is more energetically favourable at higher temperatures. The formation energies per point defect at 2750 K, close to the melting point of 2753 K [12], are calculated to be 1.47 eV for the OF⟨133⟩, 5.48 eV for the CeF⟨111⟩, 1.20 eV for the S⟨111⟩, 1.18 eV for the S⟨110⟩, and 1.40 eV for the S⟨100⟩.
These values are between 1%-24% smaller than the values calculated at 0 K. Whilst temperature has a significant impact on the magnitude of ∆A f it does not change the order of stability over the temperature range examined.
The distance between the point defects in the Frenkel and Schottky models has a significant impact, with greater distances producing a smaller ∆A f . For CeO 2 , Cooper et al [42] found oxygen vacancies placed as 3 rd nearest neighbours as the part of a Schottky defect gave a 17% smaller ∆A f than when the vacancies were placed as (1 st ) nearest neighbours. This trend is also seen for the Frenkel defects [42], such that when the interstitial-vacancy pair are placed as 2 nd nearest neighbours the ∆A f are 7% smaller than when placed as 1 st nearest neighbours.  Stratton and Tuller [45] determined a defect formation energy of 1.83 eV per point defect for the anion Frenkel defect in CeO 2 (0.1% U-doped) based on oxygen self-diffusion measurements at high temperature. While this is slightly smaller than our computed value of ∆A f = 1.96 eV per point defect for the OF⟨133⟩, our value is in agreement with the theoretical literature where the ∆A f range from 1.95-3.21 eV per point defect [19,[42][43][44][46][47][48]. As this experimental study was carried out at high temperature, the activation of the imaginary mode on the O sublattice may affect the derived energy [40], which is not accounted for in these calculations.
We predict the cation Frenkel defect to be much less favourable than the anion Frenkel, with a formation energy of 5.51 eV per point defect. Theoretical values for this defect range from 3.68-15.94 eV per point defect [42-44, 47, 48]. Using a PM method, Cooper et al [42] calculated a formation energy of 6.28 eV per point defect for a cation Frenkel defect with the interstitial-vacancy pair aligned along the ⟨111⟩ direction, which is in good agreement with our prediction.
The defect concentrations N D as a function of temperature were estimated using a Boltzmann distribution, as shown in equation (10), and are plotted in figure 3: However, if we instead use energies per disorder the anion Frenkel disorder reaches the baseline concentration at a lower temperature of 890 K compared to both the Schottky disorder (1000-1140 K) and the cation Frenkel disorder (2630 K). Keating et al [19] predicted the Schottky disorder to be present at a concentration of 10 3 cm −3 at 900 K, some two orders of magnitude higher than the anion Frenkel, based on energies per disorder and an oxygen partial pressure p(O 2 ) = 10 −10 atm. At a higher p(O 2 ) of 0.2 atm, the concentrations of Frenkel and Schottky disorders were similar to those calculated at lower p(O 2 ), i.e. the oxygen partial pressure is not expected to have a large impact on the concentrations. To compare with Keating et al [19], we predicted defect concentrations using energies per disorder at 900 K and obtained concentrations of 5.21 cm −3 for the anion Frenkel defect, between 4 × 10 −6 and 1 × 10 −2 cm −3 for the Schottky defect, and a negligible concentration of the cation Frenkel. The predicted concentration of the anion Frenkel is comparable to that from Keating et al [19] but their predicted value of the Schottky disorders is orders of magnitude higher than ours. We attributed this difference to the theoretical methodology employed in this study, leading to a smaller Schottky defect formation energy (cf table 1), together with the exponential dependence of the defect concentration on the formation energy.

Structural dynamics
The phonon dispersion of stoichiometric CeO 2 , calculated along the Γ-X-U and K-Γ-L paths through the BZ is shown together with the atom-projected phonon density of states in figure 4 and compared with inelastic neutron scattering (INS) data collected at 293 K by Clausen et al [49]. There is excellent agreement between our calculated phonon dispersion and the experimental data [49], and also with the theoretical results of Buckeridge et al [40] obtained using a similar method. At each wavevector q there are nine branches composed of three lower-frequency acoustic phonons and six higher-frequency optic phonons. The motions of the heavier Ce sublattice are the main contributor to the acoustic phonons, which split into longitudinal and transverse acoustic (LA/TA) branches. The optic phonons are comprised mainly of motions of the lighter O sublattice, and split into longitudinal and transverse optic (LO/TO) modes. The PDOS shows that Ce motion dominates at frequencies below 6 THz, whereas O motions dominate above 6 THz, in agreement with calculations performed by Malakkal et al [50]. There is a significant overlap between the acoustic and optic modes, suggesting the possibility of interactions that could limit the thermal conductivity. A comparison of our calculated frequencies to literature data is shown in table 2. Our calculations agree well with experimental measurements from Raman [51,52], IR [53,54] and optical reflectivity [55], as well as other modelling using DFT [40,56,57] and PM methods [47]. The largest differences between our calculations and experiments is in the F 1U LO mode, for which the experimental value of 17.5-17.9 THz [51-55, 58] is around 1.2 THz higher than our calculated frequency.
For this particular mode, calculations using the HSE06 hybrid functional [40] provide better agreement with experimental results, albeit at a higher computational cost.
While, the phonon dispersion of stoichiometric CeO 2 has been widely studied [19, 21, 40, 48-51, 57, 59, 60], the impact of the Frenkel and Schottky defects on the phonon spectrum has yet to be examined. Using the band unfolding methodology in Phonopy [30,31], we determined the phonon dispersions of the Frenkel and Schottky defect models with respect to the primitive cell, which are shown in figure 5 along with the Table 2. Predicted phonon frequencies (THz) of the Γ, X and L modes in CeO2 compared to experimental studies using inelastic neutron scattering (INS), infrared (IR), Raman and optical reflectivity, and to computational studies using potential models (PMs), and density-functional theory (DFT).    corresponding PDOS. In all the defective phonon dispersions, broadening is observed across the whole spectrum and is particularly prominent in the LO and TO modes around 14 and 10 THz, respectively, in the vicinity of Γ, and in the LA mode around ∼5 THz in the vicinity of X. We present an estimate of the range of frequencies at each high symmetry q point in the dispersion in table 3. A comparison of the PDOS of the stoichiometric and defective systems is given in figure S8.

IR spectra
The Γ-point phonon modes can be probed using IR and Raman spectroscopy. Pristine CeO 2 has a triply-degenerate F 2g Raman active mode at 13.5 THz, a doubly-degenerate F 1u IR active mode at 8.8 THz and a singly-degenerate F 1u IR active mode at 16.5 THz. We compare the IR spectra of stoichiometric ceria with those of the defective structures in figure 6. We have applied a nominal linewidth of 0.52 THz (17.5 cm −1 ) based on the full-width half-maximum of 0.51 THz derived for the TO bands in CeO 2 using hyper-Raman spectroscopy [58]. This linewidth is comparable to the value of 0.75 THz (25.0 cm −1 ) derived from room-temperature IR reflectivity measurements of the similar fluorite-structured ThO 2 [61]. We have also computed the linewidths of the F 1u and F 2g modes to be 0.037 and 0.105 THz, respectively, at 300 K, based on the phonon lifetimes computed as part of the thermal-conductivity calculations described below. Additional IR active modes are present in the spectra of the defective models and are associated with vibrations localised around the defect sites. The defects break the symmetry of the stoichiometric structure and result in a portion of the zone-boundary modes being folded to the Γ point, leading to a larger number of spectroscopically-active modes.
The spectra of the anion Frenkel and three Schottky defect models are similar to that of stoichiometric ceria, with the main difference being small shoulder features and less symmetric peaks. The cation Frenkel defect on the other hand introduces a pronounced splitting of both features in the IR spectrum of the stoichiometric material. Based on the simulated spectra, it may be possible to identify the presence of anion Frenkel and Schottky defects from the line-shapes and/or the presence of small shoulder features, but it would be difficult to distinguish between them, whereas the cation Frenkel defect produces a clear spectral signature and should be easily identifiable from IR measurements.

Thermal conductivity
The macroscopic thermal conductivity κ latt (T) was determined using the single-mode relaxation time approximation (SM-RTA) as described in equation (11) [35]: where the sums run over the individual phonon modes λ with modal heat capacities C λ computed from equation (12), group velocities v λ calculated using equation (13) and lifetimes τ λ from equation (14), N is the number of wavevectors q included in the summation and V is the volume of the unit cell.
where ω λ are the phonon frequencies, W λ are the corresponding eigenvectors, D (q) is the dynamical matrix at the wavevector q, and Γ λ are the phonon linewidths calculated from the imaginary part of the phonon self-energy using third-order perturbation theory [35]. A scalar average of the thermal conductivity can be computed from equation (15) by averaging the diagonal elements κ xx , κ yy , and κ zz of the κ latt tensor, which we denote κ latt .
where we have omitted the explicit temperature dependence for brevity. Figure 7 compares the calculated κ latt obtained for stoichiometric CeO 2 between 300-1500 K with other literature reports. Good agreement is obtained across the temperature range examined, in particular with equilibrium molecular dynamics (EMD) simulations and EMD using the Green-Kubo method (EMD-GK) [50]. The experimental values were obtained from laser flash analysis (LFA) on CeO 2 samples with different stoichiometry, crystallinity and density, which makes comparison somewhat difficult [9,10,16,50,62,63]. Some of the experimental values obtained using LFA are lower than our values. Our calculations are based on a single crystal of CeO 2 , whereas the experimental samples were polycrystalline and would likely contain extended defects such as grain boundaries that reduce the thermal conductivity. Furthermore, the theoretical densities of the experimental samples were less than 100%, and an increase from 95% to 100% was found to increase the thermal conductivity by 5% at 773 K [9,16]. Our calculations do however produce comparable results to those obtained using computational methods based on PMs [50,62] and DFT [21,50,59]. At 700 K we compute a κ latt of 8.4 W m −1 K −1 for stoichiometric CeO 2 , which is in good agreement with the 8.8 W m −1 K −1 calculated using EMD [50], while the experimental κ latt of 7.8 W m −1 K −1 reported at 675 K is 7% smaller than our prediction [62]. Theoretical studies based on the athermal equilibrium structure may differ from experimental measurements at finite temperature where the structure has undergone thermal expansion, which can both decrease the group velocities and enhance phonon scattering. Additionally, our calculations do not consider higher-order (e.g. four-phonon) and phonon-defect scattering, both of which would occur in the real material.
A simple boundary-scattering model was used to predict the thermal conductivity of stoichiometric CeO 2 as a function of crystal grain size at selected temperatures (figure 8), from which we find that the thermal conductivity falls sharply for grain sizes below ∼5 µm. A comparison of the thermal conductivity as a function of crystal grain size for the stoichiometric and defective CeO 2 models is provided in figure S9, supporting information.  Calculating the phonon lifetimes and thermal conductivity of the defective systems is computationally unfeasible due to the resources required to evaluate the 3 rd -order force constants. Instead, we therefore utilise the constant-relaxation time approximation (CRTA) model [36,64] in which the τ λ are replaced with a constant τ CRTA as shown in equation (16).
It should be noted that both τ CRTA and C λ are temperature dependant, and we have omitted the explicit temperature dependence for clarity. The summation in the right-hand side expression in equation (16) is calculable within the harmonic approximation, and comparison between systems allows the impact of differences in the phonon spectra and group velocities to be quantified. We carefully tested the convergence of the harmonic CRTA functions as a function of the q-point sampling for the stoichiometric and defective models (see figure S1). Comparisons of the CRTA decomposition for the primitive cell and supercell of stoichiometric CeO 2 using our selected q-point sampling meshes are given in figures S2-S5. In figure 9 we compare this harmonic function for the stoichiometric and defective systems between T = 300-1500 K. We predict a reduction of ∼84%-89% at 600 K compared to stoichiometric CeO 2 , which we attribute to the defects disrupting the regular bonding networks and reducing the group velocities. As seen in figure 9 the OF⟨133⟩ and CeF⟨111⟩ defects have the largest impact and the Schottky defects the smallest, but all are predicted to significantly reduce the κ latt compared to pristine ceria. Within the SM-RTA the phonon linewidths are calculated as the sum of the contributions from three-phonon scattering processes between the reference mode λ and pairs of interacting modes λ ′ and λ ′′ according to equation (17): where Φ λλ ′ λ ′ ′ are the three-phonon interaction strengths, which include an expression for conservation of (crystal) momentum, n λ are the occupation numbers defined in equation (18) and the δ function enforces conservation of energy.
The expression for the linewidths can be approximated according to equation (19) [35]: P λ is the averaged phonon-phonon interaction strength for the mode λ defined in equation (20): where n a is the number of atoms in the primitive cell. N 2 (q, ω) is a weighted joint two-phonon density of states (w-JDOS) that counts the number of accessible energy and momentum-conserving scattering channels for a mode λ with wavevector q λ and frequency ω λ as shown in equation (21): Separate w-JDOS function N 2 (q, ω) and N (2) 2 (q, ω) are calculated for collision (Class 1) and decay (Class 2) process respectively as shown in equations (22) and (23), where the function ∆ enforces the conservation of crystal momentum: Using this model, allows the impact of changes to the phonon spectra in the defective systems on the phonon linewidths and lifetimes to be assessed. By averaging over q, w-JDOS functions analogous to those in equations (22) and (23) as a function of frequency only can be obtained as in equation (24).
2 (ω) andN 2 (ω) computed for the defective systems to those of stoichiometric ceria at T = 700 K. To enable comparison between the systems a normalization factor of 3n a 2 was applied, which is normally folded into the P λ term in equation (19). TheN 2 (ω) suggest the defective systems would show a higher density of energy-and momentum-conserving scattering channels due to a combination of more collision channels at higher frequencies and more decay channels at lower frequencies.
Notably, there is a marked enhancement in the number of scattering pathways at acoustic-mode frequencies where modes are likely to make larger contributions to the thermal conductivity. This analysis thus suggests broader linewidths and shorter lifetimes in the defective models. Comparisons of the w-JDOS functions calculated for the stoichiometric and defective systems at additional temperatures are given in figures S6 and S7.
An approximate value of κ latt can be determined for the defective systems using equation (20) if the P λ are replaced with a constant valueP, which can be estimated from the explicit calculation performed on the Figure 11. Comparison of the predicted κ latt of the five defective models examined in this work calculated using the approximate model of the phonon linewidths in equation (19) with the averaged interaction strengths P λ replaced by a constant valueP fitted for the stoichiometric CeO2 supercell at T = 600 K.
The predicted κ latt of the five defect models obtained using this method is compared to stoichiometric CeO 2 in figure 11. This method predicts a similar reduction to the analysis of the harmonic CRTA function, with reductions of 80%-88% at 600 K. This indicates that the primary impact of the defects on the thermal transport is likely to occur through a reduction of the group velocities, with a small secondary effect on the phonon lifetimes. Our predicted values of κ latt at this temperature are 1.4 W m −1 K −1 for the anion Frenkel model, and 1.6-1.8 W m −1 K −1 for the three Schottky defect models, all of which are much lower than the 9.9 W m −1 K −1 of stoichiometric CeO 2 .

Conclusions
In this study we have applied DFT to investigate the structure, defect formation free energies, lattice dynamics and thermal transport of the intrinsic anion Frenkel, cation Frenkel and Schottky defects in CeO 2 .
We predict the relative stabilities between 0-2750 K, in terms of eV per point defect, to be Schottky > anion Frenkel > cation Frenkel. We find that the phonon contributions to the free energy reduce the defect formation free energies by 1%-24% over the range of temperatures examined.
The phonon spectra of the defective systems show significant broadening, particularly around the Γ and X wavevectors, and the reduction of symmetry due to the defects introduces additional spectral signatures into the IR spectra that could enable the defects to be identified in experimental samples.
Calculations of the lattice thermal conductivity using approximate methods suggest that these intrinsic defects can have a significant impact on the thermal transport, reducing the κ latt of pristine ceria by as much as 88% at 600 K. This occurs primarily through a reduction in the phonon group velocities, with a small secondary reduction of the phonon lifetimes. This reduction is highly significant and could have a large impact on the performance of these materials as a surrogate for nuclear fuel materials, particularly given that the large and localised energy releases associated with nuclear fission events may produce high non-equilibrium concentrations of intrinsic defects. An important avenue for future work is the prediction of the impact of other point defects on the thermal conductivity. However, this would need to include appropriate charge corrections and account for environmental conditions such as oxygen partial pressure in the case of oxygen vacancies [65,66].
Finally, we comment on the limitations and generality of our methods. Using first-principles methods such as DFT imposes a practical upper limit of one to a few hundred atoms on the size of the supercells that can be used to model defective systems, which will inevitably correspond to relatively high defect concentrations in the order of 1%. Such small supercells may also impose an artificial periodic ordering of defects on a nm length scale that may not be present in real materials. Furthermore, real materials are likely to contain complex mixtures of both intrinsic and other (extrinsic) defects, which again cannot be practically modelled using first-principles methods at present. Finally, the computational overhead of lattice-dynamics calculations imposes further limits on the models, for example making it infeasible to use hybrid functionals in place of the less-demanding GGA + U method used in this study. Some of these restrictions may be partially lifted with cheaper PMs, including, as these methods become more mature, machine-learning (ML)-derived potentials with accuracy approaching that of first-principles methods. In this context, we consider the defect models in the present study to be close to the limit of what is practically achievable using first-principles techniques, and we expect our results to form a basis for more sophisticated models in the future as the theoretical methods and computing power required to tackle them become more widely available. On that note, the methods employed in this work are general and could be used to establish the impact of defects in other material systems on the structural dynamics and thermal transport, albeit subject to the computations being technically feasible. While the latter requirement may limit applications to less complex systems, there are a number of important materials that would benefit from the microscopic insight obtained from our approach, including, for example, a number of other important oxide materials with similar structure and defect chemistry to ceria.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.