Phonon Scattering and Thermal Conductivity of Actinide Oxides with Defects

: In the present study, we examine the e ﬀ ect of point defects and ﬁssion gases on thermal transport in representative actinide oxides used in modern reactors. In particular, oxygen interstitials and Kr / Xe ﬁssion gas bubbles are of primary focus. Reverse non-equilibrium molecular dynamics is employed to investigate thermal transport in UO 2 and PuO 2 with oxygen interstitials at the defect concentrations of 0.1%, 1%, and 5%. The results show that any alteration to the lattice structures of these fuels reduce their thermal conductivities signiﬁcantly. For the largest UO 2 structure simulated in the present study, for example, 0.1% oxygen interstitials decreased the thermal conductivity by 18.6%. For the case of the e ﬀ ect of ﬁssion gas bubbles, serious modiﬁcation to phonon dispersion in oxide fuels is caused by the presence of a single ﬁssion gas bubble, resulting in a large temperature drop in their temperature proﬁles. The average interfacial thermal resistance across a ﬁssion gas bubble (comprised of 30 Kr and / or Xe atoms) is estimated to be 2.1 × 10 − 9 Km 2 / W.


Introduction
Nuclear power reactors rely on efficient thermal transport between the fuel within the core fuel rods and the liquid or gas of the primary loop. Understanding the thermal transport properties of any given nuclear fuel is paramount when considering thermal power output and reactor safety. As the fuel undergoes fission, for example, a large temperature gradient occurs, and the center of the fuel pellet reaches a temperature well over the outer layer of the fuel pellet. In this environment, the fuel pellet undergoes thermal stressing which can result in the generation of cracks. Moreover, the hot spots in fuel pellets accelerate fission gas release and can cause fuel pellet swelling. Therefore, it is critical to understand the microscopic behavior of the thermal transport in nuclear fuels and the effects of various impurities such as point defects and fission gas bubbles on the thermal transport properties of nuclear fuels.
Noble gases such as Xe and Kr dissolve poorly in the oxide fuel matrix and form fission gas bubbles. Fission gas bubbles have negative impacts on the performance of nuclear fuels in various ways. They increase the internal pressure in fuel rods and promote the pellet-cladding mechanical interaction. They also cause embrittlement in grain boundaries, degrading the fuel's mechanical quality further. In particular, fission gas bubbles alter the thermal transport properties of oxide fuels, significantly impeding the efficient distribution of thermal energies [1].
Point defects in oxide fuels are also considered to have a detrimental effect on the thermal transport in oxide fuels, since they alter the vibrational modes of the crystalline structures [2][3][4]. Point defects are understood as the interruption in the fluorite crystal structures of oxide fuels; point defects exist as single to multi-point vacancies, interstitial species, and substitutional atoms. Amongst them, oxygen 2 of 13 interstitials are major byproducts of the nuclear fuel cycle; since oxygen interstitials have negative formation energy, oxide fuels can be oxidized easily and become hyper-stoichiometric. Therefore, it is critical to analyze their effect on the thermal/mechanical behavior of oxide fuels to understand radiation damage and non-stoichiometry. The effect of oxygen interstitials on actinide oxides and their diffusion process have been important research topics for a long time. In 1978, Breitung [5] showed summarized experimental data for diffusion coefficients of oxygen interstitial in UO 2 ± x. In 1981, Kim and Olander [6] reported that oxygen interstitials and vacancies contribute significantly to oxygen diffusion in stoichiometric UO 2 . More specifically, the thermal transport properties of actinide oxides with various defects have been investigated by many researchers. For example, Lucuta et al. [7] explored point defect scattering from oxygen interstitial in hyperschiometric SIMFUEL (simulated spent nuclear fuel). Pakarinen et al. [8] reported thermal conductivity reduction in UO 2 caused by He2+ ion irradiation. However, systematic investigations on the effect of oxygen interstitials on the thermal properties of actinide oxides are rare. Recently, Watanabe et al. [9] studied the thermal conductivity of hypo-and hyper-stoichiometric UO 2 with different defect concentrations at different temperatures. Along with the formation of oxygen interstitials, oxygen vacancies can come along. Günay [10] explored how fission induces these oxygen Frenkel pairs, which are generated when an atom is dislodged, creating a vacancy defect, and then relodged into the structure as an interstitial defect site. It was reported that Frenkel pairs are produced in multiple stages as the fission dose is increased; once isolated Frenkel pairs are produced at a constant rate, newly produced displaced atoms are trapped at existing defect sites, resulting in the formation of interstitial clusters.
In the present study, we provide a microscopic understanding of the effect of oxygen interstitials and fission gases on the thermal conductivities of oxide fuels. Molecular dynamics (MD) is used to estimate the thermal conductivities of oxide fuels with different defect types and concentrations. Interfacial thermal resistances caused by fission gas bubbles are analyzed. Moreover, phonon density of states are calculated to better understand phonon scattering caused by such defects.

Reverse Non-Equilibrium Molecular Dynamics (RNEMD)
Various classical MD algorithms, such as direct non-equilibrium molecular dynamics (NEMD), reverse non-equilibrium molecular dynamics (RNEMD), equilibrium molecular dynamics (EMD), and homogeneous non-equilibrium molecular dynamics (HNEMD) have been used to explore thermal transport in solids. In the present study, RNEMD [11] is selected. As a cause and effect reversed algorithm, RNEMD imposes a heat flux in a simulation structure to induce a temperature gradient across the structure. A simulation box is divided into multiple imaginary bins, and the bin located at the center of the simulation structure is chosen to be a hot bath while the bins at the two ends are chosen to be cold baths. A heat flux is applied by exchanging energy between the hot bath and the cold baths, inducing two identical temperature gradients across the structure, as shown in Figure 1. Thermal conductivity, k, is calculated from the induced temperature gradient, dT/dx, and the imposed heat flux, q, by using the following Fourier's law.
Here, <q> is the heat flux and <dT/dx> is the temperature gradient averaged over time and space. In the present study, the force field developed by Cooper et al. [12][13][14][15] is used with LAMMPS simulator [16]. In this force field, the interatomic energy in actinide oxides takes the following potential form.
Appl. Sci. 2020, 10, 1860 The potential energy of the ith atom is determined from the combination of the pairwise interaction between the ith atom and its surrounding atoms, φ αβ (r ij ) and many body perturbations, σ β (r ij ). G α is the constant of proportionality. This potential has been used successfully to predict thermal properties of actinide oxides [12]. For example, Rahman et al. [17] demonstrated the effectiveness of this model, and the authors estimated the thermal transport properties of actinide oxides with defects using the same force field [3,4]. For the interactions between noble gas-actinide and noble gas-oxygen, parameters developed by Cooper et al. [15] are used; during their parameter fitting procedure for the gas interactions with actinide oxides, they used the gas-gas interaction developed by Tang-Toennies [18].
Here, <q> is the heat flux and <dT/dx> is the temperature gradient averaged over time and space. In the present study, the force field developed by Cooper et al. [12][13][14][15] is used with LAMMPS simulator [16]. In this force field, the interatomic energy in actinide oxides takes the following potential form.
The potential energy of the ith atom is determined from the combination of the pairwise interaction between the ith atom and its surrounding atoms, ϕαβ (rij) and many body perturbations, σβ(rij). Gα is the constant of proportionality. This potential has been used successfully to predict thermal properties of actinide oxides [12]. For example, Rahman et al. [17] demonstrated the effectiveness of this model, and the authors estimated the thermal transport properties of actinide oxides with defects using the same force field [3,4]. For the interactions between noble gas-actinide and noble gas-oxygen, parameters developed by Cooper et al. [15] are used; during their parameter fitting procedure for the gas interactions with actinide oxides, they used the gas-gas interaction developed by Tang-Toennies [18]. The following equilibration steps are taken to equilibrate the simulation structure before running RNEMD. The purpose of equilibration is to make the atoms in the simulation structure forget their initial velocities and reach a thermodynamically stable state in order to obtain thermodynamic properties that are statistically meaningful. First, velocities are randomly imposed onto each atom in the simulation structure, and the simulation structure is energy-minimized with a time step of 1 fs. Then, the simulation structure is heated up to 500 K using isenthalpic (NPH) ensemble with Langevin thermostat. Next, the simulation structure is cooled down to the target temperature (300 K) with an isothermal-isobaric (NPT) ensemble. Energy exchange using the RNEMD algorithm is performed for 2 ns with a time-step of 1 fs (femtosecond) once the equilibration is complete. The authors observed that the amount of energy exchanged between the hot bath and cold bath is maintained to be constant in all simulation structures, and temperature profiles in the left and right segments of the simulation structures reach a steady state quickly, after only 0.2 ns.

Sample Preparation for Actinide Oxides with Point Defects
The unit cell of the stoichiometric sample and the unit cells with different point defects are illustrated in Figure 2. Three different concentrations (0.1%, 0.5%, and 5%) are selected for the study. The defect concentrations are defined as the ratio of the number of defect sites to the number of atoms in the stoichiometric sample. The sizes of simulation structures show a slight difference depending on the The following equilibration steps are taken to equilibrate the simulation structure before running RNEMD. The purpose of equilibration is to make the atoms in the simulation structure forget their initial velocities and reach a thermodynamically stable state in order to obtain thermodynamic properties that are statistically meaningful. First, velocities are randomly imposed onto each atom in the simulation structure, and the simulation structure is energy-minimized with a time step of 1 fs. Then, the simulation structure is heated up to 500 K using isenthalpic (NPH) ensemble with Langevin thermostat. Next, the simulation structure is cooled down to the target temperature (300 K) with an isothermal-isobaric (NPT) ensemble. Energy exchange using the RNEMD algorithm is performed for 2 ns with a time-step of 1 fs (femtosecond) once the equilibration is complete. The authors observed that the amount of energy exchanged between the hot bath and cold bath is maintained to be constant in all simulation structures, and temperature profiles in the left and right segments of the simulation structures reach a steady state quickly, after only 0.2 ns.

Sample Preparation for Actinide Oxides with Point Defects
The unit cell of the stoichiometric sample and the unit cells with different point defects are illustrated in Figure 2. Three different concentrations (0.1%, 0.5%, and 5%) are selected for the study. The defect concentrations are defined as the ratio of the number of defect sites to the number of atoms in the stoichiometric sample. The sizes of simulation structures show a slight difference depending on the defect types and concentrations once the structures are equilibrated. Therefore, we avoided using the physical dimensions of simulation structures. Instead, the sizes of simulation structures are represented by the number of unit cells to avoid any confusion that results from the slight differences in the dimensions.
The number of unit cells in y (N y ) and z (N z ) directions are fixed to be 6, while the number of unit cells in x (N x ) is increased progressively. Since 1D heat conduction is assumed by exchanging heat energy between the hot bath and the cold baths slowly, the 1D thermal transport is not likely to be altered by the surface boundaries (x − y and x − z planes). In their previous study [4], the authors already observed that thermal transport property is not affected much by the cross sectional areas (N y × N z ). This implies that 1D heat conduction is well-maintained in the present RNEMD simulations. The same phenomenon has been identified by other researchers [19][20][21]. Therefore, the authors decided not to examine the effect of cross-sectional area on the 1D heat conduction in the present study.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 14 defect types and concentrations once the structures are equilibrated. Therefore, we avoided using the physical dimensions of simulation structures. Instead, the sizes of simulation structures are represented by the number of unit cells to avoid any confusion that results from the slight differences in the dimensions. The number of unit cells in y (Ny) and z (Nz) directions are fixed to be 6, while the number of unit cells in x (Nx) is increased progressively. Since 1D heat conduction is assumed by exchanging heat energy between the hot bath and the cold baths slowly, the 1D thermal transport is not likely to be altered by the surface boundaries (x − y and x − z planes). In their previous study [4], the authors already observed that thermal transport property is not affected much by the cross sectional areas (Ny × Nz). This implies that 1D heat conduction is well-maintained in the present RNEMD simulations. The same phenomenon has been identified by other researchers [19][20][21]. Therefore, the authors decided not to examine the effect of cross-sectional area on the 1D heat conduction in the present study.
As mentioned earlier, two similar temperature profiles are obtained from the RNEMD simulation in the left and right segments of the simulation structure. Therefore, the sample length must be defined by the half of the total simulation structure length; Nx is the number of unit cells of the half size of the simulation structure in x-direction, which is the direction of the thermal conductivity calculation. Thermal conductivities are estimated for PuO2 or UO2 with various point defect types and concentrations for seven different sample lengths in the direction of thermal transport (Nx = 10, 20, 30, 40, 50, 75, and 100)

Sample Preparation for Actinide Oxides with Fission Gases
As shown in Figure 3, two fission gas bubbles are generated in a simulation box to keep the symmetry of RNEMD. In all structures simulated, the number of unit cells in y and z directions is fixed to be 5 while the number of unit cells in x is fixed to be 100. During RNEMD, however, the middle portion is selected as the hot bath and two ends are selected as the cold baths, reducing the characteristic sample length to be half of the total simulation box size. Therefore, the characteristic sample length considered in the present paper corresponds to 50 unit cells (Nx = 50) of the UO2 cubic structure. The total numbers of atoms (Kr and/or Xe) in fission gases are fixed to be 30 while varying their compositions. Xe are selected to see the effects of the fission gas composition on thermal transport in UO2. A custom MATLAB program is used to prepare simulation structures. Once a structure with 5 × 5 × 100 unit cells of UO2 is created, spaces that correspond to 3 × 3 × 3 unit cells in the pristine UO2 structure are emptied and fission gases are inserted in the empty spaces. A periodic boundary condition is used in the direction of the thermal conductivity calculation to mimic an infinitely large structure. However, it should be noted that phonons with wavelength larger than the size of the simulation box cannot be simulated in MD As mentioned earlier, two similar temperature profiles are obtained from the RNEMD simulation in the left and right segments of the simulation structure. Therefore, the sample length must be defined by the half of the total simulation structure length; N x is the number of unit cells of the half size of the simulation structure in x-direction, which is the direction of the thermal conductivity calculation. Thermal conductivities are estimated for PuO 2 or UO 2 with various point defect types and concentrations for seven different sample lengths in the direction of thermal transport (N x = 10, 20, 30, 40, 50, 75, and 100)

Sample Preparation for Actinide Oxides with Fission Gases
As shown in Figure 3, two fission gas bubbles are generated in a simulation box to keep the symmetry of RNEMD. In all structures simulated, the number of unit cells in y and z directions is fixed to be 5 while the number of unit cells in x is fixed to be 100. During RNEMD, however, the middle portion is selected as the hot bath and two ends are selected as the cold baths, reducing the characteristic sample length to be half of the total simulation box size. Therefore, the characteristic sample length considered in the present paper corresponds to 50 unit cells (N x = 50) of the UO 2 cubic structure. The total numbers of atoms (Kr and/or Xe) in fission gases are fixed to be 30 while varying their compositions. Five different gas units are simulated in this research (Figure 3b-f); 30 Kr, 20 Kr/10 Xe, 15 Kr/15 Xe, 10 Kr/20 Xe, and 30 Xe are selected to see the effects of the fission gas composition on thermal transport in UO 2 . A custom MATLAB program is used to prepare simulation structures. Once a structure with 5 × 5 × 100 unit cells of UO 2 is created, spaces that correspond to 3 × 3 × 3 unit cells in the pristine UO 2 structure are emptied and fission gases are inserted in the empty spaces. A periodic boundary condition is used in the direction of the thermal conductivity calculation to mimic an infinitely large structure. However, it should be noted that phonons with wavelength larger than the size of the simulation box cannot be simulated in MD simulations. Since only one sample length (N x = 50) is used in the current study, the effects of sample length and phonon boundary scattering are not examined. The authors will provide a more rigorous study by simulating fission gas bubbles in samples with various sample lengths to eliminate the effect of phonon boundary scattering in the near future.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 14 simulations. Since only one sample length (Nx = 50) is used in the current study, the effects of sample length and phonon boundary scattering are not examined. The authors will provide a more rigorous study by simulating fission gas bubbles in samples with various sample lengths to eliminate the effect of phonon boundary scattering in the near future.

The Effect of Point Defects on Thermal Transport in Actinide Oxides
Representative temperature profiles after RNEMD are shown in Figure 4. A temperature gradient is estimated directly from the temperature profile after 2 ns of RNEMD using linear regression function in MATLAB. As can be seen in the temperature profile, temperature profiles are not linear near the hot bath and cold baths. This is a result of the unphysical manner during the energy exchange between the hot bath and the cold baths; the frequency of energy swap between the hot bath and the cold baths seems to be faster than the diffusion of thermal energy across imaginary bins during RNEMD. To obtain a more statistically reliable temperature gradient, linear regression is performed six times graphically on the temperature profile plots. Error bars in Figure 5 are obtained from the combination of these multiple attempts of linear regression and the intrinsic statistical errors in the linear regression function.

The Effect of Point Defects on Thermal Transport in Actinide Oxides
Representative temperature profiles after RNEMD are shown in Figure 4. A temperature gradient is estimated directly from the temperature profile after 2 ns of RNEMD using linear regression function in MATLAB. As can be seen in the temperature profile, temperature profiles are not linear near the hot bath and cold baths. This is a result of the unphysical manner during the energy exchange between the hot bath and the cold baths; the frequency of energy swap between the hot bath and the cold baths seems to be faster than the diffusion of thermal energy across imaginary bins during RNEMD. To obtain a more statistically reliable temperature gradient, linear regression is performed six times graphically on the temperature profile plots. Error bars in Figure 5 are obtained from the combination of these multiple attempts of linear regression and the intrinsic statistical errors in the linear regression function. simulations. Since only one sample length (Nx = 50) is used in the current study, the effects of sample length and phonon boundary scattering are not examined. The authors will provide a more rigorous study by simulating fission gas bubbles in samples with various sample lengths to eliminate the effect of phonon boundary scattering in the near future.

The Effect of Point Defects on Thermal Transport in Actinide Oxides
Representative temperature profiles after RNEMD are shown in Figure 4. A temperature gradient is estimated directly from the temperature profile after 2 ns of RNEMD using linear regression function in MATLAB. As can be seen in the temperature profile, temperature profiles are not linear near the hot bath and cold baths. This is a result of the unphysical manner during the energy exchange between the hot bath and the cold baths; the frequency of energy swap between the hot bath and the cold baths seems to be faster than the diffusion of thermal energy across imaginary bins during RNEMD. To obtain a more statistically reliable temperature gradient, linear regression is performed six times graphically on the temperature profile plots. Error bars in Figure 5 are obtained from the combination of these multiple attempts of linear regression and the intrinsic statistical errors in the linear regression function.   Figure 5 shows the thermal conductivities of UO 2 and PuO 2 with various point defects at room temperature. The most significant finding is that oxygen interstitials have a detrimental effect on the thermal conductivity of UO 2 and PuO 2 comparable to oxygen vacancy defects; a 0.1% interstitial defect decreases thermal conductivity for UO 2 by 18.6% (14.99 W/m-K → 12.19 W/m-K), and PuO 2 Appl. Sci. 2020, 10, 1860 6 of 13 by 17.9% (17.15 W/m-K → 14.08 W/m-K) when compared to the respective stoichiometric structures for the simulation structures with 100 unit cells in the direction of thermal conductivity calculation (N x = 100). Any perturbations to the lattice structures of crystalline solids are expected to disturb phonon transport in the structures. In this regards, Park et al. [22] showed that the thermal transport is affected by a defect site, exhibiting interfacial thermal resistance across the defect.   Figure 5 shows the thermal conductivities of UO2 and PuO2 with various point defects at room temperature. The most significant finding is that oxygen interstitials have a detrimental effect on the thermal conductivity of UO2 and PuO2 comparable to oxygen vacancy defects; a 0.1% interstitial defect decreases thermal conductivity for UO2 by 18.6% (14.99 W/m-K → 12.19 W/m-K), and PuO2 by 17.9% (17.15 W/m-K → 14.08 W/m-K) when compared to the respective stoichiometric structures for the simulation structures with 100 unit cells in the direction of thermal conductivity calculation (Nx = 100). Any perturbations to the lattice structures of crystalline solids are expected to disturb phonon transport in the structures. In this regards, Park et al. [22] showed that the thermal transport is affected by a defect site, exhibiting interfacial thermal resistance across the defect.
In all structures simulated, thermal conductivities are increased with an increase in sample length. This indicates that the thermal transport in actinide oxides is still ballistic, even when they contain point defects. When the sample is larger than its intrinsic phonon mean free path, the thermal conductivity begins to reach the plateau and finally approaches its bulk thermal conductivity. However, this ballistic thermal transport is not observed when the sample contains a large number of defects; the thermal conductivity is not increased with an increase in the sample length when the concentration of point defects is 5%. The same phenomenon can be observed in the temperature profiles in Figure 4; temperature slips near hot baths, and cold baths that are observed in the simulation structure with low (0.1%) defect concentration are not observed in the simulation structure with high (5%) defect concentration. This indicates that the thermal transport in structures with high defect concentrations is readily diffusive even when the sample size is smaller than its phonon mean free path.
According to the popular analogy of phonon transport to the kinetic theory of gas, the effective phonon mean free path, leff, is expressed using Matthiessen rule [23,24] as where lph-ph is phonon-phonon scattering length, and lb is the phonon-boundary scattering length which is the length of the simulation box. As discussed in the simulation method section, the cross-sectional area (i.e., sample lengths in y and z directions) does not affect the 1D heat conduction in the simulation In all structures simulated, thermal conductivities are increased with an increase in sample length. This indicates that the thermal transport in actinide oxides is still ballistic, even when they contain point defects. When the sample is larger than its intrinsic phonon mean free path, the thermal conductivity begins to reach the plateau and finally approaches its bulk thermal conductivity. However, this ballistic thermal transport is not observed when the sample contains a large number of defects; the thermal conductivity is not increased with an increase in the sample length when the concentration of point defects is 5%. The same phenomenon can be observed in the temperature profiles in Figure 4; temperature slips near hot baths, and cold baths that are observed in the simulation structure with low (0.1%) defect concentration are not observed in the simulation structure with high (5%) defect concentration. This indicates that the thermal transport in structures with high defect concentrations is readily diffusive even when the sample size is smaller than its phonon mean free path.
According to the popular analogy of phonon transport to the kinetic theory of gas, the effective phonon mean free path, l eff , is expressed using Matthiessen rule [23,24] as where l ph-ph is phonon-phonon scattering length, and l b is the phonon-boundary scattering length which is the length of the simulation box. As discussed in the simulation method section, the cross-sectional area (i.e., sample lengths in y and z directions) does not affect the 1D heat conduction in the simulation structure. Therefore, it is concluded that the dominant factor that limits thermal transport in actinide oxide is the sample length in the direction of 1D heat conduction (x-direction in the present study). Since l b represents the phonon scattering that arises from the simulation boundaries, it can be considered to be the sample length in x-direction. If thermal conductivity is assumed to be proportional to the effective phonon mean free path based on the kinetic theory of gas, the thermal conductivity of the sample can be expressed as This relationship indicates that the inverse of thermal conductivity, k, is a linear curve when it is plotted as a function of the inverse of a simulation structure size. Therefore, thermal conductivity for a bulk size sample can be obtained by the extrapolation of the plot; i.e., in the limit when 1/l b → 0. However, the bulk thermal conductivity for the simulation structures with high defect concentrations cannot be estimated using this traditional Matthiessen rule [23], because the thermal conductivity (or effective phonon mean free path) is independent on the sample length (or l b ) ( Figure 5).
The thermal conductivity calculation in the present study is supported by the literature. Arima et al. [25] estimated the thermal conductivity of PuO 2 to be 15.2 W/m-K. Qin et al. [26] reported 14.6 W/m-K for the thermal conductivity of UO 2 using MD. Cooper et al. [21] estimated the thermal conductivities of UO 2 and PuO 2 to be ≈17.5 W/m-K and ≈22.5 W/m-K respectively by using MD. Fink [27] calculated the thermal conductivity of UO 2 to be ≈13.5 W/m-K. Wang et al. [28] used first-principles simulations to calculate the thermal conductivities of UO 2 and PuO 2 to be ≈12.5 W/m-K and 7.5 W/m-K.
As can be seen in Figure 5, the effect of oxygen interstitials on the thermal conductivities of UO 2 and PuO 2 is similar to that of oxygen vacancies. Several researchers have obtained correlations for the thermal conductivity actinide oxide fuels as a function of oxygen-to-metal ratio using experimental results. For example, Duriez et al. [29] proposed λ = 1/(A + BT) for the thermal conductivity (λ) of PuO 2 . Here, A = 2.85x + 0.035 m-K/W and B = (−7.15x + 2.86) × 10 −4 m/W, where x is the deviation from stoichiometry (x = 2 -O/M). If this correlation is used, the thermal conductivities of PuO 2 with 1%, 2%, and 5% oxygen vacancies are calculated to be 5.0, 3.58, and 1.93 W/m-K respectively. Our results exhibit thermal conductivities that are consistent with these experimental results; for the case of PuO 2 with N x = 100, thermal conductivities are calculated to be 5.82, 3.93, and 2.16 W/m-K for 1%, 2%, and 5% oxygen vacancy defects, respectively. However, it should be noted that the direct comparison of the results obtained in the present study to the previous experimental results is difficult. In most of the experiments, researchers simply control the oxygen-to-metal ratio of oxide fuels, and thus, they cannot control the amount of oxygen vacancies and oxygen interstitials independently. For example, if the oxide fuel is missing with 10% of oxygen, it could mean that the fuel has 20% oxygen interstitials and 30% oxygen vacancies (+20% − 30% = −10%). Or, it could mean that the fuel has 10% oxygen interstitials and 20% oxygen vacancies (+10% − 20% = −10%). In both cases, the total deficit in oxygen is the same (−10%), while the actual amount of oxygen interstitials is different. Therefore, more experimental efforts are required to examine the results obtained in the present research. For example, instead of the traditional thermogravimetry method (widely used in experimentalists), more advanced methods that utilize nanoscale characterization along with nanoscale imaging should provide a breakthrough.
It should be noted that the charge imbalance is ignored throughout this research because the defect concentrations are relatively small. UO 2 can be charged with a −2 for oxygen interstitial, −4 for uranium vacancy, and +2 for oxygen vacancy. Researchers [30,31] reported that −2 charged oxygen interstitial is expected to be stable over a wide range of compositions, while +2 charged oxygen vacancy is more stable in the Fermi level close to the valence band. Therefore, further study is required to clearly understand the effect of charge imbalance on the thermal transport property of oxide fuels. For example, the effect of charge imbalance can be more systemically studied by introducing the corresponding charge while compensating for defects in the other regions of the simulation structure. Figure 6a illustrates the change in the thermal conductivity of PuO 2 as a function of defect concentration. The effect of uranium/plutonium vacancy defects is notably greater than that of oxygen vacancies and interstitials. The addition or subtraction of larger atoms is expected to cause a greater change in phonon transport in crystalline structures, because it alters the lattice vibration in the simulation structure more abruptly. The authors have identified the same results in their previous research study [3] that investigated thorium vacancies and oxygen vacancies in ThO 2 . Other researchers [22,32,33] also reported the same observations in nanostructured crystalline solids. For the case of PuO 2 with N x = 100, 0.5% plutonium vacancy shows a 61% decrease in the thermal conductivity when compared to its stoichiometric structure, while 0.5% oxygen interstitials cause a 45% decrease in the thermal conductivity when compared to its stoichiometric sample.
concentration. The effect of uranium/plutonium vacancy defects is notably greater than that of oxygen vacancies and interstitials. The addition or subtraction of larger atoms is expected to cause a greater change in phonon transport in crystalline structures, because it alters the lattice vibration in the simulation structure more abruptly. The authors have identified the same results in their previous research study [3] that investigated thorium vacancies and oxygen vacancies in ThO2. Other researchers [22,32,33] also reported the same observations in nanostructured crystalline solids. For the case of PuO2 with Nx = 100, 0.5% plutonium vacancy shows a 61% decrease in the thermal conductivity when compared to its stoichiometric structure, while 0.5% oxygen interstitials cause a 45% decrease in the thermal conductivity when compared to its stoichiometric sample. Phonon density of states (PDOS) is plotted for PuO2 with different point defects in Figure 6b. PDOS is obtained directly from MD simulations through velocity auto-correlation using the same potential field at 300K [34]. The shaded area represents the PDOS of stoichiometric PuO2. The locations of the transverse acoustic (TA) peak (near 12 meV) and longitudinal acoustic (LA) peak are consistent with the results by other researchers [35][36][37], while the location of high-frequency optical phonon peak (near 75 meV) is slightly different from the previous results. By calculating phonon spectra using the DMFT, Yin et al. [38] showed that the LA phonon mode dominates the thermal transport because of its high group velocity and small anharmonicity implied by its Gruneisen constant, while the contribution of high frequency optical phonons to the thermal conductivity of PuO2 is negligible. Therefore, one may conclude that the significant decrease in thermal conductivity by plutonium vacancy is a result from the shift of the LA peak (from 20 meV to ≈18 meV) and the reduced mode degeneracies (flattened peak); the absence of large atoms such as plutonium is expected to alter the thermal transport in actinide oxide crystals more abruptly, since they decrease the phonon group velocity because of the increased structural compliance. In the case of uranium substitution, however, the alteration to phonon spectra is not evident because the size of the plutonium atom is similar to that of the uranium atom; the thermal conductivity of PuO2 is not degraded much by uranium substitution (Figure 6a). For both cases of oxygen vacancies and oxygen interstitials, the locations of LA peaks are not shifted. Although the alteration to the TA phonon mode in Phonon density of states (PDOS) is plotted for PuO 2 with different point defects in Figure 6b. PDOS is obtained directly from MD simulations through velocity auto-correlation using the same potential field at 300K [34]. The shaded area represents the PDOS of stoichiometric PuO 2 . The locations of the transverse acoustic (TA) peak (near 12 meV) and longitudinal acoustic (LA) peak are consistent with the results by other researchers [35][36][37], while the location of high-frequency optical phonon peak (near 75 meV) is slightly different from the previous results. By calculating phonon spectra using the DMFT, Yin et al. [38] showed that the LA phonon mode dominates the thermal transport because of its high group velocity and small anharmonicity implied by its Gruneisen constant, while the contribution of high frequency optical phonons to the thermal conductivity of PuO 2 is negligible. Therefore, one may conclude that the significant decrease in thermal conductivity by plutonium vacancy is a result from the shift of the LA peak (from 20 meV to ≈18 meV) and the reduced mode degeneracies (flattened peak); the absence of large atoms such as plutonium is expected to alter the thermal transport in actinide oxide crystals more abruptly, since they decrease the phonon group velocity because of the increased structural compliance. In the case of uranium substitution, however, the alteration to phonon spectra is not evident because the size of the plutonium atom is similar to that of the uranium atom; the thermal conductivity of PuO 2 is not degraded much by uranium substitution (Figure 6a). For both cases of oxygen vacancies and oxygen interstitials, the locations of LA peaks are not shifted. Although the alteration to the TA phonon mode in PuO 2 with oxygen interstitials is more severe than PuO 2 with oxygen vacancies, it is not likely to cause a meaningful difference in their thermal transport properties, since the contribution of TA phonon mode is trivial when compared to that of LA phonon mode [38].

Temperature Effect on Thermal Transport in Actinide Oxides
To investigate the temperature effect on the thermal conductivity of actinide oxides with oxygen interstitials, RNEMD has been performed for UO 2 and PuO 2 with different defect concentrations at 100 K, 200 K, 300 K, 600 K, 900 K, and 1200 K. As can be seen in Figure 7, the thermal conductivities decrease with an increase in temperature quickly. This can be explained by Umklapp phonon scattering. In crystalline solids, phonons with large wave vectors are populated at high temperatures. These phonons with large wave vectors fall outside the first Brillouin zone when they scatter from other phonons with large wave vectors, resulting in the loss of their momentum. Umklapp phonon-phonon scattering in actinide oxides has been investigated by different researchers. Gofryk et al. [39] incorporated a Umklapp scattering term in measuring the thermal conductivity of UO 2 . Park et al. [4] showed evidence of Umklapp phonon scattering by estimating the thermal conductivity of ThO 2 at different temperatures. Matsumoto et al. [40] also showed that the Umklapp process causes a thermal conductivity reduction in actinide oxides at high temperatures. actinide oxides has been investigated by different researchers. Gofryk et al. [39] incorporated a Umklapp scattering term in measuring the thermal conductivity of UO2. Park et al. [4] showed evidence of Umklapp phonon scattering by estimating the thermal conductivity of ThO2 at different temperatures. Matsumoto et al. [40] also showed that the Umklapp process causes a thermal conductivity reduction in actinide oxides at high temperatures.
The simulation structures with large defect concentrations (5%) do not show the temperature dependency that is observed in the simulation structures with low defect concentrations (0.1%). This implies that in the simulation structures with high defect concentrations, phonons are scattered greatly from their defect sites before they are scattered by other phonons with large wave vectors. In these simulation structures with high defect concentrations, the generation of long-wavelength phonons is prohibited by their many defect sites, even at very low temperatures (100 K or 200 K). Therefore, the thermal transport properties of the structures with high defect concentrations cannot be improved by decreasing the temperature, unlike stoichiometric samples or samples with low defect concentrations.  Figure 8 shows the temperature profile for the stoichiometric UO2 and the UO2 with two fission gas bubbles. For the case of UO2 structures with point defects, temperature slips are observed at the hot bath and the two cold baths when the defect concentrations are low (Figure 4). In the UO2 structure with The simulation structures with large defect concentrations (5%) do not show the temperature dependency that is observed in the simulation structures with low defect concentrations (0.1%). This implies that in the simulation structures with high defect concentrations, phonons are scattered greatly from their defect sites before they are scattered by other phonons with large wave vectors. In these simulation structures with high defect concentrations, the generation of long-wavelength phonons is prohibited by their many defect sites, even at very low temperatures (100 K or 200 K). Therefore, the thermal transport properties of the structures with high defect concentrations cannot be improved by decreasing the temperature, unlike stoichiometric samples or samples with low defect concentrations. Figure 8 shows the temperature profile for the stoichiometric UO 2 and the UO 2 with two fission gas bubbles. For the case of UO 2 structures with point defects, temperature slips are observed at the hot bath and the two cold baths when the defect concentrations are low (Figure 4). In the UO 2 structure with fission gas bubbles, however, the temperature slip is not clearly seen at the hot bath or the cold baths, indicating the diffusive nature of the thermal transport ( Figure 8a). As stated in the previous section (Section 3.1) and the authors' previous research [3], this diffusive thermal transport is observed in actinide oxides with uniformly distributed point defects. Considering that the simulation structure in Figure 8 contains only two fission gas bubbles, it is interesting that these two fission gas bubbles (two phonon scattering sources) make the thermal transport across the entire structure highly diffusive. In another thermal transport study that simulated a nanostructure with two interfacial thermal resistance sites, Park et al. [41] observed that thermal transport still maintains a ballistic transport nature, showing temperature slips at the hot bath and the cold baths when only a few interfaces are present. As can be seen in Figure 8b, fission gas bubbles distort the lattice structure of actinide oxides and become the sources of enormous phonon scattering; they correspond to a large number of point defects that are uniformly distributed across the whole simulation structure.

The Effect of Fission Gases on Thermal Transport in Actinde Oxides
another thermal transport study that simulated a nanostructure with two interfacial thermal resistance sites, Park et al. [41] observed that thermal transport still maintains a ballistic transport nature, showing temperature slips at the hot bath and the cold baths when only a few interfaces are present. As can be seen in Figure 8b, fission gas bubbles distort the lattice structure of actinide oxides and become the sources of enormous phonon scattering; they correspond to a large number of point defects that are uniformly distributed across the whole simulation structure. Using RNEMD, thermal conductivities are estimated for pristine UO2 with Nx = 50 and UO2 with a fission gas bubble, and results are plotted in Figure 9a. A significant decrease (approximately 78% decrease) in thermal conductivity is observed in all structures with fission gas bubbles. In their previous study [2], the authors observed that 2% of oxygen vacancy in the UO2 structure decreases its thermal conductivity by 77%. Considering that the number of Kr and/or Xe atoms in the present study is only 30 atoms, and the number of point defects that correspond to 2% is 150 atoms, one may come to the conclusion that fission gas bubbles have a more disastrous effect on thermal conductivity when compared to point defects. To examine the effect of the composition in a fission gas bubble to the thermal transport in UO2, interfacial thermal resistances across fission gas bubbles are estimated for UO2 with fission gases of different Kr/Xe compositions (Figure 9a). The composition of a fission gas bubble does not seem to make a meaningful difference in interfacial thermal resistance across the fission gas bubble. Both Kr and Xe are large atoms; the atomic mass of Kr is 83.798 u and the atomic mass of Xe is 131.293 u. When these large atoms are present in the middle of a UO2 crystal structure, they are expected to alter the lattice vibration of UO2 significantly. It is interesting that the number of Kr/Xe atoms in a fission gas bubble does not alter the interfacial thermal resistance across the fission gas bubble much. In the present study, we created the same amount of empty space (that corresponds to 3 × 3 × 3 unit cells of UO2) to place a fission gas bubble. Therefore, when the number of Kr/Xe atoms in a fission gas bubble is decreased, the density of the fission gas bubble is decreased at the same time. This could increase the gap between the fission gas bubble and surrounding UO2 crystals, giving rise to further phonon scattering. Using RNEMD, thermal conductivities are estimated for pristine UO 2 with N x = 50 and UO 2 with a fission gas bubble, and results are plotted in Figure 9a. A significant decrease (approximately 78% decrease) in thermal conductivity is observed in all structures with fission gas bubbles. In their previous study [2], the authors observed that 2% of oxygen vacancy in the UO 2 structure decreases its thermal conductivity by 77%. Considering that the number of Kr and/or Xe atoms in the present study is only 30 atoms, and the number of point defects that correspond to 2% is 150 atoms, one may come to the conclusion that fission gas bubbles have a more disastrous effect on thermal conductivity when compared to point defects. To examine the effect of the composition in a fission gas bubble to the thermal transport in UO 2 , interfacial thermal resistances across fission gas bubbles are estimated for UO 2 with fission gases of different Kr/Xe compositions (Figure 9a). The composition of a fission gas bubble does not seem to make a meaningful difference in interfacial thermal resistance across the fission gas bubble. Both Kr and Xe are large atoms; the atomic mass of Kr is 83.798 u and the atomic mass of Xe is 131.293 u. When these large atoms are present in the middle of a UO 2 crystal structure, they are expected to alter the lattice vibration of UO 2 significantly. It is interesting that the number of Kr/Xe atoms in a fission gas bubble does not alter the interfacial thermal resistance across the fission gas bubble much. In the present study, we created the same amount of empty space (that corresponds to 3 × 3 × 3 unit cells of UO 2 ) to place a fission gas bubble. Therefore, when the number of Kr/Xe atoms in a fission gas bubble is decreased, the density of the fission gas bubble is decreased at the same time. This could increase the gap between the fission gas bubble and surrounding UO 2 crystals, giving rise to further phonon scattering. In their previous study [2], the authors observed a similar phenomenon for UO 2 with vacancy defects. Although the atomic mass of oxygen is much smaller than the atomic mass of uranium, the oxygen vacancies have a detrimental effect on the thermal conductivity of UO 2 , comparable to uranium vacancies.
In Figure 9b, the local PDOS near the location of a fission gas bubble and the local PDOS far from the fission gas bubble are compared to the PDOS of stoichiometric UO 2 . For the case of the local PDOS near the location of the fission gas bubble, a large amount of phonon scattering is observed. The degeneracy of LA and TA phonons is reduced significantly, and their peaks are shifted drastically; in fact, the peaks of LA and TA phonons become almost unidentifiable. In the case of the local PDOS far from the location of the fission gas bubble, however, the peak for LA phonon mode (≈ 20 meV) is not shifted. This indicates that LA phonon modes are quickly recovered after being scattered by the fission gas bubble. Since LA phonon mode is a dominant thermal energy carrier in actinide oxides, thermal transport property in actinide oxides is expected to be recovered promptly once phonons pass through the location of fission gas bubbles [38]. However, the peak for TA phonon mode (near 17 meV) shows a slight shift and reduced degeneracy, and optical phonons are altered significantly.
Yin et al. [38] calculated the contribution of TA phonon mode to the thermal conductivity of UO 2 to be approximately 4% and the contribution of optical phonon modes to the thermal conductivity to be about 6%. Therefore, thermal transport is expected to become more diffusive, even at the region far from the location of fission gas bubbles, when compared to the thermal transport in stoichiometric UO 2 . This is supported by the temperature profile in Figure 8a where no temperature slips are identified near the hot bath and the cold baths. For the case of UO 2 with a fission gas bubble of 15 Kr/15 Xe, the thermal conductivity is estimated to be still very low (3.8 W/m-K), even when only the linear portion of the temperature profile is considered for the thermal conductivity calculation, (completely ignoring the temperature jump near the location of the fission gas bubble); for the same size of stoichiometric UO 2 , the thermal conductivity is estimated to be higher than 12 W/m-K [2]. This again demonstrates a diffusive nature of thermal transport in UO 2 with a fission gas bubble.
fact, the peaks of LA and TA phonons become almost unidentifiable. In the case of the local PDOS far from the location of the fission gas bubble, however, the peak for LA phonon mode (≈ 20 meV) is not shifted. This indicates that LA phonon modes are quickly recovered after being scattered by the fission gas bubble. Since LA phonon mode is a dominant thermal energy carrier in actinide oxides, thermal transport property in actinide oxides is expected to be recovered promptly once phonons pass through the location of fission gas bubbles [38]. However, the peak for TA phonon mode (near 17 meV) shows a slight shift and reduced degeneracy, and optical phonons are altered significantly. Yin et al. [38] calculated the contribution of TA phonon mode to the thermal conductivity of UO2 to be approximately 4% and the contribution of optical phonon modes to the thermal conductivity to be about 6%. Therefore, thermal transport is expected to become more diffusive, even at the region far from the location of fission gas bubbles, when compared to the thermal transport in stoichiometric UO2. This is supported by the temperature profile in Figure 8a where no temperature slips are identified near the hot bath and the cold baths. For the case of UO2 with a fission gas bubble of 15 Kr/15 Xe, the thermal conductivity is estimated to be still very low (3.8 W/m-K), even when only the linear portion of the temperature profile is considered for the thermal conductivity calculation, (completely ignoring the temperature jump near the location of the fission gas bubble); for the same size of stoichiometric UO2, the thermal conductivity is estimated to be higher than 12 W/m-K [2]. This again demonstrates a diffusive nature of thermal transport in UO2 with a fission gas bubble.

Conclusions
In the present study, reverse non-equilibrium molecular dynamics (RNEMD) was used to explore thermal transport in two actinide oxides (UO2 and PuO2) with oxygen vacancies, uranium/plutonium vacancies, and oxygen interstitials. The effects of these defects were observed to be detrimental across all actinide oxides simulated in this research. Even the smallest percentages of these defects have a notable

Conclusions
In the present study, reverse non-equilibrium molecular dynamics (RNEMD) was used to explore thermal transport in two actinide oxides (UO 2 and PuO 2 ) with oxygen vacancies, uranium/plutonium vacancies, and oxygen interstitials. The effects of these defects were observed to be detrimental across all actinide oxides simulated in this research. Even the smallest percentages of these defects have a notable impact on thermal conductivity. The vacancies of primary atoms (uranium and plutonium) reduced the thermal conductivities of oxide fuels the most significantly; for the case of PuO 2 of Nx = 100 with 0.5% plutonium vacancy, a 61% decrease in the thermal conductivity was observed. Oxygen interstitials were found to have an almost identical effect on the thermal transport in UO 2 and PuO 2 to oxygen vacancies; 0.1% oxygen interstitials decreased the thermal conductivity of UO 2 of Nx = 100 by 18.6%, while 0.1% oxygen vacancies decreased the thermal conductivity of the same sample by 19.4%. In phonon density of states (PDOS) plots with 1% of oxygen vacancies and interstitials, a similar degree of modification to the phonon spectra was observed; the both defect types do not alter LA phonon mode seriously, which is a dominant thermal energy carrier in actinide oxides.
The thermal transport property of actinide oxides becomes highly diffusive by the presence of only one fission gas bubble. The presence of a single fission gas bubble (that is consisted of 30 Kr/Xe atoms) was found to generate the interfacial thermal resistance of 2.1 × 10 −9 K-m 2 /W. From the local PDOS plot of UO 2 with a fission gas bubble, it was observed that thermal transport is still diffusive even in the region far from the location of the fission gas bubble. It is interesting to note that changing the composition of a fission gas bubble (ratio of Kr and Xe atoms) does not alter the interfacial thermal resistance much, despite the large difference between their atomic masses.