## Abstract

Using molecular dynamics, the thermophysical properties of the (U_{x},Th_{1−x})O_{2} system have been investigated between 300 and 3600 K. The thermal dependence of lattice parameter, linear thermal expansion coefficient, enthalpy and specific heat at constant pressure is explained in terms of defect formation and diffusivity on the oxygen sublattice. Vegard's law is approximately observed for solid solution thermal expansion below 2000 K. Different deviations from Vegard's law above this temperature occur owing to the different temperatures at which the solid solutions undergo the superionic transition (2500–3300 K). Similarly, a spike in the specific heat, associated with the superionic transition, occurs at lower temperatures in solid solutions that have a high U content. Correspondingly, oxygen diffusivity is higher in pure UO_{2} than in pure ThO_{2}. Furthermore, at temperatures below the superionic transition, oxygen mobility is notably higher in solid solutions than in the end members. Enhanced diffusivity is promoted by lower oxygen-defect enthalpies in (U_{x},Th_{1−x})O_{2} solid solutions. Unlike in UO_{2} and ThO_{2}, there is considerable variety of oxygen vacancy and oxygen interstitial sites in solid solutions generating a wide range of property values. Trends in the defect enthalpies are discussed in terms of composition and the lattice parameter of (U_{x},Th_{1−x})O_{2}.

## 1. Introduction

UO_{2} has been studied extensively as the main component of conventional nuclear fuel. It is also blended with other actinide oxides, such as ThO_{2} [1] and PuO_{2} [2,3], forming mixed oxide (MOX) fuel. Alternatively, long-lived minor actinides can be separated from nuclear waste and blended with UO_{2} or MOX for transmutation in a reactor- [4] or accelerator-driven system [5,6]. Owing to its relative abundance, thorium is considered as an important candidate for MOX, whereby Th^{232} transmutates in reactor to U^{233}, which can then undergo fission. As such, it is important to understand the underlying mechanisms that govern the thermophysical and diffusion properties in mixed oxides, owing to a non-uniform cation sublattice.

For UO_{2}, there is a deviation from linear thermal expansion and a classical Debye description of the constant pressure specific heat above 1300 K [7–16]. At 2670 K (0.85T_{m}), there is a peak in specific heat owing to a pre-melting transition or superionic transition as seen in other fluorite structures [7]. Below the transition, it is not yet clear to what extent the excess specific heat and thermal expansion is driven by oxygen disorder versus electronic-defect contributions, or over what temperature ranges these effects may dominate.

A great deal of work has been carried out using atomistic simulation to study candidate actinide oxide components of nuclear fuel [17–21]. The ability of an interatomic potential to accurately reproduce the thermophysical properties of UO_{2}, such as lattice parameter, elastic constants, thermal conductivity and specific heat over a wide range of temperatures, has often been used as a key discriminator for the suitability of a potential set [17–20]. Similarly, Potashnikov *et al.* [22] compared the ability of a number of interatomic potentials to predict oxygen diffusivity in UO_{2} using molecular dynamics (MD). However, experimental data for oxygen diffusion will be influenced by the presence of point defects arising owing to materials processing conditions, and it not necessarily comparable with the perfect crystal calculations of Potashnikov [22]. For example, the enhancement of oxygen diffusivity owing to non-stoichiometry in UO_{2} has been demonstrated recently by Govers *et al.* [23] and shown experimentally by Belle [24]. Similarly, enhanced diffusivity owing to Schottky defects was also identified by Potashnikov *et al.* [22]. Therefore, it is not surprising that many simulations on perfect crystals predict lower oxygen diffusivity than experiment.

Recently, a potential set has been derived that accurately reproduces a wide range of thermomechanical and thermophysical properties for AmO_{2}, CeO_{2}, CmO_{2}, NpO_{2}, PuO_{2}, ThO_{2} and UO_{2}, between 300 and 3000 K [25]. In particular, this potential accurately represents the individual elastic constants of the actinide oxides and reproduces the Cauchy violation (C_{12}≠C_{44}) by introducing many-body interactions using the embedded atom method (EAM) [26] without the necessity for the shell model [27]. As a result, a significant improvement in the ability of empirical interatomic potentials to reproduce the bulk modulus over a large range of temperatures has been achieved. Importantly, this potential set employs the same description of oxygen–oxygen interactions throughout, enabling the simulation of actinide oxide solid solutions. Furthermore, it accurately reproduces the melting points of UO_{2} and ThO_{2} well, making it particularly suitable for investigating (U_{x},Th_{1−x})O_{2} solid solutions.

Here, we investigate, using atomistic simulation, the lattice parameter, linear coefficient of thermal expansion, enthalpy and specific heat at constant pressure for (U_{x},Th_{1−x})O_{2} between 300 and 3600 K for *x* = 0.00, 0.25, 0.50, 0.75 1.00. Furthermore, the influence of solid solution composition on oxygen-defect formation, oxygen diffusivity and the superionic transition is reported.

## 2. Methodology

MD simulations are carried out using LAMMPS [28], and the set of interatomic potentials derived previously [25].^{1} The model combines a pair potential description of each system with the many-body EAM description of Daw & Baskes [26]. As such, the potential energy, *E*_{i}, of an atom *i* with respect to all other atoms can be written as
*i* and *j*, separated by *r*_{ij}, is given by *ϕ*_{αβ}(*r*_{ij}). This has both long range electrostatic, *ϕ*_{C}(*r*_{ij}), and short range contributions. Coulombic interactions are calculated using the Ewald method [31] with the particle–particle particle–mesh (PPPM) implementation of the method being adopted within MD calculations in order to improve computational efficiency [32]. The short range contributions are described using Morse, *ϕ*_{M}(*r*_{ij}), and Buckingham, *ϕ*_{B}(*r*_{ij}), potential forms, as given by equation (2.2) [33,34]. *α* and *β* are used to label the species of atom *i* and atom *j*, respectively,
*A*_{αβ}, *ρ*_{αβ}, *C*_{αβ}, *D*_{αβ}, *γ*_{αβ} and *r*_{0} are empirical parameters that describe the pair interactions between atom *i* and atom *j*. These have been reported previously for AmO_{2}, CeO_{2}, CmO_{2}, NpO_{2}, PuO_{2}, ThO_{2} and UO_{2} [25]. For the study of solid solutions, further mixed cation–cation pair interactions (e.g. *ϕ*_{U−Th}) must also be defined. As was the case for the self cation–cation pair interactions [25], the mixed cation–cation interactions are dominated by Coulombic interactions at the separations exhibited by the fluorite structure. Hence, it is not possible to fit these parameters to experimental bulk properties. Instead, the assumptions made previously for self-interactions are extended here to mixed cation–cation pair potentials. The description of covalency predicted by the Morse potential is not required for cation–cation pairs and is therefore excluded for these interactions (e.g. *D*_{U−Th}=0 *eV*). The pre-exponential term of the Buckingham potential is the same for all cation–cation pairs and is based on the parameter reported by Grimes & Catlow [35] (e.g. *A*_{U−Th}=18600 *eV*). However, the reported self cation–cation *ρ*_{αα} parameters are scaled to cation radius [25]. Equation (2.6) [36] is used to determine *ρ*_{αβ} parameters for mixed cation–cation pairs whose values are reported in table 1,

The second term in equation (2.1) uses the EAM to introduce a many-body perturbation to the more dominant pairwise interactions. The derivation of the parameters and a description of the functional terms used in the EAM component are given in reference [25].

The thermal expansion, specific heat and oxygen diffusivity in the (U_{x},Th_{1−x})O_{2} system are investigated for the UO_{2}, (U_{0.75},Th_{0.25})O_{2}, (U_{0.5},Th_{0.5})O_{2}, (U_{0.25},Th_{0.75})O_{2} and ThO_{2} compositions. Solid solution crystal structures are created by randomly distributing U^{4+} and Th^{4+} cations on the 4a Wyckoff sites (fluorite actinide sites) throughout a supercell of 10×10×10 fluorite unit cells. These structures are equilibrated for 40 ps for temperatures between 300 and 3600 K at 25 K intervals with the thermophysical properties (lattice parameter and enthalpy) obtained from averages taken over the final 2 ps of the simulation. A 2 fs timestep is used in the NPT ensemble with Nosé–Hoover thermostat and barostat times of 0.1 and 0.5 ps, respectively. For each composition, this is repeated for 10 randomly generated structures.

The equilibrated solid solution structures are also used to determine oxygen diffusivity. The oxygen mean-squared displacement (MSD) is calculated in the NVE ensemble for 1 ns with a 1 fs timestep for a range of temperatures from 2000 to 3600 K at 100 K intervals. From this the diffusivity, *D*, is calculated using the following equation [37],
*t* is the simulation time.

The point defects present in a simulation box are counted by analysing structural information extracted from the configuration. The following procedure is carried out independently on the two sublattices present (oxygen and actinide). The average distribution of atoms around a site of the sublattice *A* is calculated. This is a three-dimensional equivalent of the partial radial distribution function:
*μ*_{A}(**r**) is proportional to the probability of finding an atom of sublattice *A* in a volume *dV* around a position **r** from a lattice site; *i* and *j* are two atoms of sublattice *A* and *N*_{A} is the number of atoms on this sublattice. The local maxima of this scalar field indicate the positions most likely to be occupied by a neighbour. When all the sites on a sublattice have the same local symmetry, these are the positions of the neighbours {**R**_{A,λ}}, with *λ*∈{1,…,*Z*^{A}}, *Z*^{A} being the connectivity of the sublattice (number of neighbours for each site). For the fluorite structure, *Z* is 6 for the oxygen sublattice (simple cubic), and 12 for the actinide sublattice (face centred cubic). The local maxima of *μ*_{A}(**r**) have a width that can be estimated from the full width at half maximum (FWHM). This is related to the fluctuations of the neighbour positions caused by thermal oscillations around the lattice sites. Thus, one can define a distance criterion *δ*_{A} based on the FWHM, so that an atom *j* is said to occupy the neighbour site *λ* around the atom *i*, if
*Z*_{i}, of the atom *i* is then the number of such valid neighbours.

Having calculated *Z* for each atom, those that have the same connectivity as their sublattice are in a perfect environment and ignored for the rest of the defect counting procedure. The vacancies are detected using the fact that several atoms would have a common missing neighbour. A virtual atom, not present in the real configuration, is added to such positions and treated as a regular atom during the remaining defect detection. The last step is repeated until no valid vacancy position is found. At this point, the atoms that have a low connectivity are either part of extended defects, which are ignored in this study, or interstitial atoms. The interstitials are characterized by a connectivity number lower than 3.

This procedure gives the number of vacancy and interstitial defects. When carried out on successive snapshots of constant-temperature MD simulations, the concentration of defects can be estimated as a function of temperature. In this work, snapshots are taken at 1 ps intervals for simulations lasting 25 ps equilibrated at temperatures ranging from 2300 to 3200 K every 100 K. The defect concentrations are then averaged over all snapshots for a given composition to obtain the thermally equilibrated defect concentrations in each solid solution.

For the energy minimization calculations of isolated oxygen vacancy and interstitial defect enthalpies, the same supercells described above are employed. Given the complex structures of these solid solutions, the perfect supercells are subjected to a rigorous energy minimization procedure to ensure they are fully relaxed. This consists of an initial minimization under constant volume conditions using a damped dynamics algorithm [38] followed by a constant pressure step using a conjugate gradient method before a final optimization step employing a steepest descent procedure with fixed lattice parameters. Once the solid solution supercells are fully optimized, point defects were introduced into the simulation supercells by either removing (vacancy) or adding (interstitial) oxygen atoms into the supercell. The defective supercells are energy minimized with the lattice parameters fixed in order to represent the dilute limit with the defect enthalpy, *dE*, calculated using
*H*_{perfect} and *H*_{defect} are the total enthalpies of the perfect and defective cells, respectively. The oxygen Frenkel enthalpies for UO_{2} and ThO_{2} are converged to within 0.1 eV for the 10×10×10 supercell compared with the fully isolated enthalpies given previously [25]; therefore, the defect enthalpies are considered to be converged with respect to system size.

For the oxygen vacancy simulations, an oxygen is removed from each of the oxygen lattice sites, in all of the 10 simulation supercells for each composition, resulting in a total of 80 000 defect simulations. Similarly, the oxygen interstitial defect enthalpy is calculated at every possible interstitial site, in all of the supercells, leading to a total of 40 000 defect simulations. This large number of simulations allows us to access the statistical distribution of the defect enthalpies arising from the random arrangement of cations on the 4a Wyckoff sites. As this approach generates a very large number of unique defect energies, the dataset has been grouped into bins of width 0.01 eV for ease of manipulation and to enable useful presentation of these results.

## 3. Results and discussion

### (a) Oxygen diffusivity

By assuming an Arrhenius relationship, equation (3.1), *D* is plotted logarithmically as a function of 1/*T*, so that, the gradient is proportional to the activation enthalpy, *H*_{a}.
*D*_{0} is the pre-exponential term, *k*_{B} is the Boltzmann constant and *T* is temperature. For each composition, *D* is averaged over all 10 randomly generated structures and plotted in figure 1*a*; error bars indicate the standard deviation. Regions of constant gradient, and thus activation enthalpy, indicate temperature regimes with a common diffusion mechanism. As in the study of Potashnikov *et al.* [22], figure 1*a* shows that the transition between the fully crystalline low temperature and superionic high temperature behaviour occurs over a range of temperature specific to each composition. Similarly, figure 2 highlights the change in activation enthalpy during the transition. It can be seen that the transition occurs at a higher temperature in ThO_{2} compared with UO_{2}. However, the addition of thorium to form (U_{0.75},Th_{0.25})O_{2} and (U_{0.5},Th_{0.5})O_{2} solid solutions does not appear to significantly increase the superionic transition temperature (and may even have reduced it for (U_{0.75},Th_{0.25})O_{2}). It is not until the (U_{0.25},Th_{0.75})O_{2} solid solution that the thorium content is significant enough to increase the superionic phase transition temperature. In table 2, the approximate range of superionic transition temperatures has been summarized. The significance of this for the high temperature thermophysical properties is discussed in §3*d*,*e*.

Figure 1*b* shows the oxygen diffusivity as a function of uranium composition for a range of temperatures. There is a clear enhancement of oxygen diffusivity below the superionic transition temperature for the solid solutions compared with pure UO_{2} and ThO_{2}. From about 2200 to 2700 K, oxygen diffusivity is greatest in (U_{0.75},Th_{0.25})O_{2}. For 2000 and 2100 K, (U_{0.5},Th_{0.5})O_{2} exhibits the highest oxygen diffusivity and (U_{0.25},Th_{0.75})O_{2} is approaching the same level as pure UO_{2}. Furthermore, if the trend for (U_{0.25},Th_{0.75})O_{2} in figure 1 continues below 2000 K, it is possible that the oxygen diffusivity may also exceed that of UO_{2}; however, further diffusivity calculations at lower temperatures must be carried out in order to confirm this prediction. The role of such solid solutions in enhancing oxygen diffusivity in mixed oxide fuels is therefore expected to be significant at reactor operating temperatures and may impact the release of fission products that occupy and migrate via the oxygen sublattice (e.g. I^{−}) [35].

### (b) Oxygen-defect concentrations

Oxygen-defect concentrations underpin the thermophysical properties of actinide oxides as well as the transition to superionicity. It is, therefore, important to calculate the defect concentrations that are generated during our high temperature simulations in order to correctly identify which defects are responsible for the high temperature phenomena. The defect concentrations are calculated using a method developed for the purpose of this study. Figure 3 shows the oxygen interstitial concentration (fraction of oxygen ions off their lattice sites) plotted logarithmically as a function of 1/*T*, such that the gradient is proportional to the defect formation enthalpy for an interstitial–vacancy pair. Below the superionic transition, the defect populations were dominated by tightly bound pairs of oxygen vacancies and interstitials, such that _{2} and 4.70 eV for ThO_{2}, corresponding to the gradients of figure 3, are therefore not equivalent to those reported previously [25] that were calculated for stable Frenkel pairs (which have a greater vacancy–interstitial separation). This spontaneous recombination of first-neighbour pairs on the oxygen sublattice has also been observed in UO_{2} using other empirical potentials [39,40] and in other similar structures, such as pyrochlores [41]. Figure 3 shows that at higher temperatures the fraction of oxygen anions in these pseudo-Frenkel positions increases. As a consequence, there is an increase in system enthalpy and lattice parameter associated with the pseudo-defect enthalpies and psuedo-defect volumes, respectively (see §3*d*,*e*). As the oxygen disorder increases during the superionic transition (table 2), the defect counting method can no longer identify a clear oxygen sublattice and, as a result, fails to identify the correct number of interstitials. This point is characterized by a sharp decrease in the gradient of figure 3 that helps identify the superionic transition, showing that the transition temperature is lower in UO_{2} than ThO_{2} and potentially even lower in (U_{0.75},Th_{0.25})O_{2}.

### (c) Oxygen point-defect enthalpies

To study the possibility of enhanced permanent oxygen-defect concentrations in the solid solutions, it is necessary to identify the point defect formation enthalpies. Unlike the pure end member systems oxygen sites in a solid solution are different and consequently, there is a wide range of defect enthalpies owing to the various environments surrounding the vacancy. A similar observation has been made for the As vacancy in In_{x}Ga_{1−x}As [42].

Figure 4 identifies the fraction of oxygen sites that lie within 0.005 eV of a given oxygen vacancy formation enthalpy. This shows that for each composition there are five peaks that correspond to the five different first nearest cation neighbour coordination of the oxygen site, with the lowest and highest enthalpy peaks coordinated by 4 and 0 uranium ions, respectively. The skew in peak heights corresponds to the solid solution composition. For example, the proportion of sites fully coordinated by uranium ions (lowest enthalpy peak) is greatest in (U_{0.75},Th_{0.25})O_{2}. Additionally, there is a shift in the peaks owing to lattice parameter, whereby, all oxygen vacancy enthalpies are shifted down for solid solutions with a greater lattice parameter (i.e. higher thorium content; see §3*d*). Therefore, the peak corresponding to fully uranium-coordinated sites is always lower for solid solutions compared with the pure UO_{2} system (shown by the vertical red line). The lattice parameter effect is most clear for the (U_{0.25},Th_{0.75})O_{2} system despite this peak being small. Similarly, oxygen interstitial enthalpies are shifted down for solid solutions with a larger lattice parameter, as illustrated in figure 5.

At lower temperatures, the oxygen vacancies with higher formation enthalpies play a proportionately lesser role in oxygen transport than at higher temperatures. This is demonstrated by enhanced diffusivity in the solid solutions that exhibit a greater proportion of oxygen defects with lower formation enthalpies (figure 1*b*). This is seen clearly for (U_{0.75},Th_{0.25})O_{2} and (U_{0.5},Th_{0.5})O_{2}, which both have a significant number of vacancy formation enthalpies below that exhibited by UO_{2}. The prediction in §3*a* that oxygen diffusivity in (U_{0.25},Th_{0.75})O_{2} may also exceed diffusion in UO_{2} below 2000 K can now be understood as a consequence of it containing the lowest enthalpy peak in figure 4. For the full oxygen Frenkel enthalpy, oxygen interstitials must also be included. Figure 5 further supports enhanced oxygen disorder as all solid solution compositions exhibit a significant number of interstitials with lower formation enthalpies than for the end members.

### (d) Thermal expansion

Figure 6 shows the increase in the lattice parameter as a function temperature for a given composition averaged over the 10 randomly generated 10×10×10 structures. Experimental data for UO_{2} [16], ThO_{2} [43] and (U_{0.55}Th_{0.45})O_{2} [44] are also included and show very good agreement with the predictions of the potential. Figure 6 illustrates a significant increase (or ‘bump’) in thermal expansion for all compositions of solid solution as well as for the pure systems at high temperature (2300–3300 K), below this temperature, Vegard's law [45] is obeyed. The ‘bump’ in lattice parameter can first be attributed to the defect volumes associated with the formation of a large number of oxygen defects and after that by the volume change owing to the superionic transition. This is more clearly demonstrated by using the first derivative of the lattice parameter with respect to temperature to calculate the linear thermal expansion coefficient, see equation (3.2).
*L*/∂*T*)_{P}, is calculated by fitting a straight line to the lattice parameter at a given temperature and the data points at ±25 K either side. Figure 7 shows the variation of linear thermal expansion coefficient as function of temperature. The peak in linear thermal expansion coefficient corresponds closely to the range of temperatures for the superionic transition for each composition (table 2). For UO_{2}, (U_{0.75}Th_{0.25})O_{2} and (U_{0.5}Th_{0.5})O_{2}, the peak is at around 2600 K, in close agreement with the experimental value for the superionic transition temperature of 2670 K for UO_{2} [7]. For (U_{0.25}Th_{0.75})O_{2} and ThO_{2}, the peak is at approximately 2700 and 2950 K, respectively. For UO_{2}, (U_{0.75}Th_{0.25})O_{2} and (U_{0.5}Th_{0.5})O_{2} a second very high temperature peak is identified that is associated with the creation of cation defects; however, as these peaks are above the UO_{2} melting point predicted by this potential [25], it is outside the regime of interest for this study.

### (e) Enthalpy and specific heat

In addition to the volume change owing to oxygen-defect formation and the superionic phase transition, there is an associated change in system enthalpy. Figure 8 shows the enthalpy increment (increase in enthalpy with respect to standard conditions) as a function of temperature (i.e. H(T)–H(298 K)) averaged over the 10 randomly generated structures for each solid solution composition. The enthalpy increment increases approximately linearly with temperature below 1500 K. Between 2300 and 3300 K, the enthalpy increment as a function of temperature increases more significantly. The first derivative of the enthalpy increment with respect to temperature is used to calculate the specific heat capacity at constant pressure using the following relationship,
*n* is the number of moles and the first derivative of enthalpy, (∂*H*/∂*T*)_{P}, is calculated by fitting a straight line to the enthalpy at a given temperature and the data points at ±25 K either side. Figure 9 indicates a gradual increase in the specific heat until around 2000 K at which point the specific heat increases more rapidly due the enthalpy required to create oxygen disorder. The peak in specific heat is commensurate with the superionic transition (table 2), which occurs close to the same temperature as the peak in the linear thermal expansion coefficient for each composition.

Despite the omission of electronic defects from these empirical calculations, figures 6–9 show that anion disorder contributes significantly to excess thermal expansion and specific heat capacity, while also indicating the peaks are commensurate with the superionic transition. However, comparison of our results with the experimental data for UO_{2} specific heat capacity [16] indicates that oxygen disorder is not sufficient to account for excess specific heat at intermediate temperatures which may be accounted for by electronic defects [12].

## 4. Conclusions

Using MD, the superionic transition in (U_{x},Th_{1−x})O_{2} is investigated for compositions where *x* equals 0.00, 0.25, 0.50, 0.75 and 1.00. This is identified by the change in activation enthalpy, and thus diffusion mechanism, for oxygen migration. It is shown that reduced oxygen-defect enthalpies in the three solid solution compositions studied here contribute to enhanced oxygen diffusivity below the superionic transition.

The creation of oxygen pseudo-Frenkel pairs and subsequently the superionic transition causes a ‘bump’ in the lattice parameter, thermal expansion coefficient, enthalpy and specific heat capacity for all (U_{x},Th_{1−x})O_{2} compositions (including end members). The onset of the superionic transition was calculated to occur at similar temperatures for UO_{2}, (U_{0.75},Th_{0.25})O_{2} and (U_{0.5},Th_{0.5})O_{2}, although it appears to occur slightly earlier in (U_{0.75},Th_{0.25})O_{2} when the oxygen activation enthalpy is considered. The superionic transition temperatures for (U_{0.25},Th_{0.75})O_{2} and ThO_{2} are higher. The change in volume owing to the creation of oxygen disorder explains the high temperature lattice expansion, whereas the latent heat required to undergo the superionic transition is responsible for a peak in the specific heat.

The enhanced low temperature-defect formation and anion diffusion in high U content solid solutions, compared with the UO_{2} end member, has implications for the mobility of fission products, such as iodine, which may be transported via the oxygen sublattice [35].

## Funding statements

Funding for M.W.D.C was provided through the EPSRC, grant no. EP/I0364400/1, and the NDA. Computational resources are due to the Imperial College High Performance Computing Service. Funding for S.T.M. was provided through RCUK, grant no. EP/K00817X/1.

## Footnotes

↵1 Supplementary material describing the use of this potential for use in GULP [29], LAMMPS [28] and DL-POLY [30] is provided at http://abulafia.mt.ic.ac.uk/potentials/actinides.

- Received May 28, 2014.
- Accepted August 1, 2014.

© 2014 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.