<EMAIL_ADDRESS>
# First-principles determination of the phonon-point defect scattering and
thermal transport due to fission products in ThO2

Linu Malakkal1 [  Ankita Katre2<EMAIL_ADDRESS>Shuxiang Zhou1 Chao
Jiang1 David H. Hurley3 Chris A. Marianetti4 Marat Khafizov5
<EMAIL_ADDRESS>1Computational Mechanics and Materials Department, Idaho
National Laboratory, Idaho Falls, Idaho 83415, USA  2Department of Scientific
Computing, Modeling and Simulation, SP Pune University, Pune-411007 India
3Idaho National Laboratory, Idaho Falls, ID, 83415, USA  4Department of
Applied Physics and Applied Mathematics, Columbia University, 500 West 120th
Street, New York, New York 10027, United States  5Department of Mechanical and
Aerospace Engineering, The Ohio State University, 201 West 19th Ave, Columbus,
Ohio 43210, United States

###### Abstract

This work presents the first-principles calculations of the lattice thermal
conductivity degradation due to point defects in thorium dioxide using an
iterative solution of the Pierels-Boltzmann transport equation. We have used
the non-perturbative Green’s function methodology to compute the phonon-point
defect scattering rates that consider the local distortion around the point
defect, including the mass difference changes, interatomic force constants,
and structural relaxation near the point defects. The point defects considered
in this work include the vacancy of thorium (ThV) and oxygen (OV),
substitutions of helium (HeTh), krypton (KrTh), zirconium (ZrTh), iodine (ITh)
and xenon (XeTh) in the thorium site, and the three different configurations
of the Schottky defects. The results of the phonon-defect scattering rate
reveal that, among all the considered intrinsic defects, the thorium vacancy
and helium substitution in the thorium site scatter the phonon most due to the
substantial changes in the force constant and structural distortions. The
scattering of phonons due to the substitutional defects unveils that the
zirconium atom scatters phonons the least, followed by xenon, iodine, krypton,
and helium. This is contrary to the intuition that the scattering strength
follows HeTh $>$ KrTh $>$ ZrTh $>$ ITh $>$ XeTh based on the mass difference.
This striking difference in the zirconium phonon scattering is due to the
local chemical environment changes. Zirconium is an electropositive element
with valency similar to thorium and, therefore, can bond with the oxygen
atoms, thus creating less force constant variance compared to iodine, an
electronegative element, noble gas helium, xenon, and krypton. These results
can serve as the benchmark for the analytical models and help the engineering-
scale modeling effort for nuclear design.

Suggested keywords

††preprint: APS/123-QED

##  I Introduction

Thorium dioxide (ThO2) is expected to play a vital role as a material for
advanced nuclear fuel cycles in future energy needs and is considered a better
and safer alternative to the currently used uranium dioxide [1]. Since the
thermal conductivity of nuclear fuel determines the efficiency and safety of a
nuclear reactor, it is imperative to understand the thermal transport of ThO2
in all possible reactor conditions [2]. It is well-known that point defects
are ubiquitous in irradiated materials and can scatter phonons, causing a
significant reduction in the material’s thermal conductivity, especially at
low temperatures or high defect concentrations [3]. Enormous efforts are being
undertaken to accurately model and ultimately control the thermal transport in
nuclear fuels. A detailed description of the new computational and
experimental tools for the same is detailed in a recent review article [2].
However, a fundamental understanding of how different fission products scatter
the phonons in nuclear fuels is still lacking. Therefore, an accurate
prediction of the phonon-point defect scattering rate by first-principles
calculations is imperative to shed light not only on the fundamental physics
of phonon transport in irradiated fuels but also in modeling the thermal
transport of ThO2 during regular operation and accident conditions in a
nuclear reactor.

The scattering of phonons by point defects is described as a perturbation to
harmonic lattice dynamics. The responsible mechanisms include the mass
difference, harmonic force constant changes, and strain field surrounding the
defect, often captured by ionic radius difference [4, 5]. So far, the widely
used description of the phonon-point defect scattering is based on the seminal
work by Klemens [6]. In this model, Klemens derived the analytical expression
for the phonon elastic scattering cross section by point defects,
dislocations, and the grain boundary using second-order perturbation theory.
However, the Klemens model is restricted to small perturbations due to mass
variance, and the equations are derived assuming a single atom unit cell and
linear phonon dispersion [7, 8, 9, 10]. Another popular mass difference model
used to describe the phonon scattering by the point defect utilizing the
realistic phonon dispersion was implemented by Tamura [11]. The main
difference between the Klemens and Tamura models is that the mass difference
term in the Tamura model is weighted by the eigenvector components
corresponding to the atom in the incident and final vibrational modes. The
Born approximation is another standard method that treats the defect in the
lowest order perturbation theory [12]. These models become questionable for
point defects with large force constant variance and structural changes [13]
(i.e., for conditions beyond the low-order perturbation regime). Hence,
defining the phonon-defect scattering by a nonperturbative method that
includes the scattering due to large variance in force constants becomes
vital. Recently, an exact calculation of the phonon-defect scattering rate
using Green’s function approach, which is nonperturbative and includes the
interatomic force constants (IFCs), has been proposed [14, 15]. To date, this
method has been successfully used to describe the phonon scattering by point
defects in various materials, such as diamond [16], boron arsenide [17],
silicon carbide [18], gallium nitride [19], graphene [20], and indium nitride
[13]. These calculations have already identified shortcomings of the
simplified perturbative methods. They have been instrumental in providing
fundamental insights into the phonon-defect scattering of the respective
materials.

In ThO2, where there is a colossal mass difference between the constituent
atoms and the possibility of sizeable interatomic force constant differences
with point defects, it becomes essential to describe the phonon scattering by
point defect with Green’s functions methodology. To our knowledge, Green’s
function-based methods have not been employed to study the thermal transport
in any actinide materials. In the literature, a few first-principles
investigations of the thermal conductivity by iteratively solving the Pierels-
Boltzmann transport equation of the defect-free ThO2 are available [21, 22,
23, 24]. In 2018, Liu et al. [21], using the generalized gradient
approximation pseudopotential, computed the second-order force constant from
the finite difference method and predicted the thermal conductivity of ThO2 at
300 K to be 12.4 W/mK, which was underpredicted by $\approx$40% in comparison
with experimental results (20.4 W/mK at 300 K) for a single crystal [25]. The
following year, using the generalized gradient approximation pseudopotential,
Malakkal et al. [22] reported a thermal conductivity of 15.4 W/mK at 300 K
(underpredicted by $\approx$25%) by computing the second-order force constant
using the density functional perturbation theory. Using the local density
approximation (LDA) pseudopotential, Jin et al. [23] predicted a thermal
conductivity of 18.6 W/mK (underpredicted by $\approx$9%) at room temperature.
Recently, Xiao et al. [24] calculated the phonon lifetimes and thermal
conductivity of ThO2 with the strongly constrained and approximately normed
functional [26] with an excellent agreement to experiments, indicating that,
over the years, researchers were able to provide a good description of the
thermal transport in the pristine ThO2 using the first-principles simulations
by refining the accuracy of calculated interatomic forces.

In practical application as nuclear fuel, a range of defects are produced in
ThO2 and from the fundamental physics perspective the role of these defects on
phonon scattering needs to be clearly understood [2]. Multiple researchers
have already applied various computational techniques to understand the effect
of point defects in the thermal transport of ThO2. For instance, Copper et al.
[27] reported the degradation of thermal conductivity due to the nonuniform
cation lattice of (UxTh1-x)O2 solid solutions, using the non-equilibrium
molecular dynamics method (NEMD), from 300 to 2000 K. Park et al. [28] used
reverse non-equilibrium molecular dynamics (r-NEMD) to investigate the effect
of vacancy and uranium substitutional defects in the thermal transport of
ThO2. Rahman et al. [29] studied the dependence of thermal conductivity on
fission-generated products (Xenon and Krypton) and vacancies in ThO2 within
300-1500 K using NEMD simulations. Recently, Jin et al. [30], using the NEMD
simulations, reported the phonon-defect scattering cross section for a range
of defects to be used with reduced order analytical models. In all these
studies, the quantum effects are ignored. The phonons follow the Boltzmann
distribution instead of the physically correct Bose-Einstein distribution [31]
and, therefore, are restricted to temperatures above the Debye temperature.
But the phonon-point defect is more dominant at low temperatures; thus, the
molecular dynamics (MD) based simulations are not ideal for describing phonon-
point defect scattering. The other challenge for the MD-based simulation is
that the results depend on the accuracy of the interatomic potential. The work
by Jin et al. [23] suggests that the empirical potential used previously needs
to be further optimized for the robust prediction of thermal conductivity of
ThO2 both in perfect crystals and in the presence of complex defects. Apart
from the MD simulations, Ryan et al. [5] studied the effect of point defects
on the thermal conductivity of ThO2 using the Klemens model for phonon
relaxation times that result from the change in mass and induced lattice
strain associated with point defects. This model computes the difference in
force constants and atomic radii using MD simulation in ThO2. Again, this
method is also limited to small perturbations and does not consider the
complete scattering to all orders.

Therefore, for the first time, the phonon scattering due to point defects and
their impact on the thermal conductivity of ThO2 is described by an ab initio
calculation utilizing the Green-function T-matrix method. In the current
study, we focus on the thermal conductivity degradation in ThO2 due to a
selected set of point defects, including the substitutions of transmutation
and radioactive decay product helium (HeTh) and fission products, such as
krypton (KrTh), zirconium (ZrTh), iodine (ITh), and xenon (XeTh) in thorium
sites, along with other point defects, such as the single vacancy of thorium
(VTh) and oxygen (VO) and three configurations of Schottky defects (SD100,
SD110, SD111).

##  II Methodology

The lattice thermal conductivity ($k_{L}$) of ThO2 with dilute point defects
sites is calculated by summing over contributions from all phonons as

$\displaystyle{k_{L}={\frac{1}{3\Omega}\sum_{j\textbf{q}}c_{j\textbf{q}}v^{2}_{j\textbf{q}}\tau_{j\textbf{q}}},}$
(1)

where $\Omega$ is unit cell volume, $c_{j\textbf{q}}$, $v_{j\textbf{q}}$,and
$\tau_{j\textbf{q}}$ are the specific heat capacity, group velocity, and
relaxation time for the phonon in mode $j\textbf{q}$, respectively, and the
$j$ and q represent the phonon branch index and wave vector, respectively. The
inverse of the phonon relaxation time $\tau^{-1}_{j\textbf{q}}$ is called the
phonon scattering rate, and the total phonon scattering rate can be written as
the sum of the contribution from different scattering mechanisms; in this
work, we have considered contributions from three major scattering phenomena,
including the intrinsic three-phonon scattering
($\tau^{-1}_{j\textbf{q},anh}$) [32], phonon-isotope scattering
($\tau^{-1}_{j\textbf{q},iso}$), and phonon-defect scattering
($\tau^{-1}_{j\textbf{q},def}$) [14]:

$\displaystyle{\tau^{-1}_{j\textbf{q}}=\tau^{-1}_{j\textbf{q},anh}+\tau^{-1}_{j\textbf{q},iso}+\tau^{-1}_{j\textbf{q},def}.}$
(2)

For brevity, the equation governing the scattering rate for the intrinsic
anharmonic scattering and the isotopic scattering rate is not provided here
and can be found in Ref. [32] and [11], respectively. However, we provide a
detailed discussion on the phonon-defect scattering. The phonon-defect
scattering is computed using Green’s function methodology using T matrix
[14].The advantages of using this approach is that, firstly, it captures the
scattering due to variance in both mass and force constants, with the latter
incorporating the structural relaxation effects; secondly, it also accounts
for complete scattering to all orders to capture the effect of resonances that
would otherwise be missing in the lower order perturbative methods. This
approach is valid under the assumption of random and dilute concentration of
defects, which means the phonon scatters independently from different defects
at a different location. The scattering cross-section $\sigma$, an essential
parameter in the phonon-defect scattering using Green’s function approach, is
explained in detail in Ref. [15] and in the supplementary materials of the
Ref. [18]. The scattering rate due to all possible elastic scatterings of
phonons processes by a random and dilute distribution is given by Equation
(3):

$\displaystyle{\tau^{-1}_{j\textbf{q},def}={\pi\chi_{def}\frac{\Omega}{V_{def}}}{\frac{1}{{\omega_{j\textbf{{q}}}}}}\sum_{j^{{}^{\prime}}\textbf{q}^{{}^{\prime}}}{{|\bra{j^{{}^{\prime}}\textbf{q}^{{}^{\prime}}}\textbf{T}\ket{j\textbf{q}}|}^{2}}}{\delta(\omega^{2}_{j^{{}^{\prime}}\textbf{q}^{{}^{\prime}}}{-}\omega^{2}_{j\textbf{q}}),}$
(3)

where $\chi_{def}$ is the number fraction of defects, $V_{def}$ is the volume
of a defected supercell, $\ket{j\textbf{q}}$ and
$\ket{j^{{}^{\prime}}\textbf{q}^{{}^{\prime}}}$ represent the incident and
scattered phonon eigenstates, respectively, $\omega$ is the phonon angular
frequency, and T is the matrix associated with the point defect scatterer. The
T matrix is defined in terms of the perturbation matrix ${(\textbf{V})}$ and
the retarded green’s function [15] for the perfect structure
${(\textbf{g}^{+})}$ as:

$\displaystyle{\textbf{T}=(\textbf{I}-\textbf{V}\textbf{g}^{+})^{-1}\textbf{V},}$
(4)

where ${(\textbf{I})}$ the identity matrix. Here, the perturbation matrix
includes the mass ${(\textbf{V}_{M})}$ and harmonic force constant
${(\textbf{V}_{K})}$ differences between the perfect and defective systems, as
shown below:

$\displaystyle{\textbf{V}=(\textbf{V}_{M}+\textbf{V}_{K}).}$ (5)

The mass term is a diagonal matrix, with nonzero matrix elements corresponding
to defect site given as
${{\textbf{V}_{M,i,i}}}$=$-{(\frac{{\text{M}^{{}^{\prime}}_{i}}-\text{M}_{i}}{\text{M}_{i}})\omega^{2}}$,
where $\text{M}^{{}^{\prime}}_{i}$ is the mass of the defect atom,
$\text{M}_{i}$ is the mass of the original atom at the $i^{th}$ site, and
$\omega$ is the angular frequency of the incoming phonon [18]. Similarly, the
presence of a point defect will also change the harmonic force constant around
the defect, and this is included in the $\textbf{V}_{\text{K}}$ matrix, which
is computed from the difference between the IFCs for the defective
$\text{K}^{{}^{\prime}}$ and perfect K structures as:

$\displaystyle{\textbf{V}_{K,i\alpha,k\beta}={\frac{{K^{{}^{\prime}}_{i\alpha,k\beta}}-K_{i\alpha,k\beta}}{\sqrt{M_{i}M_{k}}}}.}$
(6)

Here ${i}$ and ${k}$ are atom indices and ${\alpha}$ and ${\beta}$ are the
cartesian axes. Born approximation is achieved by Taylor series expansion of
Eq. 4 and approximation $\textbf{T}\approx\textbf{V}$.

It is worth noting that the ${\textbf{V}_{M}}$ is proportional to the square
of phonon frequency, hence are smaller near the zero frequency. Whereas
${\textbf{V}_{K}}$ is independent of the phonon frequency, and therefore, a
significant changes in the harmonic force constants cannot be considered small
at any frequency. Moreover, the IFCs are calculated at zero Kelvin
temperature, hence, the dependence of temperature on the scattering rate is
not included. But as the phonon defect scattering rates are frequency
dependent, and the phonon population are both frequency and temperature
dependent, the scattering rate indirectly influence the thermal conductivity
as a function of temperature [13].

In this approach, the perturbation matrix is defined in real space, and to
limit the size of the matrix, a cutoff radius ${r_{cut}}$ needs to be enforced
for the structural distortion. A large ${r_{cut}}$ corresponding to the fifth
neighbour shell in ThO2 is chosen for all the defects. Past this cutoff the
IFC changes are found to be negligible. Within ${r_{cut}}$, the changes in
IFCs up to the second nearest neighbors are considered. ${\textbf{V}_{K}}$
must obey the translational invariance rule, and when computed numerically,
this symmetry must be enforced. This can be done by adding a minor correction
to the harmonic IFCs of the perturbed system following a procedure outlined in
[18].

We simplified calculation represented by Equation (3) with the use of Optical
theorem. The total elastic scattering of phonon by defects can be efficiently
computed from the imaginary part of the diagonal element of the T-matrix, as
given as:

$\displaystyle{\tau^{-1}_{j\textbf{q},def}={{-}\chi_{def}\frac{\Omega}{V_{def}}}{\frac{1}{{\omega_{j\textbf{q}}}}}{\text{Im}\\{{{\bra{j\textbf{q}}\textbf{T}\ket{j\textbf{q}}}\\}}}.}$
(7)

When using optical theorem, Born approximation requires inclusion of higher
order terms in Taylor series expansion
$\textbf{T}\approx\textbf{V}+\textbf{V}\textbf{g}^{+}\textbf{V}$.

##  III Simulation details

All the density functional theory [33] calculations in this work were
performed using the projector-augmented-wave [34] method implemented in the
plane wave software package, the Vienna Ab initio Simulation Package (VASP)
[35], with the LDA [36] pseudopotential for the exchange and correlation. We
performed the geometry optimization of ThO2 (space group $Fm\bar{3}m$) on a
primitive unit cell by minimizing the total energy with respect to changes in
both cell parameters and atom positions using the conjugate gradient method.
We obtained the energy convergence of ThO2, using an electron wave vector grid
and plane wave energy cutoff of a 12 × 12 × 12 mesh and 550 eV, respectively.
The criteria for the electronic energy convergence were set at $10^{-8}$ eV.
The lattice parameter for the obtained relaxed structure of ThO2 at 0 K is
5.529 Å, which is consistent with the previously reported value [23] for the
LDA pseudopotentials from VASP and is comparable with the experimental lattice
constant value of 5.60 Å [37]. A 5 × 5 × 5 supercell with 375 atoms and
kpoints of 2 × 2 × 2 were used to evaluate the harmonic force constant using
the finite displacement method as implemented in PHONOPY [38]. The third-order
force constants (anharmonic force constants) were also performed on 5 × 5 × 5
supercells at the gamma point using the Thirdorder.py [39], and the force
cutoff distance was set to the fifth nearest neighboring atoms.In addition,
ThO2 being a polar material, the non-analytical contribution was considered,
and the Born charges and dielectric constant required to evaluate the non-
analytical correction [40] were calculated using the density functional
perturbation theory [41]. The phonon-phonon scattering processes are evaluated
using Fermi’s golden rule from the cubic force constants. Finally, $k_{L}$ was
calculated using the iterative solutions of the BTE as implemented in AlmaBTE
[42]. All the $k_{L}$ presented in this work are the fully iterative solutions
of the Peierls-equation. The converged $k_{L}$ values were obtained using a
cubic force constant, considering the fifth nearest-neighbor interaction and
the number of grid planes along each axis in the reciprocal space for solving
the BTE was 24 × 24 × 24.

All point defects considered in this work were created in a 5 × 5 × 5
supercell. For vacancy defects, we removed a Th and O atom from their lattice
position. We built the substitutional defects by replacing on Th atom by the
corresponding fission products. The Schottky defects were formed by removing a
unit of ThO2 in the three-bound configuration of the SD, the SD100, SD110, and
SD111 [43, 44]. Each defected supercell was relaxed by keeping the cell volume
constant. For the substitutional defects, a slight displacement of the defect
atom or its nearest neighbors is introduced before the relaxation to avoid the
system getting trapped in a saddle point of the potential energy surface. We
confirmed the dynamic stability of the relaxed structures by ensuring that the
phonon band structure did not have any modes with imaginary frequency. The
IFCs of the system with point defects were calculated using the finite-
displacement method. The space group symmetry for VTh and VO was face centered
cubic, while for the HeTh, KrTh, ZrTh, ITh, and XeTh substitution, the space
group symmetry was orthorhombic, and the Schottky defects were orthorhombic
for SD100 and SD110 and trigonal for SD111 using the spglib implemented in the
PHONOPY code. Finally, the phonon-defect scattering rates were calculated on a
uniformly spaced q-points mesh of 24 × 24 × 24 with an 18 × 18 × 18 grid for
the Green’s function, using the tetrahedron method to integrate over the
Brillouin zone. Finally, to understand the charge transfer of the substituted
defects to the neighbouring atoms, the charge density difference plot was
obtained using the VESTA software package [45] using the expression:
${\Delta\rho}_{\text{diff}}$ = ${\rho[\text{ThO}_{2}}+\text{defect}]$ \-
${\rho[\text{ThO}_{2}]}$ \- ${\rho[\text{defect}]}$ where
${\rho[\text{ThO}_{2}}+\text{defect}]$, ${\rho[\text{ThO}_{2}]}$, and
${\rho[\text{defect}]}$ are the charge density of ThO2 with point defect, the
pristine ThO2, and isolated atom corresponding to that particular defect,
respectively.

##  IV Results and Discussion

We plotted the scattering rates ($\tau^{-1}$) of each point defect type at a
1% concentration as a function of phonon frequency $(\omega)$ to provide
insights into the phonon scattering by various defects in ThO2. Fig. 1
compares the $\tau^{-1}$ due to vacancy point defects to the anharmonic
scattering rate of a pristine ThO2 at 300 K. The analysis of the phonon-defect
scattering rate suggests that the thorium vacancy scatters the phonons more
strongly than an oxygen vacancy. In vacancy defects, where the host atom is
absent, the contribution to $\tau^{-1}$ comes from the IFC differences. To
understand the extent of the force constant changes due to single vacancies of
thorium and oxygen, we computed the Frobenius norm ratio of the IFC matrix.
The values of the Frobenius norm ratio of the IFC matrix for the thorium
vacancy (0.17) and oxygen vacancy (0.05) defects suggest that a thorium
vacancy causes significantly more changes to the force constants than the
oxygen vacancy, which explains why the thorium vacancy scatters the phonons
more than the oxygen vacancy.

Figure 1:  Comparison of the $\tau^{-1}$ of the VTh (solid orange dots) and VO
(solid green dots) as a function of phonon frequency $(\omega)$ from the
Green-function T-matrix method compared to the $\tau^{-1}$ predicted by
Klemens model (anharmonic scattering [black dashed lines] and the defect
scattering [blue dashed lines]) and Born approximation for the VTh (solid cyan
dots) and the VO (solid red dots). Solid blue dots represent the anharmonic
scattering of the pristine ThO2 due to three phonon interactions at 300 K.

Additionally, Fig. 1 also compares the frequency trends suggested by Klemens
model [6] for anharmonic scattering and point defects. According to the
Klemens model, the anharmonic scattering rate follows an $\omega^{2}$
dependence (black dashed lines) and the defect scattering has a $\omega^{4}$
dependence (blue dashed lines). The latter is generally referred to as
Rayleigh scattering. Klemens model captures the anharmonic scattering in ThO2
in agreement with the ab initio predictions. In contrast, Klemens expression
for the point defect scattering, while correctly captures the low frequency
trends is unable to reproduce the high frequency phonons, especially the
optical modes which appear to have much weaker dependence on frequency. One
can argue that this is due to oversimplified description of phonon density of
states employed in the Klemens expression not valid for high frequency modes.
Nevertheless, the Tamura model, which uses the exact phonon dispersions is
unable to capture high freqeuncy trends either [11]. In fact Tamura model
makes the situation even worse, as the density of phonons increases at high
frequency, it predicts scattering rates that are even larger than one expects
from Klemens expressions [4, 5]. The high frequency overestimates of Klemens
and Tamura models is in fact due to their consideration of only mass term
which has strong frequency dependence, whereas the perturbation due to IFC is
frequency independent, suggesting that IFC term dominates at high frequencies.

Figure 2:  Comparison of the $\tau^{-1}$ of the VTh (solid orange dots) and VO
(solid green dots) as a function of phonon frequency $(\omega)$ from the
Green-function T-matrix method compared to Born approximation for the VTh
(solid cyan dots) and the VO (solid red dots). Solid blue dots represent the
anharmonic scattering of the pristine ThO2 due to three phonon interactions at
300 K.

Compared to the Green-function T-matrix methods, our calculations show that
the single-atom vacancy scattering centers are not well described by the Born
approximation (Fig 2), with a significant underestimation of the $\tau^{-1}$
at low phonon frequencies and overestimation at high phonon frequencies for
both the thorium and oxygen vacancies. It is worth noting that there are
prominent peaks observed at 2 Thz for the thorium vacancy and 3 THz for oxygen
vacancy. Further modal analysis reveals that these peaks are associated with
longitudinal acoustic (LA) wave and appear at the location where transverse
acoustic (TA) mode have large density of states. Under Born approximation the
scattering of LA modes to TA is not possible as the momentum is not conserved.
Previously such behavior in silicon carbide and gallium arsenide has been
attributed to resonant scattering [18, 46]. Green’s functions approach
considers full scattering to all orders, thus capturing the effect that leads
to enhanced scattering of phonon due to emergence of localized mode.

Apart from vacancy defects, neutron irradiation will introduce substitutional
defects. To discern the phonon scattering due to the substitutional atoms in
ThO2, we considered some common fission products such as Kr, Zr, I, and Xe as
well as He in the thorium sites. Our reason for only considering thorium sites
and not the oxygen sites is that the Th atom’s vibration governs the acoustic
modes and contributes to $\approx$70% of the thermal conductivity [22]. Hence,
defects on the Th sites, particularly those that have a large mass or IFC
variance, can degrade the thermal conductivity more significantly than those
of the light O atoms. Unlike the vacancy defects, the scattering from
substitutional atoms arises from the mass difference between the
substitutional and host atom and the IFCs variations around the substituted
site due to the structural relaxation and modified chemical environment.

Figure 3:  Scattering rates of phonons as a function of phonon frequency
$(\omega)$ due to substitutional point defects of HeTh (crimson red solid
dots), KrTh (purple solid dots), ZrTh (brown solid dots), ITh (pink solid
dots), and XeTh (gray solid dots) in thorium sites compared to the VTh (solid
orange dots) and pristine ThO2 due to three phonon interactions at 300 K
(solid blue dots).

The previous study using same methodology has shown that in general the
substitutional atoms scatter phonons less effectively than vacancies [13]. Our
ThO2 results shows that a substitution by a significantly lighter He atom
results in phonon scattering with a strength comparable to thorium vacancy, to
an extent that they overlap. The primary reason for He substitution to
strongly scatter phonons is because He as a light chemically inert atom
induces a force constant variance very similar to that of the thorium
vacancy(supplementary information [SI]: Table 1). As expected, a similar
resonant peak as seen in thorium vacancy is also observed for the helium
substitution at the same frequency.

Among the other substitutional atoms, intuitively based on mass differences,
one expects the scattering strength to be HeTh $>$ KrTh $>$ ZrTh $>$ ITh $>$
XeTh. However, contrary to this intuition, the predicted scattering strengths
of the fission products were HeTh $>$ KrTh $>$ ITh $>$ XeTh $>$ ZrTh, where
the ZrTh substitution scattered phonons the least. To understand why the ZrTh
substitution showed minimal phonon scattering, we computed the Frobenius norm
ratio of the IFC matrix (SI: Table 1) and found that the force constant
perturbations were the lowest when Zr was in thorium sites. This can be
attributed to the fact that, unlike noble gas elements (He, Xe and Kr) and a
electronegative element (I), Zr is a electropositive element, and Zr in the
cation lattice site can donate electrons to the oxygen and form bonds. The
charge transfer in the substitutional defects was quantified (SI: Fig. 10) to
provide insights into the redistribution of charge and bond strength. In the
case of a Zr substitution, the charge transfer takes place from the Zr atom
out towards the surrounding O atoms. Also, the charge accumulated around the
Zr atom is negligible since Th and Zr have similar ionic charges. This
observation clearly indicates that the local chemical modification can
significantly impact the force constant perturbations and, thus, phonon
scattering.

Finally as an example of larger defects, we investigated the role of Schottky
defects in ThO2 in scattering the phonons, as depicted in Fig. 4. Among the
three different Schottky defect configurations, the scattering strengths are
SD111 $>$ SD110 $>$ SD100, which indicates that ThO2 with an SD111 defect
scatters phonons the most.

Figure 4:  Scattering rates of phonons by the three different configurations
of the Schottky defects (SD100 [earth green solid dots], SD110 [light blue
solid dots], SD111 [Indian red solid dots]) in ThO2 compared with the phonon
scattering due to three phonon interactions at 300 K. Figure 5:  The
calculated $k_{L}$ of ThO2 with vacancy defects (VTh (0.1$\%$ and 0.25$\%$) ,
and VO (0.1$\%$)) at a lower temperature range of 100–300 K. The experiments
are from Deskins et al. [47] corresponding to ThO2 single crystals irradiated
with 2 MeV protons to a dose of 0.004 displacements per atom at room
temperature.

Using the calculated scattering rates, including the anharmonic, isotropic,
and phonon-defect scattering, the $k_{L}$ of ThO2 was computed using an
iterative solution of the BTE. Figure 5 shows the $k_{L}$ of pristine ThO2 and
with vacancy defects as a function of temperature. Thorium vacancy caused the
most significant reduction in thermal conductivity, consistent with the
observation made with the phonon-point defect scattering rate (Fig. 1). The
oxygen vacancy at 0.001 concentration caused a marginal decrease in thermal
conductivity. The theoretical results are compared with the recent
experimental measurement done on a proton-irradiated single crystal of ThO2 at
a small dose of 0.004 displacements per atom [47]. At these low dose and room
temperature irradiation conditions, the thorium vacancies and interstitial are
the most likely defects, as predicted by the rate theory calculation presented
in Dennett et al. [3]. However, current T-matrix formalism can not be applied
to interstitial defects located at octahedral sites [48], therefore to account
for this type of defects some assumptions need to be made. The scattering
cross sections of the vacancies and interstitials are similar based on the
NEMD simulations by Jin et al. [30]. Based on these arguments, we estimate
that 0.004 dpa corresponds to an equivalent concentrations of ThV to be
$\approx$ 0.25 $\%$. Using the concentration of 0.25 $\%$ VTh yields an
excellent agreement between the predicted and the experimental thermal
conductivity.

Figure 6:  The calculated $k_{L}$ of ThO2 with 0.1$\%$ fission products of
HeTh, KrTh, ZrTh, ITh, XeTh at a lower temperature range of 100–300 K.

Figure 6 shows the degradation of thermal conductivity due to the
substitutional defects. Among the considered elements, He with the largest
perturbation caused the largest reduction in thermal conductivity followed by
KrTh $>$ ITh $>$ XeTh $>$ ZrTh. Figure 7 shows the $k_{L}$ as a function of
defect concentration at 300 K. The results indicate that, for the point
defects to impact the thermal transport, the concentration has to be on the
order of 10-2 atomic $\%$. For brevity, the $k_{L}$ degradation for 1% defect
concentrations for all the vacancies and substitutional defects, the $k_{L}$
reduction due to the Schottky defects for two different concentrations of 0.1%
and 1%, and the thermal conductivity degradation at high temperatures are
provided in SI (Fig S1-S4). The results for the Schotkky defects suggest that,
among the three different configurations, the SD111 impacted the thermal
transport, followed by SD110 and SD100. All these observations shed light on
the effect of the point defects on the thermal transport of thorium dioxide
and help in understanding behavior of nuclear reactors [2].

Lastly, we should point out that experimental measurements of thermal
conductivity in irradiated ThO2 and (Th,U)O2 alloys provided evidence presence
of phonon resonant scattering as the temperature dependence deviate from
1/T-trend at low temperatures [47, 49]. Our phonon scattering rates shown in
Figs. 1 and 3 provide some indication of the resonances, which for other
materials have resulted in calculated conductivity departing from 1/T-trend
[46]. However, our calculated conductivity values reported Figs. 5 and 6 do
not reflect this. This observations suggests that further work is needed to
improve estimation of low temperature thermal conductivity.

Figure 7:  ThO2 $k_{L}$ variation with a defect concentration for VTh, VO,
HeTh, KrTh, ZrTh, ITh, XeTh, and SD111 at 300 K.

##  V Conclusion

For the first time, using a fully ab initio Green’s function T-matrix method,
we have presented an extensive analysis of the phonon-defect scattering rate
of ThO2 with point defects. The phonon-point defect scattering rates were
unified with the Peierls-Boltzmann transport equation to determine the effects
of defects on phonon thermal conductivity in ThO2. The results revealed that,
among all the intrinsic defects considered in this work, the thorium vacancy
and helium substitution in the thorium site caused the most significant
reduction in the thermal conductivity of ThO2. The similar behaviour of the
thorium vacancy and helium substitution was due to the comparable force
constant difference induced by these two defects. Also, a large peak, which is
a signature of the resonant scattering was observed in the case of thorium
vacancy and helium substitution. Among the substitutional defects, the
zirconium substitution had the smallest reduction in thermal conductivity
despite having a significant mass difference compared to iodine and xenon. The
primary reason for the relatively weaker scattering by zirconium substitution
is its similar valency to thorium, thus enabling the bond formation with
oxygen, inducing the least force constant differences. Unlike zirconium, the
chemical environment due to helium, krypton, iodine, and xenon creates a
larger force constant difference, leading to higher scattering rates. This
study revealed that the modification of the local chemical environment and the
associated structure and force variations caused by substitutional defects
should be considered while studying the phonon-defect scattering. This work
enhances our understanding of the phonon-defect scattering due to various
fission products in ThO2, which will serve as the benchmark for the analytical
models and help improve the engineering-scale modeling predictions for nuclear
reactor design.

##  VI ACKNOWLEDGMENTS

This work was supported by the Center for Thermal Energy Transport under
Irradiation (TETI), an Energy Frontier Research Center (EFRC) funded by the
U.S. Department of Energy, Office of Science, and Office of Basic Energy
Sciences. The authors also acknowledge that this research made use of the
resources of the High Performance Computing Center at Idaho National
Laboratory, which is supported by the Office of Nuclear Energy of the U.S.
Department of Energy and the Nuclear Science User Facilities under Contract
No. DE-AC07-05ID14517. This manuscript has been authored by Battelle Energy
Alliance, LLC under Contract No. DE-AC07-05ID14517 with the U.S. Department of
Energy. The United States Government retains and the publisher, by accepting
the article for publication, acknowledges that the U.S. Government retains a
nonexclusive, paid-up, irrevocable, world-wide license to publish or reproduce
the published form of this manuscript, or allow others to do so, for U.S.
Government purposes.

## References

  * iae [2005] _Thorium Fuel Cycle - Potential Benefits and Challenges_, TECDOC Series No. 1450 (INTERNATIONAL ATOMIC ENERGY AGENCY, Vienna, 2005).
  * Hurley  _et al._ [2022] D. H. Hurley, A. El-Azab, M. S. Bryan, M. W. D. Cooper, C. A. Dennett, K. Gofryk, L. He, M. Khafizov, G. H. Lander, M. E. Manley, J. M. Mann, C. A. Marianetti, K. Rickert, F. A. Selim, M. R. Tonks, and J. P. Wharry, Thermal energy transport in oxide nuclear fuel, Chemical Reviews 122, 3711 (2022).
  * Dennett  _et al._ [2021] C. A. Dennett, W. R. Deskins, M. Khafizov, Z. Hua, A. Khanolkar, K. Bawane, L. Fu, J. M. Mann, C. A. Marianetti, L. He, D. H. Hurley, and A. El-Azab, An integrated experimental and computational investigation of defect and microstructural effects on thermal transport in thorium dioxide, Acta Materialia 213, 116934 (2021).
  * Gurunathan  _et al._ [2020] R. Gurunathan, R. Hanus, M. Dylla, A. Katre, and G. J. Snyder, Analytical models of phonon–point-defect scattering, Phys. Rev. Applied 13, 034011 (2020).
  * Deskins  _et al._ [2021] W. R. Deskins, A. Hamed, T. Kumagai, C. A. Dennett, J. Peng, M. Khafizov, D. Hurley, and A. El-Azab, Thermal conductivity of ThO2: Effect of point defect disorder, Journal of Applied Physics 129, 075102 (2021).
  * Klemens [1955] P. G. Klemens, The scattering of low-frequency lattice waves by static imperfections, Proceedings of the Physical Society. Section A 68, 1113 (1955).
  * Asheghi  _et al._ [2002] M. Asheghi, K. Kurabayashi, R. Kasnavi, and K. E. Goodson, Thermal conduction in doped single-crystal silicon films, Journal of Applied Physics 91, 5079 (2002).
  * Zou  _et al._ [2002] J. Zou, D. Kotchetkov, A. A. Balandin, D. I. Florescu, and F. H. Pollak, Thermal conductivity of GaN films: Effects of impurities and dislocations, Journal of Applied Physics 92, 2534 (2002).
  * Fletcher [1987] R. Fletcher, Scattering of phonons by dislocations in potassium, Phys. Rev. B 36, 3042 (1987).
  * Malekpour  _et al._ [2016] H. Malekpour, P. Ramnani, S. Srinivasan, G. Balasubramanian, D. L. Nika, A. Mulchandani, R. K. Lake, and A. A. Balandin, Thermal conductivity of graphene with defects induced by electron beam irradiation, Nanoscale 8, 14608 (2016).
  * Tamura [1983] S. Tamura, Isotope scattering of dispersive phonons in Ge, Phys. Rev. B 27, 858 (1983).
  * Sakurai [1967] J. Sakurai,  _Advanced Quantum Mechanics_, Always learning (Pearson Education, Incorporated, 1967).
  * Polanco and Lindsay [2018a] C. A. Polanco and L. Lindsay, Thermal conductivity of InN with point defects from first principles, Phys. Rev. B 98, 014306 (2018a).
  * Mingo  _et al._ [2010] N. Mingo, K. Esfarjani, D. A. Broido, and D. A. Stewart, Cluster scattering effects on phonon conduction in graphene, Phys. Rev. B 81, 045408 (2010).
  * Kundu  _et al._ [2011] A. Kundu, N. Mingo, D. A. Broido, and D. A. Stewart, Role of light and heavy embedded nanoparticles on the thermal conductivity of SiGe alloys, Phys. Rev. B 84, 125426 (2011).
  * Katcho  _et al._ [2014] N. A. Katcho, J. Carrete, W. Li, and N. Mingo, Effect of nitrogen and vacancy defects on the thermal conductivity of diamond: An ab initio green’s function approach, Phys. Rev. B 90, 094117 (2014).
  * Protik  _et al._ [2016] N. H. Protik, J. Carrete, N. A. Katcho, N. Mingo, and D. Broido, Ab initio study of the effect of vacancies on the thermal conductivity of boron arsenide, Phys. Rev. B 94, 045207 (2016).
  * Katre  _et al._ [2017] A. Katre, J. Carrete, B. Dongre, G. K. H. Madsen, and N. Mingo, Exceptionally strong phonon scattering by B substitution in cubic SiC, Phys. Rev. Lett. 119, 075902 (2017).
  * Katre  _et al._ [2018] A. Katre, J. Carrete, T. Wang, G. K. H. Madsen, and N. Mingo, Phonon transport unveils the prevalent point defects in GaN, Phys. Rev. Materials 2, 050602 (2018).
  * Polanco and Lindsay [2018b] C. A. Polanco and L. Lindsay, Ab initio phonon point defect scattering and thermal transport in graphene, Phys. Rev. B 97, 014303 (2018b).
  * Liu  _et al._ [2018] J. Liu, Z. Dai, X. Yang, Y. Zhao, and S. Meng, Lattice thermodynamic behavior in nuclear fuel ThO2 from first principles, Journal of Nuclear Materials 511, 11 (2018), special Section on ”18th International Conference on Fusion Reactor Materials”.
  * Malakkal  _et al._ [2019] L. Malakkal, A. Prasad, E. Jossou, J. Ranasinghe, B. Szpunar, L. Bichler, and J. Szpunar, Thermal conductivity of bulk and porous ThO2: Atomistic and experimental study, Journal of Alloys and Compounds 798, 507 (2019).
  * Jin  _et al._ [2021] M. Jin, M. Khafizov, C. Jiang, S. Zhou, C. A. Marianetti, M. S. Bryan, M. E. Manley, and D. H. Hurley, Assessment of empirical interatomic potential to predict thermal conductivity in ThO2 and UO2, Journal of Physics: Condensed Matter 33, 275402 (2021).
  * Xiao  _et al._ [2022] E. Xiao, H. Ma, M. S. Bryan, L. Fu, J. M. Mann, B. Winn, D. L. Abernathy, R. P. Hermann, A. R. Khanolkar, C. A. Dennett, D. H. Hurley, M. E. Manley, and C. A. Marianetti, Validating first-principles phonon lifetimes via inelastic neutron scattering, Phys. Rev. B 106, 144310 (2022).
  * Mann  _et al._ [2010] M. Mann, D. Thompson, K. Serivalsatit, T. M. Tritt, J. Ballato, and J. Kolis, Hydrothermal growth and thermal property characterization of ThO2 single crystals, Crystal Growth & Design 10, 2146 (2010).
  * Sun  _et al._ [2015] J. Sun, A. Ruzsinszky, and J. P. Perdew, Strongly constrained and appropriately normed semilocal density functional, Phys. Rev. Lett. 115, 036402 (2015).
  * Cooper  _et al._ [2015] M. Cooper, S. Middleburgh, and R. Grimes, Modelling the thermal conductivity of UxTh1-xO2 and UxPu1-xO2, Journal of Nuclear Materials 466, 29 (2015).
  * Park  _et al._ [2018] J. Park, E. B. Farfán, K. Mitchell, A. Resnick, C. Enriquez, and T. Yee, Sensitivity of thermal transport in thorium dioxide to defects, Journal of Nuclear Materials 504, 198 (2018).
  * Rahman  _et al._ [2020] M. Rahman, B. Szpunar, and J. Szpunar, Dependence of thermal conductivity on fission-product defects and vacancy concentration in thorium dioxide, Journal of Nuclear Materials 532, 152050 (2020).
  * Jin  _et al._ [2022] M. Jin, C. A. Dennett, D. H. Hurley, and M. Khafizov, Impact of small defects and dislocation loops on phonon scattering and thermal transport in ThO2, Journal of Nuclear Materials , 153758 (2022).
  * Zhou  _et al._ [2018] Y. Zhou, Z. Fan, G. Qin, J.-Y. Yang, T. Ouyang, and M. Hu, Methodology perspective of computing thermal transport in low-dimensional materials and nanostructures: The old and the new, ACS Omega 3, 3278 (2018).
  * Zim [2001] _Electrons and Phonons: The Theory of Transport Phenomena in Solids._, Oxford University Press (Oxford University Press, University of Bristol, 2001).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B 23, 5048 (1981).
  * Idiri  _et al._ [2004] M. Idiri, T. Le Bihan, S. Heathman, and J. Rebizant, Behavior of actinide dioxides under pressure: UO2 and ThO2, Phys. Rev. B 70, 014113 (2004).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scripta Materialia 108, 1 (2015).
  * Li  _et al._ [2014] W. Li, J. Carrete, N. A. Katcho, and N. Mingo, Shengbte: A solver of the boltzmann transport equation for phonons, Computer Physics Communications 185, 1747 (2014).
  * Yi  _et al._ [2016] W. Yi, S. Shun-Li, F. Huazhi, L. Zi-Kui, and C. L. Qing, First-principles calculations of lattice dynamics and thermal properties of polar solids, npj Computational Materials 2, 16006 (2016).
  * Baroni  _et al._ [2001] S. Baroni, S. De Gironcoli, A. Dal Corso, and P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Reviews of Modern Physics 73, 515–562 (2001).
  * Carrete  _et al._ [2017] J. Carrete, B. Vermeersch, A. Katre, A. van Roekeghem, T. Wang, G. K. Madsen, and N. Mingo, almaBTE : A solver of the space–time dependent Boltzmann transport equation for phonons in structured materials, Computer Physics Communications 220, 351 (2017).
  * Iwasawa  _et al._ [2009] M. Iwasawa, T. Ohnuma, Y. Chen, Y. Kaneta, H. Geng, A. Iwase, and M. Kinoshita, First-principles study on cerium ion behavior in irradiated cerium dioxide, Journal of Nuclear Materials 393, 321 (2009).
  * Dorado  _et al._ [2009] B. Dorado, M. Freyss, and G. Martin, GGA+U study of the incorporation of iodine in uranium dioxide, The European Physical Journal B 69, 203 (2009).
  * Momma and Izumi [2008] K. Momma and F. Izumi, VESTA: a three-dimensional visualization system for electronic and structural analysis, Journal of Applied Crystallography 41, 653 (2008).
  * Kundu  _et al._ [2019] A. Kundu, F. Otte, J. Carrete, P. Erhart, W. Li, N. Mingo, and G. K. H. Madsen, Effect of local chemistry and structure on thermal transport in doped gaas, Phys. Rev. Mater. 3, 094602 (2019).
  * Deskins  _et al._ [2022] W. R. Deskins, A. Khanolkar, S. Mazumder, C. A. Dennett, K. Bawane, Z. Hua, J. Ferrigno, L. He, J. M. Mann, M. Khafizov, D. H. Hurley, and A. El-Azab, A combined theoretical-experimental investigation of thermal transport in low-dose irradiated thorium dioxide, Acta Materialia 241, 118379 (2022).
  * Jiang  _et al._ [2022] C. Jiang, L. He, C. A. Dennett, M. Khafizov, J. M. Mann, and D. H. Hurley, Unraveling small-scale defects in irradiated tho2 using kinetic monte carlo simulations, Scripta Materialia 214, 114684 (2022).
  * Hua  _et al._ [2023] Z. Hua, S. Adnan, A. R. Khanolkar, K. Rickert, D. B. Turner, T. A. Prusnick, J. M. Mann, D. H. Hurley, M. Khafizov, and C. A. Dennett, Thermal conductivity suppression in uranium-doped thorium dioxide due to phonon resonant scattering (2023), arXiv:2303.01659 [cond-mat.mtrl-sci] .