Thorium dioxide is used industrially in high temperature applications, but more insight is needed into the behavior of the material as part of a mixed-oxide (MOX) nuclear fuel, incorporating uranium. We have developed a new interatomic potential model including polarizability via a shell model, and commensurate with a prominent existing UO_{2} potential, to conduct configurational analyses and to investigate the thermophysical properties of uranium-doped ThO_{2}. Using the GULP and Site Occupancy Disorder (SOD) computational codes, we have analyzed the distribution of low concentrations of uranium in the bulk material, where we have not observed the formation of uranium clusters or the dominance of a single preferred configuration. We have calculated thermophysical properties of pure thorium dioxide and Th_{(1−x)}U_{x}O_{2} which generated values in very good agreement with experimental data.

## I. INTRODUCTION

Although the majority of nuclear power reactors operating today use uranium and/or plutonium-based fuels such as low-enriched uranium (LEU) or mixed oxide fuels (MOX), alternative nuclear fuels are increasingly considered as well.^{1} The most promising alternatives are fuels based on thorium, which is generally estimated to be between three and four times more abundant than uranium in the Earth’s crust, with large stores in India, the United States, and other areas around the world. Thorium compounds can be used to fuel nuclear reactors and they thus present an attractive alternative to the standard uranium and plutonium-MOX fuels.^{2}

Whereas thorium is a fertile material, it is not fissile and therefore a sustained conversion to the fissile ^{233}U isotope can only occur in the presence of a neutron source. Although plutonium may serve as the source, there are several advantages to blending thoria with urania that have been confirmed experimentally, including changes in the decay heat, melting point, and thermal conductivity of the material, in addition to a decrease in the release rate of fission gases despite an actual increase in the production of fission gases.^{3}

Pure thorium dioxide (thoria) is a ceramic material with a melting point of 3651 K, the highest of any known oxide.^{4} Thoria has the fluorite structure, which is shared with the dioxides of, for example, cerium, plutonium, and uranium. Unlike uranium and plutonium oxides, ThO_{2} also represents the highest oxidation state of the material, which makes it exceptionally stable in the presence of oxygen or oxygenated water. This corrosion resistance has positive implications for the use of thoria in nuclear fuels.^{3} Since a neutron source is required for its application as a nuclear fuel, and as ThO_{2} and UO_{2} are isostructural, there is significant interest in thoria-urania solid solutions and both experimental and theoretical investigations have been carried out.^{5–8} As it is useful to study thoria in the context of uranium-doping, it is advantageous to have a reliable computational method to model these materials due to the safety hazards and high costs of conducting experiments on radioactive materials.

We are not aware of any robust theoretical investigations into the different possible configurations of uranium substitution in thoria, which might elucidate the distribution of the uranium atoms serving as a neutron source in a thorium dioxide nuclear fuel, and the thermodynamic properties of these configurations. It is important to understand how uranium behaves when inserted into thoria in order to understand the behavior of potential large-scale MOX fuel rods. Since thorium has no further oxidation states beyond Th^{4+}, UO_{2} is more susceptible to oxidation than ThO_{2}, and so any uranium clustering at surfaces or grain boundaries could lead to oxidation and thereby compromise the corrosion-resistance of the material. In this paper, we present our computational investigation into the numerous independent configurations obtainable at a variety of uranium concentrations in ThO_{2}, using a new interatomic potential to describe the interaction between thorium and oxygen in this system.

In experimental test reactors using thorium fuels, typically only small concentrations, in the region of 5-15%, of uranium or plutonium are used. The role of U or Pu in these fuel blends is not as a fuel itself but rather as a neutron source for the production of ^{233}U by ^{232}Th. To model such small concentrations of uranium dopant in thoria, we need a sufficiently large system, which requires the use of affordable computational methods based on interatomic potentials. Early potentials, such as those used by Benson et al to study surface energies of thoria, were developed prior to the publication of experimental data on the bulk properties of thorium dioxide and therefore relied on estimating those properties from other ionic oxides.^{9} Later studies either did not publish the potential parameters or did not closely reproduce the available experimental data.^{10–12} While several good ThO_{2} potentials were published by Behera, et al in 2012, we have developed a thoria potential to be compatible with a leading UO_{2} potential^{13–15} and we have applied this model to the distribution of uranium in relatively large ThO_{2} systems.

## II. METHODOLOGY

### A. Energies and geometries

The energies of different configurations of the mixed solid were evaluated using the General Lattice Utility Program (GULP), version 3.4.^{16,17} GULP implements the Born model of solids, which assumes that the ions in the crystal interact via short-range repulsive interactions, longer-range attractive interactions, and long-range Coulombic interactions.^{18} The short-range energy contribution due to interactions between each particle and all others is summed within a predetermined cut-off. The contribution of the long-range electrostatic forces is summed using the Ewald sum.^{19} The electronic polarizability of ions is included via the shell model of Dick and Overhauser in which each of the polarizable ions, in this case the oxygen atoms, are represented by a core and a massless shell, connected by a spring.^{20} The polarizability of the model ion is then determined by the spring constant and the charges of the core and shell. The advantage of using the shell model is that it better reproduces the elastic and dielectric properties of ionic solids than rigid ion models, though it increases the computational cost of the simulations.^{21}

The GULP code is widely used in solid state chemistry in simulations of materials using periodic boundary conditions. We have used GULP to derive the potential parameters and compute bulk properties for ThO_{2}, with initial values taken from a suitable existing UO_{2} potential, and fitting to experimental crystallographic data. Details of the interatomic potentials used are given in section III A, equation (15). In order to make geometry predictions, the lattice energies are minimized with respect to the structural parameters, until the forces acting on the ions are all less than 0.001 eV Å^{−1}. All structures reported are the result of constant pressure energy minimizations, where not only the ionic positions but also the cell parameters are allowed to vary to find the energy minimum. The external pressure was set to zero in all calculations.

### B. Representation of the solid solutions

The Th_{1−x}U_{x}O_{2} solid solutions were represented by a symmetry-adapted ensemble of configurations in a supercell, using the methodology implemented in the SOD (Site Occupancy Disorder) program,^{22,23} which has been employed previously in the simulation of a range of mineral solid solutions. Successful applications of SOD include a study of the thermophysical properties of carbonates, investigations into the mixing thermodynamics of oxides and a simulation of cation distribution and thermodynamics in sulfides.^{24–31} This program generates the complete configurational space for each composition in a supercell of the structure, before extracting the subspace of symmetrically inequivalent configurations, for which energies and other properties are evaluated. In this case a 2x2x2 thoria supercell (Figure 1) was used as the base simulation cell, and we have worked with uranium concentrations up to x = 0.16. The computational cost of investigating higher U concentrations in a system of this size was considerable and of limited practical concern, as the amount of uranium-doping required in thoria fuels is small. However, to investigate the full range of the solid solution and derive thermodynamic properties that can be compared with experimental values, a full configurational analysis was carried out on a 1x1x2 supercell.

Once the configurational spectrum is obtained, it is possible to derive thermodynamic properties through further application of statistical mechanics. In the simplest formulation, the extent of occurrence of one particular configuration, *n*, in the disordered solid in configurational equilibrium at a temperature *T* can be described by a Boltzmann-like probability, **P**_{n}:

where E_{n} is the energy of configuration *n* and Z is the partition function given by the following equation:

The energy of the system is now calculated as a sum as follows:

Then in the reduced configurational space, a new probability can be determined for an inequivalent configuration, *m*. This is calculated from the energy of that configuration, E_{m}, and its degeneracy Ω_{m} (the number of times that the configuration is repeated in the complete configurational space):

where m = 1, …, M (M is the number of inequivalent configurations), R is the gas constant, and Z is the partition function (Equation (5)), which guarantees that the sum of all the probabilities equals 1 and also gives access to the calculation of configurational free energies and entropies.

E_{m(red)} is the reduced energy and relates to the energy E_{m} and the reduced, or degeneracy, entropy, S_{m}, via:

From this relationship, it is possible to compare the inequivalent configurations as the degeneracy and the energies are considered. For two inequivalent configurations of the same energy, the one with the highest degeneracy will then have the higher degeneracy entropy and therefore a lower reduced energy and thus higher probability (equation (4)). It is possible in this way to obtain temperature-dependent configurational thermodynamics functions and equilibrium degrees of disorder.^{32–34} In this work we have focused on the distribution of uranium in a range of concentrations in the thoria lattice and the derivation of thermodynamic properties of the solid solution.

The configurational free energy *G* of the disordered solid can be obtained directly from the partition function:

and any average observable, including the enthalpy *H* of the solution, can be estimated using configurational averaging:

Having obtained the enthalpy and free energy of the solid solution, it is useful to evaluate the enthalpy of mixing:

which can be compared with experimental calorimetric determinations, and the free energy of mixing:

which can also be found experimentally from equilibrium composition measurements.

Technically important properties like the thermal linear expansion coefficient and heat capacity can be calculated:

where *a* is the lattice parameter of the cell at a given temperature, T, and *a _{293}* is the lattice parameter at 293 K. The values of

*a*at seven temperatures were calculated directly in GULP followed by fitting of the values to a 3

^{rd}order polynomial (Table I) to enable the extrapolation to higher values. As noted, due to the computational cost of calculating the configurational thermodynamics for the full ThO

_{2}-UO

_{2}solid solution using the large 2x2x2 supercell, a 1x1x2 cell was used instead for these calculations.

Material
. | c
. _{0} | c
. _{1}×10^{5} | c
. _{2}×10^{8} | c
. _{3}×10^{12} | α
. _{293} |
---|---|---|---|---|---|

ThO _{2} | 5.60793 | 2.79036 | 1.46253 | -3.43684 | 5.608 |

UO _{2} | 5.45337 | 2.54896 | 1.45446 | -3.76733 | 5.462 |

Th_{0.87}U_{0.13}O _{2} | 5.58811 | 2.93563 | 1.23247 | -2.01114 | 5.598 |

Th_{0.75}U_{0.25}O _{2} | 5.56882 | 2.71268 | 1.70802 | -4.78603 | 5.578 |

Material
. | c
. _{0} | c
. _{1}×10^{5} | c
. _{2}×10^{8} | c
. _{3}×10^{12} | α
. _{293} |
---|---|---|---|---|---|

ThO _{2} | 5.60793 | 2.79036 | 1.46253 | -3.43684 | 5.608 |

UO _{2} | 5.45337 | 2.54896 | 1.45446 | -3.76733 | 5.462 |

Th_{0.87}U_{0.13}O _{2} | 5.58811 | 2.93563 | 1.23247 | -2.01114 | 5.598 |

Th_{0.75}U_{0.25}O _{2} | 5.56882 | 2.71268 | 1.70802 | -4.78603 | 5.578 |

GULP will also calculate certain thermodynamic properties of a system when a phonon calculation is performed. We have taken advantage of this fact to extract the specific heat capacities at constant volume of our system while computing phonons, where the constant volume heat capacity, C_{v}, is calculated from the vibrational partition function, Z_{vib}:

## III. RESULTS

### A. Development of the thoria potential

We have used the thoria structure^{35,36} in our derivation of a Th-O Buckingham potential, V_{Buck}(r_{i,j}):

where **r**_{ij} is the distance between ions i and j, and A, ρ, C, are variable parameters representing the pair-wise repulsion.

Thoria has a cubic unit cell with the space group $Fm 3 \xaf m$ (N° 225), with four cation sites in the unit cell. The least-squares fitting against structure and other properties was done in GULP using experimental data from Idiri and Macedo and co-workers.^{35,36} While there has been extensive debate in the literature as to the correct value of the bulk modulus of thorium dioxide, in this work we have used the experimental value of 198 ± 2 GPa determined by Idiri et al using energy dispersive x-ray diffraction (EDXRD)^{35} Staun Olsen et al also studied the bulk modulus of thoria with high-pressure x-ray diffraction (XRD) and synchrotron radiation, finding a similar value of 195 ± 2 GPa.^{37} The oxygen-oxygen interaction in ThO_{2} is described by the potential published by Catlow, whereas initial values of A and ρ for the Th-O interaction were taken from the Catlow U-O Buckingham potential.^{15} Although U^{4+} (1.00 Å) has a smaller ionic radius than Th^{4+} (1.05 Å), due to the similarities between urania and thoria, these values provided a good starting point for the derivation of the Th-O potential, which is introduced in Table II.^{38}

. | Charge (e)
. | . | |
---|---|---|---|

Species
. | core . | shell . | Spring constant (eV Å^{−2})
. |

Th | +4.000 | - | - |

U | +4.000 | - | - |

O | +2.400 | −4.400 | 292.98 |

Buckingham Potential | A (eV) | ρ (Å) | (CeV Å^{6}) |

Th—O_{shell} | 1281.775 | 0.3910 | 0.0 |

U—O_{shell} | 1217.800 | 0.3871 | 0.0 |

O_{shell}—O_{shell} | 22764.300 | 0.1490 | 112.2 |

. | Charge (e)
. | . | |
---|---|---|---|

Species
. | core . | shell . | Spring constant (eV Å^{−2})
. |

Th | +4.000 | - | - |

U | +4.000 | - | - |

O | +2.400 | −4.400 | 292.98 |

Buckingham Potential | A (eV) | ρ (Å) | (CeV Å^{6}) |

Th—O_{shell} | 1281.775 | 0.3910 | 0.0 |

U—O_{shell} | 1217.800 | 0.3871 | 0.0 |

O_{shell}—O_{shell} | 22764.300 | 0.1490 | 112.2 |

Table III compares structural properties of ThO_{2} calculated using the newly developed potential with available experimental data. There is of course excellent agreement with the lattice parameter, which was used in the derivation, and good agreement with the elastic properties. By using a shell model potential, we are also able to correctly model the Cauchy violation observed in fluorite structures (C_{12}≠C_{44}). This potential will now be used to study the Th_{1−x}U_{x}O_{2} solid solutions.

ThO_{2}
. | Experimental Value^{35,36}
. | Calculated Value . | % Error in Calculated Value . |
---|---|---|---|

Lattice Parameter (Å) | 5.601 | 5.6001 | 0.00178 |

Bulk Modulus (GPa) | 198 | 206.4 | 4.24 |

Shear Modulus (GPa) | 98 | 94.4 | −4.08 |

Elastic Constant C_{11} (GPa) | 367 | 389.5 | 6.14 |

Elastic Constant C_{12} (GPa) | 106 | 114.9 | 8.36 |

Elastic Constant C_{44} (GPa) | 79.7 | 65.8 | −16.72 |

ThO_{2}
. | Experimental Value^{35,36}
. | Calculated Value . | % Error in Calculated Value . |
---|---|---|---|

Lattice Parameter (Å) | 5.601 | 5.6001 | 0.00178 |

Bulk Modulus (GPa) | 198 | 206.4 | 4.24 |

Shear Modulus (GPa) | 98 | 94.4 | −4.08 |

Elastic Constant C_{11} (GPa) | 367 | 389.5 | 6.14 |

Elastic Constant C_{12} (GPa) | 106 | 114.9 | 8.36 |

Elastic Constant C_{44} (GPa) | 79.7 | 65.8 | −16.72 |

### B. Thermal Expansion of Thorium Dioxide

It is important to have a thorough understanding of the temperature-dependent properties of any nuclear fuel in order to understand and predict its behavior. For example, the thermal expansion must be understood before the fuel rods and cladding are designed.^{39} If the fuel-cladding gap is either too small or too large, failure of the fuel cladding can occur. If the gap is too small, as the fuel expands, pressure is exerted on the cladding and deformation of the fuel can occur. If the gap is too large, the fuel may vibrate excessively and so damage the cladding.^{40}

We have calculated the value of the thermal linear expansion coefficient, α, at 6.4 × 10^{−6} K^{−1} at 300 K. This value is in line with the experimental values reported by Marples (7.3 × 10^{−6} K^{−1}), Kempter and Elliot (8.2 × 10^{−6} K^{−1}), and, at 293K, Yamashita et al (8.43 × 10^{−6} K^{−1})^{41–43} While our thermal expansion coefficient does underestimate the values found experimentally, this could be explained by the small difference in the initial lattice parameters used in the experimental systems and in our system. We have used a value for the lattice parameter, a, of ThO_{2} of 5.6 Å in our simulations, whereas Kempter, using XRD measurements, found an initial a value of 5.597 Å and Yamashita one of 5.58 Å, also using XRD. As ThO_{2} does not undergo much expansion with temperature, it is possible that by using a larger initial lattice parameter, our calculations failed to account for a small amount of expansion and so underestimated the coefficient of thermal expansion. It should also be noted that Kempter reported impurities in his original samples of ThO_{2} and these could have had a non-negligible effect on the experimental values reported.

A computational study based on Density Functional Theory (DFT) by Lu et al reports values at higher temperatures, i.e. α = 33–37 × 10^{−6} K^{−1}, an order of magnitude above the values determined experimentally and that are estimated by our simulation.^{44} Martin et al have used a rigid ion potential in molecular dynamics simulations to calculate an α of 9.2 × 10^{−6} K^{−1} at 1500 K, which is in the same range as our value and the experimental values but this method does not account for the polarizability of the ions or the contribution of phonons in their simulation^{7} Similarly, using a Born-Mayer-Huggins potential with a partially ionic model, Ma et al have calculated an α at 1500 K of 11.9 × 10^{−6} K^{−1}.^{45} These calculated values at 1500 K are lower than the experimental values reported at 1200 K of Yamashita et al and that at 1273 K of Kempter and Elliot. These findings indicate that interatomic potentials used to determine the coefficient of linear expansion tend to underestimate the values but are of the same order of magnitude as experiment and a distinct improvement on the values obtained from DFT for this material.

### C. Configurational Variations in U-Substituted Thorium Dioxide

The equilibrium geometries and energies of the Th_{(1−x)}U_{x}O configurations were calculated both with and without considering vibrational effects. The SOD code identified a number of independent configurations for the uranium-doped thoria systems under investigation (Table IV). The number of inequivalent configurations increases significantly with each additional substitution until a U mole fraction of x = 0.5 is reached in the 1x1x2 supercell.

Supercell . | Th_{(1−x)}U_{x}O_{2}
. | Mole fraction of Uranium, x . | Total Number of Configurations . | Total number of Inequivalent Configurations . |
---|---|---|---|---|

Th_{7}U_{1}O_{16} | 0.125 | 8 | 1 | |

Th_{6}U_{2}O_{16} | 0.25 | 28 | 4 | |

Th_{5}U_{3}O_{16} | 0.375 | 56 | 4 | |

1x1x2 | Th_{4}U_{4}O_{16} | 0.5 | 70 | 8 |

Th_{3}U_{5}O_{16} | 0.625 | 56 | 4 | |

Th_{2}U_{6}O_{16} | 0.75 | 28 | 4 | |

Th_{1}U_{7}O_{16} | 0.875 | 8 | 1 | |

Th_{31}U_{1}O_{64} | 0.031 | 32 | 1 | |

Th_{30}U_{2}O_{64} | 0.063 | 496 | 5 | |

Th_{29}U_{3}O_{64} | 0.094 | 4960 | 14 | |

2x2x2 | Th_{28}U_{4}O_{64} | 0.125 | 35960 | 71 |

Th_{27}U_{5}O_{64} | 0.156 | 201376 | 223 | |

Th_{26}U_{6}O_{64} | 0.188 | 906192 | 874 | |

Th_{25}U_{7}O_{64} | 0.219 | 3365856 | 2706 |

Supercell . | Th_{(1−x)}U_{x}O_{2}
. | Mole fraction of Uranium, x . | Total Number of Configurations . | Total number of Inequivalent Configurations . |
---|---|---|---|---|

Th_{7}U_{1}O_{16} | 0.125 | 8 | 1 | |

Th_{6}U_{2}O_{16} | 0.25 | 28 | 4 | |

Th_{5}U_{3}O_{16} | 0.375 | 56 | 4 | |

1x1x2 | Th_{4}U_{4}O_{16} | 0.5 | 70 | 8 |

Th_{3}U_{5}O_{16} | 0.625 | 56 | 4 | |

Th_{2}U_{6}O_{16} | 0.75 | 28 | 4 | |

Th_{1}U_{7}O_{16} | 0.875 | 8 | 1 | |

Th_{31}U_{1}O_{64} | 0.031 | 32 | 1 | |

Th_{30}U_{2}O_{64} | 0.063 | 496 | 5 | |

Th_{29}U_{3}O_{64} | 0.094 | 4960 | 14 | |

2x2x2 | Th_{28}U_{4}O_{64} | 0.125 | 35960 | 71 |

Th_{27}U_{5}O_{64} | 0.156 | 201376 | 223 | |

Th_{26}U_{6}O_{64} | 0.188 | 906192 | 874 | |

Th_{25}U_{7}O_{64} | 0.219 | 3365856 | 2706 |

The variation of the lattice parameter, a, with the mole fraction, x, in the ThO_{2}–UO_{2} solid solutions modeled follows Vegard’s law (Figure 3), in good agreement with the experimental measurements of Hubert et al and Banerjee et al^{5,46} We can further compare our models with these data by examining the cation-anion and cation-cation distances Th–O, U–O, Th–Th, and U–U (Figure 4 and Figure 5).

Before we consider the distances in the full ThO_{2}–UO_{2} solid solution, it is useful to examine the effect of a single uranium atom substitution in the thoria supercells. While the 1x1x2 supercells is of practical use in modeling the full solid solution, we observe that our values of the interatomic distances at x = 0.13 are higher than the experimental values at this mole fraction of uranium (Figure 4 and Figure 5). To investigate why this might be, we have investigated the effects of a single uranium substitution in both the 1x1x2 and 2x2x2 supercells (Figure 2). Of course a single substitution yields different mole fractions of uranium in the cells, where x = 0.13 in a 1x1x2 and 0.03 in a 2x2x2 supercell. However, we can see from the interatomic distances in Table V, that small changes are introduced to the lattice even at the lowest mole fraction of uranium studied, x = 0.03. While in this cell, the introduction of one U atom does not significantly affect the Th-Th or Th-O distances, the Th-U and U-O distances are smaller than the Th-Th and Th-O distances being replaced. Moving to the 1x1x2 cell, with one U substitution, all four of the possible cation-cation and anion-anion distances are affected by the presence of a uranium atom. Due to the shape of the 1x1x2 cell, a larger percentage of the oxygen atoms in the lattice are affected by the presence of uranium than in the 2x2x2 cell. These oxygen atoms are moved closer to uranium than they would be in the pure ThO_{2} lattice, thus affecting the other interatomic distances in the small supercell, resulting in an exaggerated effect when a single substitution is introduced into a small supercell. In a larger supercell with an equivalent mole fraction of uranium, there are many available configurations. We have calculated the interatomic distances in a 2x2x2 supercell for x = 0.13 as an average of all 71 inequivalent configurations (Table V). These distances reproduce more closely the experimental values of Hubert et al compared to the distances in the singly substituted 1x1x2 cell. Therefore, we consider that the overestimation of the interatomic distances at x = 0.13 in Figure 4 and Figure 5 is due to the effect of a single uranium substitution in a small simulation supercell. However, the 1x1x2 cell is able to provide a good representation of the solid solution at higher levels of uranium substitution, when the concentration is more balanced.

. | Th-Th . | Th-U . | Th-O . | U-O . |
---|---|---|---|---|

ThO_{2} | 3.961 | — | 2.425 | — |

Th_{7} U_{1} O_{16} | 3.978 | 3.970 | 2.437 | 2.422 |

Th_{31} U_{1} O_{64} | 3.960 | 3.948 | 2.426 | 2.389 |

Th_{28} U_{4} O_{64} | 3.952 | 3.940 | 2.423 | 2.386 |

. | Th-Th . | Th-U . | Th-O . | U-O . |
---|---|---|---|---|

ThO_{2} | 3.961 | — | 2.425 | — |

Th_{7} U_{1} O_{16} | 3.978 | 3.970 | 2.437 | 2.422 |

Th_{31} U_{1} O_{64} | 3.960 | 3.948 | 2.426 | 2.389 |

Th_{28} U_{4} O_{64} | 3.952 | 3.940 | 2.423 | 2.386 |

As noted, to obtain the interatomic distances, we have taken the average values over all the inequivalent configurations. Overall, our Th-Th distances are slightly overestimated compared to those distances found by Hubert et al at low concentrations of U, and slightly underestimated at high U concentrations. The U–U distances are more consistent but are underestimated across the full range of solid solutions. Hubert et al did not publish data on Th-U distances for us to use for comparison, but we have found in our calculations that the Th-U distances are greater than those of U-U but less than Th-Th for each value of x. In absolute terms all of the calculated cation-cation distances are extremely close to experiment and consistent with the fact that U is the smaller cation and that its incorporation into thoria creates stress in the lattice to maintain the four-fold coordination of the oxygen in the fluorite structure. The U–U separation is larger at low values of *x*. As typically only small amounts of uranium are added to a thorium-uranium MOX fuel, the larger U–U separation at these concentrations indicates that uranium does not cluster within the thoria lattice. The Th-O and U-O distances show similar behavior to the Th-Th distances, as expected from the smaller ionic radius of uranium compared to thorium.

To analyze the Boltzmann probability distributions, it is now useful to consider the larger 2x2x2 supercell. By using the larger cell, we can not only analyze the system at the low concentrations of uranium that are of particular interest in MOX fuels, but we can also more clearly observe any changes in the lattice that may occur on the addition of uranium. In Figure 6, with a mole fraction of uranium of x = 0.06, we show that configurations 4 and 5 together are clearly dominant, with probabilities of 38 and 41%, respectively. The corresponding structures may be seen in Figure 8. As the concentration of uranium and the total number of inequivalent configurations increases for Th_{29}U_{3}O_{64}, there are no longer just one or two dominant configurations as illustrated in Figure 7.

To further analyze the relationship between the structures of the inequivalent configurations and the probabilities that a given configuration will occur at equilibrium, it is useful to examine the energies of each symmetrically inequivalent configuration. In Table VI, the energies of the independent configurations obtained from GULP are listed relative to the lowest energy configuration of the Th_{30}U_{2}O_{64} (uranium mole fraction of x = 0.06) system. As shown in Table VI, the energies of the different configurations do not vary much, indicating that the U distribution can be expected to be fairly random. Following the simple arguments by McLean^{47} for metals and confirmed for ionic materials,^{48} that dopant incorporation is based on the reduction of elastic strain, the smaller U cation should be easily incorporated in the thoria lattice, which is confirmed by the lack of site preference shown by the U cations in the thoria lattice.

Configuration . | Degeneracy (Ω) . | E_{m} eV
. | E_{m(red)} eV
. |
---|---|---|---|

1 | 48 | 0.0147 | 0.0667 |

2 | 48 | 0.0187 | 0.0707 |

3 | 16 | 0 | 0.0993 |

4 | 192 | 0.0101 | 0.0024 |

5 | 192 | 0.0077 | 0 |

Configuration . | Degeneracy (Ω) . | E_{m} eV
. | E_{m(red)} eV
. |
---|---|---|---|

1 | 48 | 0.0147 | 0.0667 |

2 | 48 | 0.0187 | 0.0707 |

3 | 16 | 0 | 0.0993 |

4 | 192 | 0.0101 | 0.0024 |

5 | 192 | 0.0077 | 0 |

However it is important to consider the degeneracy, Ω, of the system, which here is the number of symmetrically equivalent configurations that can be represented by one inequivalent configuration. The energy of the degenerate system, E_{m(red)}, considers the degeneracy of the system and therefore differs from the energy of the corresponding structure in the full configurational space (equations (6) and (7)). Therefore the probability of finding any one configuration must consider the energy of the degenerate system (E_{m(red)}) and not simply the energy of the configuration alone.

For Th_{30}U_{2}O_{64}, as noted in Table IV, there are 496 initial configurations which, in the reduced configurational space, yield five inequivalent configurations. The degeneracy and the relative energy of each configuration with and without considering the degeneracy (E_{m(red)} and E_{m}, respectively) are presented in Table VI. Configuration 3, which is the lowest probability configuration seen in Figure 6, is the structure with the lowest energy originally (E_{m} = 0). However, as this configuration only occurs 16 times in total, compared to the Ω values of 48 for configurations 1 and 2 and 192 for configurations 4 and 5, configuration 3 has the highest E_{m(red)} and therefore the lowest probability.

In the two dominant configurations at x = 0.063 and the lowest energy structure of x = 0.094 (Figure 9), there is a distortion of the lattice. This same distortion was observed in 2x2x2 supercells at the other concentrations considered, and can be explained by the changed cation-cation and cation-anion distances in the lower energy configurations as compared to the original lattice positions and the highest energy (lowest probability) structures. Due to the U-U and U-O distances being shorter than the corresponding Th-Th and Th-O distances, as uranium substitutions are made in the cell, oxygen atoms move closer to the U atom. This moves some of the nearby Th atoms away, as the distance between Th and oxygen atoms is preferably longer. These movements create local distortions in those regions of the cell, particularly where two uranium atoms are close together, as observed in Figure 8a and b and Figure 9. When the energies of the Th_{29}U_{3}O_{64} (U mole fraction of x = 0.094) are considered along with the degeneracy of each independent configuration (Table VII), Configuration 6 has the highest relative energy, both including the degeneracy of the configuration and without. The lowest energy configuration is Configuration 10, but the configurations with the lowest degenerate energy, Ed, are configurations numbered 13, 12, and 9. These three structures each occur 768 times in the full configurational space, the highest degeneracy of any of the 14 independent configurations. These configurations are therefore also the three highest probability configurations, with probabilities of 17%, 16%, and 14% respectively. However, as Figure 7. illustrates, no single configuration is dominant at a uranium mole fraction of x = 0.094.

Configuration . | Degeneracy (Ω) . | E_{rel} eV
. | Ed_{rel} eV
. |
---|---|---|---|

1 | 96 | 0.0304 | 0.1122 |

2 | 96 | 0.0158 | 0.0977 |

3 | 384 | 0.0173 | 0.0394 |

4 | 384 | 0.0123 | 0.0344 |

5 | 384 | 0.0147 | 0.0368 |

6 | 32 | 0.0388 | 0.168 |

7 | 192 | 0.0212 | 0.0732 |

8 | 192 | 0.0165 | 0.0685 |

9 | 768 | 0.0188 | 0.0111 |

10 | 384 | 0 | 0.0221 |

11 | 256 | 0.0122 | 0.0517 |

12 | 768 | 0.0102 | 0.0024 |

13 | 768 | 0.0078 | 0 |

14 | 256 | 0.0054 | 0.0449 |

Configuration . | Degeneracy (Ω) . | E_{rel} eV
. | Ed_{rel} eV
. |
---|---|---|---|

1 | 96 | 0.0304 | 0.1122 |

2 | 96 | 0.0158 | 0.0977 |

3 | 384 | 0.0173 | 0.0394 |

4 | 384 | 0.0123 | 0.0344 |

5 | 384 | 0.0147 | 0.0368 |

6 | 32 | 0.0388 | 0.168 |

7 | 192 | 0.0212 | 0.0732 |

8 | 192 | 0.0165 | 0.0685 |

9 | 768 | 0.0188 | 0.0111 |

10 | 384 | 0 | 0.0221 |

11 | 256 | 0.0122 | 0.0517 |

12 | 768 | 0.0102 | 0.0024 |

13 | 768 | 0.0078 | 0 |

14 | 256 | 0.0054 | 0.0449 |

For all concentrations of uranium investigated within the 2x2x2 supercell, there is no single dominant configuration. Above x = 0.063, there are no longer even two dominant structures and instead several configurations share nearly equal probabilities and none of these probabilities exceed 17%. We do find that in the highest probability (lowest energy) configurations at all values of x, the U–U distances are smaller than these distances in the lowest probability (highest energy) configurations, but since there is a large number of possible independent configurations available in the system and no single one is dominant, including the lowest energy configuration, we predict that uranium atoms will be distributed throughout the cell without forming clusters.

### D. Thermodynamic Properties of Th_{(1−x)}U_{x}O_{2}

From the average configurational enthalpies computed with the SOD code, using equation (5). we have calculated the enthalpies of mixing (ΔH_{mix}) at several temperatures for a range of solid solutions of the 1x1x2 supercell system of Th_{(1−x)}U_{x}O_{2}. By using a smaller supercell than that evaluated for the configurational energies and probability, the number of inequivalent configurations and hence computational cost was reduced and the full solid solution could be calculated. All values of ΔH_{mix} are positive as shown in Figure 10, indicating that mixing of the two solids is endothermic. For an ideal solid solution, by definition ΔH_{mix} = 0, whereas Th_{(1−x)}U_{x}O_{2} deviates from an ideal solid solution, with a maximum ΔH_{mix} of approximately 2.5 kJmol^{−1} near the midpoint of the solid solution range. From the values of ΔH_{mix}, we have calculated the values of ΔG_{mix} (Figure 11), to see if mixing is spontaneous at temperatures above 298 K. To determine the vibrational contributions to the free energy of the system, the k-points were checked in GULP until the energy converged. We have used a k-points mesh of 2x4x4 to calculate the phonon frequencies of Th_{(1−x)}U_{x}O_{2}. When vibrational effects are included at 600 K (Figure 12), ΔG_{mix} is lowered, showing the importance of including vibrational contributions to ΔG_{mix} to obtain more accurate values. On close examination of the mixing curves in Figure 12, a miscibility gaps is observed. However, the gap is only 0.1 kJmol^{−1} without consideration of vibrational effects, increasing to just 0.05 kJmol^{−1} when vibrational contributions are included. As this miscibility gap is so small in energetic terms, we do not expect to observe the formation of two separate domains at 600 K. Additionally, a slight shift towards the lowest values of x is observed in the ΔG_{mix} curves in both Figure 11 and Figure 12 with the more thorium-rich Th_{(1−x)}U_{x}O_{2} systems having a lower ΔG_{mix} than uranium-rich Th_{(1−x)}U_{x}O_{2}.

Understanding the heat capacity and thermal expansion is important in any nuclear fuel material in order to predict the temperature of the fuel in the reactor. We have calculated these properties for the full range of the ThO_{2}-UO_{2} solid solution in the 1x1x2 supercell and compared our values to those found experimentally or using other computational methodologies where available. The results of the constant-volume heat capacity calculations for the pure materials and the solid solutions are presented in Figure 13 and agree closely experiment. As the C_{v} values of ThO_{2} we calculated are in line with the DFT values and have used a less compute-intensive methodology than the value found by Lu et al^{44} at all temperatures studied here, we expect that this methodology will be useful in determining C_{v} of even lower concentrations of uranium in the ThO_{2}-UO_{2} solid solution, or other nuclear materials.

The coefficient of thermal expansion, α, of pure thoria has been discussed in section III B but in MOX fuels it is important to understand the effect of both cations on the thermal expansion. Generation IV reactors are typically designed to withstand temperatures up to 1873 K, so we have calculated α from room temperature (298 K) to 1900 K to represent typical operating temperatures of this fuel type. In Figure 14, we present the results of Th_{(1−x)}U_{x}O_{2} where x = 0.125, 0.25, and 0.38 compared to ThO_{2}. The values of α are in good agreement with the pure material at all temperatures calculated. There is particularly close agreement at low temperatures, with deviation from the pure material increasing above approximately 1000 K.

## IV. CONCLUSIONS

We have developed a new interatomic shell-model potential that can be used to accurately model the structure and properties of pure ThO_{2}. This potential can be combined with a leading UO_{2} potential to investigate uranium-doped thorium dioxide. In this study we have employed these potentials in a study of the distribution of uranium atoms in thoria supercells and we have closely reproduced available experimental data for the lattice parameters, coefficient of thermal expansion, and constant volume heat capacity.

We have shown that uranium atoms do not cluster in the bulk material but are instead distributed throughout the cell. In a Th_{(1−x)}U_{x}O_{2} solid solution the enthalpy of mixing is endothermic across the full range of the solution, but while the Gibbs free energy of mixing is favorable above 600 K, when vibrational effects are included in the calculations, the free energies of mixing show a very small miscibility gap. However, as the energy penalty is so small, we do not predict the formation of separate domains at high or low concentrations of uranium in the thoria lattice.

## ACKNOWLEDGMENTS

Via our membership of the UK’s HPC Materials Chemistry Consortium, which is funded by EPSRC (**EP/L000202**), this work made use of the facilities of HECToR and ARCHER, the UK’s national high-performance computing service, which is funded by the Office of Science and Technology through EPSRC’s High End Computing Programme. AS gratefully acknowledges funding from the Molecular Modeling and Materials Science Centre for Doctoral Training at UCL, and NHdL thanks the Royal Society for an Industry Fellowship and AWE for funding.