Structural dynamics of Schottky and Frenkel defects in ThO 2 : a density-functional theory study †

Thorium dioxide (ThO 2 ) is a promising alternative to mixed-oxide nuclear fuels due to its longer fuel cycle and resistance to proliferation. Understanding the thermal properties, in particular the thermal conductivity, under reactor conditions is critical to the success of any candidate fuel material. ThO 2 has a higher thermal conductivity and thus a lower operating temperature than other fuel systems. However, the presence of defects in real materials directly in ﬂ uences the structural dynamics and physical properties, and the impact of defects on the properties of ThO 2 is largely unexplored. We have employed density-functional theory calculations to study the structure and energetics of the intrinsic Schottky and Frenkel defects in ThO 2 and their impact on the thermophysical properties. We identify the anion Frenkel defect to be the most stable, and we identify characteristic spectral signatures of the defects in the phonon dispersions and infrared spectra. We further employ two approximate models to assess the impact of the defects on the thermal transport and ﬁ nd that both types of defect are predicted signi ﬁ cantly to reduce the thermal conductivity. The methodology we present facilitates the prediction of the thermophysical and transport properties of defective materials at an atomistic level, and should be readily transferrable to existing and in-development nuclear fuel systems.


I. Introduction
As fossil fuels are phased out in favour of renewable energy to limit global warming, 1 nuclear power plays an increasingly critical role in meeting worldwide energy demand.Thorium-based fuels are ideal candidates to address this challenge, due to the potential safety and economic benets compared to more traditional actinide-dioxide fuels such as UO 2 .ThO 2 can be used in mixed-oxide (MOX) fuels with uranium or plutonium to improve reactor safety, as its higher thermal conductivity results in a lower operating temperature and prevents overheating. 2,3ThO 2 also reduces the risk of proliferation and the amount of highly radioactive spent fuel waste generated.Indeed, due to its large domestic supply of thorium (up to a third of global reserves), 3 India has positioned itself to take full advantage of these benets by aiming to use ThO 2 in advanced heavy water reactors (AHWRs). 4 complete understanding of the thermophysical properties of ThO 2 is critical to controlling fuel performance at the high temperatures and pressures in nuclear reactors.An important consideration here is the defect chemistry.The energy release associated with nuclear events in reactors can generate large nonequilibrium concentrations of defects, and these can have a direct inuence on the physical properties.Since experimental measurements on ThO 2 are challenging, computational modelling using rst-principles methods such as density-functional theory (DFT), in conjunction with modern high-performance computing (HPC), plays a key role in contemporary nuclear research, particularly when atomistic or nanoscale details are required. 5nder ambient conditions ThO 2 is non-magnetic with a cubic structure in the Fm 3m space group (no.][8][9][10] A phase transition to an orthorhombic Pnma structure (no.62) begins at actinide oxides is controlled by heat transport through phonons. 12,13The thermal conductivity is however generally low for oxide materials due to the limited transport through the F 2g mode, despite its large mode group velocity. 12There is also relatively little known about how defects might impact the thermal transport, although a recent study of Deskins et.al. 14 predicts a decrease in the thermal conductivity of ThO 2 from point defects including oxygen and thorium vacancies.This has also been observed in other work on ThO 2 (ref.8][19] While a variety of point defects and defect clusters can be found in uorite structures, Schottky and Frenkel defects are simple, charge-neutral defects.These will both be produced under reactor operating conditions, and are likely to have an impact on the thermal conductivity, 15,19 but unravelling this at the nanoscale requires computational approaches since disentangling the effects of specic defects using experimental approaches is challenging and currently limited.
In this work, we use DFT modelling to perform a comprehensive study of the intrinsic defects in ThO 2 and to assess their impact on the structural dynamics.Whereas the structural dynamics of stoichiometric ThO 2 are well documented, [20][21][22][23][24] little is currently known of defective ThO 2 due to the computational challenges inherent in studying defective systems.We model the phonon dispersion and density of states of stoichiometric ThO 2 and compare these to models with intrinsic Schottky and Frenkel defects.This allows us to determine key spectral signatures for identifying these defects in material samples.We also quantify the impact of the lattice dynamics on the defect formation thermodynamics.Finally, we estimate the potential impact of the defects on thermal transport, which is an important material property both for reactor design and the storage of waste material.

II. Methods
DFT calculations were performed using the Vienna Ab initio Simulation Package (VASP) code. 25 96-atom 2 Â 2 Â 2 supercell expansion of the 12-atom ThO 2 conventional unit cell was constructed and used to build models of three intrinsic defects, namely cation (Th) Frenkel, anion (O) Frenkel and Schottky defects, with different relationships between the vacancy/interstitial sites (Fig. 1).For comparison with the literature and to verify that our choice of cell size was appropriate we also modelled the defects in 3 Â 3 Â 3 supercells with 324 atoms.
The spin-polarised simulations used a plane-wave basis set with the electron-ion interactions described using projector augmented-wave (PAW) pseudopotentials 26,27 with frozen [Xe] and [He] cores for Th and O respectively.The plane-wave cutoff was set to 500 eV.0][31] Spin-orbit coupling was not included as this has been shown to have a comparatively small effect in ThO 2 . 32The Brillouin zone of the 12-atom ThO 2 conventional cell was sampled using a G-centred Monkhorst-Pack k-point grid 33 with 4 Â 4 Â 4 subdivisions, which was proportionally reduced to a 2 Â 2 Â 2 sampling for the 2 Â 2 Â 2 and 3 Â 3 Â 3 supercells.Pristine ThO 2 was geometry optimised at constant pressure to tolerances of 10 À8 eV and 10 À3 eV A À1 on the electronic energy and ionic forces respectively.The defective models were optimised with the cell shape and volume xed to that of the equivalent supercell of the geometry optimized pristine ThO 2 .
To model the lattice dynamics, the 2 nd -order harmonic interatomic force constants (IFCs) were obtained using the nite-displacement approach implemented in the Phonopy code 34 with VASP as the force calculator.The displacements were generated in the 96-atom cells used for the stoichiometric and defective models.The atom-projected phonon density of states (PDOS) curves for the models were evaluated on G-centred grids of 10 Â 10 Â 10 phonon wavevectors (qpoints).For the stoichiometric model, a transformation matrix was used to convert to the three-atom primitive cell to evaluate the phonon dispersion.For the defective models, we made use of the band unfolding methodology 35,36 implemented in the Phonopy code to project the dispersions onto the ThO 2 primitive cell for comparison to the stoichiometric structure.A non-analytical correction was applied using Born effective charges and high-frequency dielectric constants computed with the density-functional perturbation theory (DFPT) routines in VASP. 37To evaluate the IR and Raman activities of the G-point phonon modes, we used the Phonopyspectroscopy project, 38 which is documented in ref. 39.Brief details of the methodology used for simulating the IR and Raman spectra are given in the ESI.† The Phono3py code 40 was used to compute the thermal conductivity of stoichiometric ThO 2 within the single-mode relaxation-time approximation by calculating the 3 rd -order IFCs using the 96-atom supercell.The phonon group velocities, lifetimes and thermal conductivity were deemed to be converged with a 48 Â 48 Â 48 sampling mesh for the ThO 2 primitive cell.As described in the text, we then estimated the thermal conductivity of the 96-atom defect models using one of two approximate models. 41In principle, a q-point mesh with 15 Â 15 Â 15 subdivisions should be adequate for sampling the 32Â smaller Brillouin of the supercells.However, we were only able to use a smaller mesh with 10 Â 10 Â 10 subdivisions due to memory constraints.As discussed in the text below, we estimate that this leads to an error in the calculated properties on the order of 20%.Finally, we also generated for the supercell models two-phonon weighted joint density of states (w-JDOS) functions, using the same 10 Â 10 Â 10 mesh.
We note that our chosen supercell size may have an impact on the structural dynamics calculations, i.e. the calculated frequencies, simulated phonon spectra, and predicted thermal transport.However, the substantial workload means that it is computability infeasible to perform these calculations on the larger 3 Â 3 Â 3 defect supercell models, so we are unable to quantify the impact as we do for the defectformation energies.

A. Structures of the stoichiometric and defective systems
Our calculated lattice constant for stoichiometric ThO 2 (Fig. 1A) is 5.5605 Å, compared to the experimental range of 5.5921-5.604Å. 6,[46][47][48][49][50] Standard (athermal) DFT does not take into account the effect of temperature.An experimental lattice parameter of 5.5828 Å at 0 K can be extrapolated from the empirical function of Tyagi and Mathews, 51 which is still higher than our calculated value.This indicates that the underestimation of the lattice parameter may be due to the DFT functional used.
The Schottky defect involves the removal of a complete ThO 2 formula unit to create two oxygen vacancies and a Th vacancy Following Govers et al. 45 we considered only the rst nearestneighbour positions for the vacancies.We investigated three potential arrangements where the two V O are arranged along the h111i, h110i and h100i directions with respect to the central V 0000 Th , which we denote Sh111i, Sh110i and Sh100i respectively (Fig. 1B).
The cation and anion Frenkel defects are a type of intrinsic disorder and involve the formation of pairs of vacancies and interstitials according to: To prevent spontaneous recombination during the geometry optimisation, the interstitial and vacancy must not be adjacent.As shown in Fig. 1C, we therefore created O Frenkel defects with the O 00 i and V O arranged as 2 nd -nearest neighbours along the h111i and h133i directions.These congurations, which we refer to as OFh111i and OFh133i, are equivalent to those examined in other literature. 42,43,45In our Th Frenkel models (Fig. 1D), the Th i and V 0000 Th are arranged along the h111i and h201i directions, which we denote ThFh111i and ThFh201i respectively.The ThFh201i has previously been studied, 42,43,45 but as discussed below our calculations predict the ThFh111i conguration to be energetically more stable.

B. Defect formation energies
The calculated defect formation energies are presented and compared to other computational studies in Table 1.
Experimental studies of the O Frenkel disorder indicate defect formation enthalpies in the range of 2.28-4.68eV from diffusion measurements using the isotope-exchange method, 52 and a value of 4.2 eV from tting high-temperature heat capacity data. 53These are a good match to our calculated formation energy of 4.42 eV for the most stable OFh133i defect.The higher stability of the OFh133i defect compared to the OFh111i also agrees with potential-model (PM) calculations using larger 4 Â 4 Â 4 42 and 5 Â 5 Â 5 supercells. 44Ando et al. concluded that O Frenkel and Schottky defects would need to have energies of 2.82 eV and 4.2 eV, respectively in order to explain the selfdiffusion behaviour measured in a ThO 2 single crystal, which indicates that the O transport is primarily via Frenkel defects. 54his is the only experimental value we could nd for the Schottky defect, but it is reassuring that it has been calculated to be higher in energy than the O Frenkel defect, thus supporting our DFT values (Table 1).
The Schottky defect has been studied computationally using DFT 43 and PMs using the shell, 55 rigid ion 42,56,57 and embedded atom 44 formalisms.Using a 2 Â 2 Â 2 supercell, we predict defect formation energies between 5.16-5.73eV, in the order Sh110i < Sh111i < Sh100i, in agreement with free-energy calculations on ThO 2 at 1000 K and 0.2 atm 43 and similar calculations on UO 2 . 58Although the energy differences between the congurations are marginal, the relative stabilities are changed when using larger supercells such that our 3 Â 3 Â 3 supercell models predict a stability order of Sh111i < Sh110i < Sh100i.PM calculations using larger supercells predict the same stability trends and small energy differences as our DFT calculations on larger cells.We attribute the differences between cell sizes to better screening of the coulombic repulsion between positively charged oxygen vacancies in the larger models.
For the Th Frenkel disorder, we nd that the ThFh111i conguration is more stable than the ThFh201i defect examined in other literature, 45 although the formation energy of 12.57 eV makes this defect considerably less favourable than the Schottky and O Frenkel defects.
The trend in stability of the different Frenkel pair congurations appears to be related to the distance between the vacancy and interstitial: the more stable ThFh111i and OFh133i congurations have interstitial-vacancy distances of 4.48 and 5.87 Å, respectively, whereas the less stable ThFh201i and OFh111i congurations have larger separations of 5.96 and 6.88 Å.We therefore conclude that the closer the interstitial and vacancy sites are, the more stable the defect pair is, likely because shorter distances optimise the attractive coulombic interactions between the oppositely-charged point defects.
56,59 The potential model studies in ref. 42 and 44 predict a qualitative stability ordering of OF < S < ThF, in agreement with our results.While the trend is unchanged when the formation energies are corrected to the innite-dilution limit, the correction predicts a greater stabilisation of the Schottky defect relative to the O Frenkel.We infer from this that at the high defect concentrations associated with small supercell models the interactions between defects may have a stabilising effect.
To account for thermal effects, we also calculated defectformation free energies DA f from 0 K up to the melting temperature of ThO 2 (3630 K). 60 The Helmholtz free energies of formation, DA f , of the Schottky and Frenkel defects are calculated as: Table 1 Formation energies (eV) for the Schottky and Th/O Frenkel defects examined in this study compared to other theoretical literature.The ThF N and OF N structures denote defect energies corrected to the infinite dilution limit.The size of the supercell used in the calculations is given in the first row of the table.Comparison studies are identified as being performed using DFT, as in the present work, or using classical potential models (PM).Where available, athermal energies are also compared to free energies computed at where A Stoich , A Schottky and A Frenkel are the free energies of the supercell models of stoichiometric and defective ThO 2 , and A ThO 2 is the free energy of a single ThO 2 formula unit.All free energies are temperature dependent.The Helmholtz energies in eqn (4)/( 5) are calculated based on the phonon contributions to the free energy A vib (T), computed as: Here Z vib (T) is the vibrational partition function at temperature T, k B is the Boltzmann constant, the sums run over phonon modes with frequency u qj indexed by a wavevector q and a band index j, and the normalisation factor N is the number of unit cells in the crystal (equivalent to the number of q in the summation).The internal energy and entropy of the phonon system, U vib (T) and S vib (T), are given by: The total Helmholtz energy of the system, A(T), can then be computed by adding the lattice internal energy U 0 , which we equate to the DFT total energy: Under the assumption of zero pressure, the Gibbs free energy G(T) and the Helmholtz energy A(T) can be equated.
However, as our calculation are performed at constant volume there may be a residual cell pressure in the defect models, and so we strictly work with A(T).To determine the A vib (T) we use one of two methods.In method 1, we consider frequencies u qj sampled at wavevectors q across the Brillouin zone.In Method 2, we consider only the u j at the zone centre q ¼ G.While an approximation, this simplied model is widely used for treating large unit cells where further supercell expansions to calculate frequencies at q s G without interpolation are not practical.
Fig. 2 shows the defect-formation free energies of the Sh100i, Sh110i, Sh111i, ThFh111i and OFh133i defects as a function of temperature calculated using method 1.The formation energies decrease with temperature, which implies that it is easier to generate defects at elevated T. This is in line with experimental observations, 52 which show that the generation of O Frenkel defects at higher temperatures leads to enhanced oxygen diffusion. 21Comparing the formation energies at T ¼ 0 and 3630 K predicts decreases of 0.85-1.02eV for the Schottky defects, 1.33 eV for the OFh133i defect, and 0.30 eV for ThFh111i defect.We therefore conclude that temperature effects are signicant, but do not affect the overall stability trend.Method 2, using three different models for the contributions from the acoustic modes, predicts a small $0.25-0.5 eV difference in the formation energies at T ¼ 3630 K (see ESI †).For both methods the trend in stability ordering remains the same, with the defect formation energies in the order of OF < S < ThF.
To understand how the changes in the defect formation energies impact the defect concentrations N D , we estimated the concentrations using the Boltzmann distribution: where N is density of "sites".Our defect formation energies are relative to the supercell, so N ¼ 1/V SC ¼ 7.10 Â 10 20 cm À3 .We note that we do not account for degeneracies due to the equivalent O/Th sites in the supercell, which a crude estimate suggests could potentially increase the concentrations by between 10-100Â.The temperature dependence of the concentrations of the three Schottky and the OFh133i and ThFh111i defects are shown in Fig. 3. Taking a concentration of N D ¼ 1 cm À3 as a threshold, we predict signicant concentrations of OF to form at a much lower T ¼ 1030 K than the more energetically unstable ThF defect (T ¼ 3000 K).

C. Phonon dispersion and phonon density of states
Since the majority of the thermal transport in ThO 2 is through the phonons, an accurate description of the structural dynamics is crucial.The calculated phonon dispersion and atom-Fig.2 Defect-formation free energies for the OFh133i, ThFh111i, Sh100i, Sh110i and Sh111i defects as a function of temperature calculated using method 1.
projected density of states are shown in Fig. 4 and compared to data from inelastic neutron-scattering (INS) measurements. 21he ThO 2 primitive cell contains three atoms and the dispersion thus has three acoustic and six optic branches at each wavevector q.The three lower-frequency acoustic modes split into longitudinal (LA) and transverse (TA) modes and are mainly attributed to motions of the heavier Th sublattice.The six higher-frequency optic modes again divide into longitudinal and transverse (LO/TO) modes and are predominantly associated with motion of the lighter O sublattice.The general features of the spectrum are common to UO 2 and PuO 2 and have been demonstrated in calculations using DFT 61 and potential models, 62 and in experimental time-of-ight INS measurements. 61,63 list of modes and frequencies at each of the high-symmetry q in the Fm 3m Brillouin zone point is provided in  and O).These assignments are supported by previous work. 70The predicted frequencies are also in excellent agreement with INS measurements collected at 273 K, with an average 6% and maximum 13% difference.The discrepancies may in part be attributed to the neglect of temperature effects (i.e.thermal expansion and higher-order anharmonicity) in our calculations. 24he overlap between the LA and TO branches around q ¼ X between $6.2-6.6THz was predicted for ThO 2 previously using LDA calculations, 23 and has also been predicted to occur in UO 2 and PuO 2 using LDA with dynamical mean-eld theory (DMFT). 12his feature is noteworthy because the acoustic modes tend to be the majority heat carriers, and the overlap with low-lying optic modes may facilitate enhanced phonon scattering and a reduced thermal conductivity in these materials. 23For UO 2 , the Debye approximation suggests that the F 2g (LO) mode only contributes 0.11 W m À1 K À1 along the LO branch from G to X, compared to 1.68 W m À1 K À1 along the LA branch. 12Park et al. also found that the balance of the thermal transport shis towards higherfrequency modes with shorter wavelengths at elevated temperature, resulting in reduced thermal transport. 73hile the phonon spectrum of stoichiometric ThO 2 is well documented in the literature, [20][21][22][23][24]70 this is the rst study to explore the impact of defects. Thephonon dispersions of the defective structures, obtained in terms of the ThO 2 primitive cell using band unfolding, 35,36 are shown together with the calculated PDOS in Fig. 5.In Table 3, we present estimates of the range of frequencies spanned by modes at the highsymmetry q for comparison with Table 2.
The defective structures all show some degree of broadening across the phonon dispersion, including notable changes to the acoustic modes that are most pronounced along G-X direction around 2 THz.For the Frenkel defects (Fig. 5 D and E) the dispersions show notable mode soening at X around 6.6 THz.This suggests that the defects affect the vibrations of the O sublattice and have similar effects to the temperature-induced soening that occurs to this mode. 24Ghosh et al. 71 have linked reductions in the frequencies of the X-point phonon modes at 3.5 and 6.2 THz in ThO 2 to an increase in lattice strain, and have suggested that large strains lead to the phonon mode at 6.2 THz becoming imaginary.This effect has also been observed in CeO 2 , 74 where the appearance of an imaginary mode under thermal expansion at high temperature predicates a superionic transition to a disordered oxygen sublattice.Taken together with our results, this provides a potential mechanism for how O Frenkel defects could enhance oxygen diffusion.We note however that while a mode soening is visible in Fig. 5E we do not observe an imaginary mode, which may be because we do not take into account temperature effects.
A handful of recent studies have investigated the effect of point defects on the phonon scattering in ThO 2 , UO 2 and PuO 2 . 17,18,75For ThO 2 , non-equilibrium molecular dynamics simulations 75 found that Th and O vacancies affected the lowerfrequency modes with the largest contributions to the thermal Fig. 3 Predicted defect concentrations as a function of temperature for the Sh100i, Sh110i, Sh111i, ThFh111i and OFh133i defects using the calculated defect formation free energies.transport.Th vacancies were predicted to have a larger impact, with 0.1% Th and O vacancies predicted to induce 62 and 20% reductions in the thermal conductivity respectively.Comparing the PDOS curves of the defective models in Fig. 5 with the PDOS curve of pristine ThO 2 reveals an enhanced contribution from O at lower frequencies, suggesting that the defects may facilitate enhanced coupling between the O and Th sublattices.There are also some subtle changes to the DOS of the acoustic modes below $2.5 THz, which may indicate changes in the number of phonon scattering pathways.

D. Infrared spectra
The G-point optic modes can be investigated using infrared (IR) and Raman spectroscopy.For stoichiometric ThO 2 there is a triply-degenerate Raman-active F 2g mode at 13.6 THz, a doubly-degenerate IR-active F 1u TO mode at 8.5 THz, and a singly-degenerate IR-active F 1u LO mode at 17.1 THz.Fig. 6 compares simulated IR spectra of pristine ThO 2 to the predicted spectra of the defective models.To allow for comparison between the spectra, and for compatibility with existing literature, we applied the nominal linewidth of 0.75 THz (including instrumental broadening) estimated from room-temperature IR reectivity measurements. 69We note however that our calculated intrinsic linewidths at T ¼ 300 K are considerably sharper, at 0.034 and 0.098 THz for the F 1u and F 2g modes respectively.We attribute this in part to the difficulties inherent in extracting accurate linewidths from experimental spectra.However, we also note that the linewidths we calculate are the intrinsic linewidths due to three-phonon scattering.This approximation neglects higher-order scattering and additional linewidth contributions from phenomena such as phonon-defect scattering that may be present in experimental samples.We also consider the structure at 0 K, and thermal expansion at nite temperature would usually narrow the phonon density of states and enhance the phonon scattering and linewidths.
The calculated frequencies of the F 1u TO and LO modes at 8.5 and 17.1 THz compare well to single-crystal measurements performed by Knight et al. 68  The defects break the symmetry of stoichiometric ThO 2 and result in some of the zone-boundary modes being "folded" to the G-point.This in turn leads to additional IR-and Ramanactive modes in the defective systems, which are oen associated with localised vibrations of atoms around the defect sites.As can be seen in Fig. 6B-F this results in more featured IR spectra.Inspection of the eigenvectors of the IR-active modes indicates that the vibrations are somewhat complex with contributions from both the Th and O sublattices.
Due to the much higher computational cost of calculating Raman activities, we were unable to predict Raman spectra for the defective models, but for completeness we show a simulated Raman spectrum of stoichiometric ThO 2 in Fig. S2.† The spectrum features a single peak at 13.6 THz from the F 2g mode, in line with other computational studies, 16,20,22,71 but underestimated compared to the experimental literature. 9,64-67

E. Thermal conductivity
Using the single-mode relaxation time approximation (RTA) model, the macroscopic thermal conductivity k latt (T) is obtained as a sum of contributions from individual phonon modes l according to: 40 where C l are the modal heat capacities, v l are the group velocities, s l are the lifetimes, V is the volume of the unit cell, and N is the number of wavevectors q included in the summation.The C l are given by: where u l are the phonon frequencies.The v l are given by: where the W l are the phonon eigenvectors and D(q) is the dynamical matrix.Finally, the s l are obtained from the phonon linewidths G l as: The G l are obtained as the imaginary part of the phonon selfenergy using third-order perturbation theory, and the method is documented in detail in ref. 40.Computing the linewidths requires the third-order force constants, and obtaining these represents a substantial computational workload.
Since the k latt are tensor quantities, we discuss the scalar average of the diagonal elements, which we denote k latt in the following, i.e.
We have omitted the explicit temperature dependence for brevity.The k latt are representative of what would be measured for a pressed pellet or polycrystal.
The thermal conductivity of ThO 2 is well studied in the literature using DFT 23,24,70 and other theoretical ).Given the scatter in the experiments, we suspect our prediction overestimates the measurements because it is difficult to prepare high-quality bulk samples that are free from defects such as grain boundariesthese defects limit the transport through modes with long mean free paths and thus reduce the overall k latt .In keeping with this, Mann et.al. 86 obtained a room-temperature k latt between 18-20 W m À1 K À1 for single crystal ThO 2 , which is higher than previous experimental results and closer to our prediction.
Our results are also consistent with other studies using DFT 23,24 and potential models. 42,70At 300 K, the values of 20.6 and 17.9 W m À1 K À1 calculated using lattice dynamics coupled with the Boltzmann transport equation (LD-BTE) 42 and equilibrium molecular dynamics (EMD) 70 respectively compare well with our value of 20.8 W m À1 K À1 .The DFT study in ref. 23 used a similar technical setup to the present work and reported a value of 9.3 W m À1 K À1 at 600 K, which compares well to our prediction of 10 W m À1 K À1 .However, we note that the choice of functional appears to have a signicant impact, as the calculations in ref. 24 using the WC functional predicted a somewhat lower 7.7 W m À1 K À1 .
To investigate the effect of crystal grain size on the thermal conductivity we employed a simple boundary-scattering model to limit the phonon mean free paths.As shown in Fig. 8, we nd that beyond a grain size of 5 mm the thermal conductivity Table 3 Predicted frequency ranges (THz) of the G, X and L modes in the defective ThO 2 models considered in this work, estimated from the unfolded dispersions in Fig. 5, for comparison with the data for stoichiometric ThO 2 in Table 2 Mode Sh100i Sh110i Sh111i OFh133i ThFh111i plateaus, in agreement with the ndings of Nakamura and Machida. 23At elevated temperature, the phonon mean free paths decrease and grain-boundary effects become less prominent.Although most experimental studies do not report a grain size, those that do indicate grain sizes in the range of 1-6 mm, 5,8,75 in which case we would not expect the grain size to have a signicant impact on the k latt .
Computing the third-order force constants to calculate the phonon lifetimes of the defective structures would be prohibitively expensive and is therefore not feasible, so we explored instead the approximate models developed in ref. 41.We rst consider the constant relaxation-time approximation (CRTA) model, 41,87 where the s l are replaced with a constant relaxation time, s CRTA according to: We omit the temperature dependence of the terms for brevity but note that both the s CRTA and the summand are temperature dependent, the latter by virtue of the C l term (c.f.eqn (11)/( 12)).
Since the terms in the summand are calculable from the harmonic force constants, this model allows for quantitative assessment of the impact of the defects on the k latt due to changes in the harmonic phonon frequencies and group velocities.Remarkably, we nd that the defects reduce the harmonic term in eqn ( 16) by between 83.2 and 90.4% compared to stoichiometric ThO 2 at 600 K (Fig. 9).The Sh111i defect is predicted by this model to have the smallest impact, whereas the two Frenkel defects have the largest impact.Given that the k l in eqn (11) are proportional to the square of the group velocities, and eqn ( 12) is a shallow function of temperature, we can attribute these changes to disruption of the bonding network reducing the group velocities.
We next consider the approximate model for the phonon linewidths G l outlined in ref. 40.In the perturbative treatment the G l are obtained as a sum of three-phonon scattering processes between phonon triplets consisting of a reference mode l and two interacting modes l 0 and l 00 : À n l 00 Þ½dðu þ u l 0 À u l 00 Þ À dðu À u l 0 À u l 00 Þg (17)   In this equation the F ll 0 l 00 are the three-phonon interaction strengths, the n l are the occupation numbers given by, and the function d enforces the conservation of energy.An approximate expression for the linewidths Gl can be written: 40 In this expression P l is the averaged phonon-phonon interaction strengths for the mode l given by: where n a is the number of atoms in the primitive cell and there are 3n a bands at each q.N 2 (q,u) is a weighted joint two-phonon density of states (wJDOS) that counts the number of energy and momentum-conserving scattering channels available to a mode l with wavevector q l and frequency u l : N 2 (q,u) ¼ N (1) 2 (q,u)+N (2) 2 (q,u) Fig. 10 Comparison of the weighted joint density of states (wJDOS) functions defined in eqn (21) for stoichiometric ThO 2 (grey) and the Sh110i (blue), Sh100i (purple), Sh111i (red), ThFh111i (green) and OFh133i (orange) defect models.The first column (A/D/G/J/M) compares the N (1) 2 for collision processes, the second column (B/E/H/K/N) compares the N (2) 2 for decay processes, and the third column (C/F/I/L/O) compares the sum N 2 .
The two components N (1) 2 (q,u) and N (2) 2 (q,u) represent collision and decay (Class 1/Class 2) processes, respectively, and are dened as: where the function D enforces the conservation of (crystal) momentum.This model serves two purposes.Firstly, the wJDOS functions can be calculated using the harmonic frequencies and enable a partial assessment of how changes in the phonon spectrum are likely to impact on the phonon linewidths and lifetimes.This is most easily done by averaging over q to obtain a function of frequency only: Secondly, if assume the P l in eqn (20) can be replaced with a constant value P, e.g.tted for stoichiometric ThO 2 , we can also predict an approximate k latt . 41he averaged wJDOS functions computed for ThO 2 and the defective systems are compared in Fig. 10.To allow for comparison between systems, we have normalized the functions by (3n a ) 2 , which is folded into the P l term in eqn (19).Peaks in the N 2 (u) at low frequencies are likely to indicate an enhanced broadening of the acoustic modes, which as noted above typically contribute the most to the thermal conductivity, leading to shorter lifetimes and suppressing the heat transport.Fig. 10 shows that the defects generally produce small increases in the N 2 (u) at low frequencies, and thus suggests that changes to the phonon frequency spectrum may play a role in reducing the thermal conductivity.
To apply the approximate model in ref. 41 a value of P ¼ 9.722 Â 10 13 was determined such that for the 96-atom stoichiometric supercell and a 10 Â 10 Â 10 q-point sampling the k latt for bulk ThO 2 at T ¼ 600 K is recovered.Although the F ll 0 l 00 in eqn ( 17) are formally temperature independent, in practice the limitations of using a single value results in a weak temperature dependence, and so 600 K was chosen as a compromise.The k latt of pristine ThO 2 calculated using eqn (11), for the primitive cell, is compared to that obtained using the supercell model with the tted P in Fig. S3.† We note that the 10 Â 10 Â 10 sampling mesh produces a harmonic function in eqn ( 16) that is $20% smaller than for the primitive cell, indicating that the n l are not fully converged with respect to the sampling mesh (Fig. S4 †).We found that the harmonic function for the supercell was slow to converge with respect to the sampling mesh, and we were unable practically to increase it due to memory constraints.This is a technical limitation of the current implementation in Phonopy.However, we would expect a similar relative error for the defect supercells, and we also note that the tted P may partially compensate for this issue.
Fig. 11 shows the thermal conductivity of the defective models predicted using this approximate model.Interestingly, we nd the same 83.2-90.4% reductions in the k latt at 600 K as predicted by considering the harmonic function in eqn (16).This suggests that the small changes in the numbers of energy-conserving scattering channels suggested by the wJDOS functions in Fig. 10 do not have a material impact on the phonon lifetimes.While perhaps counter-intuitive, it is possible that the averaged functions in Fig. 10 will mask larger variation in the N 2 (q,u) of individual wavevectors (i.e.increases and decreases in the modal k l at different q may compensate one another).As noted above, our relatively small defect supercells represent a signicantly higher defect concentration than is likely to be achieved under thermal equilibrium (Fig. 3), and so the reduction in the k latt is unlikely to be as prominent in real materials.On the other hand, the signicant energy release associated with nuclear decay events in working reactors is likely to produce regions with much larger non-equilibrium local concentrations of defects, and our predictions suggest these may form "hotspots" where heat cannot effectively be dissipated.This could have several important consequences, including cracking of the fuel rods as a result of accumulated thermal stresses and swelling due to accelerated release of ssion gases. 17

IV. Conclusions
In summary, we have performed a comprehensive study of the structure, energetics, lattice dynamics and thermal transport of stoichiometric ThO 2 and models of the intrinsic Schottky and O/Th Frenkel defects.
First-principles DFT calculations predict that the O Frenkel defect is the most energetically favourable.We have identied a more stable conguration of the Th Frenkel defect, not found in previous studies, with the interstitial and vacancy arranged along the h111i rather than h201i direction.We have also investigated the phonon contributions to the thermodynamic free energy and found that these terms may selectively stabilise the defects relative to pristine ThO 2 and thereby reduce the defect formation energies at elevated temperature.
Our lattice-dynamics calculations demonstrate that the presence of defects leads to a signicant broadening of the phonon dispersion and introduces additional IR-active modes that may serve as useful spectral signatures for identifying the defects in real materials using routine spectroscopy.
Finally, we also applied two approximate models to assess the impact of the defects on the thermal transport.We predict that the defects can reduce the (local) heat transport by as much as $90% by disrupting the regular bonding network of the stoichiometric material and reducing the phonon group velocities.Interestingly, we also nd that, within the approximation of a constant phonon-phonon interaction strength, the changes in the phonon spectra due to the presence of defects do not lead to signicant changes to the phonon lifetimes.While verifying the predicted reduction in k latt would require explicit computation of the third-order force constants, it would be somewhat surprising if the defects led to so much weaker phonon-phonon interactions as to counteract the large reduction in the group velocities.
Overall, our modelling indicates that the presence of defects can have a signicant impact on the thermal transport, and this is particularly signicant in operating reactors where nuclear events are likely to lead to large non-equilibrium local defect concentrations.On the other hand, our results suggest that the impact on the thermal conductivity can be minimised if the defect concentrations can be controlled, which would improve the safety and performance proles of the fuel.If the thermal performance of an alternative fuel material like ThO 2 could be made tolerant to the changes in its structure, i.e. to the formation of defects, then it is reasonable to suggest that this may be benecial from a safety perspective compared to UO 2 . 88urther research in this direction will require reliable means to identify and quantify defects in real materials and to assess their impact on the physical properties; modelling can play an important role in this, for example, as we have shown, by identifying spectral signatures that could be measured using routine spectroscopic techniques.

Fig. 1
Fig. 1 Intrinsic ThO 2 defects investigated in this work: (A) pristine ThO 2 with the Th and O atoms in green and pale pink respectively; (B) three Schottky defect models ðV 0000 Th þ 2V O Þ with the V 0000 Th and V O shown in orange and yellow and arranged along the h111i, h110i and h100i directions; 42-45 (C) two O Frenkel defect models with the O 00 i and V O (blue/yellow) arranged along the h111i and h133i directions; and (D) two Th Frenkel defect models with the Th i and V 0000 Th (purple/orange) arranged along the h111i and h201i directions.

Fig. 4
Fig. 4 Calculated phonon dispersion and atom-projected density of states of stoichiometric ThO 2 .Inelastic neutron-scattering (INS) measurements from ref. 21 are overlaid as red dots.

Fig. 5
Fig. 5 Phonon dispersion curves and atom-projected density of states for the Sh100i, Sh110i and Sh111i (A-C), ThFh111i (D) and OFh133i (E) defect models.The spectral functions in the dispersion curves are shown on a logarithmic scale with the phonon dispersion of pristine ThO 2 overlaid as a yellow line.Each PDOS shows the total DOS in grey with the contributions from Th and O overlaid as blue and orange lines.

Fig. 6
Fig. 6 Simulated IR spectra of stoichiometric ThO 2 (A) and models with Th (B) and O (C) Frenkel defects and one of three Schottky defect configurations (D-F).Each spectrum is generated with the nominal linewidth of 0.75 THz derived by Axe et al.69

Fig. 7
Fig. 7 Comparison of the k latt calculated for stoichiometric ThO 2 using the RTA model in eqn (11) with selected experimental and theoretical studies from the literature.

Fig. 8
Fig.8Thermal conductivity of stoichiometric ThO 2 at T ¼ 300, 600, 900, 1200 and 1500 K as a function of the crystal-grain size, predicted using a boundary-scattering model.

Fig. 9
Fig. 9 Comparison of the harmonic function in eqn (16) calculated for the 96-atom stoichiometric ThO 2 supercell, three configurations of Schottky defects, and the lowest-energy O and Th Frenkel defects.

Fig. 11
Fig.11Comparison of the k latt of the five defective models considered in this work calculated using the constant P l model with the P ̃fitted for the stoichiometric ThO 2 supercell at T ¼ 600 K.

Table 2
and64ompared to INS,21,64Raman9,64-67 and IR68,69 measurements.At the G point, stoichiometric ThO 2 has a doubly-degenerate TO mode at 8.5 THz (alternate motion of Th and O atoms), a triply-degenerate Raman-active mode at 13.6 THz (alternate motion of O), and a singly-degenerate LO mode at 17.1 THz (alternate motion of Th

Table 2
72lculated phonon frequencies (THz) of the modes at G, X and L in ThO 2 compared to data from other computational studies and experimental inelastic neutron-scattering (INS), Raman and IR measurements.Since predicted frequencies are often sensitive to the cell volume,72we also show the predicted lattice constants from the computational studies a o[ Å] Open Access Article.Published on 07 January 2022.Downloaded on 9/16/2023 6:25:48 AM.This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.