First-principles study of UO$_2$ lattice thermal-conductivity: A simple description

TL;DR

This study models the lattice thermal conductivity of bulk UO₂ in its high-temperature paramagnetic state using a simplified first-principles approach, revealing dominant phonon modes and impurity effects with results aligning well with experimental data.

Contribution

It introduces a straightforward first-principles method to estimate UO₂'s thermal conductivity without complex magnetic or Hubbard-U corrections, emphasizing phonon lifetime analysis.

Findings

TA phonon branches dominate thermal conductivity

Optical phonons have shorter lifetimes and less contribution

Impurities slightly reduce thermal conductivity

Abstract

Modeling the high-$T$ paramagnetic state of bulk UO$_2$ by a non-spin-polarized calculation and neglecting the Hubbard-U correction for the $f$ electrons in U atoms, the lattice thermal conductivity of bulk UO$_2$ is investigated by the exact solution of the Boltzmann transport equation for the steady-state phonon distribution function. The results show that TA branches corresponding to U-atoms vibrations have the largest lifetimes and therefore have dominant role in thermal conductivity, while the optical branches corresponding mainly to O-atoms vibrations have the shortest lifetimes. Using this simple model, our results for the thermal conductivity show a very good agreement with the experiments. The calculations are repeated for bulk UO$_2$ with different U-235 concentrations of 3\%, 5\%, 7\%, and 20\%, and the results show a small decrease of thermal conductivity which arise from scattering of phonons by impurities.

Tables

Table 1. Table 1: Zero-pressure lattice constant, in Å italic-Å \AA , as a function of k 𝑘 k -mesh and plane-wave energy cutoff, in Ry, for wavefunction expansion. The cutoff for density expansion is taken as 10 times of E c subscript 𝐸 𝑐 E_{c} . The chosen parameters for k 𝑘 k -mesh and E c subscript 𝐸 𝑐 E_{c} are 8 × 8 × 8 8 8 8 8\times 8\times 8 and 60 Ry, respectively; and the corresponding lattice constant is 5.3968 Å italic-Å \AA .

			$E_c$
$k$-mesh	40	50	60	70	80
$2\times 2\times 2$	5.4007	5.4011	5.4010	5.4014	5.4015
$3\times 3\times 3$	5.3993	5.3994	5.3996	5.3999	5.3999
$4\times 4\times 4$	5.3937	5.3937	5.3939	5.3943	5.3942
$5\times 5\times 5$	5.3979	5.3985	5.3982	5.3987	5.3987
$6\times 6\times 6$	5.3956	5.3957	5.3959	5.3962	5.3965
$7\times 7\times 7$	5.3965	5.3965	5.3963	5.3967	5.3966
$8\times 8\times 8$	5.3965	5.3968	5.3968	5.3970	5.3969
$10\times 10\times 10$	5.3965	5.3963	5.3964	5.3971	5.3966

Table 2. Table 2: Cubic cutoff [# of shells/distance ( Å italic-Å \AA )] in terms of number of neighboring shells/distance, the resulting number of displaced structures, and the computed thermal conductivity.

cubic cutoff	# of disp. struct.	$\kappa$ (W/mK)
1/3.5	44	7.8728
2/4.0	108	7.2686
3/5.0	176	7.1736
4/5.4	208	7.2347

Equations

$$\nabla T\cdot {\mathbf{v}}_{\lambda }\frac{\partial f_{\lambda }}{\partial T}={\frac{\partial f_{\lambda }}{\partial t}|}_{\mathrm{scattering}},$$

$${\omega }_{\lambda }\pm {\omega }_{{\lambda }^{\prime }}={\omega }_{{\lambda }^{\prime \prime }},$$

$$\mathbf{q}\pm {\mathbf{q}}^{\prime }={\mathbf{q}}^{\prime \prime }+\mathbf{Q},$$

$${\tau }_{\lambda }^{(\alpha )}={\tau }_{\lambda }^0(1+{\Delta }_{\lambda }^{\alpha }),$$

$$\frac{1}{{\tau }_{\lambda }^0}=\frac{1}{N}\left(\sum _{{\lambda }^{\prime }{\lambda }^{\prime \prime }}^{+}{\Gamma }_{\lambda {\lambda }^{\prime }{\lambda }^{\prime \prime }}^{+}+\frac{1}{2}\sum _{{\lambda }^{\prime }{\lambda }^{\prime \prime }}^{-}{\Gamma }_{\lambda {\lambda }^{\prime }{\lambda }^{\prime \prime }}^{-}+\sum _{{\lambda }^{\prime }}{\Gamma }_{\lambda {\lambda }^{\prime }}\right),$$

$${\Delta }_{\lambda }^{\alpha }=\frac{1}{N}(\sum _{{\lambda }^{\prime }{\lambda }^{\prime \prime }}^{+}{\Gamma }_{\lambda {\lambda }^{\prime }{\lambda }^{\prime \prime }}^{+}({\xi }_{\lambda {\lambda }^{\prime \prime }}^{\alpha }{\tau }_{{\lambda }^{\prime \prime }}^{(\alpha )}-{\xi }_{\lambda {\lambda }^{\prime }}^{\alpha }{\tau }_{{\lambda }^{\prime }}^{(\alpha )})$$

$$+\frac{1}{2}\sum _{{\lambda }^{\prime }{\lambda }^{\prime \prime }}^{-}{\Gamma }_{\lambda {\lambda }^{\prime }{\lambda }^{\prime \prime }}^{-}({\xi }_{\lambda {\lambda }^{\prime \prime }}^{\alpha }{\tau }_{{\lambda }^{\prime \prime }}^{(\alpha )}+{\xi }_{\lambda {\lambda }^{\prime }}^{\alpha }{\tau }_{{\lambda }^{\prime }}^{(\alpha )})$$

$$+\sum _{{\lambda }^{\prime }}{\Gamma }_{\lambda {\lambda }^{\prime }}{\xi }_{\lambda {\lambda }^{\prime }}^{\alpha }{\tau }_{{\lambda }^{\prime }}^{(\alpha )})$$

$${\kappa }^{\alpha \beta }=\frac{1}{V{k}_B{T}^2}\sum _{\lambda }(\hslash {\omega }_{\lambda }{)}^2{f}_{\lambda }^0(f_{\lambda }^0+1)v_{\lambda }^{\alpha }{v}_{\lambda }^{\beta }{\tau }_{\lambda }^{(\beta )}.$$

Full text

Uranium dioxide; Lattice thermal conductivity; Boltzmann transport equation; Acoustic branch; Optical branch; Phonon lifetime; Density-functional perturbation theory

I Introduction

One of the vastly used fuels in nuclear power reactors is uranium-dioxide, UO2. The performance of a nuclear fuel is highly correlated to its thermal conductivity, and therefore studying the thermal conductivity of nuclear fuel and understanding the mechanisms behind it, is one of the most active fields of research in nuclear industry. Due to the fission processes of the uranium atoms, large amounts of heat are released, leading to large temperature gradients throughout the fuel rod. Having an efficient thermal conductivity, the generated temperature-gradients in the fuel are immediately balanced and the heat is easily extracted by the coolant so that the fuel system does not meet any safety problems because of any temperature-increase or thermal stresses.

Experiments have determined the crystal structure of UO2 as a 3k-order antiferromagnet (AFM) at $T<30^{\circ}$K, and paramagnetic at higher temperatures [1, 2]. The uranium atoms sit on the sites of an FCC structure with lattice constant $a=5.47\AA$, while the oxygen atoms are positioned at sites with $Pa\bar{3}$ symmetry[3]. Uranium-dioxide is electrically an insulator material (the so-called Mott insulator), and because of the localized partially-filled $f$-orbitals on U atoms, it is a strongly-correlated electron system. Theoretical description of electronic properties of such a system by ordinary density-functional theory (DFT)[4, 5] approximations usually lead to incorrect metallic behavior, and for a correct prediction, one has to somehow take into account the “localized” behavior of the $f$ electrons in U atoms. Two methods that are commonly used for this purpose, are the DFT+U[6, 7] and using orbital-dependent hybrid functional for the exchange-correlation (XC) part of the energy[8]. Since in the bulk UO2 there is no delocalized electrons, the thermal properties depends only on lattice dynamics of the system. The properties of lattice dynamics are characterized by phonon dispersion, phonon density of states (PhDOS), and anharmonic behaviors.

The lattice dynamics and thermal conductivity of bulk UO2 have been studied by many researchers both experimentally and theoretically. For example, Dolling and coworkers[9], using the neutron inelastic-scattering technique, were the first ones that experimentally determined the phonon dispersion and density-of-states curves of UO2; Godfrey and coworkers[10] have used a radial heat flow technique and measured the thermal conductivity of polycrystalline UO2 in the range -57∘ to 1100∘C; Goldsmith and coworker[11] using the laser-flash method have measured the thermal conductivity of porous stoichiometric and hyper-stoichiometric uranium dioxide in the temperature range of 670-1270∘K; Fink and coworkers in two separate works[12, 13] have reviewed the experimental data on thermodynamic and transport properties of solid and liquid UO2, and by analyzing the data have obtained consistent equations for the thermophysical properties. On the other hand, in theoretical investigations, Motoyama and coworkers[14] resorting to nonequilibrium classical molecular dynamics (MD) have calculated the thermal conductivity of UO2 pellet; Yamada et al[15], employing the partially ionic model within their MD simulations, have studied the molar specific heat and the thermal conductivity of UO2; Arima and coworkers[16] using the Born–Mayer–Huggins interatomic potential with the partially ionic model, have performed equilibrium MD simulations and investigated the thermal properties of UO2 and PuO2 between 300 and 2000∘K; Kaur et al[17] within DFT+U method in generalized gradient approximation (GGA)[18] and applying an external hydrostatic pressure of 7 GPa, have optimized the geometry of UO2 crystal in AFM 1k-order configuration, and then using the resulted lattice parameters have studied the thermal properties using the density functional perturbation theory (DFPT)[19] in conjunction with the quasi-harmonic approximation (QHA); Pang and coworkers[20], in a joint experimental and theoretical investigation, using high-resolution inelastic neutron scattering, have measured phonon lifetimes and dispersion of UO2 at 295 and 1200∘K, and then analyzing the calculated result of thermal conductivity within the relaxation-time approximation (RTA), have concluded that longitudinal optical (LO1) branch of phonons carries the largest amount of heat; Resnick and coworkers[21] have used an MD simulation to study the thermal transport of plutonium dioxide and uranium dioxide with point defects; Torres and coworker[22] had performed the solution of BTE on top of a DFT+U calculation with 1k-order AFM configuration for UO2, but because of the diversity of their work, they had not provided a satisfactory analysis of the problem.

One of the ways to investigate the heat transport in solids is using the Boltzmann transport equation (BTE). The earlier BTE studies on lattice thermal conductivity of UO2 were based on the relaxation-time approximation (RTA). Although the RTA gives good results in cases where the “umklapp” scattering processes are dominant, the exact solution of the BTE is needed whenever the “normal” processes dominate. In this work, we have obtained the exact solution of BTE for phonons’ distribution from which the lattice thermal conductivity is calculated. Our results show a very good agreement with the experiments.

In nature, uranium is found as U-238 (99.2739%), U-235 (0.7198%), and a very small amount of U-234 (0.0050%). Since in different nuclear fuels the relative abundances of U-238 and U-235 are different, the scattering rates by “impurities” differs from one fuel to other and this affects the thermal conductivity of the fuel. We have therefore repeated the procedure for the cases U-235 (3%), (5%), (7%), and (20%), and the results showed a small decrease of thermal conductivity by increasing the impurity level.

The organization of this paper is as follows. In Section II we present the computational details; in Section III the results are presented and discussed; Section IV concludes this work. Finally, in A the convergency issues of our calculations are detailed.

II Computational details

Since the simulation of the high-$T$ paramagnetic state of UO2 system with randomly oriented magnetic moments requires a very large supercell, it is computationally very prohibitive and it is common to model it with a 1k-order AFM. However, in this work, we model the paramagnetic state of bulk UO2 by a non-spin-polarized calculation, and neglect the Hubbard-U correction for the localized orbitals in the study of lattice dynamics. It has already been observed that although not using the Hubbard correction in the DFT calculations gives incorrect metallic ground state for UO2, but surprisingly the phonon properties are comparable with experiment [23].

II.1 Geometry and Harmonic lattice dynamics

For the electronic structure calculations, we have used the first-principles DFT method as implemented in the Quantum-ESPRESSO code package[24] and used Ultra-soft pseudopotentials with the rev-PBE XC-functional[25]. The reason for choosing the rev-PBE method is that it gives the best lattice constant for this system. Comparing the results for functionals LDA[5], PBE[18], PW91[26], and rev-PBE; we obtained the values of 5.27, 5.36, 5.36, and 5.39 Åfor the lattice constants, respectively. The kinetic-energy cutoffs for the expansion of the Kohn-Sham (KS)[5] orbitals and charge densities were chosen as 60 and 600 Ry, respectively. For the integrations over the Brillouin-zone (BZ), a $\Gamma$-centered $8\times 8\times 8$ grid were used. The optimization of the geometry was performed with a maximum $10^{-5}$Ry/au of residual force on each atom. The convergency tests with respect to the parameters are detailed in A.

To calculate the second-order (2nd) interatomic force constants (IFCs) and phonon frequencies, we have employed the DFPT method implemented in QE package for the optimized geometry of 3-atom primitive cell using a $\Gamma$-centered $12\times 12\times 12$ uniform grid in the reciprocal space (The convergency test is presented in A). To ensure the translational invariance of the symmetrized dynamical matrix, the acoustic sum-rule was applied. We have obtained the phonon frequencies by diagonalizing the calculated dynamical matrix. The effect of different atomic masses of the uranium isotopes on the phonon frequencies was checked and no meaningful differences obtained, and therefore, in subsequent calculations all phonon frequencies were calculated with the mass of U-238 isotope.

II.2 Thermal conductivity

The lattice thermal conductivity of bulk UO2 is calculated by solving the linearized-BTE[27, 28] for the steady-state phonon distribution function $f_{\lambda}$:

[TABLE]

where the left hand side of the equation corresponds to the phonon diffusion due to temperature gradient, and the term in the right hand side is the time rate of change of phonon distribution due to all allowed scattering processes. Here, ${\bf v}_{\lambda}$ is the group velocity of phonon in mode $\lambda$ and $\lambda\equiv(s,{\bf q})$ with $s$ and ${\bf q}$ being the phonon branch index and wave vector in reciprocal space, respectively. In the linearized-BTE, the distribution $f_{\lambda}$ differs from the equilibrium Bose-Einstein distribution function $f_{\lambda}^{0}=1/(\exp(\beta\hbar\omega_{\lambda})-1)$ by a linear term in $\nabla T$. Here $\beta=1/k_{B}T$ and $\omega_{\lambda}$ is the phonon frequency in mode $\lambda$. The solution of the BTE is done using a full iterative algorithm[29, 30, 31, 32] employing ShengBTE code package[33]. To calculate the thermal conductivity, the third-order (3rd) IFCs were computed up to the second shell of neighbors in a $3\times 3\times 3$ supercell with a $\Gamma$-centered $2\times 2\times 2$ uniform mesh of $k$ points. The 3rd-IFCs corresponding to the displacements of $(i,j,k)$ atoms along directions $(\alpha,\beta,\gamma)$ were computed using a three-point finite-difference method. To ensure the translational invariance, we have imposed the constraint $\sum_{i}\Phi_{\alpha\beta\gamma}(i,j,k)=0$ according to the prescription given by Esfarjani et al.[34] and Li et al.[35].

For the scattering contributions, which is treated within perturbation theory, we have considered scattering by isotopes[36, 37] and all the three-phonon processes satisfying the energy and quasi-momentum conservation:

[TABLE]

where ${\bf Q}$ is a reciprocal lattice vector. The linearized-BTE may be recasted in a form expressed in terms of a set of coupled equations for the phonon lifetimes, $\tau^{(\alpha)}_{\lambda}$ as[29, 30, 31, 32]:

[TABLE]

where $\tau^{0}_{\lambda}$, the phonon lifetime in the single-mode relaxation-time approximation (RTA), is defined by:

[TABLE]

in which $N$ is the number of unit cells, and $\tau^{(\alpha)}_{\lambda}$ corresponds to phonon modes propagating in $\alpha$ direction. The “+” and “–” symbols denote the sums are over two different types (”combination” and “decay” processes, respectively) of three-phonon processes defined by Eq.(2) and Eq.(3) including both “normal” (${\bf Q}=0$) and “umklapp” (${\bf Q}\neq 0$) processes. The quantities $\Gamma^{\pm}_{{\lambda}{\lambda^{\prime}}{\lambda^{\prime\prime}}}$ and $\Gamma_{{\lambda}{\lambda^{\prime}}}$ are respectively the three-phonon and phonon-impurity scattering rates. The quantity $\Delta^{\alpha}_{\lambda}$ in Eq.(4) is defined by:

[TABLE]

where $\xi^{\alpha}_{{\lambda}{\lambda^{\prime}}}=v^{\alpha}_{\lambda^{\prime}}\omega_{\lambda^{\prime}}/v^{\alpha}_{\lambda}\omega_{\lambda}$. The solution of Eq.(4) is achieved by the iterative process of $\tau_{\lambda}^{(\alpha),(n)}=\tau_{\lambda}^{0}(1+\Delta_{\lambda}^{\alpha,(n-1)})$ in which $\Delta_{\lambda}^{\alpha,(0)}=0$ and $\tau_{\lambda}^{(\alpha),(0)}=\tau_{\lambda}^{0}$. When the self-consistent $\tau_{\lambda}^{(\alpha)}$ is determined, the thermal conductivity tensor is calculated from:

[TABLE]

Having 2nd- and 3rd-IFCs at hand, the iteration of Eq.(4) is started by calculating the scattering rates $\Gamma^{\pm}_{{\lambda}{\lambda^{\prime}}{\lambda^{\prime\prime}}}$ and $\Gamma_{{\lambda}{\lambda^{\prime}}}$, and $\tau^{0}_{\lambda}$ from Eq.(5), and continued by computing $\Delta^{\alpha}_{{\lambda}}$ from Eq.(II.2). At each step of iteration, the conductivity is calculated from Eq.(7) and the iteration continues until the relative change of the conductivity norm $\sqrt{\sum_{\alpha}\sum_{\beta}|\kappa^{\alpha\beta}|^{2}}$ is less than $10^{-5}$. To calculate the thermal conductivity form Eq.(7), we have used an appropriate BZ sampling density of $10\times 10\times 10$ and Gaussian smearing of 0.1, which gives almost equal conductivity values obtained using $30\times 30\times 30$ and Gaussian smearing of 1.0.

III Results and Discussions

III.1 Geometry and Harmonic lattice dynamics

Using a non-spin-polarized calculation, we have first fully optimized the geometry of the UO2 primitive unit cell with space group 225 and obtained the optimized value of lattice constant for the FCC structure as $5.396\AA$.

Using the optimized geometry lattice parameter, the 2nd-IFCs were calculated employing the DFPT method. To this end, we first performed an scf calculation using the optimized lattice parameter with kinetic energy cutoffs of 60 and 600 Ry for the wavefunction and density expansion in terms of plane waves with a $\Gamma$-centered $8\times 8\times 8$ mesh in $k$ space with a tighter convergency threshold of $10^{-12}$ Ry. Using the resulting wavefunctions and eigenvalues as unperturbed quantities, the inhomogeneous KS equations in DFPT were solved in $q$-space self-consistently for the potential and density variations and thereafter the dynamical matrix is calculated using the atomic mass of U-238 isotope. The threshold for self-consistency of the potential variation was taken as $10^{-14}$. Diagonalizing the dynamical matrix, the phonon frequencies were obtained. In Fig. 1 the experimental and the calculated phonon dispersion curves along the high-symmetry directions in the BZ as well as the PhDOS are shown.

In Fig. 1(a), among the nine calculated branches, the lowest three belong to transverse acoustic (TA) and longitudinal acoustic (LA) that correspond to the vibrations of heavy U atoms. The higher frequency branches belong to the transverse optical (TO) and longitudinal optical (LO) modes that correspond to vibrations of lighter O atoms. As is seen, the agreement between experiment and our calculation is excellent for the acoustic branches. However, because we have not included the non-analytical terms in the dynamical matrices in our simplified model, near the $\Gamma$ point in the $M-\Gamma-X$ and $\Gamma-L$ paths we see no splittings of LO-TO branches. As we will see below, these fine structures does not spoil the calculated lattice thermal conductivity which is a sum over all branches.

In Fig. 1(b), the plotted PhDOS shows that the states in acoustic branches belong to the vibrations of the heavy U atoms whereas the states of optical branches are comprised from the vibrations of the lighter O atoms. It should be mentioned that the calculated phonon frequencies using the atomic mass of U-235 isotope did not result in a meaningful differences in the phonon dispersion, and therefore, we use the same phonon frequencies obtained from U-238 in the calculations of thermal conductivities of bulk UO2 with different isotope abundances.

Before starting the discussion on thermal conductivity, as was mentioned earlier, in reality the bulk UO2 solid is electronically an insulator[8] with an experimental gap of 2.1 eV. To obtain the exact phonon dispersion curve, some people[20] have used the sophisticated DFT+U method. However, in our simple model, the ”exact dispersion” showing LO-TO splitting is obtained just by applying the non-analytic correction to the dispersion of Fig 1(a). In the correction, the dielectric constant was chosen as 5.2, and the Born charges of U and O as 4.7 and -2.35 units. The result is shown in Fig 2. As is seen from Fig 2(a), the non-analytic term correction applied on the second-order force constants in our simple model, could reproduce the LO-TO splitting around the $\Gamma$ point which leads to excellent agreement with experimental results.

III.2 Thermal conductivity

From Eq.(7),to calculate the thermal conductivity we need phonon frequencies, equilibrium distribution function, group velocities, and phonon lifetimes. All needed quantities but phonon lifetimes were determined from our harmonic calculations. To calculate the phonon lifetimes we have solved the equation Eq.(4) by iteration until self-consistency. The calculated phonon lifetimes at $300^{\circ}K$ are shown in Fig. 3.

As seen from Fig. 3, the TA branches have the largest lifetimes and therefore have dominant role in thermal conductivity. On the other hand, we expect that the optical branches having the shortest lifetimes (corresponding to the vibrations of O atoms) to have smaller contributions in the total conductivity. This result contradicts the result reported in Ref.[20] which claims the optical phonons have the largest contributions in thermal conductivity. However, the experimental results reported in Ref. [20] that show the dominant role is played by optical phonons should be verified in another theoretical study using DFT+U or other methods.

In Fig. 4, we have shown the lifetimes of the phonons at two temperatures of $300$ and $1000^{\circ}K$. It is evidently seen that the phonon lifetimes decrease with temperature.

In Fig. 5, the calculated thermal conductivities in both RTA and full-iterative schemes are compared with experimental data. As is seen, the agreement between full-iterative results and experiment is very good at relatively lower temperatures. This is because, in experiment at lower temperatures the four- and higher- phonon processes are not activated and the dominant contribution comes from the three-phonon processes which is consistent with our calculations.

In Fig. 6, the thermal conductivity is resolved into the elemental contributions and as is seen, the contributions from U atoms are dominant at all temperatures. This is consistent with the PhDOS plot of Fig. 1(b) in which all contributions of acoustic branches comes from the U atoms.

Finally, in Fig. 7, we have plotted the thermal conductivities for fuels with different abundances of U-235 isotope. Inspecting the values for different abundances, we observe that with increasing the concentration of U-235, the thermal conductivity decreases, which is due to the increasing the scattering rates due to impurities, $\Gamma_{{\lambda}{\lambda^{\prime}}}$. This argument is confirmed by looking at Fig. 8. However, the decrease of thermal conductivity is not significant because, the phonon lifetime, which is determined from the combination of anharmonic and isotopic scattering rates through Matthiessen’s rule[38], does not change significantly.

IV Conclusions

In this work the lattice thermal conductivity of bulk UO2 is studied by the exact solution of the BTE equation for the steady-state phonon distribution function. In this lattice-thermal-conductivity study, we have modeled the high-$T$ paramagnetic state of bulk UO2 by a non-spin-polarized calculation, and neglected the Hubbard-U correction for the localized f electrons in the U atoms. The computed phonon spectra showed that among the nine branches, the lowest three belonging TA and LA, corresponding to U atoms, were in excellent agreement with experimental data. However, because in the simplified model we had not included the non-analytical terms in the dynamical matrices, no LO-TO splittings were observed near the $\Gamma$ point in the $M-\Gamma-X$ and $\Gamma-L$ paths in phonon spectra. Even with existing these differences in the experimental phonon spectra and our calculated results, the computed thermal conductivity was in good agreement with experiment because, the sum over phonon branches in Eq.(7) somehow washes out the detailed information of the phonon spectra. Our calculated thermal conductivity showed a small deviation from the experiments which was due to the fact that we had only taken into account the three-phonon processes while in reality the higher-order heat-carrying processes come into play at higher temperatures. The calculated phonon lifetimes at 300${}^{\circ}K$ showed that the TA branches had the largest lifetimes and therefore had dominant role in thermal conductivity. On the other hand, the optical branches having the shortest lifetimes (corresponding to the vibrations of O atoms) had smaller contributions in the total conductivity. Due to the fact that nuclear power reactors may use fuels with different relative concentrations of U-235 and U-238, we had also repeated the thermal-conductivity calculations for the cases with U-235 concentrations of 3%, 5%, 7%, and 20%, and observed a small decrease of thermal conductivity by increasing the impurity level. This fact was explained to be due to the increase of the phonon scattering rates from the impurity atoms (U-235).

Acknowledgement

This work is part of research program in School of Physics and Accelerators, NSTRI, AEOI.

Appendix A Convergency tests

In the course of thermal conductivity calculations, at each step we have to ensure the convergency of the calculated quantities with respect to relevant parameters. These steps include the DFT calculations, harmonic lattice-dynamic calculations, 3rd-IFCs calculations, and the solution of the BTE. First we ensure the convergency of lattice constant with respect to basis-set kinetic-energy cutoffs and $k$-grid for structural optimization. The results are tabulated in Table 1:

We have chosen a $\Gamma$-centered $8\times 8\times 8$, 60 Ry, and 600 Ry for the $k$-mesh, wavefunction cutoff, and density cutoff, respectively for our DFT calculations and obtained the converged value of $5.3968\AA$ for lattice constant.

In the DFPT calculations, we have started with a $q$-grid of $8\times 8\times 8$ and calculated the phonon frequencies and the 2nd-IFCs. To calculate the 3rd-IFCs, we have used the supercell approach. In this step we determine the appropriate supercell size, magnitude of atomic displacement, and the number of neighboring shells. The appropriate supercell size is intimately connected to the range of forces[39]. We have displaced one atom in a $3\times 3\times 3$ supercell with 81 atoms from its equilibrium position by 0.2 Bohr. Then, we determined the magnitude of forces acting on each of 81 atoms in the supercell and plotted with respect to the distance from the displaced atom (See Fig. 8).

As shown in Fig. 9, the forces drop to $\sim 0.01$ of the maximum force at a distance of $4\AA$ which encloses up to second neighbor atomic shells. This estimate indicates that the adopted supercell size of $3\times 3\times 3$ is appropriate for our 3rd-IFCs calculations. Using this force cutoff as an initial guess for the cubic cutoff, we have then determined the appropriate magnitude for the atomic displacements. To this end, considering all the symmetries of the system, we have generated sets of 108 displacements sufficient for calculating 3rd-IFCs for each magnitude of 0.005, 0.01, 0.03, 0.05, and 0.1$\AA$. Testing the thermal conductivities obtained from these trial 3rd-IFCs (See Fig. 10), we found that the appropriate magnitude of $0.03\AA$ gives stable value for the thermal conductivity at $300^{\circ}K$.

Next, we have tested the convergency of thermal conductivity with respect to the number of interaction shells in 3rd-IFCs calculation. Due to significant increasing of the number of displaced structures (108 compared to 208), but no significant gain in accuracy of thermal conductivity, we adopted up to second neighboring shell in our calculations (See Table 2).

Finally, we have tested the convergency of thermal conductivity with respect to the $q$-mesh in the DFPT calculations, and found that the $q$-mesh of $12\times 12\times 12$ leads to converged result.

Data availability

The raw or processed data required to reproduce these results can be shared with anybody interested upon sending an email to M. Payami.