Predicting point defect equilibria across oxide hetero-interfaces: model system of ZrO 2 /Cr 2 O 3 †

We present a multi-scale approach to predict equilibrium defect concentrations across oxide/oxide hetero-interfaces. There are three factors that need to be taken into account simultaneously for computing defect redistribution around the hetero-interfaces: the variation of local bonding environment at the interface as epitomized in defect segregation energies, the band oﬀset at the interface, and the equilibration of the chemical potentials of species and electrons via ionic and electronic drift-diﬀusion fluxes. By including these three factors from the level of first principles calculation, we build a continuum model for defect redistribution by concurrent solution of Poisson’s equation for the electrostatic potential and the steady-state equilibrium drift-diﬀusion equation for each defect. This model solves for and preserves the continuity of the electric displacement field throughout the interfacial core zone and the extended space charge zones. We implement this computational framework to a model hetero-interface between the monoclinic zirconium oxide, m-ZrO 2 , and the chromium oxide Cr 2 O 3 . This interface forms upon the oxidation of zirconium alloys containing chromium secondary phase particles. The model explains the beneficial eﬀect of the oxidized Cr particles on the corrosion and hydrogen resistance of Zr alloys. Under oxygen rich conditions, the ZrO 2 /Cr 2 O 3 heterojunction depletes the oxygen vacancies and the sum of electrons and holes in the extended space charge zone in ZrO 2 . This reduces the transport of oxygen and electrons thorough ZrO 2 and slows down the metal oxidation rate. The enrichment of free electrons in the space charge zone is expected to decrease the hydrogen uptake through ZrO 2 . Moreover, our analysis provides a clear anatomy of the components of interfacial electric properties; a zero-Kelvin defect-free contribution and a finite temperature defect contribution. The thorough analytical and numerical treatment presented here quantifies the rich coupling between defect chemistry, thermodynamics and electrostatics which can be used to design and control oxide hetero-interfaces.


I. Introduction
Novel phenomena at oxide hetero-interfaces have been of great interest for both fundamental studies 1 and device applications. 24][5][6] Understanding the redistribution of defects across such interfaces serves as the foundation for studying reaction and diffusion processes in related systems.However, addressing the electronic and ionic structures of the interfacial zone is not a trivial task.The nature of bonds between ions and chemical composition in the oxide crystal could be significantly altered near interfaces. 7hree factors, the band offset, chemical potential difference and segregation energy, rendered in Fig. 1 as (a), (b) and (c), need to be addressed comprehensively and at the same level of theory in order to self-consistently describe the redistribution of charged defects at a hetero-interface.General theory of these three factors have been described in previous works, including defect chemistry at interfaces of ionic materials [8][9][10][11] and semiconductors, [12][13][14] band offset, 15 bulk defect chemistry, 16 and electric field discontinuity. 17Here we provide a summary of the involved theory.First (Fig. 1(a)), before considering any electronic or ionic defect, the discontinuities of valence band DE v and conduction band DE c immediately at the interface need to be calculated.Theoretically, this band discontinuity corresponds to the large electrostatic potential gradient at ideally abrupt defect-free interfaces.This discontinuity is localized in the few atomic layers in the immediate vicinity of the interface and is independent of external perturbations such as dilute doping and applied electric bias.Second (Fig. 1

(b)), charged species
This journal is © the Owner Societies 2017 have different electrochemical potential in equilibrium bulk materials.When two materials are put into contact, electronic and ionic defects start to diffuse, driven by this chemical potential difference.This leads to build-up of local electric field and in turn drives drift flux in the opposite direction.9][20] Third (Fig. 1(c)), atoms near an interface are in a distinct bonding environment compared to atoms in bulk.Defect formation energy will change due to the interfacial lattice relaxation and chemical composition variation.This adds the segregation energy, E seg D q , to defect formation energy of ions in the vicinity of the interface.This region where the segregation energy is non-negligible is termed as the core zone.
Each of these three factors represented in Fig. 1 have been recognized in literature, but were usually treated separately with different modeling methods.[34][35][36] In the rest of this section, we provide more details on these numerical implementations from literature.Predicting band offsets (factor (a)) from density functional theory (DFT) calculations has been of interest since decades ago and a number of theories have been developed. 15One commonly accepted scheme was first proposed by Van de Walle and Martin 21 and applied in a number of studies, 22,23 which aligns the band positions to the electrostatic potential across the heterojunction.These studies consider defectfree interfaces and only evaluate the band offsets.The electrostatic potential, f, built by the drift-diffusion mechanism (factor (b)) can extend in space for several hundreds to thousands of nanometers, and is usually studied with finite-element methods in the context of semiconductor engineering.The drift-diffusion model is welldeveloped for semiconductor heterojunctions, [12][13][14]32,33 where only electrons and holes are considered as the mobile species. In tht case, there is no need to consider segregation energies since free holes and electrons are assumed to be unaffected by the changes in the bonding at the interface.In such model, only the space-charge potential causes changes in holes and electrons at the interface.
When it comes to ionic or mixed electronic-ionic conductors that are used in electrochemical devices, the ionic defect contribution cannot be neglected any more.Ionic and electronic drift-diffusion models based on the chemical potential difference (factor (b)) have been built for mixed electronic-ionic pn homojunctions made of p-type and n-type doped SrTiO 3 . 24Similar continuum level studies have also been implemented for a metal/oxide interface 29 and a gas/oxide interface 26 based on first-principle predictions.However, in these two systems the segregation energy term (factor (c)) is still missing from the consideration of the charged point defect redistribution.In fact in the case of the gas/oxide interface, there is the segregation of defects at the interface within the oxide side of the interface, and there is also the adsorption of charged species from the gas side.While ref. 26 chose to treat the adsorption of charged gas species only, ref. 18 chose to treat segregation of defects within the oxide only.A complete treatment of segregation/adsorption phenomena at gas/oxide interfaces starting from first-principles all the way to continuum level treatment is still missing as well, and the recent work on hydrated BaZrO 3 surface 27 presents a very encouraging attempt in that direction.
The defect segregation energy E seg D q , which have discrete values in each atomic layer close to the interface, can be evaluated using atomistic simulation methods 25,27,28,30,31 and are explicitly accounted for in this work.In fact these segregation energies are responsible for the formation of the core zone at the interface which exhibit very different defect concentrations in comparison with the extended space charge zone.
Attempts have been made to establish interfacial defect chemistry theoretically. 9,11,37These theories consider equilibrium defect reactions at the interface as boundary condition.Counting defect reactions corresponds to a canonical formalism in which defects and in particular electronic defects are tracked particle by particle.Such a canonical formalism, while tractable if one knows a priori the important defect reactions, is not in the spirit of the predictive power of first-principle modeling. 16lternatively, a grand canonical formalism that treats individual point defects and electronic defects such that each of them is in equilibrium with a reservoir is much more compatible with atomistic simulations. 16So far, first-principles-based models have only been developed for symmetric grain boundaries in last few years, [34][35][36] in which the segregation energy is included self-consistently into the drift-diffusion model.A symmetric grain boundary is a simplified scenario where there is no band offset on the two sides of the interface and the electrostatic profiles are symmetric across the boundary.Moreover, the electric field in the core zone was not explicitly included into the numerical treatment of grain boundaries.
In this work, we implement a generalized model of heterointerfaces where no symmetry is assumed and the continuity of the displacement field 17 is adhered to throughout the system by including explicitly the interface core into the solution.To our knowledge this is the first attempt to simultaneously include these three factors listed in Fig. 1 completely from ab initio calculations, without a priori assumption of dominant defect species.On the continuum modeling level, we report for the first time the solution of Poisson's equation throughout the core zone and the space-charge zone by using defect formation energy profiles of step-function behavior, while preserving the continuity of the electric displacement field throughout the entire space including the core zone.Finally, by using this ability we deconvolute the contributions to the interfacial electrostatic profiles arising from the 0 K-perfect and the finite temperature defects effects.This model and the explored concepts are applicable to any type of interface system.
Defect equilibria in bulk provide the boundary conditions needed for the interfacial defect redistribution.In previous work we assessed the defect equilibria in bulk monoclinic ZrO 2 . 38The key idea is to capture the defect structures and formation energies with density functional theory (DFT) calculations and feed it to the thermodynamic model where the total Gibbs free energy of the defective crystal is minimized. 16owever, additional difficulty to tackle is the asymmetry of the hetero-interface.We performed DFT calculations to obtain band alignment and defect formation (segregation) energies at the interface and in bulk of both sides of the heterojunction (described below).Using these results as input, we are able to solve for the electrostatic potential generated by defect redistribution and defect concentrations across the space-charge region and in the interfacial core layer.
We selected ZrO 2 /Cr 2 O 3 hetero-interface as our model system which arises in the utilization of zirconium alloys as nuclear cladding material in light water reactors. 39The oxidation of zirconium alloys by the reactor coolant (light water) could be a severe design-limiting issue under aggressive conditions inside nuclear reactors.Zirconium alloys contain additions of tin, iron, nickel and chromium which tend to precipitate in the form of intermetallics.Upon oxidation, both the host zirconium metal and the intermetallic precipitates form their respective oxides.The net outcome is multiple interfacial zones at the boundary between ZrO 2 and oxides of the precipitates. 40,41xygen diffusion kinetics in ZrO 2 , which is a major consideration in designing oxidation-resistant zirconium alloys, can be significantly altered near the interface between ZrO 2 and the precipitate metal oxide.In this work we focus on the interface between monoclinic ZrO 2 and hexagonal Cr 2 O 3 .Our results show that throughout all temperatures and oxygen partial pressures ZrO 2 exhibits a downward band bending and Cr 2 O 3 exhibits and upward one.The core and space charge zones of ZrO 2 are negatively charged, whereas those of Cr 2 O 3 are positively charged.However, the predominant defect in each zone varies depending on temperature, oxygen partial pressure and the spatial location with respect to the junction indicating a rich coupling between defect chemistry, thermodynamics and electrostatics.In particular in the core and space charge zones of ZrO 2 and under oxygen rich conditions we observe a depletion of doubly charged oxygen vacancies and a similar depletion in the sum of electron and holes.This depletion indicates a reduction in electronic and oxygen transport near Cr 2 O 3 particles dispersed in ZrO 2 .This, in turn, provides an explanation for the corrosion resistance of Zr alloyed with Cr.Finally we discuss how the potential drop, interface dipole, and interface electric field are made out of two components; a 0 K defect-free contribution and a finite temperature defect contribution.

II. Methods
In this section we present a comprehensive methodology to model hetero-interfaces.First, we give a brief review of establishing point defect equilibria in the bulk of an oxide.Second, we describe how to construct a stable interface and obtain band offsets by aligning bulk band positions to a reference potential.
Third, we introduce our continuum model for the equilibrium space-charge zone.Finally, we demonstrate how we explicitly treat defects in the interfacial, core zone.

A. Bulk defect equilibria
We adopt the framework established in our previous work 50 to predict intrinsic point defect equilibria for bulk ZrO 2 and Cr 2 O 3 .The results for bulk monoclinic ZrO 2 were presented in ref. 38 and the detailed analysis for Cr 2 O 3 defect chemistry is to be published elsewhere.The formation energy of a defect D q whose charge is q is computed by: Here E DFT defective is the DFT energy of the defective cell, E DFT perfect is the DFT energy of the perfect cell.m i denotes the chemical potential of each element involved in the cell, and DN i the relative number of each type of atoms with respect to the perfect cell (e.g.DN Zr = À1, DN O = 0 for zirconium vacancies).E VBM is the valence band maximum of the perfect cell.m e is the Fermi level of electrons relative to the valence band maximum.E MP is the first-order Makov-Payne correction 51 for the finitesize supercell errors in charged defect calculations.
The chemical potential of oxygen as a function of temperature T and oxygen partial pressure P O 2 is calculated by: where E O 2 DFT is the DFT-calculated energy of an oxygen molecule.E over is the correction term for O 2 overbinding error created by semilocal functionals, 52 which has a value of 1.22 eV for the parameters used in this work. 38 Fermi level is determined from the charge neutrality condition imposed on point defects D q , electrons n c , and holes p v X where in this work [D q ] is expressed using the Boltzmanndistribution Here n D is the number of possible sites for the defect D per unit formula.Eqn ( 5) is an approximate form of the expression for defect concentration sought by minimizing global Gibbs free energy. 29They are equivalent when E f D,q c k B T. Similarly, the concentrations of electrons and holes, n c and p v are also functions of Fermi level and are computed by integrating the electronic density of states g(E) as follows: Here E CBM is the edge of the conduction band minimum.

B. Interface band alignment
1. Interface construction.We selected oxygen-terminated m-ZrO 2 (001) surface 54 and Cr 2 O 3 (01 À 12) surface 55 to construct a coherent interface.Both surfaces are non-polar.This choice was made by surveying all low-energy surface combinations and minimizing the interface misfit.We built a supercell from these two surfaces with 13 atomic layers of ZrO 2 and 10 atomic layers of Cr 2 O 3 , as shown in Fig. 2(c).The as-built slab model was then fully relaxed with a 4 Â 4 Â 1 k-points grid until the forces on ions were less than 0.02 eV Å À1 , allowing cell shape and volume to change.e comment on the method of constructing heterointerfaces with the following: a stable interface should be built from two low-energy surfaces which have matching lattice parameters.It should be noted that the lowest-energy interface doesn't necessarily arise from two lowest-energy surfaces. 57The most rigorous method is to test all surface combinations and search for the lowest interface energy.Some of these combinations might require an extremely large supercell if the two surfaces have different periodicity.In this work m-ZrO 2 (001)/ Cr 2 O 3 (01 ) is a good combination because both surfaces have rectangular cell and the mismatch is small.Additionally, neither surface is polar, so it avoids electrostatic complications with the periodic slab model. 58To confirm the stability of the constructed interface, we also performed molecular dynamics simulation on the slab model (details in ESI, † Section S1), which shows that the interface retains its phases and structure.
2. Band position.Much interest has been devoted to predicting band offsets from density functional theory for the last few decades and a number of theories have been developed. 15he key of such predictions is to align the band positions of the interface components sought from bulk calculations to a common reference level.Here we adopt the scheme proposed by Van de Walle and Martin 21 and further applied in a number of studies, 23,59 which uses the electrostatic potential across the heterojunction as the reference level.
The electrostatic potential obtained from DFT calculation of the interface under consideration is shown in Fig. 3(a).To get the macroscopic average, electrostatic potential is first averaged over the plane parallel to the interface to get V(z), then averaged over one atomic layer along the z direction. 60Þ ¼ Here o(z) is the filter function which represents the periodicity of the crystalline material.Since the lattice parameters along z axis for the two materials, c 1 and c 2 , are mismatched, we use the double filtering method 60 where o(z) is defined by and is zero otherwise.The valence-band offset is determined by We would like to clarify here the concept of band offset.As shown in Fig. 3(a), the potential in the perfect slab calculation becomes constant within two atomic layers away from the interface.This potential should not be confused with the space-charge layer potential.The equilibrium space-charge zone forms upon redistribution of various electronic and ionic charge carriers.This redistribution is a finite temperature excitation which cannot be captured directly by a zero-Kelvin DFT calculation.The band offset computed here resulted solely from the localized charge density redistribution due to the abrupt change in bonding environment at the interface.In the next section we will superimpose the defect redistribution on the defect-free heterojunction results and solve for the space-charge zone with a continuum-level model.The band discontinuities will serve as an input.The relationship between defect-free model and defective model will be further discussed at the end of this paper.

C. Equilibrium of space-charge zone
From the above calculation, we obtained the electrostatic profile of the perfect interface.Now we use these results for the perfect interface as the neutral background and superimpose the defect equilibria on it.The redistribution of defects will lead to the formation of an extended space-charge layer around the interface and a localized core layer right at the interface.The core layer is discussed in Section D below.In the space-charge layer, defect formation energies are different from bulk by a term of electrostatic potential f scl .Adopting a one-dimensional continuum model, f scl is governed by Poisson's equation e 0 and e are the permittivity of vacuum and relative permittivity of the host material, respectively.r(z) is the total charge concentration including electrons, holes and all types of charged point defects.
Defect concentrations in equilibrated heterogeneous system change with the electrostatic potential.In dilute limit where Boltzmann-like distribution is satisfied, they can be expressed with respect to their bulk value as where Df scl (z) is the electrostatic potential in the space charge layer referenced to its bulk value.Eqn ( 11) is solved with the following boundary conditions z equals 0 and L represent the starting and end points of the modeled heterojunction, respectively.L needs to be large enough to retrieve bulk value far away from the interface.At the interface, the boundary conditions are the continuity of f scl and the continuity of electric displacement field D. In this case, continuity of D is equivalent to the global charge neutrality condition A detailed derivation of this equivalence is in ESI, † Section S2.
The zero point of f scl is arbitrary.Here we take the convention of setting f scl | z=0 to zero.Built-in voltage f built-in ¼ f scl j z¼L Àf scl j z¼0 will be determined upon completing the solution of Poisson equation throughout the space charge zones and the core zones.A detailed description of the boundary conditions used when core zones are added is provided in ESI, † Section S2.
We point out here that using eqn (13) saves the trouble of solving defect concentrations from drift-diffusion equation coupled with Poisson equation for the electrostatic potential.The simplified form in eqn ( 13) is valid when E f D;q ðzÞ ) k B T for all defects.When the formation energies start to approach k B T, their concentrations need to be described by Fermi-Dirac distribution. 29In the latter case no analytical solution for the drift-diffusion equation exists and it must be solved numerically coupled with Poisson's equation.A detailed derivation is provided in ESI, † Section S3.It is worth recalling that Poisson equation in which the charge density stems from defects described by Boltzmann-like expression under equilibrium 13,61 as in this work is sometimes also termed as the Poisson-Boltzmann equation. 8,18,20

D. Core layer
The method presented above is already sufficient to obtain meaningful information of the space-charge region, such as the built-in potential and space-charge length.However, there still exists an important term we have ignored thus far.Defects close to interface are influenced not only by the electrostatic potential, but also by the varying chemical bonding environment in the interfacial layer.This effect is rather localized in the first few atomic layers next to the interface, compared to the space-charge zone which could extend to hundreds of nanometers spatially.To evaluate this term, we define the atomic layer in the vicinity of the interface as the core layer and treat defects in the core layer explicitly.This definition is conventional in studies of grain boundary 34 and dislocations. 62owever, in these existing models the core layer is not solved consistently with the space-charge layer, but treated as a homogeneous zero-electric-field zone.As a consequence of this inconsistency the electric displacement field is not continuous across the interface although it should be. 17In this section we discuss how defects in the core layer are simulated and included in the continuum model.
The width of the core layer is taken to be equal to the lattice parameter in the z direction.We assume that outside the core layer, and except for the electrostatic effect, defect formation energies are not altered by the non-uniformities of the interface and remain in their bulk value.Defect formation energies in the core layer are obtained from DFT calculation to accurately evaluate the effect of structural change near interface.The cell used in Section B is not sufficient to simulate interfacial defects.A larger cell is required to model isolated, noninteracting point defects.We performed DFT calculations for defects in the core layer with a supercell doubled in the x and y direction and 2 Â 2 Â 1 k-points grid.Expression for defect formation energy is the same as eqn (1) except now all values are obtained from the interface calculations.
Following the arguments in Section C, it is straightforward to extend eqn (13) to the core layer and the defect concentration expression in that case reads is the difference in formation energy between defects in bulk and core layer, excluding the contribution of electrostatic potential.Values of E seg Dq for selected defects are listed in Table 2.This value reflects the change of formation energy due to interfacial bonding environment variation and is independent of temperature and oxygen partial pressure.
Poisson equation is solved explicitly inside the core zone subject to the boundary conditions described in Section C and subject to the global charge neutrality as described in eqn (16).A detailed solution scheme is in ESI, † Section S2.It is this explicit solution of Poisson equation inside the core zone that preserves the continuity of the displacement field throughout the heterojunction.Poisson equation is discretized onto a grid with mesh spacing of 1 Å outside core layer and 0.1 Å inside core layer and numerically solved with Newton's method.More details of the numerical implementation are discussed in ESI, † Section S2.

III. Results
We divide this section into four parts.First we analyze the electrostatic profiles resulted by the redistribution of charged defects across the interface.In particular we focus on the nonuniformities of the core layer and explicitly inspect the electric field and charge concentrations in this region.In the second part we take a close look at the spatially-resolved concentration profiles of the dominant defects and discuss how the existence of such interfaces influences the corrosion of zirconium alloys.In the third part we elucidate the general implications that this work provides on studying hetero-interface by making a comparison between the contribution of the defect-free interface and defect redistribution.We present all the results obtained at 1200 K for a range of oxygen pressures.Results obtained for 800 K are in the Supporting Information, presenting the same trends as shown here.Last, we list possible future improvements for the computational method.

A. Electrostatic profiles
In this section we represent electrostatic profiles under a selected thermodynamic condition (1200 K, P O 2 = 10 À5 atm) to demonstrate the features of defect redistribution across the interface.Fig. 4 shows the electrostatic potential, the band bending diagram, and the formation energy of V O across this hetero-interface.Under these conditions, the total built-in potential is f built-in = 1.98 V.At the interface, the band discontinuities from core/ space-charge boundary in ZrO   From Fig. 4 we can see that the electrostatic potential is not smooth at the interface, implying a discontinuity of the electric field.To further examine the interfacial electric field, we zoom on the interfacial zone in Fig. 5. Several features are captured in this plot.First, the electric displacement field D = ee 0 E (Fig. 5(b)) is continuous throughout, as required by Poisson's equation.At ZrO 2 /Cr 2 O 3 boundary, there is a discontinuity in both the charge density and the dielectric constant, this results in a discontinuity in the electric field at z = 0.In addition, at the boundaries core/space charge within each oxide, we assume e is constant and not affected by the electric fields.However, because of the abrupt change in charge density (Fig. 5(c)) at these boundaries, the first derivative of electric field (its divergence) is discontinuous.
Under the above mentioned thermodynamic conditions and because of the downward band bending of ZrO 2 and upward band bending of Cr 2 O 3 , ZrO 2 side is negatively charged and Cr 2 O 3 side is positively charged.The core layer of ZrO 2 is more heavily charged in comparison with its space charge layer by several orders of magnitude.This is mainly caused by the segregation of Cr 0 0 Zr and V 0 0 0 0 Zr to the core layer.On the other hand, the core layer of Cr 2 O 3 is positively charged due to the abundance of valence band holes and Zr Cr .However, contrary to ZrO 2 , the charge density in the core layer of Cr 2 O 3 is less than the charge density in it space charge zone.Details of the defect chemistry will be discussed in the following section.
It is common in literature to treat the core layer as a homogeneous region where the electrostatic potential is flat and defect concentrations have constant values, excluding it from the continuum solution.This is a reasonable approximation because the core is a single atomic layer.4][65][66][67][68] However, this approximation did not capture the magnitude of electric field in the core region.There are some encouraging attempts of treating core layer explicitly in literature.9][70] Mebane et al. took another approach by treating the non-zero segregation energy in the core zone as a boundary condition.This approach bypasses the complication of solving the Poisson's equation with a stepfunction-like formation energy profile and does not account for the large electric field that arises in this core zone. 71,72g. 5

View Article Online
The localized field is important in understanding localized processes at the interface.Recent studies have both measured experimentally 73 and predicted theoretically 74,75 the existence of such large electric fields at a water/oxide interface, an oxide surface, and oxide heterostructures, and their impact on the field-dependent properties such as local bonding interaction, dielectric constant, and 2D electron gas density.When these properties are of interest, accurate prediction of local electric field accounting for both perfect crystal contribution and defect contribution would be necessary.Here we address these two contributions for ZrO 2 /Cr 2 O 3 system.
In this study, we explicitly included the variations of the potential and field inside core layer into the continuum model.Thus, Poisson's equation is adhered to and the displacement field D is continuous throughout the core layer and the spacecharge layer.Based on our results, the potential drop in the core layer is on the order of 10 À3 -10 À1 V, as such defects concentrations do not deviate much from the constantpotential approximation.The advantage of our method, however, is that we can now quantitatively evaluate the electric field at the interface caused by defect redistribution.

B. Defect concentrations
In this section we look at the calculated defect concentrations across ZrO 2 /Cr 2 O 3 interface and deduce how Cr 2 O 3 grown on the Cr secondary phase particles may affect the degradation processes in ZrO 2 passive film formed on zirconium alloys which are used in nuclear reactors.The alloying of Fe and Cr to the zirconium alloys significantly enhances their corrosion resistance. 41Zircaloy-4, for example, contains 0.05-0.15%Cr and 0.12-0.18%Fe.These alloying elements have very low solubility in zirconium metal and precipitate as intermetallic secondary phase particles (SPPs) Zr(Cr,Fe) 2 .These intermetallic precipitates are incorporated in the oxide layer as it starts to grow, and later, these incorporated precipitates get oxidized.It has been observed that Cr 2 O 3 will grow as an outer shell as Zr(Cr,Fe) 2 precipitates are oxidized and thus, form an interface with ZrO 2 host. 40The existence of these interfaces will impact the transport of oxygen, hydrogen, and electrons inside the ZrO 2 layer and thus, affect the corrosion rate of the underlying zirconium alloy.To quantify this interfacial effect, defect concentrations in the space-charge and core layers are needed.
In Fig. 6(a) we show the defect concentrations of undoped bulk m-ZrO 2 at 1200 K and different oxygen partial pressures.ZrO 2 exhibits oxygen deficiency at low P O 2 with positively-charged oxygen vacancies compensating electrons.The dominant type of oxygen vacancy is V O at P O 2 o 10 À15 atm and then switches to V O at higher oxygen partial pressures.ZrO 2 becomes an intrinsic semiconductor at P O 2 4 10 À7 atm where electronic defects (holes and electrons) predominate.The results here differ from our previously reported result for bulk m-ZrO 2 38 because the vibrational free energy is excluded here in order to be consistent with the results for the interface, where vibrational free energy calculations are prohibitively too expensive.When ZrO 2 is brought into contact with Cr 2 O 3 , chromium ions diffuse into ZrO 2 and introduce chromium substitutional defects.We assume the diffusion length L D of defects are longer than the space-charge zone width L scl , thus complete equilibrium is reached.This journal is © the Owner Societies 2017 compensated by oxygen vacancy at low P O 2 (dominant charge states for the vacancy changes from +1 to +2 as P O 2 gets higher).
At much higher P O 2 Cr 0 Zr defects are compensated by holes.Dissolution of chromium into ZrO 2 extends the regime of P O 2 over which the oxygen vacancy predominates.Approaching to the interface in the space-charge layer in Fig. 6(c), positively-charged defects are depleted and negatively-charged defects are accumulated.Cr 0 Zr and electrons predominate in this zone, leading to a negatively-charged space-charge layer.In the core zone, Cr Zr defects are enriched significantly, as a sign of intermixing between ZrO 2 and Cr 2 O 3 .Multiple oxidation states (0 to +3) of chromium exist in the core zone with concentrations much appreciable than the rest of ZrO 2 thanks to their segregation energy in the core zone.Moreover V 0 0 0 0 Zr and O 0 0 i start to populate the core zone of ZrO 2 at high P O 2 , showing that the segregation energies in the core zone not only affect extrinsic defects, but intrinsic ones, as well.By tracing the dominant defect in the core zone as a function of P O 2 we observe a switch from doubly charged Cr 0 0 Zr when P O 2 o 10 À7 atm to quadruply charged V 0 0 0 0 Zr when P O 2 4 10 À7 atm.This shows the strong coupling between thermodynamic conditions, defect chemistry and the electrostatics of oxide heterojunctions with implications for the possibility of experimentally controlling these junctions via controlling growth conditions.
Similar analysis can be done on the Cr 2 O 3 side, as shown in Fig. 7.In undoped bulk Cr 2 O 3 , there is no significant concentration of ionic defects (Fig. 7(a)) and only electronic defects predominate.This is consistent with the conclusions of the conductivity experiments conducted on undoped Cr 2 O 3 . 76Far away from the interface such that local charge neutrality is maintained but within the diffusion length of Zr defects, Zr Cr and Zr Â Cr (Zr in +4 and +3 oxidation states, respectively) start to show up with appreciable concentrations.Getting closer to the interface, the bands of Cr 2 O 3 bend up, and the concentration of the positively-charged Zr Cr is increased.Inside the core layer, oxygen vacancies tend to segregate but the predomination in this zone remains for Zr Cr and holes.The core zone of Cr 2 O 3 maintains the predomination of singly positively charged defect throughout all oxygen partial pressures considered here at 1200 K but with a switching from Zr Cr to holes at P O 2 = 10 À18 atm.This again brings up the rich coupling between thermodynamics, chemistry and electrostatics in oxide heterointerfaces.
For symmetric grain boundaries in oxides it is customary to postulate that the core zone has a charge opposing and compensating that of the space charge zone. 34,36The symmetry of the grain boundary justifies this postulate.At a heterointerface, however, the symmetry is broken and in the case of the ZrO 2 /Cr 2 O 3 system studied here we found that the core and the space charge have similar charge (both are negative in the case of ZrO 2 and both are positive in the case of Cr 2 O 3 ).When the space charge and the core zone possess similar charges, the electric compensation occurs on the other side of the junction (the other oxide).In fact this charge compensation scenario in heterojunctions suggests the possibility of the presence of a thermodynamic point of zero charge in analogue with what happens when an oxide is immersed in an aqueous solution. 77hat is a point in which the two sides of the junction are uncharged.We systematically varied temperature (between 400 K to 1400 K) and oxygen partial pressures (between 10 À30 atm to 1 atm) to search for such a neutral point in the ZrO 2 / Cr 2 O 3 system but did not observe it.Its absence in this system does not preclude its existence in other oxide hetero-interfaces.We turn now to the implications of this analysis on zirconium alloys in nuclear reactors.The growth rate of zirconium oxide layer is thought to be limited by either electron transport or oxygen diffusion in ZrO 2 . 78Fig. 8 shows the concentrations of free electrons, holes and V O throughout the heterojunction at various oxygen partial pressures at 1200 K.By looking at the concentration distribution of oxygen-related defects and electronic defects, we can analyze how electron transport and oxygen diffusion are changed in the space-charge layer of ZrO 2 compared to its bulk due to the change in the concentration of these defects in the space-charge layer.In Fig. 8(a), because the bands bend down on the ZrO 2 side, the concentration of electrons is increased in the space charge layer compared to bulk, and subsequently holes are decreased as well in this layer.This enhancement of electron concentration has a favorable effect on preventing hydrogen, by the corrosion reaction, from getting into ZrO 2 . 38,79In addition, the sum of electron and hole concentrations is plotted in Fig. 8(b).If we assume that the mobilities of electrons and holes are on the same order of magnitude, this quantity should give us an estimation of the change in total electronic conductivity in the space-charge layer.Fig. 8(b) shows that when P O 2 o 10 À12 atm, the space charge layer of ZrO 2 is in its n-type regime.In this case, the total electronic conductivity in the space-charge zone is increased.On the other hand, when ZrO 2 becomes p-type at higher P O 2 , the total electronic conductivity is decreased in the space-charge layer in comparison to bulk due the downward band bending.If electron transport is the rate limiting step of the corrosion reaction, then the oxidized Cr precipitates lying in the oxygen rich, and hence the p-type part of ZrO 2 , closer to the interface with water coolant in the nuclear reactor, have a desirable of effect of slowing down the electronic conductivity and consequently reducing the corrosion rate.In fact, the ZrO 2 oxide layer grown is exposed to a gradient in oxygen chemical potential between an oxygen poor zone at the interface with the underlying Zr alloy and an oxygen rich zone at the interface with water.As such, even if the electronic conductivity is enhanced around the precipitates in the oxygen-poor part of ZrO 2 (n-type), it is inevitable that this conductivity is reduced around the precipitates in the oxygen-rich part (p-type).
Cr 2 O 3 has low oxygen diffusivity (B10 À16 cm 2 s À1 at 1200 K, 0.1 atm 80 ) compared to monoclinic ZrO 2 (B10 À13 cm 2 s À1 at 1200 K, 0.4 atm 81 ), which indicate that the in-precipitate oxygen diffusion should be negligible, especially that the concentration of V O in Cr 2 O 3 is several orders of magnitude lower than in ZrO 2 as shown in Fig. 8(c).Since the precipitates are discrete particles that do not form a continuous path through which oxygen ions can diffuse, we focus here on oxygen diffusion in the matrix, ZrO 2 .The extended space charge zones around the oxidized Cr precipitates can form a significant fraction of the matrix with possibilities for space charge percolation.Fig. 8(c) shows that oxygen vacancies are strongly depleted throughout the space-charge layer and in particular in the core zone on ZrO 2 side, which indicates that oxygen diffusion in ZrO 2 is suppressed around the Cr oxidized precipitates.This observation suggests that if oxygen diffusion is the rate limiting step of the corrosion reaction, then the presence of oxidized Cr precipitates has also a desirable effect of slowing down oxygen diffusion in the extended space charge zones around them.While the conclusions here are based on the defect profiles evaluated at 1200 K, they are still applicable to lower temperatures closer to the actual reactor operation temperature (see ESI, † Section S4 for calculated profiles at 800 K).Modeling the space charge zone at lower temperatures becomes more challenging in the dilute limit because the difference in the spatial extent of the core layer and the space charge layer becomes increasingly large.It is worth noting that the above analysis considers chromium as the only alloying element in ZrO 2 , and Zr as the only alloying element in Cr 2 O 3 .In actual zirconium alloys used by the nuclear industry, other alloying elements dissolving into ZrO 2 can also impact its defect chemistry and corrosion resistance. 38,82Fe alone (zircaloy-4) or Ni and Fe (zircaloy-2) are key alloying elements in Zr alloys that contain Cr as well.In addition, hydrogen from the aqueous environment has a nonnegligible contribution to the defect chemistry.Recently we have established the defect equilibria in M-ZrO 2 when it is co-doped with hydrogen and 3d transition metals. 38Future work will build on the analysis presented here and our previous work, to quantify the combined effect of these doping elements on the interfacial defect distributions and electrical properties in oxide hetero-interfaces including ZrO 2 .
C. Potential drop, electric field and interface dipole: finite temperature vs. 0 K In this study, we have obtained two sets of solutions of Poisson's equation: (1) from first-principle calculations on the defect-free interface at 0 K, and (2) from continuum level modeling of defect equilibria, where ionic and electronic defects redistribute due to chemical and electrical forces and reach thermal equilibrium at finite-temperature.Therefore, our multi-scale method allows us to separate the contributions of the perfect crystal (0 K) and the defect formation and segregation (finite T) on the potential drop, the electric field, and the interface dipole obtained from solving Poisson's equation.These three quantities are all experimentally measurable 83 and are important for understanding the physical and chemical properties of hetero-structures. 84,85In this section we compare the origins and consequences of these two contributions.The summary of the comparison is shown in Table 3.
Fig. 9 shows a comparison between the macroscopic averages of the solution of Poisson's equation of the defectfree interface and the solution of Poisson's equation for the charged defects at T = 1200 K and P O 2 = 10 À5 atm.The potential drop across perfect interface DV, as shown in Fig. 9(a), is À3.8 V.This potential drop results from the forming of polar bonds at the interface.Electrons distribute unevenly on the two sides of the interface and a localized interface dipole is built.Notice from Fig. 9 that the electrostatic potential becomes flat, and electric field and charge concentration of the defect-free interface become zero already at one atomic layer away from the interface.Because of the highly localized nature of this potential drop, electric field and charge concentrations have very high values right at the defect-free interface.In the system under consideration, the maximum electric field is 4.9 Â 10 9 V m À1 and the interface dipole density is 3.0 Â 10 11 e m À1 .
In contrast, the effect of the defects is extended in space (hence the denotation space charge).As shown in Fig. 9, the length of space-charge layer is several hundred nanometers.The potential drop is 2.0 V across the space-charge layer, with only 0.2 V drop in the core zone.The maximum electric field is À4.3 Â 10 8 V m À1 and the dipole density is 1.2 Â 10 9 e m À1 , which are order-of-magnitude lower than the perfect interface contribution.The two contributions are, however, mathematically additive as required by the superposition principle. 86hese two effects coexist in every interface between two semiconductors (or insulators) and can both influence the interface properties.However, because of the difference in spatial scales of these two effects, either could be dominating at different spatial positions around the interface.When electrostatic profiles are discussed in literature, these two effects are seldom distinguished or contrasted.We believe that the highly localized nature of the field and the dipole arising from the defect-free interface gives rise to the catalytic activity of exposed or barely hidden interfaces. 87It can also influence the transport processes right at the interface (parallel or across).On the other hand, the field and dipole generated by defects at finite temperatures influence the transport processes over an extended region around the interface.

D. Possible future improvements for the computational method
There are several improvements that could be made to further refine our predictions.First, Fermi-Dirac distribution could be used instead of Boltzmann distribution in calculating defect concentrations.This requires including explicit solution of the drift-diffusion equations for all defects.This is particularly important when the formation energies start approaching k B T. Second, we did not consider finite temperature effects in this work.To accurately evaluate defect concentrations at various temperatures, formation energy should be replaced by formation free energy, by adding the contribution of phonons and the electronic entropy. 50Third, it has been reported in literature that the formation energy of a point defect with respect to its distance to the interface should be an oscillating function, which gradually approaches its bulk value after a few atomic layers. 25,30,88In this work we only treat the first atomic layer as the core layer, which is the computationally affordable approximation for evaluating segregation energies for all types of defects in the system at this time.In principle we could obtain spatiallyresolved defect formation energy profiles from first-principles calculations by inserting defects in each atomic layer and add as many core layers as needed into the continuum model.Fourth, we treated electrons and holes as free carriers in bulk, space-charge layer and the core zone and thus their concentrations were calculated by applying Fermi-Dirac distribution to the electronic density of states.However, it is known that these electronic defects can self-trap in polarizable oxides forming small polarons which is indeed the case for holes in monoclinic ZrO 2 . 89Small polarons should be treated as point defects when calculating their concentrations.The self-trapping energy for holes in monoclinic ZrO 2 is only 0.13 eV 89 indicating that the results presented here remain qualitatively correct even if self-trapping is accounted for.However, it is of interest to explore the differences in self-trapping energies between bulk and the core region in future work.Finally, more accurate electronic structure methods, such as hybrid functionals 22 and G 0 W 0 , 90 might be needed to correctly predict band alignment at the interface and electronic structures of the defects.Additional modeling difficulties could arise while going to more complicated oxide systems, but the framework provided here addresses defect chemistry at hetero-interface in a generalized way and could potentially be applied to other oxide interfaces and surfaces.Lastly, results for temperatures lower than 800 K are not presented because the space charge layer length becomes too large (B100 microns at 600 K), requiring to use a very coarse mesh to do the calculations, and also because these really large space charge zones are beyond the thickness of the corrosion passive films.In reality, the space charge length could be significantly reduced by other alloying elements that exist in the system, which is not considered in this work.

IV. Conclusion
In this paper, we implemented a thorough framework to predict the equilibrium defect concentrations across a hetero-interface starting from first principles-based calculations.Our method takes into account both long-range Coulomb interaction in the space charge layer and the localized effects of band discontinuity and bonding environment variation, allowing for a comprehensive description of defect redistribution across hetero-interfaces.
The addition of core layers explicitly into the solution provides a closer view on the electric field and defect chemistry in the vicinity of the interface.Our model features a consistent treatment of Poisson's equation throughout the core and space charge layers preserving the continuity of the displacement field, whereas prior studies treat the core layer as a field-free homogeneous zone.Our results suggest that the variation in the formation energy of defects can be up to several eVs in the core region, causing a localized change in the electrostatic potential on the order of 0.1 V inside the core zone.Consequently, a giant electric field arises in this core zone and hence it deserves explicit modeling when studying oxide interfaces.
Our selected model interface, ZrO 2 /Cr 2 O 3 , is of interest to understand and improve the performance of zirconium alloys as nuclear cladding materials.This interface arises from the oxidation of the host Zr and the Cr secondary phase particles embedded inside the host.Here we clarified the beneficial effect of this interface in improving the corrosion and hydrogen resistance of Zr alloys.In particular, in the extended space charge layer around these interfaces and under oxygen rich conditions, we predict a decrease in the total electronic defect concentrations and in the oxygen vacancy concentration.Both of these results should decrease the metal oxidation kinetics.In spite of the reduction in the sum of electronic charge carriers under oxygen rich conditions, we found an increase of free electron concentration in the ZrO 2 space-charge zone compared to its bulk for all oxygen partial pressures.This effect should promote proton reduction and hydrogen evolution at the surface, and thus, suppress hydrogen pickup into the alloy.Both effects are beneficial towards protecting the integrity of zirconium alloys.The model and the analyses presented here provide a utility to understand the effect of alloying elements on the corrosion resistance of zirconium alloys when these alloying elements tend to precipitate as secondary phase particles and oxidize within the ZrO 2 passive layer.
This systematic study also facilitates a comparison of two contributions to the interfacial electrostatic properties; a 0 K contribution originating from the formation of a perfectly abrupt interface, and a finite temperature excitation due to the equilibration of electronic and ionic defects.Finally, we discussed several improvements to enhance the predictive power of the herein presented model.Although here we presented the application of the model to the problem of zirconium alloy corrosion resistance, we believe that it is broadly applicable to oxide interfaces arising in other important technologies such as lithium batteries, fuel cells, and computer memories.

Fig. 1
Fig. 1 Schematic illustration of the comprehensive framework adopted in this work for predicting point defect equilibria across hetero-interfaces (ref.8-17 in the Introduction section).(a) Band offset of the ideal interface.(b) Any species has different chemical potentials, m 1 and m 2 , when bulk defect equilibria is reached on each side of the interface.This difference drives diffusion flux as the two materials are put into contact and lead to the formation of a space-charge zone.Upon equilibrium, Fermi level and the chemical potential of each element should be constant throughout the heterojunction, and the electrochemical potential of each type of ionic defect should be constant on the oxide side where it exists.(c) Defect segregation or depletion due to interfacial structural and compositional change.After assessing these three factors at the atomistic level by density functional theory (DFT), these are provided as input into the steady state drift-diffusion model.In this work, the steady state drift diffusion model is solved in a self-consistent way to generate (d) the electrostatic profiles in the space-charge and core zones, and (e) the defect redistribution profiles across the interface with explicit consideration for the interfacial core zone.When equilibrium is reached, the diffusion flux driven by concentration gradient rc and drift flux driven by electrostatic potential gradient rf cancel out.Core regions extend for only several atomic layers in the immediate vicinity of the interface and are drawn with enlarged scale and light colors in (d) and (e).
chemical potential for O 2 gas at reference pressure P 0 = 0.1 MPa.53 k B is the Boltzmann constant.Once m O is known, the chemical potential of cation elements are calculated by m M T; P O 2 where M x O y represents the monoclinic phase of ZrO 2 for zirconium and Cr 2 O 3 for chromium.E MxOy DFT is the DFT-computed energy of the oxide material.

Fig. 2
Fig. 2 (a) m-ZrO 2 (001) and (b) Cr 2 O 3 (01 À 12) surface cross sections and (c) relaxed interface.Oxygen ions are drawn small and red, Cr ions are medium-sized and blue, Zr ions are largest and green.This figure is generated with visualization software VESTA.56

E ZrO 2 v
and E Cr 2 O 3 v are the valence band energies with respect to the local electrostatic potential obtained from two separate bulk calculations.D % V ZrO2-Cr2O3 is the potential difference of the bulk-like regions of the interface, where flat potentials are reached.The aligned band positions are shown in Fig. 3(b).

Table 1 Fig. 3
Fig. 3 (a) Planar (solid black line) and macroscopic average (dashed red line) of electrostatic potential across the ZrO 2 /Cr 2 O 3 interface obtained from DFT calculation.D % V ZrO2-Cr2O3 = 3.77 V is the potential difference on both sides of the interface where flat potentials are reached.(b) Aligned band positions.Band gaps of ZrO 2 and Cr 2 O 3 from DFT calculations are 3.49 eV and 3.38 eV, respectively.
2 to core/space-charge boundary in Cr 2 O 3 are DE v = 1.19 eV and DE c = 1.08 eV, as sought from the method described in Section IIB.Bands bend down on the ZrO 2 side and up on the Cr 2 O 3 side.As such, the positively charged defect species are depleted and negatively charged defects accumulate in both the space-charge and core zones of ZrO 2 .The converse occurs on Cr 2 O 3 side.The presence of the core layer adds a potential drop of 0.05 V on the ZrO 2 side and 0.10 V on Cr 2 O 3 side.Therefore, right at the interface (core/core) both band discontinuities are increased by 0.05 eV due to defect segregation in the core layer.We plot the formation energy of V O in Fig. 4(c) with the circle and square explicitly showing values in the core layer on ZrO 2 and Cr 2 O 3 side, respectively.Outside the core layer (in the space charge zone), E f V O differs from bulk value by twice the electrostatic potential change as it

Fig. 4
Fig. 4 (a) Electrostatic potential, (b) band edges positions for valence band E v and conduction band E c , and the Fermi level E f , (c) formation energy of doubly charged oxygen vacancy calculated at 1200 K and P O 2 = 10 À5 atm.The interface is placed at z = 0 with ZrO 2 on the left side and Cr 2 O 3 on the right side.In (c), the blue circle (ZrO 2 side) and red square (Cr 2 O 3 side) show explicitly the formation energy in core layers on each side.
Fig. 5 (a) Electrostatic potential f, (b) the magnitude of the electric field E and electric displacement field D and (c) charge concentration r close to the interface.Regions with green and blue background color represent core layer on ZrO 2 and Cr 2 O 3 side, respectively.

Fig. 6
Fig. 6 Defect concentrations at 1200 K and different oxygen partial pressures in ZrO 2 (a) in bulk, (b) far from the space-charge zone but within the reach of Cr defects, (c) outside the core layer within the space-charge zone, and (d) in the core layer.

Fig. 7
Fig. 7 Defect concentrations in Cr 2 O 3 at 1200 K and different oxygen partial pressures (a) in bulk, (b) far from the space-charge zone but within the reach of Zr defects, (c) outside the core layer within the space-charge zone, and (d) in the core layer.

Fig. 8
Fig. 8 Spatially-resolved concentrations of (a) electrons, (b) sum of electrons and holes and (c) V O as a function of oxygen partial pressure at T = 1200 K. Circles and squares in panel (c) represent the concentrations in the core layer on the ZrO 2 and Cr 2 O 3 sides, respectively.

Table 2
Defect formation energy differences (in eV) between bulk and core layer OÀ2.84

Table 3
Potential drop across the interface, maximum electric field and dipole moment density calculated for the two sets of solution of Poisson's equation: the perfect interface from 0 K DFT calculations, and the defect contribution from the continuum model.The defective interface results are calculated at T = 1200 K and P O 2 = 10 À5 atm.The dipole density p is evaluated using p ¼