Effect of electric field on migration of defects in oxides : Vacancies and interstitials in bulk MgO

Al-Moatasem El-Sayed,1,2,* Matthew B. Watkins,3,† Tibor Grasser,2,‡ and Alexander L. Shluger1,§ 1Department of Physics and Astronomy and London Centre for Nanotechnology, University College London, Gower Street, London, WC1E 6BT, United Kingdom 2Institute for Microelectronics, Technische Universität Wien, A-1040 Vienna, Austria 3School of Mathematics and Physics, University of Lincoln, Brayford Pool, Lincoln LN6 7TS, United Kingdom


I. INTRODUCTION
Understanding the effect of electric fields on the migration of defects and ions is important for elucidating mechanisms of degradation of dielectric films in microelectronic devices, diffusion of ions in nanoionics and in electrocatalysis.In particular, metal oxides used as gate dielectrics in microelectronic devices experience prolonged exposure to strong electric fields.Defects are implicated in a number of reliability issues in electronic devices, such as dielectric breakdown, bias temperature instability, 1 / f noise, random telegraph noise, and stress-induced leakage current [1,2].They also play a central role in resistive switching mechanisms in oxide films and it is thought that applying a bias could lead to the creation of new defects [2][3][4][5][6].Understanding the mechanisms of defect creation and migration processes is required for better design of microelectronic devices as well as improving their reliability.However, the effect of bias application on defects in oxides is still poorly understood at the atomistic level.
The effect of applied electric fields is accounted for in phenomenological models of devices through parameters that modify the energetics of charge transfer and defect diffusion processes.In the early studies of the formation of thin metaloxide films, Cabrera and Mott described interstitial oxygen migration as a simple ion hopping mechanism whose barrier is lowered by 1  2 qaF in the direction of the field, where q is the net charge of the ion, a is the distance it migrates, and F is the strength of the applied field [7].McPherson et al. went on to develop a model known as the thermochemical E model, where the activation enthalpy of a defect process is modified by the application of an external electric field [8,9].More formally, this is expressed as: H F = H 0 − pF , where H F is the enthalpy change under an applied field of strength, F , p is a quantity known as the effective dipole moment, and H 0 is the enthalpy change when no field is applied.This model has been used to explain field dependent processes and dielectric breakdown in SiO 2 and related phenomena in other materials [10].Further molecular calculations have been applied to implicate a Si-O bond breakage process in time-dependent dielectric breakdown [11,12].However, these calculations were largely done on SiO 4 clusters using classical potentials whose parameters were derived from ab initio calculations.Furthermore, the effective dipole moments were extracted empirically from the valence and symmetry of the defect [10].Effective dipole moments were linked to the equilibrium geometry of the structure and it is, therefore, unclear how this notion should be interpreted in a dynamic process.
Electric fields are directly included in a growing number of atomistic simulations of defect processes.A sawtooth potential with a constant slope throughout the 2D periodic slab model has been used in density functional theory (DFT) calculations of migration of oxygen vacancies at the anatase surface [13].A similar approach [14] has been used to study the effect of electric field on the Cu migration in SiO 2 [15] and O vacancy migration in amorphous InGaZnO [16] using a slab model.The advent of the modern theory of polarization [17][18][19][20] has allowed one to include an applied external field in 3D periodic DFT calculations of materials and calculate properties derived from the dipole moment and polarization of periodic systems [21,22].Such simulations have helped to understand how water dissociation and proton conduction is affected by an applied field [23,24] and even to shed light on the renowned Miller experiments [25].Being able to calculate the changes in polarization has even allowed for a rigorous definition of the oxidation states of ions [26].
In solids, the thermodynamics of neutral oxygen vacancies in binary oxides has been analyzed in an electric field [27].It has been shown that the Gibbs free energy of formation of neutral vacancies can be significantly lowered by the field due to the work of polarization and the ease with which the two electrons trapped at the neutral vacancy can be polarized.The Berry phase approach within the modern theory of polarization was employed to calculate polarization properties of SrTiO 3 with the spontaneous polarization defined as the difference in polarization between the polar and nonpolar reference states [28].Polarization switching barriers were calculated along the minimum energy path between two polarization states.Minor alterations to the formalism allow the usage of constant electric displacement rather than the electric field, extending its application to open or closed circuit conditions [29].Despite this progress, detailed microscopic investigations of how adiabatic barriers for defect migration in the bulk of a solid are affected by the field are still lacking.Such calculations are important to test numerous phenomenological theories of the effect of electric field on the degradation of dielectric films, on the behavior of ions in nanoionics and in electrocatalysis.
In this paper, we present ab initio calculations of the effect of electric field strength on the barriers for oxygen vacancy and interstitial ion migration in numerous charge states in the bulk of MgO.The high symmetry of cubic MgO makes the analysis of the results more transparent.However, our main conclusions should be valid for more complex oxides where migration of oxygen vacancies and interstitial ions is commonplace.We show how the changes in the dipole moment of the system along the migration trajectory define the magnitude of the barrier change under an applied field.We find that barriers for neutral oxygen vacancy and interstitial migration change negligibly even in a strong field of 10 MV cm −1 whereas charged defects experience much bigger changes.We also show how the application of a strong field changes the profile of the adiabatic barrier for migration.Our results quantify changes in the dipole moment along defect migration processes and link it to how the barrier is affected by an applied field.This opens the door to quantification of applied fields' effects on defect processes in other technologically relevant and more complex oxides.

A. Computational details
The calculations for crystalline MgO were carried out in a 3 × 3 × 3 supercell containing 216 atoms using the PBE functional [30].All calculations were performed using the CP2K code [31] at the point and employed a mixed Gaussian and plane-waves basis set [32] to represent the electrons in the system.Test calculations performed with increasing realspace cell volumes of MgO demonstrated that the electronic structure truly converges at the 216 atom cell size.Periodic boundary conditions were enforced in all calculations apart from calculations of molecular polarizabilities presented in the Appendix, where fixed boundary conditions were employed.Double-ζ Gaussian basis sets optimized for molecules and condensed phase systems using the method of VandeVondele [33] were employed for Mg and O ions.The molecules used for optimizing the basis sets came from the work by Weigend [34].These basis sets were used in conjunction with the Goedecker, Teter, and Hutter pseudopotentials [35].The plane wave cutoff used in these calculations was set to 8844 eV (650 Ry).All geometry optimizations were performed using the BFGS optimizer to minimize forces on atoms to within 37 pN (2.3 × 10 −2 eV Å−1 ).Migration barriers at zero field were calculated using the climbing image nudged elastic band (NEB) method [36].Linear interpolation was used to generate three frames to be optimized, with each of the frames connected by a spring with a force constant of 19.5 eV Å−2 .We ran benchmark calculations with increasing numbers of NEB frames and found that only three frames were necessary to obtain a good description of the migration barrier for the defects considered in this paper.We note that all the migration barriers obtained under zero field are highly symmetric and the optimization resulted in trajectories with equal spacings between the frames.The spacings between the climbing image NEB frames of the trajectories (calculated as the norm of the displacements) of all three charge states of the O vacancy differed by up to 0.03 Å at most.Therefore in the Results section, properties involving the migration trajectories are plotted against the frame number (rather than, e.g., O ion displacement), which are equidistant from each other.
To calculate the effect of an electric field on the migration barriers, we used the following approach.Using the trajectories obtained from the NEB calculations at zero field, we applied a linear interpolation between the three frames to obtain eleven interpolated frames.We then applied an electric field with eleven different strengths along the direction of the defects' migration and a single point calculation was carried out for each frame along the interpolated migration trajectory.Ionic cores were fixed in these calculations in their interpolated positions as full geometry relaxation would have made the calculations prohibitively expensive.However, test NEB calculations on the F 2+ center at 10 MV cm −1 field found that the positions of the ions varied negligibly with respect to those at zero field.

B. Including electric field
We use the modern theory of polarization and the Berry phase operator methods to include an electric field in defect calculations [17][18][19][20][21][22].In this approach, the electric polarization of the material is calculated as the Berry phase of the periodic components of the Bloch states.An electric enthalpy functional can be defined in terms of the polarization as where H E [ρ; F] is the electric enthalpy functional, E KS [ρ; F] is the ground state Kohn-Sham energy obtained from a standard DFT calculation, is the cell volume, P[ρ; F] is the polarization of the system, and F is the electric field [21,22].The electric enthalpy is a functional of the electron density ρ and depends parametrically on the electric field.After this section, we will suppress these dependencies for notational convenience.The polarization is calculated from the Berry phase of the Kohn-Sham states as implemented in Ref. [37].Once the polarization has been calculated, it is trivial to extract the ab initio dipole moment as μ[ρ; F] = P[ρ; F].This method is designed for 3D periodic systems, making it ideal for studying processes in bulk solids.It allows one to explicitly include the electric field in the self-consistently obtained electron density.When using this approach, which we will refer to as the Berry phase method, the polarization is defined to within one modulo of the polarization quantum.This is due to the phase being defined to within 2π .To compare the results, one has to make sure that they belong to the same branch of polarization (see discussion in Ref. [20]).All results presented in this paper have been analyzed and corrected by the polarization quantum, 2π R, where R is the lattice vector, to ensure that they do indeed belong to the same polarization branch.
To test our use of the implementation of this method in CP2K [37], we calculated the isotropic polarizabilities of several small molecules using the finite field technique [38,39].The field was applied both using the Berry phase method and via including a linear change in the potential in the Hamiltonian of a finite system.The Berry phase method was also used to calculate the high frequency and static permittivity of bulk MgO.The results of these calculations are given in the Appendix and demonstrate good agreement between the two methods and with experimental data.
The defect migration barriers presented in this paper are calculated with static, homogeneous electric fields applied to the system.We use a range of field strengths up to 10 MV cm −1 .The latter value is comparable to breakdown fields characteristic of many oxides.In the calculations for charged systems, a neutralizing background was included in order to compensate for the net charge of the cell and prevent calculations of the Hartree energy from diverging.The dipole moment of charged systems is implicitly dependent on the origin chosen to calculate it from.The dipole moments of charged systems discussed in the Results section are adjusted using the defect configuration under zero field as a reference point so that the lowest value under zero field equals zero.

C. Analysis of results using simple model
We note that the calculations carried out using the Berry phase method are self-consistent and take into account changes in the electron density induced by the electric field.To compare the results with simple classical models of the type presented in Ref. [7] we use the following approximations.For moderate field strengths, we assume that the electron density does not change upon application of the field.The energy of a charge distribution under an applied field can then be described as where E F and E 0 are the energies of the charge distribution with and without an applied field, respectively, μ is the dipole moment of the charge distribution without the applied field, and F is the applied field.The last term in this equation is the work of polarization.Now let us consider the initial and transition states of a defect undergoing migration.The energy difference between the initial and transition states at zero field is the barrier to migration, E 0 B .The energy difference between these two states, herein referred to as 1 and 2, upon applying a field is then equal to where the energy at the transition state under zero field and its associated dipole moment are E 0 2 and μ 2 and those of the initial configuration, also under zero field, are E 0 1 and μ 1 , respectively.Note that the term E 0 2 − E 0 1 corresponds to the energy barrier for the process, E 0 B , at zero field.We can simplify the final term in equation ( 3) by defining an effective dipole moment, μ eff , to be μ 2 − μ 1 .To estimate how the barrier energy has changed due to an applied field we subtract the barrier under zero field term, E 0 B , to give In this approximation, the effective dipole moment at zero field determines the strength of the barrier change in a field.We will use an effective dipole moment, μ eff , at zero field together with Eq. ( 4) to compare the barrier changes in the field calculated using this approximation and the ab initio dipole moments calculated using the Berry phase method.

III. RESULTS OF CALCULATIONS
The geometry of a 3 × 3 × 3 supercell of MgO was optimized in the neutral charge state using the PBE functional.The obtained Mg-O interatomic distance of 2.10 Å is in good agreement with the experimental value of 2.11 Å [40].The calculated band gap of 4.65 eV is, however, considerably lower than the experimental value of 7.9 eV [41] but to be expected when using the PBE functional.The remainder of the paper reports the results of calculations for the oxygen vacancy and interstitial defects in MgO.For each defect, we have calculated their properties and migration barriers in three different charge states.The dipole moment changes along the migration pathways have been calculated.Finally, for each defect, we show how the application of a field affects the calculated barrier.

A. Oxygen vacancies
To study oxygen vacancy migration barriers under an applied electric field, one O atom was removed from the cell and the geometry was optimized in three charge states: neutral, positive, and double positive.These oxygen vacancies are the canonical color center defects in MgO: the F 0 , F + , and F 2+ centers [42].Introducing the vacancy in its neutral charge state changes the system's atomic structure negligibly.However, the O and Mg ions displace toward and away from the vacancy, respectively, in the F + and F 2+ centers.The displacements are largest in the F 2+ center, as is to be expected.Each center creates Kohn-Sham (KS) states in the band gap.The F 0 center induces a doubly occupied KS state at 2.4 eV above the valence band edge.The F + center introduces singly occupied and unoccupied KS states located at 1.77 and 2.99 eV above the valence band edge, respectively.The F 2+ center induces a doubly unoccupied state at 1.77 eV above the valence band.These properties do not change if the calculations are carried out in larger cells.
The migration barriers for the F centers calculated using the NEB method at zero field are shown in Table I along with the results of previous studies.The calculated barriers of 4.50, 3.54, and 2.24 eV for the F 0 , F + , and F 2+ centers, respectively, are in fair agreement with the results of Mulroue et al. [43] who also used DFT to calculate the F center migration barriers.
Qualitatively, the migration pathway is the same for all of the F centers: An O atom moves along the 110 direction into the vacancy.In further discussion, the electric field is applied along that direction, i.e., along the direction of vacancy migration.Figure 1 shows how the dipole moments of the three vacancy charge states change along the migration pathway.The dipole moments are extracted from the polarization as μ = P.The x axis is the displacement of the vacancy along the 110 direction in terms of the trajectories' frames.There are eleven lines on each graph corresponding to the field strength ranging from zero up to 10 MV cm −1 with a 1 MV cm −1 interval.Each data point in the graph was calculated as a single-point calculation, that is, without optimizing the positions of the ionic cores.The black, bottom line in each panel is the dipole moment with no applied field and the reddest line, which is always at the top, corresponds to a field of 10 MV cm −1 .The panels show the F 0 , F + , and F 2+ centers from top to bottom.We note that a negative dipole moment in this context means a dipole moment aligned 180 • with respect to the direction of the applied field while a positive one is aligned along the field's direction.
Let us first consider the results for the F 0 center, shown in the top panel of Fig. 1.The starting dipole moment is zero due to the high symmetry of the vacancy and increases as an oxygen atom displaces from the neighboring site towards the The panels show the F 0 , F + , and F 2+ centers from top to bottom.The black, bottom line in each panel is the dipole change with no applied field.The red lines are the increasing fields up to 10 MV cm −1 .The dipole moment under no applied field is always the lowest line in each graph, with the dipole moments increasing with applied field strength and the top-most, highest dipole moments stemming from the highest applied field.The insets of the F + and F 2+ centers show the change in dipole moment under no applied field without the change due to the movement of excess charge.vacancy (see Fig. 2).However, it then decreases back to zero at the transition point at frame number 5 on the graph, which is exactly halfway between the two O sites.The dipole moment then decreases to become negative and increases to return to zero once the vacancy has migrated.The dipole moment change is antisymmetric around the transition point.
We note that the F 0 center migration is a complex process involving the displacement of the O ion and redistribution of electron density between the two vacancies.One can obtain an intuitive, albeit qualitative, understanding of why the dipole moment changes in this manner by analyzing the change in electron density along the migration pathway.Figure 2 shows how the electron density of the crystal with the vacancy changes along the migration trajectory with the bulk MgO electron density subtracted at each frame from the initial to the transition state.The numbering of the frames corresponds to the numbers on the x axis of Fig. 1, where frame 0 is the initial vacancy configuration.The electron density around the vacancy can clearly be seen on the left of frame 0. It is symmetric and sits exactly at a lattice site and, therefore, has no dipole moment.As the vacancy migration begins and it moves through frame 1 to 3, the structure and electron density becomes asymmetric and a dipole begins to form around the migrating O atom.This dipole reaches its maximum at frame 3 and then decreases to zero at frame 5.One can qualitatively see in frame 5 that the system is once again symmetric and that all dipoles cancel each other out with the electron density evenly distributed between vacancies.A similar change in dipole moment occurs from frame 5 through to 10, although the dipole now points in the opposite direction.We note that the migration process studied here is adiabatic; that is, the electron transfers gradually from one vacancy site to the other and there is no electron hopping.The observed change in the dipole moment results from the interplay between the ionic displacements and the electron transfer but only the change in the dipole moment has physical meaning.
The F + and F 2+ centers' electron densities along the migration pathway evolve similarly to the F 0 center shown in Fig. 2.However, the net excess charges of these centers play a strong role in their total dipole moment change, as can be seen in Fig. 1.When a charged vacancy is translated by a distance r, the dipole moment changes by q r, as is to be expected for a charged system.As the absolute value of the dipole moment for a charged system depends on the origin chosen, it only makes sense to look at the change in dipole moment along the migration pathway.For convenience, we set the absolute values of the dipole moments of the F + and F 2+ centers so that the lowest value under zero field is also zero.For both the F + and F 2+ centers, the lowest dipole moment at zero field occurs at the initial point along the trajectory and its dipole moment is set to zero.In all cases of charged F center migration, the process remains adiabatic and the displacement r is dominated by that of the O nucleus involved in the migration.
Figure 1 shows that for the F + center, the dipole moment initially increases and then decreases towards the transition point.However, the dipole moments of the initial, transition, and final configurations (frames 0, 5, and 10) of the F + center are no longer equal and the dipole moment is higher in the final configuration than the initial.In the F 2+ center, the oscillatory behavior disappears completely and the dipole moment increases linearly along the migration pathway.Clearly for the F 2+ center, the movement of net charge dominates the change in the dipole moment.
One can effectively separate the total change in dipole moment in the charged systems into two components: one due to the polarization induced by the migration of the defect and the other depending on the movement of charge.The second component changes linearly along the migration pathway with its gradient being equal to the amount of charge moving, making it contribute strongly to the dipole moment changes of the F + and F 2+ centers.By removing this component from the total dipole moment change one is left with the change due to the system's polarization.The remaining dipole moment change, calculated at zero field, is shown in the insets of the dipole moment plots for the F + and F 2+ centers in Fig. 1 and is negligibly different for other field strengths.One can see that this component behaves similarly for all three charge states, displaying an oscillatory behavior.The peak to peak magnitude of dipole moment change due to polarization decreases from the F 0 to F + and then the shape of the dipole moment changes for the F 2+ center.We note that the dipole moment of the F 2+ center at the transition state in the inset of Fig. 1 is not zero, although it should be as the transition state has inversion symmetry.We believe this shift is due to numerical errors of the calculated NEB pathway which results in the transition state being slightly asymmetric; the ions at the transition state lie 4.52 × 10 −5 Å on average off their symmetrical positions.This noise could perhaps be removed by the use of tighter NEB convergence criteria or bigger cells, both of which come at a higher computational expense.
The mechanism of vacancy migration in all three charge states is qualitatively very similar.The major difference between the configurations is the change in net charge.For the F + center, the effect of the net charge on the vacancy is to reduce the magnitude of dipole moment change due to O migration.For the F 2+ center, the increase in charge is high enough that it actually flips the dipole moment resulting in the dipole moment change shown in the inset for the F 2+ center in Fig. 1.These changes in dipole moment are a balance between the change in electronic structure along the migration pathway and the amount of charge in the vacancy.
Figure 3(a) shows how the electric enthalpy of the F 2+ center changes along its migration pathway.One can see that it is symmetric at zero field and that the enthalpy difference between the initial and final configurations increases as a field is applied.Despite the initial and final configurations being the same structurally and chemically, there is an enthalpy difference between them in a field which stems from their dipole moment difference (see Fig. 1).We note that even at fields as high as 10 MV cm −1 the transition point remains very close to the middle of the trajectory for all of the F centers.
Figure 3(b) shows how the barriers for all three types of F centers change as the field gets stronger.The forward barrier is defined as the enthalpy difference between the initial and the transition state, while the reverse barrier is the enthalpy difference between the transition and the final states.The latter corresponds to the vacancy migrating back to its original position against the field.The black lines are the forward barriers and the dashed red lines are the reverse barriers for each F center.The top lines at around 4.5 eV belong to the F 0 center.One can clearly see that its migration barrier changes only by a small amount (0.05 eV) even at a field of 10 MV cm −1 .For the F + and F 2+ centers the barrier changes are bigger and reach 0.23 and 0.35 eV, respectively, at 10 MV cm −1 .Curiously, the barriers for these centers are lower to begin with, so these changes are relatively more pronounced and comprise 6.4 and 15.5%, respectively.These changes can be explained by the dipole moments shown in Fig. 1.Equation ( 4) describes the change in barrier under an applied field in terms of the effective dipole moment μ eff .One can clearly see in Fig. 2 that the dipole moment at the transition and initial states of the F 0 center are almost zero.Hence the effective dipole is almost zero and the barrier change is negligibly small.The effective dipole moments for the F + and F 2+ centers are much higher.The extracted values of μ eff are 0.08, 7.85, and 14.67 Debye, respectively.We note that μ eff for the F 0 center should be exactly zero.However, the numerical methods we have used in these calculations introduce a small error in both the position of the transition state and the calculated dipole moments.Plugging the effective dipole moments into equation ( 4) results in barrier changes of 0.002, 0.18, and 0.34 eV at 10 MV cm −1 .This estimate agrees very well with the calculated barrier changes and gives us confidence that μ eff is indeed a good indicator of the strength of the reduction of the barrier of a defect migration process under an applied field.

B. Interstitial oxygen
The oxygen interstitial in MgO was also modeled in three charge states: neutral, negative, and doubly negative.An extra O atom was initially placed in the center of an MgO channel and 1 Å from a site O atom in all three charge states and the geometry of each structure was then optimized using spin-polarized calculations.The lowest energy configuration in the neutral and singly negatively charged states turns out to be a dumbbell configuration; that is, a lattice O forms a bond with the interstitial O and the center of this nascent bond is located at the original lattice O site.The The migration barriers for interstitial oxygen migration at zero field were calculated using NEB.Using the configurations described above, a band of frames was set up for each charge state using linear interpolation.The lowest energy migration pathways were found to follow the so-called interstitialcy mechanism: The migrating O ion displaces towards the nearest O neighbor regardless of the charge state, displacing the neighbor in the process to take its place and pushing the neighbor into an interstitial position (see also Fig. 5).This is in contrast to an interstitial or void crossing mechanism, which was found to be much higher in energy.The interstitialcy mechanism for O 0 I center follows the 111 direction, while this mechanism for the O − I and O 2− I centers follows the 110 direction.The migration barriers at zero field were found to be 1.06, 0.13, and 0.33 eV for the neutral, negative, and doubly negatively charged states, respectively.The results are shown in Table II along with results from the literature.The values of the O 0 I and O 2− I barriers agree rather well with the literature results, particularly with those obtained in Ref. [43].However, the barrier for the O − I center is higher in comparison to the literature.Electric fields were then applied along the direction of the migration of the interstitial oxygen species.The dipole moment changes along the interstitial migration pathways are shown in Fig. 4. The dipole moments for the charged interstitial configurations were set so that the lowest value under no applied field is zero, as for the F centers.For the O − I and O 2− I centers, this occurs at the end of the trajectory; therefore, the dipole moment at zero field at the end of their trajectories is Figure 5 shows the electron densities of the interstitial centers only at their transition states.The transition state for the O 0 I center is symmetric and has no dipole moment, leading to the graph shown in I are asymmetric at their transition states.This broken symmetry leads to the associated dipole moment changes being qualitatively different from the F centers and the O 0 I center.We now turn to the calculated migration barriers for the O interstitial migration under an applied field.Figure 6(a) shows how the electric enthalpy of the O 2− I center changes during migration.Interestingly, one can see that as the field is applied, the transition state shifts to the right side of the graph.This is most pronounced at 10 MV cm −1 , where the transition state has shifted from frame 5 to 7.5 compared to the zero field.The O 2− I is the only configuration that shows such a pronounced behavior.Transition states for both the O 0 I and O − I centers remain near the same location at all the fields we applied.

IV. DISCUSSION AND CONCLUSIONS
We have investigated the effect of an applied electric field on the migration barrier of defects in MgO using DFT calculations employing the modern theory of polarization implemented in CP2K.We considered only the canonical defects-oxygen vacancies and oxygen interstitial ions-in three charge states.Electric fields of different strengths were applied in the direction of oxygen ion displacements.The considered processes are nontrivial as they include both the ionic motion and the electron density redistribution and are typical of similar processes in other oxides.However, the migration processes in all cases were adiabatic with electron density continuously following the motion of nuclei.It would be of interest to investigate field effects on migration processes involving longer distances between equilibrium configurations, e.g., O vacancy migration in BaO or polaron hopping in amorphous oxides, where electron transfer may prove nonadiabatic.We showed how the dipole moment of the system in different charge states changes as defects migrate.The dipole moment change along the migration pathway can be broken down into two components, one being the change in polarization as the electronic structure changes along the migration pathway while the second comes from the migration of the net excess charge.The effective dipole moment μ eff determines how the applied external field affects the migration barrier.We find that the dipole moment changes due solely to the change in defect's electronic structure along the migration pathway affects the migration barriers very little when a field is applied.However, a reaction leading to charge separation could prove to have a significantly different effective dipole moment than those shown in this paper.The strongest contribution to the migration barrier change upon field application comes from the movement of net excess charge.Naturally, we find that applying a field to the neutral defects affects their barriers the least while the systems with the highest charge, doubly charged in this paper, are affected more significantly.This is most pronounced for the O 2− I which has a barrier change of up to 0.54 eV, which corresponds to a 163% change, upon application of a 10 MV cm −1 field.
Our self-consistent DFT calculations include complex changes in the defect electronic structure both along the migration path and due to the field application.Remarkably, by calculating the dipole moment change μ eff at zero field and multiplying it by the field strength, one can obtain an estimate of the barrier change in excellent agreement with the DFT calculated values.Using this methodology allows one to estimate how the migration barrier will be affected without having to apply an electric field.Recent calculations for the effect of electric field on defect formation in TiO 2 led to similar conclusions [48].
It is illuminating to compare our results and methods to the literature.As mentioned in the Introduction, Cabrera and Mott [7] adjusted the barrier of their interstitial migration process by 1  2 qaF .For the neutral defect, this change will always be zero, while for the F + and F 2+ centers this approach gives a barrier change of 0.17 and 0.33 eV at 10 MV cm −1 .Similarly, for the neutral interstitial center, the barrier change will always be zero, while for the O − I it is 0.12 eV and for the O 2− I it is 0.29 eV, again at 10 MV cm −1 .The estimation based on a point charge interaction with the electric field gives remarkable agreement with experiment at these fields.This is to be expected for migration process of highly ionic species.However, this situation is likely to change for a migration process where the charge is more delocalized or separated.An additional factor that affects these estimations is the position of the transition state.In the Cabrera-Mott estimation, the transition state is always halfway along the trajectory as a field is applied, hence the factor of a half.However, our results for the O 2− I center clearly demonstrate that the position of the transition state can be significantly shifted by the field with respect to the zero-field configuration.The situation may prove even more complicated in more complex and less symmetric systems than those considered here.It would be interesting to investigate whether by modifying Cabrera and Mott's initial estimation by a more accurate reflection of the location of the transition state one could improve its predictive capability across more complex defect processes.Furthermore, it would also be interesting to make full NEB calculations on those systems whose trajectories deviate strongly as a field is applied.
In McPherson's thermochemical E model [8,9], the barrier is reduced by the interaction of the so-called effective dipole moment with the field.This quantity was estimated by McPherson and others using molecular dipole moments [10].It is unclear how one could have extracted an effective dipole moment using these rules for a system with a vacancy.Here, we define our effective dipole moment to be the difference between 3.32 [52] the dipole moments of the initial and transition states.This effective dipole moment can be extracted from a calculation with no field applied and this works rather well for the fields we have looked at.Curiously, all the effective dipole moments extracted here are far lower than values mentioned in the literature.This is true even when comparing our effective dipole moments in charged systems to those that are ascribed to neutral bonds [12,49,50].Another important issue concerns the effect of the relative alignment of the field and ionic displacement on the barrier reduction.In this work, we only considered fields oriented along the direction of ionic displacement where the effect of the field is the strongest.In the simplest case of highly symmetric systems, such as those considered here, the strongest effect occurs when the field is aligned parallel to the direction of migration and there is no effect when the field is in the perpendicular direction.However, the situation is much more complicated when complex materials and migration paths are involved and these issues will be considered in a separate publication.
To summarize, these results shed some light on effects of electric field on the mechanisms of defect processes in dielectrics at the atomistic level.Further applications of the Berry phase method to more complex materials and defect processes will help to assess the applicability of phenomenological models and elucidate linear and nonlinear effects of field application in degradation of microelectronic devices, electrocatalysis, batteries, and other applications.
Field strengths of ±10 MV cm −1 were applied along each of the coordinate axes and component ii of the polarizability tensor was calculated as The isotropic polarizability was then obtained as Table III lists the calculated isotropic polarizabilities as well as the experimental values.It can be seen that both the Berry phase and Hamiltonian methods give reasonable agreement with experiment.Although the calculated polarizabilities could be improved upon by using more extensive basis sets, it is beyond the focus of this paper.

Permittivity of MgO
A field of up to 10 MV cm −1 was applied to bulk MgO and its energies and polarization were evaluated.Two types of calculations were run: one where the atoms remained in fixed positions and another where the atoms were allowed to relax.Figure 7(a) shows how the polarization (defined as the dipole moment per unit volume) changes as a field is applied.The solid line shows the polarization obtained in the calculations with relaxed atomic positions while the dashed line shows that for the fixed atoms.One can see that a linear relationship holds in both cases, as is to be expected for a linear, simple dielectric such as MgO.The permittivity can be extracted from this graph as the change in polarization over the change in the field with the fixed and relaxed atom positions, corresponding to the high frequency and static permittivities, respectively.The calculated static and high-frequency permittivities at 8.76 and 2.65 are slightly lower than the experimental values of 9.34 and 3.01 [53].

25 FIG. 1 .
FIG.1.Dipole moments in Debye along the F center migration pathway.The panels show the F 0 , F + , and F 2+ centers from top to bottom.The black, bottom line in each panel is the dipole change with no applied field.The red lines are the increasing fields up to 10 MV cm −1 .The dipole moment under no applied field is always the lowest line in each graph, with the dipole moments increasing with applied field strength and the top-most, highest dipole moments stemming from the highest applied field.The insets of the F + and F 2+ centers show the change in dipole moment under no applied field without the change due to the movement of excess charge.

FIG. 2 .
FIG.2.Change in electron density of the F 0 center along its migration pathway.The red and light pink spheres are O and Mg, respectively.The blue and red translucent profiles depict an isosurface of the electron density change, plotted with an isovalue of 0.15 e Å−3 .The electron density is calculated as that of the vacancy minus that of the bulk.The number in the upper left corners of the panels correspond to their position along the migration pathway.This is the number found on the x axis of Fig.1.

FIG. 3 .
FIG. 3. (a)Electric enthalpy change along the migration pathway of the F 2+ center as an electric field of up to 10 MV cm −1 is applied.The black line, at the bottom on the left side and top on the right, is the electric enthalpy with no field applied and the red line at the top on the left side and bottom on the right is the highest applied field.(b) Values of the migration barriers of F centers in MgO as an electric field is applied.The black lines are the forward barriers while the dashed red lines are the reverse barriers.Each defect can be seen as a pair of lines with a common origin from 0 MV cm −1 .
O 2− I center, in contrast, occupies the void configuration.The O 0 I center has no occupied states in the band gap; however, there are unoccupied states just below the bottom of the conduction band.The O − I center introduces a singly occupied KS state into the band gap at 1.2 eV below the bottom of the conduction band.The O 2− I center introduces a doubly occupied KS state located at 3.4 eV below the bottom of the conduction band.The extra electrons introduced into the O − I and O 2− I states are strongly localized on the interstitial O.Although the interstitial O displaces a single lattice anion rather strongly in the dumbbell configuration of the O 0 I center, the surrounding ions displace very little.The ionic displacements induced by the O − I and O 2− I centers are much stronger, with the cations surrounding the O 2− I displacing 0.50 Å towards the interstitial.

5 FIG. 4 .
FIG. 4. Dipole moments in Debye along the oxygen interstitial migration pathway.The panels show the O 0 I , O − I , and O 2− I centers from top to bottom.The black, bottom line in each panel is the dipole change with no applied field.The red lines are the increasing fields up to 10 MV cm −1 .As with the F centers, the lowest dipole moment lines are the dipole moments without an applied field, while the highest line corresponds to the highest applied field.The insets of the O − I and O 2− I centers shows the change in dipole moment under no applied field without the movement of excess charge.

FIG. 5 .
FIG. 5. Electron densities of the transition states of oxygen interstitial migration in MgO.From left to right, they are the O 0 I , O − I , and O 2− I

Figure 4 .
However, the O − I and O 2−

Figure 6 (
b) shows barrier changes for all three O interstitial charge states.The forward and reverse barriers are defined in the same way as in Fig. 3. Similar to the F center, the magnitude of barrier change increases as the charge increases.The barrier changes of the O 0 I and O − I centers are 0.02 eV and 0.04 eV at 10 MV cm −1 , corresponding to 1.8% and 34% changes.The barrier changes for the O 2− I center show the strongest asymmetry.The forward barrier at 10 MV cm −1 increases by 0.52 eV while the reverse barrier decreases by 0.28 eV.This asymmetry is due to the transition state of the O 2− I changing position along the migration pathway, as can be seen in Fig. 6(a).Interestingly, due to the smaller barriers to interstitial migration, applying a field of 10 MV cm −1 almost entirely eliminates the reverse barrier for the O 2− I center.We use equation (4) and the dipole moments calculated for the interstitial centers at zero field to estimate the barrier change.The effective dipole moments are 0.52 and 3.92 Debye for the O 0 I and O − I centers, making their barriers to change by 0.012 and 0.09 eV, respectively.However, the nonzero value of μ eff for O 0 I at zero field is due to numerical errors.Similar to the F centers, the barrier change estimations are rather good.Due to the O 2− I center's change in transition state location, the effective dipole moment for the forward and reverse barriers are quite different.We extract the forward and reverse μ eff values as 23.4 and 10.5 Debye, which lead to barrier change estimates of 0.54 and 0.24 eV, respectively.Again, the barrier change estimates from the effective dipole moments are in excellent agreement with the calculated values.

FIG. 6 .
FIG. 6. Migration barriers of oxygen interstitial ions in MgO as an electric field is applied.The black lines are the forward barriers while the red lines are the reverse barriers.Each defect can be seen as a pair of lines.The top lines are the O 0 I , the middle lines are the O − I , and the bottom lines are the O 2− I centers.

Figure 7 (
b) shows the change in energy of the system as the field increases.Two lines are shown in the figure: One corresponds to the calculated DFT energies and the other is calculated as E field = E 0 − μ( F) • F, with E 0 being the energy of the system with no applied field.Both lines agree extremely well with each other, effectively lying one on top of the other.

TABLE I .
Calculated and previously reported values of F-center migration in eV.

TABLE II .
Calculated and previously reported values of interstitial oxygen migration in eV.

TABLE III .
Isotropic polarizabilities of small molecules in Å3 calculated using the finite field technique.The external field is included as either a linear potential change in the Hamiltonian (Ham.) or using the Berry phase approach (B.Phase).The errors of the calculated values compared to experiment are shown as a percentage in the brackets.