Near Surface Stoichiometry in UO

The mechanisms of oxygen stoichiometry variation in UO 2 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 UO 2 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 UO 2 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 UO 2 . By using the derived expression, it is observed that, under poor oxygen conditions, the stoichiometry of near surface is switched from hyperstoichiometric at 300K with a depth around 3 nm to near-stoichiometric at 1000K and hypostoichiometric at 2000K. Furthermore, at very poor oxygen concentrations and high temperatures, our results also suggest that the bulk of the UO 2 prefers to be hypostoichiometric, although the surface is near-stoichiometric.


Introduction
Uranium dioxide (UO 2 ), 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 [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15], 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 UO 2 under radiation [2,8].In particular, the mechanisms of oxygen release and uptake by UO 2 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 UO 2 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 UO 2 .The free energies of possible surface compositions and geometries can be explored by using DFT calculations, enabling the identification of the lowestenergy 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 UO 2 .The derived expression will enable us to investigate the near surface stoichiometry in UO 2 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.

Calculation Method and Methodology
2.1.First Principles Calculations.Quantum mechanics calculations based on DFT were carried out using the frozencore 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 (2 × 2 × 1) 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 2 × 3 along [100] × [011] 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 UO 2 without any empirical parameters.However, traditional DFT calculations predict UO 2 to be a metal, which is contradictory to the fact that UO 2 is a Mott insulator [21,22].To account for the strong on-site Coulomb repulsion among the localized U 5 electrons, a value of 3.96 eV for U eff 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 U = 4.50 eV and  = 0.54 eV, which were determined in UO 2 by Yamazaki and Kotani [23] based on the analysis of X-ray photoemission spectra.To avoid local energy minima for finding low-energy 5 orbital occupations in UO 2 , the procedure of U-ramping method is used [24,25].

Chemical Potential and Stoichiometry
of UO 2+ .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  = 0 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 O 2 molecule as a function of temperature  and partial pressure  using standard thermodynamic expressions, which is defined as Table 1: The chemical potential of O as a function of the temperature.The entropy and enthalpy changes are taken from the JANAF thermochemical tables at  0 = 0.1 atm [16].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 O 2 molecule, similar to the defect formation in bulk system [21].
The stoichiometry  of UO 2+ , assuming that contributions due to U defects (vacancies and interstitials) are negligible, can be expressed based on the point-defect model (PDM) [26] as ] .

Results and Discussion
The difference of Bader charge (Δ =  −  bulk ) 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 (/ 0 , in which  is the distance to the surface and  0 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 (/ 0 ) 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].

Defective (110) Surface Relaxation.
Since the surface is the interface in contact with atomic reservoirs such as oxygen gas (O 2 ), O 2 may be uptaken/absorbed onto or released from the UO 2 surface to form oxygen defects depending on 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.

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 UO 2 , 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 UO 2 .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, U 4+ may be reduced to U 3+ .However, the reduced U 4+ cations are adjacent to the O vacancy within same (110) plane.This result is different compared to the case of CeO 2 in which the reduced Ce 4+ 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 O I prefers to segregate to surface rather than staying inside of the bulk to form UO 2+ , 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 spinpolarized 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 Table 3: Fitting parameters for defect formation energy as a function of distance to the surface.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  2 , similar to the case for the grain boundary [35].Thus the formation energy of oxygen point defects may be written as Here  0 ,  1 , and  0 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.

Stoichiometry
from Surface to Bulk.The above DFT calculations are calculated at  = 0 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 UO 2 by varying the temperature.As shown in Figure 5, when UO 2 is in contact with air, which has an oxygen partial pressure ( O 2 ) 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 UO 2+ 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 UO 2 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 UO 2 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  O 2 = 10 −5 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 UO 2 prefers to be hypostoichiometric, although the surface is near-stoichiometric.
We note that the stoichiometry  in UO 2 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 UO 2+ 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 UO 2 .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 totalenergy 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  O and O I formation as shown in Figure 4. Geometry and periodicity effects are also contributing to errors.

Conclusions
In this paper, we studied the near (110) surface relaxation and stoichiometry in UO 2 .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 UO 2 .By using the derived expression, we found that when UO 2 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 nearstoichiometric and hypostoichiometric, depending on the temperature.Furthermore, at very poor oxygen and very high temperature, our results also suggested that the bulk of the UO 2 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.

Figure 1 :
Figure 1: The slab structure of (a) pristine in the direction perpendicular to the (110) surface, in which each layer is the stoichiometric UO 2 ; (b) O vacancy on the defective (110) plane corresponding to UO 1.833 ; and (c) O interstitial on the defective (110) plane corresponding to UO 2.167 .O atoms are shown in red color, U atoms in grey, and defects in green.The bottommost fixed bilayer is marked as layer of 0 and 1.

Figure 2 :
Figure 2: The average difference of Bader charge (Δ =  −  bulk ) for U and O between different layer and bulk ( bulk = 2.54 and −1.27 for U and O, resp.).

Figure 3 :
Figure 3: Calculated local electrostatic potential versus normalized distance of (110) surface.The integer number along the horizontal axis is the atomic layer being the bulk-truncated position and 11 is the topmost surface layer.We used -averaged local electrostatic potential.The relative position of relaxed surface layers is determined by the minimum of the local electrostatic potential.E-Fermi, O I -Fermi and  O -Fermi are the Fermi levels of the corresponding surface (bulk-truncated, O interstitial and O vacancy) calculations.

Figure 4 :
Figure 4: O defect formation energies for near surface and in the bulk UO 2 with LDA, as well as the fitting curve (a) O vacancy and (b) O interstitial, corresponding defect formation energies in the bulk UO 2 .
O 2 is the total energy of O 2 molecule from the standard DFT calculation.The reference states of  O 2 (,  0 ) =    −  at different temperatures can be obtained from experimental specific heat[16]data at  0 = 0.1 atm and those values are listed in Table where

Table 2 :
Staring and final average layer position of (110) surface.