The Conundrum of Relaxation Volumes in First-Principles Calculations of Charged Defects in UO 2

: The defect relaxation volumes obtained from density-functional theory (DFT) calculations of charged vacancies and interstitials are much larger than their neutral counterparts, seemingly unphysically large. We focus on UO 2 as our primary material of interest, but also consider Si and GaAs to reveal the generality of our results. In this work, we investigate the possible reasons for this and revisit the methods that address the calculation of charged defects in periodic DFT. We probe the dependence of the proposed energy corrections to charged defect formation energies on relaxation volumes and ﬁnd that corrections such as potential alignment remain ambiguous with regards to its contribution to the charged defect relaxation volume. We also investigate the volume for the net neutral defect reactions comprising individual charged defects, and ﬁnd that the aggregate formation volumes have reasonable magnitudes. This work highlights the issue that, as is well-known for defect formation energies, the defect formation volumes depend on the choice of reservoir. We show that considering the change in volume of the electron reservoir in the formation reaction of the charged defects, analogous to how volumes of atoms are accounted for in defect formation volumes, can renormalize the formation volumes of charged defects such that they are comparable to neutral defects. This approach enables the description of the elastic properties of isolated charged defects within an overall neutral material.


Introduction
Charged defects in non-metals, such as semiconductors and insulators, play a central role in determining electronic and electrical behavior of materials for semiconductor and battery applications [1]. Point defects are also very relevant to actinide materials such as UO 2 under extreme conditions of temperature, pressure, and irradiation, as they play a major role in accommodating stoichiometry variations [2]. In addition, defects induce lattice distortions in the perfect crystal, leading to changes in volume which can, at least in principle, be measured experimentally [3,4] and compared to theoretical calculations [4][5][6][7][8][9]. These lattice distortions are particularly important in determining the interaction of defects with external strains, dislocations, and grain boundaries [10,11], and play a central role in the evolution of microstructures in materials.
The defect formation volume is an intrinsic property of the defect; theory for computing the volume change due to point defects gained attention in the early 1980s with a focus on point defects in ionic crystals [5,6]. It is important to note that experimental measurements of the volume change associated with point defects yield the defect-averaged response of a mass and charge balanced system [3][4][5][6]; as a result there is little understanding of the volume change or lattice distortion associated with individual charged defects. Investigations of the formation volume of charged point defects in UO 2 [9], Si [12][13][14][15], ZnO [16], In 2 O 3 [17], CeO 2 [18], SrTiO 3 [4,19] indicate that the defect formation volume is a strong function of the charge state of the defect, and as electrons are added (or removed) to the system the formation volume increases (or decreases) linearly with charge.
Defect formation and relaxation volumes are closely related concepts. The defect relaxation volume is simply the change in volume of a computational cell due to relaxation of the crystal lattice around the defect site. By contrast, the defect formation volume is calculated as the change in the volume associated with the formation reaction for the defect; it, therefore, accounts for the reference volume associated with the reservoir to or from which atom(s) are moved in order to create the defect in the simulation. We describe defect relaxation and formation volume and the differences between the two in more detail in Section 2.2.
The theory of relaxation volume and lattice relaxation around the defect site can be used to obtain the so-called defect elastic dipole tensor [20,21], which relates the atomic structure of the point defect with its elastic distortion of the lattice. The elastic dipole tensor [7,20,22], G ij , captures the symmetry of the point defect and is given as: where C ijkl is the elastic constant tensor, and ∆Ω r kl is the relaxation volume tensor of the defect [23]. The elastic dipole tensor describes the change in formation energy of the defect, to first order, in the presence of external strains [10]. Therefore, computation of the defect relaxation volume tensor quantifies the lattice distortion around the defect and thus the elastic interaction with other defects and strain fields [11,24].
An accurate description of the dipole tensor is critical for developing models of transport in complex microstructures. Real materials contain dislocations, grain boundaries, and other defect structures that can be described by elastic strain fields. The dipole tensor provides a means for relating changes in energy of a given point defect as its position varies in that strain field. Thus, the dipole tensor provides a link between atomistic calculations of defects in an ideal bulk environment and in the more realistic microstructure of the material during application, allowing one to up-scale results to more applied studies, such as irradiation-induced creep which relies on the treatment of defects interacting with the existing microstructure. To establish this link, one must be able to calculate quantitatively accurate values for the dipole tensor or, related to this, the relaxation volume of defects.
Building on our previous work on point defect behavior in the presence of external strain [11], here we find the relaxation volumes of individual charged point defects in uranium dioxide (UO 2 ) to be much larger than those of the corresponding neutral defects. Furthermore, when we calculate the corresponding elastic dipole tensor, the predicted changes in energy of the defect under elastic strain fields are extremely large; indeed, they are so large that they seem to be unphysical. As the dipole tensor is directly related to defect relaxation volumes, this result motivates us to examine this issue in greater depth. The issue of large DFT relaxation volumes of charged defects as well as discrepancies between theory and experiments on the topic has been recognized in the past [4,9,14,15,25]. However, we argue that the issue is not completely resolved. We suggest that by considering the change in volume of the electron reservoir in the formation reaction of the charged defects, we can renormalize the formation volumes of charged defects, and this approach represents one possible route for resolving the issue, in addition to using a net neutral defect reaction model to determine relaxation volumes.
In this work we focus on the issue of the relaxation volume of charged defects, focusing on the actinide material UO 2 , both because of its complex electronic structure and because its fluorite crystal structure is common to many technologically important materials. Defect calculations [26][27][28][29] in UO 2 have established that defects such as Frenkel pairs and Schottky trios and fission product complexes are likely to consist of individual components that are charged. The individual volume change depends not only on the size mismatch of the defect species, but also on the more complex charge transfer mechanisms occurring in the presence of the defect. In addition to UO 2 , we present results and comparisons for relaxation volumes of charged defects in Si and GaAs, for the purpose of validating against existing data in these systems and uncovering behavior that applies to charged systems more generally.
The paper is arranged as follows. Section 2.1 presents the density-functional theory (DFT) methodology used in calculating defect properties. In Section 2.2, we discuss the definitions of defect formation and relaxation volumes. Various schemes to correct the total energy of charged supercells simulated using DFT and their relation to defect relaxation volumes are investigated in Section 3.1. Our results on the details of the lattice distortion around individual charged and neutral defects are presented in Section 3.2. In Section 3.3, we discuss the defect reaction approach for calculating the defect relaxation and formation volumes for complete defect formation reactions. In Section 4, we discuss the significance of our results and present our conclusions in Section 5.

Density Functional Theory
Here, we summarize the details of the first-principles based DFT calculations used in this study. The DFT calculations are performed with the Vienna Ab Initio Simulation Package (VASP) [30,31] using the projector-augmented-wave (PAW) method [32,33]. For UO 2 the electron exchange and correlation potential is described by the local density approximation (LDA) [34] and the Hubbard +U model is used for the strongly-correlated U 5f electrons [35]. To be in accordance with our earlier studies [11,36], we employ the LDA functional for UO 2 , and the U and J values are set to U = 4.5 eV and J = 0.51 eV; i.e., U − J = 3.99 eV. We adopt the Perdew-Burke-Ernzerhof generalized gradient approximation (PBE-GGA) functional [37] for silicon and LDA for GaAs to be consistent with reference literature [38,39]. The energy cut-offs employed are 500 eV for UO 2 , 520 eV for Si and 275 eV for GaAs. Defect properties in UO 2 are calculated using 2 × 2 × 2 and 2 × 2 × 3 supercells with 2 × 2 × 2 and 2 × 2 × 1 k-point meshes, respectively. Defect calculations are performed at both constant volume and fully relaxed (both in cell volume, shape, and atomic coordinates) conditions. Defects in silicon are simulated in a 2 × 2 × 2 and 3 × 3 × 3 size supercells with k-point meshes of 4 × 4 × 4 and 3 × 3 × 3, respectively, while those in GaAs are simulated in a 3 × 3 × 3 supercell with a k-point mesh of 2 × 2 × 2. The dielectric tensors are important for the correction schemes for the defect formation energies of charged systems and are calculated using density functional perturbation theory [40,41] as implemented in VASP. The computed bulk properties of interest for UO 2 , Si and GaAs are summarized in Table 1.
We perform convergence tests of the final total energy with respect to the energy cut-off and k-point mesh sizes. The convergence tolerance for forces on each ion is 0.01 eV/Å, resulting in total energies converged to less than 10 −5 eV. To ensure there are no issues associated with the change in plane wave basis and k-point density due to changing cell size, we also calculate the defect relaxation volume of point defects for higher energy cut-off and denser k-point mesh; we find that the relaxation volumes of neutral and charged defects are converged to 0.3-1.8 Å 3 (which is within 5% of the relaxation volume) depending on the type of defect. Energy vs. volume calculations of defects are performed starting from the converged wave functions for the minimized defect structure. This is done to avoid the system becoming stuck in the metastable energy state associated with the possible configurations of the localized uranium f-electrons [42].
Finally, to verify that our results are consistent across DFT codes, we performed calculations of the relaxation volume for defects in GaAs and Si with CASTEP [43], again using GGA-PBE with "on the fly" Table 1. Calculated lattice parameters, a 0 in Å, elastic constants, C ij , and bulk modulus, B, in GPa and the static (electronic + ionic contribution) dielectric constant, ε, from density-functional theory (DFT) calculations in UO 2 (local spin density approximation; LSDA + U), Si (GGA) and GaAs (local density approximation, LDA). Our results are indicated in italics while literature values are provided for comparison.

Defect Formation and Relaxation Volumes
The formation volume of a defect ∆Ω f is given by the change in Gibbs free energy, ∆G f (p, T), for the defect formation reaction at constant pressure and temperature via [5]: where ∆E f denotes the change in the internal energy, ∆Ω f the formation volume, and ∆S the change in entropy due to the formation of the defect at constant pressure and temperature. The defect formation volume is the difference in volume across the defect formation reaction and is related to the Gibbs free energy of defect formation by: At T = 0 K, the change in Gibbs free energy, ∆G f (p, 0) reduces to the change in enthalpy, ∆H f . However, at higher temperatures, the entropy contribution can also become relevant [52]. The defect formation enthalpy, ∆H f , from DFT calculations for a defect X, in charge state q is given as [52,53]: where, E DFT tot (X, q) is the total energy of the DFT supercell with the charged defect, E DFT tot (bulk) is the total energy of the defect-free bulk supercell, and n i and µ i are the number and chemical potential of each element contributing to the defect, ε v is the energy of the valence band maximum (VBM) calculated in the defect free bulk supercell [49], and ε F is the Fermi energy with respect to the VBM. The term q(ε v + ε F ) is the chemical potential of electrons, i.e., it accounts for exchange of electrons with the reservoir. The ∆Ω f contribution to the enthalpy in Equation (4), is usually very small at ambient conditions (~10 −5 eV per defect at p = 1 atm) and often neglected. The terms E Image corr and ∆V are the charged image interaction and potential alignment corrections, respectively. These corrections are needed to correct for finite-size artifacts within the DFT supercell approach to calculate the defect formation energies. The expression for computing formation volume from DFT calculations is often given as [23,54]: Here, ∆Ω r is the relaxation volume, and Ω(defect) and Ω(bulk) are the volumes of the supercell with a defect X in charge state q, and the perfect bulk supercell, respectively, both relaxed to zero stress. The + and − signs in Equation (5) denote vacancy and interstitial defects, respectively. The defect relaxation volume can be understood as the change in volume due to relaxation of the crystal lattice around the defect site, or the direct volume change in the system upon introducing the defect in a constant pressure calculation. By contrast, the formation volume also accounts for the reference volume (or the choice of reservoir) associated with the origin or placement of the species comprising the defect, Ω 0 . The relaxation and formation volume are actually tensor quantities, e.g., ∆Ω r ij is given as [23]: where ε ij is the strain tensor (strain between supercells computed from DFT calculations) due to the defect. This definition holds for any anisotropy, shape or size of the supercell. The magnitude of the volume change is then given as, ∆Ω r = 1/3 × tr(∆Ω r ij ), and is consistent with that computed using Equation (6). This description of the relaxation volume accounts for changes in cell shape upon introduction of the defect. It is important to note that the relaxation volume is a finite converged value even for very large simulation cell size. Finally, the elastic dipole tensor G ij is related to ∆Ω r ij (Equation (1)) and thus, in order to describe the response of charged defects to strain, it is imperative to be able to accurately determine the relaxation volumes.

Volume of Individual Species Comprising the Defect
As described in Equation (5), computation of the defect formation volume needs to take into account the reference volume, Ω 0 , of the individual atomic species constituting the defect reaction. For a single-component system like silicon, Ω 0 is equal to the cohesive volume of an atom in bulk silicon. However, for a multicomponent system, Ω 0 is given by the partial molar volume, which measures the variation of volume as chemical components are added or removed from the reservoir during the defect reaction. The choice of reservoir for partial molar volume is exactly analogous to the atomic chemical potential in the definition of defect formation energy. In a multicomponent system, defect reactions can be of two types: (1) Stoichiometric reactions comprising of defects like the Schottky trio. In this case, Ω 0 is the volume per formula unit of the multicomponent system or, analogous to the monoatomic case, the cohesive volume of a formula unit. For example, in UO 2 this is computed to be Ω UO 2 = 40.3 Å 3 (40.9 Å 3 in experiment [55]) and in GaAs it is Ω GaAs = 44.5 Å 3 ; (2) Non-stoichiometric reactions comprised of defects such as vacancies and interstitials. Here Ω 0 is the volume per unit of the reference.  [59]. From this, we derive the reference volume for O atom to be 22.6 Å 3 from Ω UO 2 and Ω U 3 O 8 as discussed in detail in Section 3.3. In all of these cases, the reference explicitly depends on the experimental conditions and is not something directly computable from electronic structure methods. That is, once a reference is chosen, the relevant thermodynamics properties can be calculated, but the reference itself is not something that can be determined solely from the calculations; it must be asserted or assumed.
We use Equations (5) and (6) to calculate the defect formation and relaxation volumes. In compounds this calculation is complicated by the fact that the stoichiometry of the system changes when individual point defects such as vacancies or interstitials are formed. For charged defects, the issue is further exacerbated as the electrons or holes added to the defect system also contribute to both the defect formation energy and volume. While the contribution of net charge to the defect formation energy has been well studied in the past and several correction terms have been developed, the contribution of the net charge to the defect formation volume has only recently been recognized [14,15]. Indeed, in those seminal works, the authors have shown that defect formation volumes are large, and propose various approaches for considering artificial effects of net charge on formation volume. However, there are still important and unanswered questions about the magnitude of the predicted formation volumes of charged defects as calculated within DFT, and the physical reasonableness of those volumes. In this work, motivated by the work of those authors, we further their efforts by detailed examination of how approaches proposed for correcting energy impact volumes. We conclude that to appropriately compute defect relaxation volumes requires accounting for the volume associated with the charge reservoir just as is done for the energies.
In the following discussion, we explicitly focus on the defect relaxation rather than the formation volume as the latter can be obtained by properly accounting for Ω 0 . Here we show that the contribution of the change in charge to defect relaxation volumes is surprisingly large. By studying point defects in three different systems, we estimate for a given system the contribution of the charge to the defect relaxation volume.

Results
The computed values of the relaxation volumes of point defects in both neutral and fully charged states are listed in Table 2. The lattice relaxation and its symmetry around a defect are quantified by the relaxation volume tensor, ∆Ω r ij (Equation (7)), which in turn determines the defect dipole tensor, G ij (Equation (1)), listed in Table 3. The relaxation volumes are relatively large for charged defects compared to their neutral counterpart and, interestingly, the sign of the relaxation volume is sensitive to the charge state of the defect. The relaxation volumes are computed using Equation (6) (when relaxing cell shape, volume and ionic positions), and Equation (7) (when keeping volume and cell shape fixed, relaxing only ionic positions). The agreement in computed relaxation volumes calculated through the two approaches is between 0.5-2 Å 3 ; the differences in total energies is 50-100 meV, a similar range to that reported by Freysoldt et al. [52].
As discussed above, the dipole tensor is used to compute the change in defect formation energy upon interaction with strains in the material [11,24]. We find that the changes in the formation energy differ significantly between the charged and neutral states for a given defect. For example, for a fully charged (2−) oxygen interstitial under 1% compressive strain the change in formation energy is about 277.4 meV, which is 6 times larger than the −46.4 meV change in energy for a neutral oxygen interstitial under the same strain state. Similarly, the segregation energy of a fully charged oxygen interstitial in UO 2 under the strain field of edge dislocation [60] is computed to be~−8.3 eV, much larger than thẽ −1.2 eV segregation energy of the neutral oxygen interstitial. The differences in formation energies [29] and the analyzed local lattice distortions (see Section 3.2) between the neutral and charged defects are, however, not so large. Therefore, this extremely large change in segregation energies (which depend on the value of the relaxation volumes) between a charged and a neutral defect interacting with the same strain field seems unphysical and suggests that something is amiss with the volumes calculated for the charged defects.  As discussed above, the dipole tensor is used to compute the change in defect formation energy upon interaction with strains in the material [11,24]. We find that the changes in the formation energy differ significantly between the charged and neutral states for a given defect. For example, for a fully charged (2−) oxygen interstitial under 1% compressive strain the change in formation energy is about 277.4 meV, which is 6 times larger than the −46.4 meV change in energy for a neutral oxygen interstitial under the same strain state. Similarly, the segregation energy of a fully charged oxygen interstitial in UO2 under the strain field of edge dislocation [60] is computed to be ~−8.3 eV, much larger than the ~−1.2 eV segregation energy of the neutral oxygen interstitial. The differences in formation energies [29] and the analyzed local lattice distortions (see Section 3.2) between the neutral and charged defects are, however, not so large. Therefore, this extremely large change in segregation energies (which depend on the value of the relaxation volumes) between a charged and a neutral defect interacting with the same strain field seems unphysical and suggests that something is amiss with the volumes calculated for the charged defects.  Figure 1 shows the dependence of the defect relaxation volume on the simulation cell size for charged and neutral vacancies. These results indicate that the large charged defect relaxation volumes are not a consequence of the finite system size limitation of DFT. In the next section, we investigate in detail the dependence of relaxation volumes on the various finite size charge corrections schemes to charge supercell calculations, which have been suspected [14] to have large volume dependence.  Figure 1 shows the dependence of the defect relaxation volume on the simulation cell size for charged and neutral vacancies. These results indicate that the large charged defect relaxation volumes are not a consequence of the finite system size limitation of DFT. In the next section, we investigate in detail the dependence of relaxation volumes on the various finite size charge corrections schemes to charge supercell calculations, which have been suspected [14] to have large volume dependence.

Finite Size Correction Schemes to Energy for Charged Supercell
As discussed in the Section 2.2, the defect formation volume follows thermodynamically from the defect formation energy (Equation (3)). Therefore, it is important to investigate any possible connection between the corrections to the defect formation energy of charged defects and their defect formation volumes.
Charged supercell calculations in DFT require an artificial background charge density (opposite in sign to the net charge of the supercell) to ensure global charge neutrality and convergence of the total energy of the charged system within the supercell approach. The genesis of this background charge is related to the DFT convention of setting the average electrostatic potential over the entire volume of the supercell to zero in order to evaluate the electrostatic Coulomb energy [61]. We note that the same convention is typically applied to calculations of charged defects using ionic potentials. The total energy of charged supercells must account for the chemical potential of the electrons, q(ε v + ε F ) in Equation (4), and needs to be corrected for finite size errors associated with the periodic supercell approach [39,52,61,62]. There are two specific effects. The first is the interaction between a charged defect and its periodically repeated charged images, adding a spurious electrostatic energy; this is referred to as the image charge correction, E Image corr . The second is the electrostatic potential between the defect and bulk system that needs to be aligned in order to accurately define the reference energy (chemical potential) of electrons, added or removed in charged supercell calculations. This is termed as the potential alignment correction, q∆V. We adopt the method by Freysoldt et al. [39] to correct for finite size errors as it has been shown to achieve good convergence [49] for defect formation energies with system size.
In the next subsections, we systematically investigate the relation between these finite-size energy corrections and the charged defect relaxation volumes. Any change in total energy that depends on the volume of the simulation cell affects the computed forces, pressure/stress tensors, and relaxation volumes in DFT simulations. Therefore, both volume and stresses need to be corrected in the same manner as the total energy.

Image Charge Correction
The image charge correction is needed because, in periodic supercell calculations, a charged defect (in the presence of the homogeneous background charge density) has long-range electrostatic interactions with its periodic images. If the dimensions of the unit cell were infinite, this electrostatic interaction would be zero. However, in periodic DFT calculations the supercell dimensions are finite; therefore we need to correct the total energy of the charged system to account for this spurious Coulomb interaction. This correction is referred to as the monopole correction by Leslie and Gillan [63], and is designed to transform the electrostatic energy of a periodic lattice of point charges in a neutralizing background into the electrostatic energy of single point charge. As such, it is not unique to DFT, but occurs whenever the energetics of a system with net charge is considered. The correction term is given as [39,63]: where α is the Madelung constant, Ω is the volume of the bulk supercell in which the defect is created, and ε represents the electronic and ionic low frequency dielectric constant of the material, computed from DFT calculations (Table 1). In our calculations, we use the average value [38,64], 1 3 tr ij of the dielectric tensor as the effective value for dielectric constant in Equation (8).
The energy correction given by Equation (8) is computed directly from VASP, and the corresponding correction to the pressure ∆P image corr can either be calculated directly from the slope of the image-charge correction (∆E Image corr ) as a function of volume or via the relation as derived in Reference [13]: with the corresponding volume correction given as, The computed image charge corrections to energy, pressure and volume are summarized in Table 4 for charged defects in UO 2 , GaAs and Si. The corrections to volume are very small compared to the corrections in the total energy and to the corresponding defect relaxation volumes: ∆Ω image corr~0 .01-0.02 ∆Ω r , for the three systems studied here. Thus, the correction to the relaxation volumes due to the image charges is only a 1%-2% effect.  16.5 Our computed value of ∆Ω image corr for the charged (2+) Si vacancy is 0.02 Ω Si , which is similar to that reported by Bruneval and co-workers [14]. We note that the correction to the relaxation volume and pressure (or stresses) decreases with the increasing size of the simulation supercell and is independent of the sign of the charge on the defect, depending only on the magnitude of the charge. For example, for the U vacancy, the energy and volume corrections are 0.82 eV, and 0.2 Å 3 in a 3 × 3 × 3 supercell, compared to 1.23 eV, and 0.32 Å 3 in 2 × 2 × 2 supercell.
The above changes in the pressure and volume are computed based on two assumptions. The first assumption is that the value of the bulk modulus, B, does not change significantly due to the presence of defects. The second assumption is that the charge is localized at the defect site within the simulated supercell size. Our computed value of the bulk modulus of UO 2 with a U vacancy (in a 144 atom supercell) is 203.3 GPa compared to the bulk value of 220.5 GPa. This difference gives an error of 0.01-0.03 Å 3 in the estimation of ∆Ω image corr . We also computed the difference in planar averaged charge density of the charged defect (averaged over planes of atoms) from its bulk structure, as shown in Figure 2a, for the charged uranium vacancy in UO 2 . The charge density is fairly localized at the defect site and saturates to the screened [65] charge value of 1 − 1 q Ω , as indicated by the dashed red line in Figure 2a.

Alignment Correction
Among the various procedures [39,52,62,66] outlined in the literature to align the average electrostatic potential in the defect supercell with that in the bulk supercell, we use the scheme provided by Freysoldt et al. [39]. The key advantage of this scheme (referred to in the literature as the FNV scheme [49]) is that the long-range 1/r potential is removed before the alignment is determined for the remaining short-range potential. If the size of the simulation cell is big enough to fully localize the defect charge density within a sub-region of the supercell, the short range potential reaches a plateau (Figure 2b) far from the defect, which yields the alignment as a well-defined quantity [52].
first assumption is that the value of the bulk modulus, B, does not change significantly due to the presence of defects. The second assumption is that the charge is localized at the defect site within the simulated supercell size. Our computed value of the bulk modulus of UO2 with a U vacancy (in a 144 atom supercell) is 203.3 GPa compared to the bulk value of 220.5 GPa. This difference gives an error of 0.01-0.03 Å 3 in the estimation of ∆ . We also computed the difference in planar averaged charge density of the charged defect (averaged over planes of atoms) from its bulk structure, as shown in Figure   is shown in green. The z-axis distance is measured from the location of the defect. The potential falls off again at a distance of~16 Å as that is the size of the cell and where the periodic image of the defect is located. Results are shown for a 144 atom UO 2 supercell with volume of 1936.76 Å 3 .
The short-range potential away from the defect site is obtained by taking the difference between the electrostatic potential from a model charge density V model and the DFT electrostatic potential from the charged defect V DFT q and the perfect bulk V DFT b : (11) The model charge density is the charge distribution used to capture the localized defect charge (as shown in Figure 2a). In the FNV alignment correction scheme, we use a Gaussian charge distribution for the model charge density, with width equal to 1 Bohr. The DFT electrostatic potential (blue curve), and the potential due to model charge density (red curve) are shown in Figure 2b. The difference is the short-range electrostatic potential (green curve) which reaches a plateau away from the defect site (centered at origin), thereby giving the alignment correction ∆V (highlighted as a dashed line in the plots). The computed value of ∆V for V 4− U is −0.17 eV, and it is computed in a similar manner for all other charged defect studied. For defects with large atomic relaxations the charge density and potential can show strong variation close to the defect site, resulting in uncertainty in extracting the exact value of the alignment correction. Recent reviews [38,49] of charged defect correction schemes highlight similar issues with relaxed geometries. While this source of error in potential alignment correction needs to be quantified, the results from Reference [46] suggest that the changes are relatively small, i.e., within 40-100 meV, and unlikely to make a significant change to the defect's volume dependence.
To obtain the corrections to pressure and volume due to potential alignment correction, we compute the defect formation energy with and without the alignment correction for a set of isotropically-scaled volumes. Figure 3 shows the defect formation energy as function of volume; from the fit to the energy-volume data we compute the volume corresponding to the energy minima, calculated with two approaches. In approach (1), the potential alignment is computed between defect and bulk, with volume change (isotropic-scaling) applied to both the defect and the bulk supercell, such that the defect and bulk have the same volume. In approach (2), the potential alignment is computed such that volume change is applied only to the defect supercell and bulk is kept fixed to its pristine volume. For approach (2) we have employed the LZ scheme [62] because the FNV scheme cannot be applied when the defect and bulk have different volume. That is, the difference in the two approaches is the volume at which the bulk is kept when computing the potential alignment correction of the defective cell.  Comparing the volumes corresponding to the minimum defect formation energies with and without the alignment correction gives an estimate of the correction to the relaxation volume due to the alignment correction. Based on approach (1), top panels in Figure 3, we find that this correction is small. For in a 144 atom supercell, the volume correction is 1.09 Å 3 (0.03 ΩUO2), while that for in a 216 atom supercell it is about 0.81 Å 3 (0.02 ΩGaAs). However, based on approach (2), bottom panels in Figure 3, the correction is very significant; for the volume correction is −40.33 Å 3 (−1.001 ΩUO2), while that for it is −51.2 Å 3 (−1.2 ΩGaAs). This gives final relaxation volumes of to be 0.05 Å 3 and that of to be −19.9 Å 3 . Approach (1) computes the alignment correction by aligning the common reference (average electrostatic potential) in the defect and bulk supercells, as volumes are kept same, as in References [36,59]. Based on the analytical expression derived using FNV scheme by Komsa et al. [49], the potential alignment follows a 1/Ω (where Ω is supercell volume) dependence, and converges much faster (Ω −2 ) than the image charge correction (Ω −4/3 ) with increasing cell size. This suggests that the potential alignment has a weak volume dependence, which falls in line with the results obtained using approach (1), as shown in Figure 3a.
The large volume change as calculated from approach (2) can be broken down into two independent components. As shown in Figure 4, the first component (very small) is actually the same as approach (1), while the second (much larger) component aligns strained bulk (in the same charge state as the defect) with the pristine bulk. Clearly, the volume change of the second component is as  Comparing the volumes corresponding to the minimum defect formation energies with and without the alignment correction gives an estimate of the correction to the relaxation volume due to the alignment correction. Based on approach (1), top panels in Figure 3, we find that this correction is small. For V 4− U in a 144 atom supercell, the volume correction is 1.09 Å 3 (0.03 Ω UO2 ), while that for V 3− Ga in a 216 atom supercell it is about 0.81 Å 3 (0.02 Ω GaAs ). However, based on approach (2), bottom panels in Figure 3, the correction is very significant; for V 4− U the volume correction is −40.33 2 Ω GaAs ). This gives final relaxation volumes of V 4− U to be 0.05 Å 3 and that of V 3− Ga to be −19.9 Å 3 . Approach (1) computes the alignment correction by aligning the common reference (average electrostatic potential) in the defect and bulk supercells, as volumes are kept same, as in References [36,59]. Based on the analytical expression derived using FNV scheme by Komsa et al. [49], the potential alignment follows a 1/Ω (where Ω is supercell volume) dependence, and converges much faster (Ω −2 ) than the image charge correction (Ω −4/3 ) with increasing cell size. This suggests that the potential alignment has a weak volume dependence, which falls in line with the results obtained using approach (1), as shown in Figure 3a.
The large volume change as calculated from approach (2) can be broken down into two independent components. As shown in Figure 4, the first component (very small) is actually the same as approach (1), while the second (much larger) component aligns strained bulk (in the same charge state as the defect) with the pristine bulk. Clearly, the volume change of the second component is as large as that of the strained defect with respect to the pristine bulk, as shown for both charged uranium and gallium vacancies. This suggests that one probably does not need to introduce a defect to calculate the volume dependence of the potential alignment correction, instead it can, to a good approximation, be computed by just using the strained bulk (in same charge state as the defect). The correction to the volume ∆Ω PA corr , based on just the dominant bulk contribution per unit charge (the last term in Figure 4), for the systems studied are calculated as −9.5 Å 3 , −17.0 Å 3 and −19.9 Å 3 for UO 2 , GaAs and Si, respectively. Results from approach (2) (Figure 3, bottom panels) agree with the conclusions from Reference [13], which reported that the alignment correction gives a large contribution to the charged defects relaxation volumes.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 12 of 24 large as that of the strained defect with respect to the pristine bulk, as shown for both charged uranium and gallium vacancies. This suggests that one probably does not need to introduce a defect to calculate the volume dependence of the potential alignment correction, instead it can, to a good approximation, be computed by just using the strained bulk (in same charge state as the defect). The correction to the volume ∆ , based on just the dominant bulk contribution per unit charge (the last term in Figure 4), for the systems studied are calculated as −9.5 Å 3 , −17.0 Å 3 and −19.9 Å 3 for UO2, GaAs and Si, respectively. Results from approach (2) (Figure 3, bottom panels) agree with the conclusions from Reference [13], which reported that the alignment correction gives a large contribution to the charged defects relaxation volumes. Our results and the defect literature thus present a difficulty in identifying the correct approach and resolving the ambiguity about the volume dependence determined from potential alignment correction. Schemes addressing the potential alignment correction assume localized defects and aim at determining the change in average electrostatic potential, perturbed by the presence of a defect, by keeping the volume in bulk and defect calculations the same. In Figure 4, we show that the origins of the large volume change due to potential alignment is related to the charged bulk and different references (average electrostatic potential) because of the comparison between the strained bulk and the pristine bulk. That is, the large correction in the volume induced by introducing the charged defect is not due to the atomic defect (or very little), but is almost entirely due to the effect of the charge in the system. Here, we have treated potential alignment using the conventional schemes, based on References [36,59] and at present we cannot determine which of the two approaches in Our results and the defect literature thus present a difficulty in identifying the correct approach and resolving the ambiguity about the volume dependence determined from potential alignment correction. Schemes addressing the potential alignment correction assume localized defects and aim at determining the change in average electrostatic potential, perturbed by the presence of a defect, by keeping the volume in bulk and defect calculations the same. In Figure 4, we show that the origins of the large volume change due to potential alignment is related to the charged bulk and different references (average electrostatic potential) because of the comparison between the strained bulk and the pristine bulk. That is, the large correction in the volume induced by introducing the charged defect is not due to the atomic defect (or very little), but is almost entirely due to the effect of the charge in the system. Here, we have treated potential alignment using the conventional schemes, based on References [36,59] and at present we cannot determine which of the two approaches in Figure 3 is correct, if either. We discuss the volume contribution of charged bulk and potential alignment results explicitly in the discussion section.

Structure of Charged Defects
In this section, we investigate in detail the atomic structure of charged defects, to gain insight into the differences in local lattice relaxation around the neutral and charged defects, and thus understand the origin of the anomalously large relaxation volumes of charged defects.

Lattice Relaxation around an Electron and Hole
An electron or hole is simulated in DFT calculation by changing the total number of electrons by +1 or −1, in the optimized bulk structures. The relaxation volumes of electrons and holes in bulk Si, GaAs, and UO 2 , obtained after relaxing (p = 0) the charged bulk structures, are summarized in Table 5. Addition and removal of an electron results in expansion and contraction of the lattice, respectively. These volume changes (~12 to 22 Å 3 or 0.3 to 0.7Ω 0 , depending on the system) are large and result in about 1% change in their respective lattice parameters for the current supercell sizes. There is not much data in the literature on the relaxation volumes of electrons and holes to compare our results to. However, a few DFT studies [12,25] on defects in Si report the relaxation volumes of an electron of 0.68Ω Si to 0.72Ω Si and a hole of −0.78Ω Si , where Ω Si is the volume per Si atom in a perfect supercell.
Our computed values for the relaxation volume of an electron and hole in Si are in a similar range: 0.74Ω Si and −0.84Ω Si , respectively. A DFT calculation of an electron or hole generally yields a delocalized electron, as it represents the exchange of electron or hole with the bulk conduction band or the valence band edge, unless the added (or removed) extra charge has a preference to localize on a specific atomic site. To verify this, we performed Bader charge analysis [67,68] and looked at the difference in charge density (planar averaged) between the electron (or hole) and bulk for both UO 2 and GaAs.
For GaAs, the difference in charge density ( Figure 5) between the charged and uncharged systems is periodic across the As and Ga atomic planes normal to the z-axis, with no evidence of charge localization. In contrast, for UO 2 , the difference in charge density is not uniform ( To gauge the local atomic relaxations due to the presence of an excess electron or hole, we analyze the change in distances of the neighboring atoms for all the sites in the perfect and defective (electron or hole) structure. Figure 7 shows the relative change in nearest neighbor distances ∆d d 0 for all the sites in the defect structure, plotted as the function of distance (d 0 ) of the given site from the defect. The change in the neighbor distance (∆d = d − d 0 ) is calculated by taking the average difference of distances for each set of neighbors at a given distance from the defect site between the bulk and defect structures.  In UO2 with a partially localized extra electron on a U site ( or U ) (Figure 7a), the local defect structure is complex. Neighboring atoms for any given site in the structure show outward relaxation, independent of the position of the site from the defect. For those sites at a distance of ~5 Å and beyond from the defect, the relative change in neighboring distance saturates to a constant value. This constant value is approximately equal to the fractional volume change . due to the addition of an electron in bulk UO2, which is about ~0.006 (or 0.6%). That is, it is what would be expected from an isotropic expansion of the cell. However, for the sites that are close to the defect (distance < 5 Å) the change in neighboring distances is much larger than 0.6%. Similarly, for a partially localized hole on U site (ℎ • or UU δ+ ) (Figure 7b), the structure exhibits a large inward relaxation of the lattice close to the defect site, and a constant inward relaxation far away from the defect site, equal to the fractional volume change (−0.6%) due to removal of an electron from bulk UO2. This is expected in UO2 since charges are known to form small polarons (known both experimentally and using higher levels of theory like DFT+U or hybrid-DFT), which by definition carry local lattice distortions [70]. By contrast, for GaAs and Si the electron (or hole) is completely delocalized and does not reside on any specific atomic site; therefore, the change in the neighboring distance is independent of the site considered: an electron (or a hole) induces a constant expansion in the lattice that is uniform over  In UO2 with a partially localized extra electron on a U site ( or U ) (Figure 7a), the local defect structure is complex. Neighboring atoms for any given site in the structure show outward relaxation, independent of the position of the site from the defect. For those sites at a distance of ~5 Å and beyond from the defect, the relative change in neighboring distance saturates to a constant value. This constant value is approximately equal to the fractional volume change . due to the addition of an electron in bulk UO2, which is about ~0.006 (or 0.6%). That is, it is what would be expected from an isotropic expansion of the cell. However, for the sites that are close to the defect (distance < 5 Å) the change in neighboring distances is much larger than 0.6%. Similarly, for a partially localized hole on U site (ℎ • or UU δ+ ) (Figure 7b), the structure exhibits a large inward relaxation of the lattice close to the defect site, and a constant inward relaxation far away from the defect site, equal to the fractional volume change (−0.6%) due to removal of an electron from bulk UO2. This is expected in UO2 since charges are known to form small polarons (known both experimentally and using higher levels of theory like DFT+U or hybrid-DFT), which by definition carry local lattice distortions [70]. By contrast, for GaAs and Si the electron (or hole) is completely delocalized and does not reside on any specific atomic site; therefore, the change in the neighboring distance is independent of the site considered: an electron (or a hole) induces a constant expansion in the lattice that is uniform over   These calculations show that the addition or removal of an electron to or from a bulk lattice can result in delocalized or localized charge depending on the system. It is intriguing that independent of the behavior (localized or delocalized) of the charge, there is an expansion (or contraction) of the bulk lattice, resulting in significant relaxation volumes for electrons and holes determined within the DFT framework. Whether such an effect is real or a behavior particular to DFT calculations is not clear at present. However, a detailed examination like this of electrons and holes in bulk structures helps us to establish the dependence of relaxation volume and corresponding elastic dipole tensor purely on the behavior of the charged bulk system.

U Vacancy in UO2
The U vacancy in the fluorite structure resides at the center of the cube formed by the first eight neighboring oxygen atoms. The difference between neutral and charged supercell calculations lies in In UO 2 with a partially localized extra electron on a U site (e or U δ− U ) (Figure 7a), the local defect structure is complex. Neighboring atoms for any given site in the structure show outward relaxation, independent of the position of the site from the defect. For those sites at a distance of~5 Å and beyond from the defect, the relative change in neighboring distance saturates to a constant value. This constant value is approximately equal to the fractional volume change ∆Ω r Ω(bulk) = 11.52 1936.76 due to the addition of an electron in bulk UO 2 , which is about~0.006 (or 0.6%). That is, it is what would be expected from an isotropic expansion of the cell. However, for the sites that are close to the defect (distance < 5 Å) the change in neighboring distances is much larger than 0.6%. Similarly, for a partially localized hole on U site (h. or U U δ+ ) (Figure 7b), the structure exhibits a large inward relaxation of the lattice close to the defect site, and a constant inward relaxation far away from the defect site, equal to the fractional volume change (−0.6%) due to removal of an electron from bulk UO 2 . This is expected in UO 2 since charges are known to form small polarons (known both experimentally and using higher levels of theory like DFT+U or hybrid-DFT), which by definition carry local lattice distortions [70]. By contrast, for GaAs and Si the electron (or hole) is completely delocalized and does not reside on any specific atomic site; therefore, the change in the neighboring distance is independent of the site considered: an electron (or a hole) induces a constant expansion in the lattice that is uniform over all atoms. Structure analysis of an electron or hole in GaAs (Figure 7c,d) shows that the change in neighboring distances between the defect and bulk structures is approximately a constant value. The structure with an extra electron shows a homogeneous expansion of~0.49% (Figure 7c) while that with an extra hole shows a homogeneous contraction of~0.40% of the lattice (Figure 7d), resulting in roughly ±1% change respectively, in the lattice parameter.
These calculations show that the addition or removal of an electron to or from a bulk lattice can result in delocalized or localized charge depending on the system. It is intriguing that independent of the behavior (localized or delocalized) of the charge, there is an expansion (or contraction) of the bulk lattice, resulting in significant relaxation volumes for electrons and holes determined within the DFT framework. Whether such an effect is real or a behavior particular to DFT calculations is not clear at present. However, a detailed examination like this of electrons and holes in bulk structures helps us to establish the dependence of relaxation volume and corresponding elastic dipole tensor purely on the behavior of the charged bulk system.

U Vacancy in UO 2
The U vacancy in the fluorite structure resides at the center of the cube formed by the first eight neighboring oxygen atoms. The difference between neutral and charged supercell calculations lies in the presence or absence of local charge compensation mechanisms around the defect site [27,28]. We examine these differences by explicitly analyzing the difference in charge density and Bader charge in neutral and charged supercell calculations. The neutral uranium vacancy is charge compensated by four U sites (h. or U U δ+ ) with Bader charges of 2.79-2.81 (much larger than Bader charge of 2.54 on U site in bulk UO 2 ) located at 2nd (3.86 Å) and 4th (5.44 Å) neighboring distances from the defect site. However, for the charged (4−) uranium vacancy there is no local charge compensation from neighboring U sites.
The changes in nearest neighbor distances (Figure 8) for the both the neutral and charged U vacancy clearly reveal the differences in the short-and long-range atomic relaxation around the defect site. The lattice distortion around the neutral vacancy extends to a relatively larger distance (~6-7 Å) compared to the charged vacancy (~4-5 Å), because of the presence of neighboring U U δ+ , which are absent for the charged vacancy. For both the neutral and charged supercell, the relative change in neighboring distances saturates away from the defect site. The fractional volume change for neutral and charged U vacancy are −0.46% and 2.1%, respectively, in a 144-atom supercell. Analysis in a larger 324-atom supercell yields similar trends in the local and long-range defect structure of neutral and charged U vacancies.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 16 of 24 the presence or absence of local charge compensation mechanisms around the defect site [27,28]. We examine these differences by explicitly analyzing the difference in charge density and Bader charge in neutral and charged supercell calculations. The neutral uranium vacancy is charge compensated by four U sites (ℎ • or UU δ+ ) with Bader charges of 2.79-2.81 (much larger than Bader charge of 2.54 on U site in bulk UO2) located at 2nd (3.86 Å) and 4th (5.44 Å) neighboring distances from the defect site. However, for the charged (4−) uranium vacancy there is no local charge compensation from neighboring U sites.
The changes in nearest neighbor distances (Figure 8) for the both the neutral and charged U vacancy clearly reveal the differences in the short-and long-range atomic relaxation around the defect site. The lattice distortion around the neutral vacancy extends to a relatively larger distance (~6-7 Å) compared to the charged vacancy (~ 4-5 Å), because of the presence of neighboring UU δ+ , which are absent for the charged vacancy. For both the neutral and charged supercell, the relative change in neighboring distances saturates away from the defect site. The fractional volume change for neutral and charged U vacancy are −0.46% and 2.1%, respectively, in a 144-atom supercell. Analysis in a larger 324-atom supercell yields similar trends in the local and long-range defect structure of neutral and charged U vacancies. The dashed lines mark the % fraction volume change of the entire supercell, relative to a defect-free neutral system, to show the overall expansion or contraction of the lattice in the respective cases.
In corresponding investigations of charged vacancies in GaAs and Si we find that vacancies show short-ranged local distortion around the localized defect sites. However, changes in nearest neighbor distances away from the localized defect site are manifestations of the overall volume change, which depend on the charge state of the defect.

Defect Reactions versus Individual Defects
In the preceding discussion, we have focused on the properties of isolated point defects. However, in reality, charged defects are never created in isolation, but by means of a series of reactions involving multiple defects; this is typically described via Kröger-Vink notation [71]. These reactions obey strict conservation rules of mass, charge, and site (maintenance of constant ratio of cationic to anionic sites in the structure) balance. There have been several theoretical studies [9,27,28,72] that develop a point defect reaction model for UO2 in order to predict defect concentration, formation energies, lattice constant, and transport properties under thermodynamic conditions prevailing in experiments [26]. In this section, we employ a similar approach based on a point defect model to compute defect reaction volumes of the net reaction, as opposed to the individual constituents. We illustrate this by using the example of the oxygen interstitial and the uranium vacancy, which are more likely to form in UO2 under oxidizing/transmutation conditions. The dashed lines mark the % fraction volume change of the entire supercell, relative to a defect-free neutral system, to show the overall expansion or contraction of the lattice in the respective cases.
In corresponding investigations of charged vacancies in GaAs and Si we find that vacancies show short-ranged local distortion around the localized defect sites. However, changes in nearest neighbor distances away from the localized defect site are manifestations of the overall volume change, which depend on the charge state of the defect.

Defect Reactions versus Individual Defects
In the preceding discussion, we have focused on the properties of isolated point defects. However, in reality, charged defects are never created in isolation, but by means of a series of reactions involving multiple defects; this is typically described via Kröger-Vink notation [71]. These reactions obey strict conservation rules of mass, charge, and site (maintenance of constant ratio of cationic to anionic sites in the structure) balance. There have been several theoretical studies [9,27,28,72] that develop a point defect reaction model for UO 2 in order to predict defect concentration, formation energies, lattice constant, and transport properties under thermodynamic conditions prevailing in experiments [26].
In this section, we employ a similar approach based on a point defect model to compute defect reaction volumes of the net reaction, as opposed to the individual constituents. We illustrate this by using the example of the oxygen interstitial and the uranium vacancy, which are more likely to form in UO 2 under oxidizing/transmutation conditions.
The reaction for incorporating oxygen at an interstitial site is given as where V × i is a vacant octahedral interstitial site, in the fluorite structure, O i is the oxygen ion at the interstitial site, and h. represents an U 5+ ion (hole) not bounded to the oxygen interstitial O i . The volume change, ∆Ω rxn , for the reaction (Equation (12)) is then given as the sum of the individual relaxation volumes (as listed in Table 2, without any corrections to volumes), ∆Ω r , obtained from separate DFT supercell calculations of the charged (2−) oxygen interstitial and the U 5+ ion, The reaction volume ∆Ω rxn O i for the oxygen interstitial reaction is computed to be −25.8 Å 3 , implying contraction of the lattice due to formation of this overall neutral collection of defects (note that, in this calculation, the interstitial and the holes are assumed to be infinitely far apart). The sum of the relaxation volumes of charged oxygen interstitial and holes is −3.2 Å 3 , similar to the relaxation volume, ∆Ω r O × i of a neutral oxygen interstitial (−3.38 Å 3 ) calculated within a DFT supercell, which is not surprising if the neutral O interstitial can be understood as a charged oxygen interstitial with two bound holes.
Similarly, a defect reaction for the formation of a uranium vacancy, is given as: where point defects V U and h. are the fully charged uranium and U 5+ ion (hole), respectively, and (UO 2 ) SRG represents the corresponding lattice atoms present at a site of repeatable growth (SRG) such as a surface, dislocation or grain boundary. The volume change for the reaction (Equation (14)) is then given as: involving individual charged-defect relaxation volumes ∆Ω r V U and ∆Ω r (h.), calculated again from separate supercell calculations. The overall computed reaction volume ∆Ω rxn V U is −10.8 Å 3 , which again corresponds to contraction of the lattice due to the formation of uranium vacancies. These findings are in good agreement with the experimental results [73] and with theoretical prediction [74] for the contraction of the UO 2 lattice due to the formation of the hyper-stoichiometric phase UO 2+x (with x < 0.5). The sum of the relaxation volumes of individual charge defects (first two terms in Equation (15)) is about −6.26 Å 3 , similar but not exactly equal to that predicted by the relaxation of neutral U vacancy (∆Ω r V × U = −9.07 Å 3 ) in a DFT supercell calculation. Our results comparing the relaxation volume of the neutral defects with the summed relaxation volumes of the charged defects calculated from separate supercells calculations are shown in Table 6. There is good agreement with regards to the sign of the lattice relaxation (expansion or contraction) due to the formation of an overall neutral defect reaction with the two approaches. The difference in the predicted magnitude can be attributed to the interaction (both elastic and electrostatic) between the individual charged defects. Thus, while the relaxation volumes of the charged defects have large magnitudes, they to a great extent cancel when summed, leading to values that are in reasonable agreement with corresponding neutral defects. Similarly, while defect reactions, such as those in Equations (12) and (14), contain components corresponding to charged defects, once these are summed in a net neutral reaction, the large volumes of the charged defects cancel, leading to relatively small volumes for the reaction itself. In a recent study by Bruneval et al. [9], similar conclusions were also reached, where authors showed that, using a net neutral defect reactions model, relaxation volumes can reproduce the experimental observed lattice-constant variation with stoichiometry. Table 6. Comparison of the defect relaxation volume (Å 3 ) of neutral individual defects computed within a single (one) supercell (labeled OS) or via a sum of defect relaxation volumes of individual charge defects computed in separate supercells.

Discussion
We have found that the magnitude of the dipole tensors of individual charged defects leads to seemingly unphysical behavior. For example, combined with elasticity theory, the dipole tensors calculated here lead to segregation energies of charge defects to dislocations that are extremely large.
In an effort to understand the origin of this behavior for charged point defects, we have applied volume correction schemes that are analogous to those proposed for correcting the energetics of charged defects.
Our results indicate that the corresponding changes to defect volumes (and thus to dipole tensors) can either be relatively very small or quite large depending on how the corrections are applied, presenting an uncertainty as to which, if either, approach is correct.
In contrast, if one considers a full defect reaction and the associated relaxation volume, while the individual contributions correspond to charged defects and thus have large relaxation volumes, the volume change of the overall reactions are physically reasonable and are similar to those of neutral defects. Thus, the large relaxation volumes of negative defects are essentially canceled out by large, but oppositely signed, relaxation volumes of positive defects. This suggests that one cannot think of a charged defect in isolation but only in the context of the overall defect reactions that produce the defect in the bulk phase. However, in numerous fields, we are accustomed to thinking of the properties of isolated charged defects as entities that migrate and interact in the material. For example, in fast ion conductors, we consider the transport of charged carriers in the material as specific entities. Thus, it seems reasonable that one should be able to describe the elastic response of the medium due to a charged defect and, thus, determine how it interacts with, for example, dislocations. Therefore, the above discussion suggests that, perhaps, the problem lies with the reference state reservoir. The defect formation volume depends on the choice of reservoir, just as does the defect formation energy. Considering a neutral defect, the atoms removed or added from the perfect bulk material have to be placed in a reservoir and the volume of both the material and the reservoir changes. Analogously, for charged defects, the electrons or holes added to the defect system also contribute to both the defect formation energy and volume. While the contribution to the defect formation energy has been well studied and several correction terms have been developed, the contribution of the change in charge to the volume of the reservoir has not been considered. From our results for three different materials systems, we can estimate the contribution of the net charge to the defect formation volume, as shown in Figure 9. We see a linear relationship between the electron count and the defect relaxation volume in all three systems, similar to that observed in previous DFT studies [12][13][14]. Defect formation volumes, as defined in Equation (5), require Ω 0 to be added or subtracted to the defect relaxation volume to account for the volume of the individual species added or removed from the reservoir to maintain overall mass balance. However, in Equation (5), the contribution of the change in charge to the volume of the electron reservoir is missing. The change in the volume of the bulk with charge can be attributed to the volume contribution per electron as coming from the change in volume of the electron reservoir, ∆Ω re f . Therefore, a complete definition of defect formation volume should be given as: Appl. Sci. 2019, 9, x FOR PEER REVIEW 19 of 24 ∆ is determined by the slope of the bulk relaxation volume as a function of charge. Values of ∆ are system dependent and for UO2, GaAs, and Si we find −11.6 Å 3 , −23.6 Å 3 and −16.14 Å 3 per unit charge, respectively. For Si, Windl and Daw [13] reported a value of about −15 Å 3 from the slope of the change in volume of the bulk with net charge, similar to the value we have obtained. One can then recompute defect formation volumes by accounting for ∆ in all three systems, as summarized in Table 7. Reaction volumes are not affected by ∆ as the reaction is overall charge neutral and, hence, the change in volume of the electronic reservoir due to charge cancels out exactly.
In our discussion so far we have assumed negligible contribution to volume from the potential alignment correction, consistent with the results described in Figure 3a. However, as discussed in Section 3.1.2, the volume correction from potential alignment can be significant, depending on the adopted approach. As mentioned above, the difference in the two potential alignment approaches ( Figure 4) comes down to how the volume of the reference bulk calculation is treated. That is, when considering the volume change induced by a charged defect and the associated corrections, do we compare to the original neutral bulk or a bulk cell expanded to the same volume as the defective state? Figure 4 indicates that the major contribution to the correction is associated with adding charge to the bulk. Comparing the change in volume of the electron reservoir ∆ from Figure 9 with the volume correction from potential alignment correction ∆ , we find that both are large (c.f. Figure  4) and significantly affect the final defect relaxation volumes. Corrections based on potential alignment with respect to pristine bulk ( Figure 3b) agree with the conclusions made in Reference [13], where potential alignment correction was shown to have a huge volume contribution to the charged defects relaxation volumes and likely provide a different viewpoint of correcting charged defect relaxation volumes.

5.75
The charged defect elastic dipole tensor is theoretically computed from its relaxation volume ∆Ω re f is determined by the slope of the bulk relaxation volume as a function of charge. Values of ∆Ω re f are system dependent and for UO 2 , GaAs, and Si we find −11.6 Å 3 , −23.6 Å 3 and −16.14 Å 3 per unit charge, respectively. For Si, Windl and Daw [13] reported a value of about −15 Å 3 from the slope of the change in volume of the bulk with net charge, similar to the value we have obtained. One can then recompute defect formation volumes by accounting for ∆Ω re f in all three systems, as summarized in Table 7. Reaction volumes are not affected by ∆Ω re f as the reaction is overall charge neutral and, hence, the change in volume of the electronic reservoir due to charge cancels out exactly. In our discussion so far we have assumed negligible contribution to volume from the potential alignment correction, consistent with the results described in Figure 3a. However, as discussed in Section 3.1.2, the volume correction from potential alignment can be significant, depending on the adopted approach. As mentioned above, the difference in the two potential alignment approaches ( Figure 4) comes down to how the volume of the reference bulk calculation is treated. That is, when considering the volume change induced by a charged defect and the associated corrections, do we compare to the original neutral bulk or a bulk cell expanded to the same volume as the defective state? Figure 4 indicates that the major contribution to the correction is associated with adding charge to the bulk. Comparing the change in volume of the electron reservoir ∆Ω re f from Figure 9 with the volume correction from potential alignment correction ∆Ω PA corr , we find that both are large (c.f. Figure 4) and significantly affect the final defect relaxation volumes. Corrections based on potential alignment with respect to pristine bulk ( Figure 3b) agree with the conclusions made in Reference [13], where potential alignment correction was shown to have a huge volume contribution to the charged defects relaxation volumes and likely provide a different viewpoint of correcting charged defect relaxation volumes. The charged defect elastic dipole tensor is theoretically computed from its relaxation volume (Equation (1)), and experimental works often report it as chemical strain [7,22,75], estimated by measuring the change in the lattice parameter as a function of deviation from the stoichiometric chemical formula. The apparent question that comes next is should the charged defect relaxation volume be corrected to account for the volume associated with the electron/hole reservoir.
The change in the volume of the bulk with charge (volume associated with the electron/hole reservoir) as given by ∆Ω re f can be understood as the lattice perturbation associated with excess electron and/or holes in the bulk, as described in Section 3.2.1. We reason that charged defects are individual entities that are present in an overall charge neutral system. Both their elastic dipole tensor and relaxation volumes are meaningful when defined with respect to the electron/hole reservoir in order to keep the overall system charge neutral, as one would expect when measuring these defect properties from experiments. Traditionally, when one calculates the elastic dipole tensor and relaxation volume of a charged defect, it is relative to the neutral bulk. Perhaps the correct reference for a charged atomic defect in a lattice is not the neutral bulk, but the charged bulk, so that the defect represents a perturbation of the charged lattice only, and not the neutral lattice.
The elastic dipole tensor G ij for a charged defect when defined in an overall neutral (with respect to the electron or hole reservoir) but locally charged lattice is then given as: where, ∆Ω r kl − q∆Ω re f represents the corrected relaxation volume. As listed in Table 8, the corrected relaxation volumes of charged defects are more reasonable in magnitude when compared to relaxation volumes for both the neutral individual defects and those of defect reactions. For example, the corrected relaxation volume of −3.28 Å 3 for the charged (2−) oxygen interstitial (overall neutral i.e., the oxygen interstitial with respect to doubly charged bulk UO 2 ) compares very well with the volume of −3.20 Å 3 of the reaction involving charged (2−) oxygen interstitial with two holes. Similar comparisons for other point defects in UO 2 are also very good, because as can be seen in Figure 9, the slopes of change in volume of defects with charge are almost parallel to the slope of change in volume of bulk with charge. However, for Si and GaAs we observe that different defects have different slopes, suggesting that the effect of electrons and holes present in different defects leads to different strains and volume changes. We also acknowledge that the difference in the dependence of the change in volume on charge between different charged defects can originate due to the deep (localized) vs. the shallow (delocalized) nature of the defect [15]. However, in this work we do not assess the impact of deep vs. shallow nature of the defect on defect relaxation volumes because in both cases the relaxation volumes show strong sensitivity to the charged state. Table 8. Comparison of the corrected defect relaxation volume [∆Ω r − q∆Ω re f ] (Å 3 ) of neutral and charged individual defects, with the sum of defect relaxation volumes of individual charge defects computed in separate supercells under the defect reaction approach. Some of these values appeared in Table 6 but are reproduced here for convenience. Computing the dipole tensor (Equation (17)) of charged defects from the corrected relaxation volumes does result in physically reasonable magnitudes in the change in defect formation energy of the charged point defect upon interaction with external strains. Accounting for the change in volume of the electron reservoir with net charge in the definition of the defect formation volume and elastic dipole tensor allows us to reconcile these two perspectives: that defect behavior not only has meaning via reactions but that isolated defects in an overall neutral lattice (which requires the reservoir) also have unique properties that determine how they behave within the lattice.

Defect Reaction
Finally, returning to Figure 4 and the associated discussion, there is an interesting parallel in the concept of the charge reservoir highlighted in Figure 9 and the corrections described in Figure 4. In Figure 4, we decomposed the potential alignment correction into two parts, one due to the net charge introduced into the cell and another due to the defect itself. On the other hand, as highlighted by Figure 9 and Table 8, the largest component of the absolute relaxation volume is due to the net charge, not the defect. That is, the biggest contributor to both the relaxation volume and the corrections to the relaxation volume is the net charge, not the atomic relaxations due to the atomic defect. This highlights the importance of a physically meaningful description of how net charge impacts the volume of a system. Thus, the conclusion seems clear: to calculate the volume of a charged defect, one must be careful in how the source of the net charge is accounted for, in the same way one does for energy. This implies that the relaxation volumes associated with the defect depend on that reference or reservoir. We note that this choice is less important when considering defect reactions as the volume of the reservoir largely cancels out in these cases, unless one is exchanging charge with the environment.

Conclusions
To conclude, we have shown that point defect relaxation volumes, and the associated defect elastic dipole tensors which measure the elastic distortion of the bulk lattice caused by the presence of the defect, are unexpectedly large for charged defects in UO 2 . Defect relaxation volumes have a finite measure even for infinitely large system sizes, and we show that the large relaxation volume of charged defects is not due to finite-size artifacts within the DFT supercell approach, as previously suspected. Relaxation volumes using an overall neutral defect reaction model based on DFT calculations of isolated charged defects agree reasonably well in terms of the overall nature of lattice relaxation and magnitude when compared with the individual neutral defects simulated within a single DFT supercell. We conclude that, similar to the defect formation energies, the defect formation volumes also depend on the choice of reservoir. We show that by taking into account the contribution of change in charge to the volume of the electron reservoir, the recomputed charged defect formation volumes and elastic dipole tensors are reasonable in magnitude, and that it is possible to describe the elastic properties of isolated charge defects within the material, not only in the context of the overall defect reactions that produce the defect but also as individual species in an overall neutral lattice. This provides a route for using atomistic studies to determine a key property-the dipole tensor-for use in more applied studies that account for the interaction of defects with microstructure, including irradiation-induced creep, long-term radiation damage evolution, and mass transport in complex microstructures.