Abstract

The mechanisms of oxygen stoichiometry variation in UO2 at different temperature and oxygen partial pressure are important for understanding the dynamics of microstructure in these crystals. However, very limited experimental studies have been performed to understand the atomic structure of UO2 near surface and defect effects of near surface on stoichiometry in which the system can exchange atoms with the external reservoir. In this study, the near (110) surface relaxation and stoichiometry in UO2 have been studied with density functional theory (DFT) calculations. On the basis of the point-defect model (PDM), a general expression for the near surface stoichiometric variation is derived by using DFT total-energy calculations and atomistic thermodynamics, in an attempt to pin down the mechanisms of oxygen exchange between the gas environment and defected UO2. By using the derived expression, it is observed that, under poor oxygen conditions, the stoichiometry of near surface is switched from hyperstoichiometric at 300 K with a depth around 3 nm to near-stoichiometric at 1000 K and hypostoichiometric at 2000 K. Furthermore, at very poor oxygen concentrations and high temperatures, our results also suggest that the bulk of the UO2 prefers to be hypostoichiometric, although the surface is near-stoichiometric.

1. Introduction

Uranium dioxide (UO2), a ubiquitous fluorite-structured ceramic oxide, is by far the most commonly used nuclear fuel for nuclear power production. Its physical, chemical, and thermodynamic properties such as electronic structure, heat capacity, and thermal conductivity, as well as behavior under extreme conditions of high temperature, pressure, and radiation dose, have been studied extensively for over half a century both experimentally and theoretically [115], especially during the last decade. However, experimental data has still remained insufficient as it is difficult to isolate unit mechanistic changes under reactor operating conditions [9]. Of particular interest is microstructural evolution and stoichiometry change in near surface UO2 under radiation [2, 8]. In particular, the mechanisms of oxygen release and uptake by UO2 crystals containing a high density of point defects are important for understanding the dynamics of defects and microstructure in these crystals. However, only limited experimental studies [7, 8, 10] have been performed to understand the atomic structure of UO2 surface and defect effects of near surface on stoichiometry, corresponding to an opened region in which the surface is contacting with external reservoirs such as oxygen gas or air, and the system can exchange atoms with the reservoirs.

Computer modeling and simulation can be complementary to experiments to explore materials properties under different conditions. However, the surface properties have not been widely studied in modeling and simulation. In this work, density functional theory (DFT) calculations and atomistic thermodynamics are used to investigate the fundamental questions surrounding the specific mechanisms that arise in off-stoichiometric UO2. The free energies of possible surface compositions and geometries can be explored by using DFT calculations, enabling the identification of the lowest-energy structure for a given condition of the thermodynamic reservoirs for the atoms and electrons. Combining DFT total-energy calculations and atomistic thermodynamics, and on the basis of the point-defect model (PDM), a general expression for the near surface stoichiometric variation is derived, in an attempt to pin down the mechanisms of oxygen exchange between the gas environment and defected UO2. The derived expression will enable us to investigate the near surface stoichiometry in UO2 as a function of temperature and oxygen partial pressure, as well as enabling the comparison with the available experimental observations. The transition or variation of stoichiometry from the surface into the bulk will also be presented in this study. This study will focus on the (110) surface, as it is charge neutral without polarity contribution from dipole or quadrupole moments.

2. Calculation Method and Methodology

2.1. First Principles Calculations

Quantum mechanics calculations based on DFT were carried out using the frozen-core all-electron projector-augmented-wave (PAW) method [17, 18] as implemented in VASP (Vienna Ab-intio Simulation Package) [19, 20]. Standard library PAW-LDA pseudopotentials for “U” and “O” were used, which give good results for surface energy and structural parameters. Periodic boundary conditions (PBC) were used in all calculations. The materials were fully relaxed until the minimized energy levels were reached. To ensure accurate results, a plane-wave cutoff energy of 500 eV was used for all calculations. Gaussian smearing was used with a smearing parameter of 0.20 eV. The Brillouin-zone integrations were performed using the Monkhorst-Pack grids and employed () for all (110) surfaces with/without defects.

The supercell of clean (110) surface contained 216 atoms with 12 layers, and each layer had 6 U atoms (by constructing along directions) and 12 O atoms, respectively. All atoms in the bottommost two layers (marked as layers 0 and 1 in Figure 1) were fixed and the thickness of vacuum was ~21 Å. The defective surfaces were generated by removing or inserting one oxygen atom from or to each of remaining 10 layers (as layer 2 to layer 11 in Figure 1). The initial position of atoms was set up to their bulk-truncated position and the volume of the supercells (including the vacuum gap) was fixed in all the surface calculations.

In general, DFT has the ability to obtain the physical properties of UO2 without any empirical parameters. However, traditional DFT calculations predict UO2 to be a metal, which is contradictory to the fact that UO2 is a Mott insulator [21, 22]. To account for the strong on-site Coulomb repulsion among the localized U electrons, a value of 3.96 eV for Ueff was used, which is aligned with the Hubbard model-based local density approximation (LDA + U) technique proposed by Dudarev and coauthors [22]. This value of effective U corresponds to  eV and  eV, which were determined in UO2 by Yamazaki and Kotani [23] based on the analysis of X-ray photoemission spectra. To avoid local energy minima for finding low-energy orbital occupations in UO2, the procedure of U-ramping method is used [24, 25].

2.2. Chemical Potential and Stoichiometry of

The energy to form a defect at a given temperature and pressure is the Gibbs energy. However, in theoretical calculations of DFT, the energies such as oxygen defect formation are typically calculated at  K. The oxygen defect formation energy for the surface in equilibrium with atomic reservoirs such as oxygen gas can be related to the chemical potential of the O2 molecule as a function of temperature and partial pressure using standard thermodynamic expressions, which is defined as where is the total energy of O2 molecule from the standard DFT calculation. The reference states of at different temperatures can be obtained from experimental specific heat [16] data at  atm and those values are listed in Table 1. Thus the oxygen defect formation energies at different temperatures and pressures near the surface can be obtained from 0 K ab initio values and integration of the chemical potential of the O2 molecule, similar to the defect formation in bulk system [21].

The stoichiometry of , assuming that contributions due to U defects (vacancies and interstitials) are negligible, can be expressed based on the point-defect model (PDM) [26] as

3. Results and Discussion

3.1. (110) Surface Properties and Relaxation
3.1.1. Pristine (110) Surface Properties

Knowledge of the charge variation and atomic arrangement near the surface region are a prerequisite to understanding of the properties of surfaces. In the framework of formal oxidation states, U and O occur as 4+ and 2−, respectively. Average atomic charge values may be obtained from population analyses from the DFT calculations, and the corresponding Bader charges in the bulk of UO2 () are 2.54 and −1.27 for U and O, respectively. The partially charge states are consistent with the fact that the UO2 is a partially ionic and partially covalent system. Furthermore, such partial charges also provide the foundation for using informal charges in the ionic potentials such as Basak potential [27] in the molecular dynamics simulations.

The difference of Bader charge () for U and O between different layer and bulk is showed in Figure 2. It is not surprising that the charges of atoms near surface deviate from their bulk values, because the atoms near the surface of a crystal have fewer neighbors than those in the bulk. As shown in Figure 2, charges are not fixed near surfaces, although the magnitude of charge deviation is not significant for the charge neutral (110) surface. It is also clearly indicated that fixed charge potential used in molecular dynamics simulations would not be able to capture this feature.

It is interesting to note that, under the topmost layer of the (110) surface, U prefers to donate more electron to the neighbor O, as shown in Figure 2. However, it is unexpected that in the fixed first and 11th layers the charge for U and O is about 2.45 and −1.23, respectively. Thus U donates less charge to O in the first layer than those in bulk and layers underneath. This result suggests that the charge accumulation/depletion close to the surface is more localized, which is attested by the pDOS analysis. In fact the surface is complemented by electronic reconstruction by a rearrangement of the -electrons of U and rehybridization of U-O at/near the surface to alleviate the surface energy, which is similar to the charge transfer mechanism on the oxidation state of transition metals in insulators [28].

According to the electronic reconstruction model, an electric field (the gradient of electrostatic potential) across the surface differing from the corresponding bulk system is required to facilitate the process of charge transfer. The surface relaxation involves the displacement of atoms perpendicular to the surface from one layer to another and hence a microscopically tailored electrostatic potential. Thus, it is highly desirable to direct mapping of the local electrostatic potential distribution. Figure 3 shows the plot of the calculated local electrostatic potential versus normalized distance (, in which is the distance to the surface and the corresponding bulk-truncated interlayer distance) of (110) surface. It is clear that the local electrostatic potential determined by the wavefunctions may have a larger weight at the surface than in the bulk as evidenced in Figure 3.

Considering the displacement perpendicular to the surface, the profile of the relaxation of the (110) surface calculated with LDA + U is also displayed in Figure 3. In this figure, normalized distances () are used to allow a fair comparison between surface and bulk. The integer number along the horizontal axis stands for the atomic layers at their bulk-truncated position. The normalized position of relaxed surface layers is determined by the minimum of the local electrostatic potential and the spacing of the bilayer near surface is the distance of the corresponding two minima. It is expected that atoms at or near clean surface generally do not maintain their bulk-truncated position. As listed in Table 2 and shown in Figure 3, the relaxation of the clean (110) surface will lead to the outward movement of the surface by about 0.241 Å. Meanwhile the intralayer spacing of topmost bilayer is reduced by ~0.019 Å and the bilayer below is increased by 0.05 Å. Subsequent interlayer spacing of bilayers also pair in a similar pattern. This pattern of alternating relaxation is also predicted by Tasker [29] using an ionic model, although the potential used is not local environmental dependent. Such oscillations are probably related to the rearrangement of the -electrons of U and rehybridization of U-O at/near the surface inducing charge transfer or redistribution rather than jellium Friedel oscillations [30, 31].

3.1.2. Defective (110) Surface Relaxation

Since the surface is the interface in contact with atomic reservoirs such as oxygen gas (O2), O2 may be uptaken/absorbed onto or released from the UO2 surface to form oxygen defects depending on the environmental temperature and pressure. These O defects, either vacancies or interstitials, near or at the surface will also affect the surface relaxation and the electrostatic potential. Figure 3 also shows the plot of the local electrostatic potential versus surface relaxation perpendicular to the surface. As expected, for O defects, either vacancies or interstitials locating near surface, the profile of the electrostatic potential and the pattern of alternating relaxation of defected surface are very similar with those of the clean (110) surface.

However, we also find that the outward movement of the (110) surface is strongly suppressed by oxygen defects of vacancies and interstitials near surface in comparison to the clean surface. Such suppression of defective surfaces is probably due to the fundamentally different kind of electronic reconstruction as mentioned above. The implication of the shrinkage as defects introduced into the surface will increase the volume of the void and may affect the thermal and mechanical properties of fuels to some extent.

3.2. Near (110) Surface Defects
3.2.1. Near (110) Surface Defects Formation Energies

In order to understand and control the fundamental solid state processes and stoichiometric variation that govern the surface of UO2, values for defect formation energetics near surface are required. Assuming that contributions are due to U defects (vacancies and interstitials), the O defect formation energies for near surface with LDA + U are displayed in Figure 4, as well as the corresponding defect formation energies in the bulk UO2. As shown in Figure 4(a), the defect formation energies for O vacancy near surface are 4 eV, larger than that in bulk (~3.78 eV). The higher O vacancy formation energy near surface indicates the O vacancies may prefer to stay away from the surface. It is interesting to note that, upon the removal of one oxygen atom, U4+ may be reduced to U3+. However, the reduced U4+ cations are adjacent to the O vacancy within same (110) plane. This result is different compared to the case of CeO2 in which the reduced Ce4+ cations are not the one adjacent to the O vacancy [32], though both materials have same fluorite structure. Meanwhile, the defect formation energies for O interstitial near surface are below −0.5 eV, lower than the corresponding one in bulk (~0.02 eV) as shown in Figure 4(b). Thus OI prefers to segregate to surface rather than staying inside of the bulk to form , which is in agreement with experimental observations [2, 33]. Note that, for the surface related calculations, spin-polarized effects are not considered, since it is extremely difficult to converge the defective surface with spin-polarized effects considered. For self-consistence and canceling systematic errors, spin-polarized effects are not considered for the corresponding defect calculations in bulk. Consequently, the formation energies for oxygen, for example, vacancy in the bulk is lower than our precious spin-polarized DFT result (~5.1 eV) [21] and recent DFT result (~4.5 eV) [34] as well. However, the relatively defective formation energies between bulk and surface will not be affected qualitatively with/without spin-polarized effects accounted.

We note that, in our LDA + U calculations, the supercells contain only 12 layers; thus the valid data are for the topmost 6 layers due to the symmetry of the slab. It is also unlikely that the surface defect formation energies will be converged to those in bulk under the current LDA + U calculations due to the finite thickness of slab. In order to map out the transition from surface to bulk and the surface influence depth, it is convenient to use an analytic equation with a quadratic dependence on the inverse of the distance as , similar to the case for the grain boundary [35]. Thus the formation energy of oxygen point defects may be written as Here , , and are fitting parameters listed in Table 3. is the distance to the surface. is the formation energy. The resulting transition curves are also shown in Figure 4 and the surface influence may reach about 100 Å from this prediction.

3.2.2. Stoichiometry from Surface to Bulk

The above DFT calculations are calculated at  K. However, combing DFT results with atomistic thermodynamics as defined in (2) and the linkage of defect formation energy from surface to bulk as (3) can be used to investigate the near surface stoichiometry and the transition of stoichiometry from the surface into the bulk. Figure 5 shows the modeling results of the near (110) surface stoichiometry in UO2 by varying the temperature. As shown in Figure 5, when UO2 is in contact with air, which has an oxygen partial pressure of 0.21 atm, the stoichiometry near (110) surface is generally hyperstoichiometric. However, the degree of hyperstoichiometry will be significantly affected by temperature. The penetration depth of off-stoichiometry of will be reduced from ~6 nm at 300 K, to ~1 nm at 1000 K, to near-stoichiometry at high temperature, that is 2000 K. We anticipate that other low index surfaces with polarity, such as (100) and (111) surfaces, will also have the similar temperature dependence of the penetration depth. In fact, such a variation in the surface chemistry of UO2 has been measured with atom probe tomography. The detailed description of experimental investigation will be presented elsewhere [33].

We further investigate how the varying oxygen partial pressure affects the near surface stoichiometry. The modeling results of the near surface stoichiometry in UO2 as the oxygen partial pressure varied are also plotted in Figure 5. Under poor oxygen conditions, such as the oxygen partial pressure being reduced to  atm, the stoichiometry of near surface is switched from hyperstoichiometric at 300 K with a depth around 3 nm to near-stoichiometric at 1000 K and hypostoichiometric at 2000 K. Furthermore, at very poor oxygen concentrations and high temperatures, our results also suggest that the bulk of the UO2 prefers to be hypostoichiometric, although the surface is near-stoichiometric.

We note that the stoichiometry in UO2 is related to the point-defect formation energies via the analytic equation (2) based on the PDM. However, an important word of caution on PDM is necessary. The basic assumptions of the PDM are that the point defects are diluted and there are no interactions between any two defects. However, in our DFT calculations each clean (110) layer only has 12 O atoms. Thus it is unlikely that the finite size effect can be screened out due to the limited size of each layer, which leads to defects interacting with its images due to the PBC applied. To reduce or avoid the various finite size effects including atomic layer and slab thickness, a much larger supercell is required. However, such requirement of a larger supercell is very difficult to carry out with the current DFT methods. Alternatively, to develop environmental dependence and charge transfer such as COMB potentials [36] for atomistic MD simulations may be an effective solution. Furthermore, due to the limited size of each layer, the stoichiometry in is more likely limited to the range of −0.167 to 0.167. The interpolation for the large deviations (>0.167 or −0.167) from the stoichiometry from (2) may just indicate the coexistence of multiple defects in one simulation layer.

It is clear that first principles calculations can be very useful for pinning down the mechanisms of oxygen exchange between the gas environment and defected UO2. But even then, there are various types of errors or uncertainties in the DFT calculations. Besides the contribution of various finite size effects to possible errors on the DFT total-energy calculations, the comparison of bulk and slab can also produce noncancelling systematic errors. Meanwhile the poor performance of density functionals in regions of rapidly varying density, that is, surfaces, is another important error that makes the self-consistent DFT calculations very difficult to converge [37], which leads to the oscillation of the data for and OI formation as shown in Figure 4. Geometry and periodicity effects are also contributing to errors.

4. Conclusions

In this paper, we studied the near (110) surface relaxation and stoichiometry in UO2. We showed that the build-up local electrostatic potential across the surface induces the charge transfer near surface and surface expansion. We also found that the outward movement of the (110) surface is strongly suppressed by oxygen defects of vacancies and interstitials near the surface in comparison to a clean surface.

Based on the combining DFT total-energy calculations and atomistic thermodynamics and on the basis of the PDM, a general and useful expression for the near surface stoichiometric variation was derived. The derived expression enables us to pin down the mechanisms of oxygen exchange between the gas environment and defected UO2. By using the derived expression, we found that when UO2 is in contact with air, the stoichiometry of near (110) surface is generally hyperstoichiometric. However, the hyperstoichiometric variation will be significantly affected by temperature. We also found that, under poor oxygen conditions, the stoichiometry of near surface is switched from hyperstoichiometric to near-stoichiometric and hypostoichiometric, depending on the temperature. Furthermore, at very poor oxygen and very high temperature, our results also suggested that the bulk of the UO2 prefers to be hypostoichiometric, although the surface is near-stoichiometric. We also found that the influence of surface on the off-stoichiometric depth is affected by the environmental temperature and oxygen partial pressure.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by the Center for Materials Science of Nuclear Fuel, an Energy Frontier Research Center (EFRC) funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award no. FWP 1356. The computations were performed using resources of the Idaho National Laboratory (INL) High Performance Computing facilities.