A density functional theory study of uranium-doped thoria and uranium adatoms on the major surfaces of thorium dioxide

(cid:1) Uranium substitution in ThO 2 is found to increase the covalent nature of the ionic bonding. (cid:1) The (111), (110), and (100) surfaces of ThO 2 are studied and the particle morphology is proposed. (cid:1) STM images of the (111), (110), and (100) surfaces of ThO 2 are simulated. (cid:1) Uranium adsorption on the major surfaces of ThO 2 is studied. Thorium dioxide is of signi ﬁ cant research interest for its use as a nuclear fuel, particularly as part of mixed oxide fuels. We present the results of a density functional theory (DFT) study of uranium-substituted thorium dioxide, where we found that increasing levels of uranium substitution increases the covalent nature of the bonding in the bulk ThO 2 crystal. Three low Miller index surfaces have been simulated and we propose the Wulff morphology for a ThO 2 particle and STM images for the (100), (110), and (111) surfaces studied in this work. We have also calculated the adsorption of a uranium atom and the U adatom is found to absorb strongly on all three surfaces, with particular preference for the less stable (100) and (110) surfaces, thus providing a route to the incorporation of uranium into a growing thoria particle. © 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license


Introduction
Thorium dioxide is a ceramic material that has industrial applications owing to its high melting point (3300 C) and its strong insulating properties [1]. It is of particular research interest owing to its applications in nuclear reactors for power generation. Currently, reactors capable of using ThO 2 fuels include the very popular pressurized water reactors (PWRs), although most PWRs presently use uranium-based fuels. Several of the new reactor designs, the so-called Generation IV reactors, are also capable of using thorium-based fuels [2]. While thorium fuels can take the form of either the oxide or the fluoride salt, current reactor technology is heavily biased in favor of mixed oxide reactors, where, uranium and plutonium oxides are used in the form of mixed oxide fuels (MOX).
As ThO 2 , PuO 2 , and UO 2 all share the fluorite structure (space group 225), ThO 2 is a good candidate for MOX fuel reactors and offers many safety benefits. As thorium is not a fissile element, a sustained nuclear decay reaction can only occur in the presence of a neutron source. Both plutonium and uranium produce neutrons and thus are suitable in thorium-based fuels. Once the neutron source is depleted, thorium ceases to undergo further radioactive decay, thereby forestalling any reactor meltdowns. Additionally, thorium dioxide is considered to be proliferation resistant due to the extreme difficulty of isolating any potential weapons-grade nuclear material from its waste products [3,4]. Reprocessed waste uranium or plutonium from conventional MOX reactors can be used, thus reducing nuclear waste stores [5,6].
The presence of uranium in thorium dioxide can arise in two ways, either through the deliberate addition of uranium to thoria in a MOX fuel or as an impurity. For example, the 232 Th isotope is transmuted via the addition of a neutron and a series of beta decays to the fissile 233 U isotope, which undergoes the fission required for nuclear power generation (Scheme 1) [7].
In previous work, we have used a statistical mechanics approach to investigate the thermophysical properties of ThO 2 and ThO 2-eUO 2 solid solutions [8], such as thermal expansion and the enthalpy of mixing. Now we wish to investigate the specific local effects of uranium substitution in the thoria bulk and when adsorbed on the surfaces. By using Density Functional Theory (DFT) calculations, we can explore the effect of uranium dopants on the nature of the atomic bonding and the charge distribution in thorium dioxide. Any changes in the bonding interactions between the component atoms in the system will affect the mechanical properties, including the bulk modulus or the thermal expansion. These properties are important to understand in any material and are particularly critical to predicting the fuel behavior and integrity in nuclear reactors. Fuel behavior, including expansion and operating temperature must be considered in the design of any nuclear reactor. We have explored the density of states (DOS) of 2 Â 2 Â 2 supercells of the pure thoria compound and a select number of doped systems. To investigate the effects of the U in the electron cloud of the material we have obtained the DOS curve for very low Several theoretical studies have investigated the surface chemistry of ThO 2 [9e11]. An early study by Benson et al. calculated the surface energies of ThO 2 and UO 2 surfaces via two different means based on experimental estimates of the compressibility of the two materials, although the authors admitted more reliable data was needed for accurate calculations [9]. Skomurski et al. simulated the adsorption of the molecular and dissociated forms of water and oxygen on the (111) surface, finding a preference for molecular water on a clean surface, although this preference decreases as adsorbate coverage increases on the surface [10]. In this work we have used DFT to calculate the surface energies of the three low Miller index surfaces of thoria and to predict scanning tunneling microscopy (STM) images of these surfaces.

Simulation of bulk thoria
ThO 2 (Fig. 1) has the fluorite structure with a lattice parameter of 5.6 Å [12]. In this work, we have generated 1 Â 1 Â 2 and 2 Â 2 Â 2 supercell structures of thoria, containing 24 and 96 atoms respectively. We have placed between 1 and 3 uranium atoms in each 2 Â 2 Â 2 supercell to simulate low-level uranium doping. These supercells were then optimized using spin-polarized density functional theory (DFT) calculations implemented within the Vienna Ab Initio Simulation Package (VASP) [13e16]. The exchangecorrelation functional of Perdew, Burke, and Ernzerhof (PBE) [17] was used with the projector augmented wave (PAW) method [18,19] to describe the core electrons and the interaction between the atomic cores and the electrons in the valence shells. For oxygen, the 1s electrons were kept frozen during the simulations and for thorium and uranium, electrons were frozen up to and including the 4f and 5d orbitals. A kinetic energy cut-off of 600 eV was used for the plane wave basis set expansion. A gamma-centered Monkhorst-Pack k-points mesh of 6 Â 6 Â 6 was used after both were carefully checked for convergence.
We have not considered spin-orbit coupling in our calculations, following previous works, where l-s coupling was found to be important for the correct reproduction of the lattice parameters, but not for the ground state energy of the bulk UO 2 material [20]. Additionally, another study determined that the orbital contribution to the magnetic moment was significant for Co ions in the bulk and sub-surface layers of CoO, but that it vanished for the atoms at the surface layer [21]. Thus, we expect that the l-s coupling can be ignored in our calculations of U adatoms without affecting significantly the accuracy of the results. However, more research into the value of the orbital contribution to the magnetic moment is needed with regard to uranium atoms in the bulk, sub-surface and surface of uranium-bearing compounds or oxides (such as UO 2 ), but a study like that is beyond the scope of this work.
To correct for the strong on-site Coulomb repulsion arising from the presence of f electrons in the actinides and to reduce hybridization with the oxygen atoms, a Hubbard Hamiltonian following the approach of Dudarev et al. was used with the Generalized Gradient Approximation (GGA þ U) where: with r s representing the density matrix, whereas U and J here are the elements of the spherically averaged matrix of the screened Coulomb interaction [17,22,23]. For the purpose of this study, U and J are not treated separately and instead the effective U parameter is considered, where: The transmutation of the most common thorium isotope into the fissile 233 U isotope.
The value of the Hubbard parameter, U eff , used for uranium was 4 eV, a commonly used value in previous DFT studies, including those by Shi et al. and Devey et al., which modeled the optical and elastic properties of UO 2 in good agreement with experiment [24,25]. For Th, we investigated many different values of U eff using both the local-density approximation þ Hubbard U correction (LDAþU) and two exchange-correlation GGAþU functionals. However, we found that any value of U eff used decreased the accuracy of the structural prediction as a whole, while offering no improvement in the density of states as compared to previous studies of the density of states of pure thoria found in the literature, when employing the GGA approximation in VASP. Fig. 2 shows that a U eff value of 0 eV is preferable for ThO 2 for the modeling of the structural parameters and that of the three exchange-correlation functionals tested, LDA significantly underestimates the experimental lattice parameter of 5.6001 Å [12] as compared to GGA functionals.
Wang et al. [26] successfully modeled ThO 2 using GGA and Lu et al. [27] similarly concluded that no U correction was needed for ThO 2 , which was also the case in DFT studies of other thorium containing compounds, e.g. ThN which was modeled using GGAþU by Lu et al. [28]. All of these results agree with our own conclusion that no additional structural accuracy resulted from the addition of the Hubbard U parameter.
Using the Site Occupancy Disorder (SOD) code [29], we have generated a number of inequivalent configurations for each all possible atomic substitutions in the crystal lattice for a selected composition, e.g. U for Th in this work. Next, symmetry operators are used to reduce the configurational space to a number of inequivalent or independent configurations. These independent configurations represent all of the symmetrically equivalent configurations found in the full configurational space. In our previous work [8], we have calculated the lowest energy configuration for each value of x in Th (8Àx) U x O 16 . We have used these lowest energy structures in our simulations in the present study.

Simulation of ThO 2 surfaces
Surfaces were generated from the geometry optimized bulk unit cell structure using METADISE [30]. This approach is based on the work of Tasker on ionic crystal surfaces, where the surface is considered to be a slab of stacked planes [31,32]. Tasker described three types of surface plane stacking sequences. Each plane in a type 1 surface has anions and cations in their stoichiometric ratios and therefore there is no dipole. Type 2 surfaces consist of planes A and B, stacked in repeating units of A/B/A, such that the net charge in each repeat unit is zero and no net dipole is found perpendicular to the surface. For example, if plane A has a charge of À1 and plane B has a charge of þ2, a surface slab stacking sequence of A/B/A/A/B/ A will be non-dipolar. The final surface type, type 3, is an alternating pattern of oppositely charged planes, thus generating a dipole moment. By moving half of the ions of the same charge from the top to the bottom of the slab, a type 3 surface may be reconstructed to contain no net dipole perpendicular to the plane of the surface. It is important to avoid dipoles, as the surface energies will diverge in the presence of a dipole perpendicular to the surface planes and dipolar surfaces are not found experimentally without the occurrence of surface reconstruction or other mechanisms to cancel the dipole [33,34]. METADISE has been used previously to create reconstructed non-dipolar surfaces of many different oxides and sulfides [30,35e38].
The natural cleavage plane of fluorite and fluorite-type materials is the anion-terminated (111) surface ( Fig. 3), which is the most stable [31,32]. In this work, however, we model not only the Oterminated (111) surface orientation but also the O-terminated (100) and (110) surfaces for comparison ( Fig. 4 and Fig. 5). Surface slabs comprising of nine atomic layers were used for the (111) and (100) surfaces and a slab of six layers was used to model the (110) surface, where the (100) slab contained 48 atoms and the (110) and (111) slabs contained 72 atoms each. A 6 Â 6 Â 1 grid of k-points was used for all surface simulations. The convergence of the energy with respect to the vacuum size was carefully checked in VASP, with other settings kept the same as for the bulk thoria simulations. A vacuum layer of 0.15 nm was sufficient to isolate one slab and the adsorbate, where applicable, from its image in the periodic calculations. Whereas for bulk calculations, all ions were allowed to fully relax, for the surfaces, only the ions in the top layers of the slab were relaxed, while the rest were held fixed at their bulk relaxed positions. For the (111) and (100) surfaces, the top five layers were relaxed and for the (110) surface, the top three layers were relaxed. Although all surfaces were intrinsically non-dipolar, following the approach by Tasker [31], we have applied dipole corrections to our surface calculations to aid convergence during the geometry optimizations [39,40]. We have used the Bader [41] charge analysis code of Henkelman et al.
[42e44] to analyze the charge distribution in our systems.

Calculation of surface energies
We can calculate the surface energies of the (111), (110), and (100) surfaces from the surface slab energies. The surface energy represents the excess energy of the surface as compared to the original bulk material and thus serves as a measure of the stability of the crystal surface relative to the bulk. The cleavage of the bulk to form the surface reduces the coordination of those atoms located in the surface slab. We can calculate the surface energy before relaxation, g u , for a stoichiometric and symmetric surface, with the equation: where E slab,u is the energy of the unrelaxed surface slab, E bulk is the energy of a bulk system comprised of the same number of atoms as the surface slab, and A is the surface area of the top and bottom of the slab, which is why A is doubled in the equation. surface slabs were then allowed to relax while the bottom part of the slab was kept fixed. This allows us to derive the relaxed surface energy, g u , given by the equation: where E slab,r is the total energy of the relaxed surface slab, where only the top layers of the simulation slab have been allowed to relax. The relaxed surface energy will be lower than the unrelaxed surface, as all surface atoms will be in their optimized positions. It is also useful to consider surface energies in the context of the energy required to cleave the bulk material to form two equivalent surfaces. The cleavage energy, E c is related to the relaxed surface energy by the equation; The difference in energies between the relaxed and unrelaxed surface energies represents the extent of relaxation that each surface undergoes, i.e. it indicates how much the atomic positions have changed during the relaxation calculations, which can be expressed as a percentage via the equation: We have used the relaxed surface energies to calculate the equilibrium morphology using Wulff's Theorem [45], which states that the distance between the surface and the center of the particle is proportional to the surface energy. For equilibrium morphologies, the Gibbs-Wulff theorem applies, where the crystal particle at thermodynamic equilibrium will form the shape that minimizes the Gibbs free energy of the surface [46]. This approach has been used successfully to predict crystal morphologies for a range of materials including sulfides [37,47], oxalates [48], and oxides [49e52].

Uranium adatom surfaces
To calculate the energy of adsorption of a uranium adatom, we placed a single uranium atom above the ThO 2 surface slab and performed a geometry optimization calculation. The initial unrelaxed positions of U may be seen in Figs. 6e8. The adsorption energy (E ads ) for one atom of the U adsorbate can be calculated using the equation: where E slabþU is the energy of the slab with the uranium adatom determined from a single point calculation, E slab is the energy of the surface slab without an adsorbate, and E U is the energy of a uranium atom in isolation. We have chosen to use U in isolation rather than in the bulk due to the well-known complications and lack of consensus [53,54] in accurately calculating the energy of metallic uranium with DFT.

Simulation of the STM images
We have calculated the STM images of the three surfaces, which in the absence of experimental STM images of the surfaces of thorium dioxide will provide insight into the surface structures and will allow a comparison with any future experiments on thoria surfaces.
We have used the STM simulation program HIVE to predict STM images from calculations of the partial charge density in VASP [55]. By applying the Tersoff-Hamann model, which considers the local density of states (LDOS) to be proportional to the tunneling current at a given position of the STM tip [56], HIVE integrates the calculated partial charge density from the Fermi energy to a sample bias (À2.5 eV in the present study). The STM tip is assumed to be an infinitely small point that varies in height above the surface in such a way that the charge density remains constant. A constant charge density represents a constant LDOS, and so follows the Tersoff-Hamann formalism. The HIVE methodology has been successfully implemented previously on the surfaces of germanium [57], gold [58], copper [59], and iron oxide [36] surfaces to simulate STM topographies, which showed good agreement with experimental STM images.

Structural effects of U-doping in the ThO 2 bulk
We have plotted the density of states (DOS) for each 2 Â 2 Â 2 supercell system to examine the effects of uranium doping on the electronic structure of thorium dioxide. While small changes to the geometry of the cell are visible in the relaxed structures of the Th (32Àx) U x O 64 systems under investigation (Fig. 9), changes to the electronic structure are clearer when the density of states (DOS) is examined. On the addition of two uranium atoms to the cell (Fig. 9b), the cationecation distance is affected and some of the thorium atoms have moved away from the defect region. The U 4þ cation is smaller than Th 4þ , with ionic radii of 0.100 nm and 0.105 nm, respectively [60]. The introduction of the smaller uranium cations to the system creates strain in the lattice, as the TheU distance is smaller by an average of 0.001 nm than the TheTh distance. When 1 U substitution is made, the TheTh distance is 0.3971 nm and the new TheU distance is 0.3965 nm. On the addition of a second U atom in the cell, the TheU distance decreases to 0.3957 nm and the TheTh distance decreases slightly to 0.3968 nm. A third uranium atom decreases both cationecation distances further to a TheTh distance of 0.3967 nm and the TheU distance becomes 0.3954 nm.
The DOS plot for a 2 Â 2 Â 2 supercell of ThO 2 , shown in Fig. 10, was found to underestimate the experimental band gap of 6 eV [61], with a calculated value of 4.12 eV, which, however, matches the band gap reported by other DFT studies of ThO 2 in the literature [62e64]. The conduction band is dominated by the Th 5f electrons while the valence band is predominately comprised of O 2p states. Uranium substitution affects the total DOS plot as seen in Figs. 11e13. While the band gap remains the same, there is a shift of the oxygen valence to a lower energy level from the pure ThO 2 DOS and a small peak appears at the Fermi level. This shift indicates an increase in the covalent character of the cationeanion interaction compared to the pure compound, which can be further investigated by calculating the charge density distribution for these systems (Fig. 12).

Charge density of the Th 8Àx U x O 16 bulk
We have completed a Bader charge analysis of the 1 Â 1 Â 2 supercell for the doped systems Table 1. The charges for pure ThO 2 are Th þ2.729 e À and O À1.365 e À , less than the full theoretical charges of Th þ4 e À and O À2 e À , which is common in DFT calculations. In the pure material, the valence charge is predominately located on the oxygen atoms, in agreement with the analysis of ThO 2 conducted by Wang et al. [63]. On the addition of uranium to the system, some of the charge density is now distributed to the uranium ions. As the number of U substitutions increases, there is

Thoria surface energies
The calculated surface energies before and after relaxation of the (111), (110), and (100) surfaces are presented in Table 2

Morphology of ThO 2 particles
From the relaxed surface energies, we have predicted the Wulff morphology of a thorium dioxide particle at thermodynamic equilibrium. We have used the three surfaces studied at a temperature of 0 K and the predicted particle is therefore limited to only these three surfaces. However, as the other possible surfaces of ThO 2 are expected to be higher in energy than the low miller index planes [32], we do not expect these to influence the morphology significantly. We find that the lowest energy (111) surface facets completely dominate the morphology as seen in Fig. 14, and the (110) and (100) surface orientations do not appear in the particle morphology. Our predicted morphology is in good agreement with other equilibrium morphology studies of fluorite systems, including UO 2 [50] and CeO 2 [51]. The only experimental ThO 2 morphology study of which we are aware is the 1980 study by Kawabuchi et al. of the influence of plasma currents on the growth morphology of thoria crystals [49]. They found that under the influence of the current, thoria formed rounded crystals and not faceted crystals. However, the experimental stoichiometry of the plasma experiment was ThO 1.64 and extreme growth conditions were used, both of which may have contributed to the final crystal morphology.

Uranium adsorption on ThO 2 surfaces
To calculate the adsorption energy (E ads ) of a uranium atom on the thoria surfaces, we have placed a single atom of uranium near the relaxed surface and allowed it to move in all directions to find the optimal position (Figs. 15e17). After geometry optimization, the uranium atom is located almost directly over a thorium atom in the second surface layer on the (111) and (100) surfaces, which are both oxygen-terminated. On the (110) surface, which has alternating thorium and oxygen atoms in the surface layer, the uranium has moved to a bridging position and is not located directly over any of the atoms in the top layer of the surface. The TheO, UeO, and OeO interatomic distances are listed in Table 3, before and after addition of the uranium adatom to the relaxed surfaces. The TheO distance in the (100) surface is increased in the presence of the U adatom as the oxygen atoms are drawn away from the thorium atoms, as seen by the decreased OeO distance when compared to the pristine surface. In the (110) and (111) surfaces, the interatomic distances change less from the initial relaxed surface than in the (100) surface. The energy of adsorption (Table 4) is most favorable for the (100) surface and the extent of atomic movement in this surface may be why U adsorption on this surface is highly favored. E ads is negative for every surface, meaning that the surface adsorption of a U adatom is favorable compared to the isolated atom on all three surfaces studied. The U adatom binds very strongly to the (100) and (110) surfaces with similar energies, whereas E ads for the (111) surface is less than half the value. However, even though the (111) surfaces has the smallest E ads at À1.89 eV, this value is still highly negative and uranium adsorption on this surface is highly feasible.
We have analyzed the effect of the U adatom on the charges of the ions. The results of the Bader charge analysis are presented in     surface and the classic hexagonal structure of the (111) surface can be seen in Fig. 19. The alternating Th and O ions are visible on the (110) surface (Fig. 20).

Conclusions
We have analyzed available DFT and DFTþU methodologies to determine which reproduces most accurately the electronic structure and physical characteristics of pure ThO 2 . We have subsequently used GGA with PBE functionals, without the Hubbard U correction, to calculate the density of states and charge density of thoria. However, on the addition of uranium to the thoria system, the Hubbard U correction to the exchange-correlation must be applied via the GGAþU framework. We have used a U eff of 4 eV for uranium, the value common to many other DFTþU studies of uranium. The calculated DOS plots for uranium-substituted thorium dioxide show an increasing covalent character in the system, a finding supported by the calculated Bader charges. Pure ThO 2 is ionic, with the electron density localized on the anions. As the uranium content in the system increases, some of the charge Table 3 Interatomic distances in the relaxed surfaces before and after addition of a uranium adatom.   density is shifted to the uranium cations. We have modeled three low Miller index surfaces of thorium dioxide: the (100), (110), and (111) planes. The (111) surface is found to be the lowest in energy and is therefore the most stable surface; it also undergoes the least amount of relaxation from the bulk termination as compared to the two other surfaces. The Wulff morphology of a ThO 2 particle at equilibrium, calculated from the relaxed surface energies, shows a particle completely bounded by the (111) surface orientation. From the calculated adsorption energies of a U adatom on the three relaxed surfaces, we expect U to adsorb onto all three surfaces, with particularly strong adsorption on the (100) and (100) surfaces, thus providing a route to the incorporation of U into the growing crystal.