First Principles Analysis of the Stability and Diffusion of Oxygen Vacancies in Metal Oxides

0031-9007= Oxygen vacancies in metal oxides are known to determine their chemistry and physics. The properties of neutral oxygen vacancies in metal oxides of increasing complexity (MgO, CaO, -Al2O3, and ZnO) have been studied using density functional theory. Vacancy formation energies, vacancy-vacancy interaction, and the barriers for vacancy migration are determined and rationalized in terms of the ionicity, the Madelung potential, and lattice relaxation. It is found that the Madelung potential controls the oxygen vacancy properties of highly ionic oxides whereas a more complex picture arises for covalent ZnO.

Metal oxides have attracted a renewed interest because of their relevance to applications as different as high-T c superconductivity, heterogeneous catalysis, photoelectrolysis, and biocompatibility [1].Nevertheless, many of these properties do only appear when the material is either doped or nonstoichiometric.In particular, the role of oxygen vacancies in metal oxides is crucial in the performance of solid state fuel cells and controls the degradation of damaged insulators.The fundamental understanding of the properties of oxygen vacancies is very scattered, having been deeply investigated only for MgO [2,3]; a unified view of their properties on different oxides is lacking.Among the fundamental parameters controlling the physics and chemistry of oxygen vacancies, their stability and diffusion have a prominent role; both are relevant to the description of vacancy agglomeration [4], dispersion of metals on oxides [5], or the release of instabilities created by polar terminations [6].Depending on the ambient gases, some oxides lose part of their oxygen content [7,8], while others are very stable even against very harsh conditions.In the latter case, vacancies can be generated under neutron radiation or highly energetic ions, cutting procedures, or synthesizing these solids under nonstoichiometric conditions.
The aim of the present study is to provide a systematic study of the stability and mobility of oxygen vacancies in simple metal oxides and to unravel the main physical mechanisms governing these processes.Previous simulations and experiments have shown that ionic radius and the charge of the ions are relevant quantities to investigate impurity stability and diffusion in ionic systems [9].A series of oxides of increasing chemical complexity-MgO, CaO, -Al 2 O 3 , and ZnO -has been selected and the energy needed to generate an oxygen vacancy, the relative stability of different structures corresponding to divacancies, and the diffusion barriers obtained from density functional theory (DFT) calculations.To avoid long-range interactions between neutral vacancies [10], large supercells have been considered corresponding to 1% O-vacancy concentration.Those are 3 3 3 for the alkaline-earth oxides, 4 4 3 for ZnO, and 3 3 1 for -Al 2 O 3 .In all cases the minimum distance between vacancies exceeds 12 A so the interaction among vacancies is assumed to be negligible.This assumption is supported by recent STM measurements for the rutile (110) surface [10] which show that, in the worst case, at 12 A the repulsion between oxygen vacancies is less than 0.02 eV.Because of the higher coordination of the atoms in the bulk, screening effects are larger than on the surface and more likely prevent long-range interactions.Test calculations varying the number of atoms allowed to relax around a given vacancy have been performed in order to reach convergence in the geometric relaxations and the total energy.The exchange and correlation contribution to the total energy was obtained with the Perdew-Wang 91 functional [11].The frozen-core-like projector augmented-wave pseudopotentials have been used to represent the inner electrons [12].One-electron states have been expanded in a plane wave basis with a 415 eV kinetic energy cutoff, and a minimum 0:15 A ÿ1 k-point density in each direction.All calculations have been carried out by with the Vienna ab initio simulation package (VASP) code [13][14][15].With this approach the unit cell parameters and internal coordinates for our set of oxides, represented in Fig. 1, are in good agreement with experimental results, the average error being 1%.An important feature of oxide properties is the O charge.This has been evaluated by integrating the electronic density  in the corresponding topological O basin of the regular structure [16].MgO, CaO, and -Al 2 O 3 are highly ionic, the charge of the O atoms being close to ÿ2e, while only ÿ0:4e is found for the oxygen charge in ZnO.
The formation energies for isolated oxygen vacancies have been determined with respect to the ground state of the O atom and show a strong dependence on the nature of the oxide (Table I).For ionic oxides, the formation energy is almost constant, 10 eV, similar to previous estimates for MgO [17] and -Al 2 O 3 [18].The higher covalent character of the ZnO is manifested in lower formation energy per O vacancy, 7:0 eV, which agrees with LDA estimates [19,20].The presence of the vacancy distorts the neighboring lattice region to different extents.For MgO and CaO, lattice relaxation around the oxygen vacancy is small; neighboring vacancy cations move outwards by 0.01 and 0:04 A for Mg and Ca, respectively, while second neighbors are immobile.The pair of electrons left in the vacancy is less strongly bound for CaO and more delocalized than for MgO due to the lower Madelung field, thus resulting in a somewhat larger relaxation.The case of -Al 2 O 3 is slightly different since both Al and O ions are found as nearest neighbors [18].An asymmetric distortion is observed for the cations; two Al move inwards by 0:10 A while the other Al relax outwards by 0:02 A, and the O atoms relax by a similar amount towards the vacancy.Nevertheless, the next shells of atoms surrounding the vacancy do not relax at all as found for the other highly ionic oxides.For oxygen vacancies in ZnO, the convergence with the number of optimized shells of atoms is slow; a minimum of 16 atoms should be included in the optimization.However, the formation energy is 0.04 eV higher than the one corresponding to optimization of the next shell of atoms (41 centers).In this case, the defect nearest neighbors move inwards by 9% (0:18 A) leading to the formation of a Zn 4 colloid, where the Zn-Zn distance is 3:25 A and induces a noticeable strain in the next Zn-O bonds, reduced when O atoms can follow the Zn displacements.
The electronic structure of single isolated vacancies depends substantially on the nature of the oxide considered.Impurity states in the density of states arising from the presence of a vacancy can be seen for MgO, CaO, and -Al 2 O 3 ; the shape of these levels indicates that the electrons are trapped in the vacancy left upon removal of one O atom.In ZnO, the presence of vacancies leads to a reduction in the band gap without evidence of individual impurity states, therefore no significant amount of charge is left at the vacancy site.For a 1% oxygen vacancy concentration, the cohesive energy, U, of the MgO, CaO, and -Al 2 O 3 highly ionic oxides does not significantly differ from that of the perfect crystal.Therefore, the vacancy formation energy, E f , can be computed by combining two Born-Haber cycles, corresponding to the regular system and to the O-deficient oxide (see Fig. 2).Since U can be assumed to be similar for both systems, the only differential contributions to E f are those corresponding to the interelectronic repulsion, E e-e , due to the electron pair confinement in the cavity left upon oxygen removal, and the extra oxygen electron affinities needed to form the regular bulk.These two terms do only involve properties related to oxygen (electron affinity and ionic radius) and that is reason behind the similarity of the E f values in ionic oxides.
Experimental barriers for oxygen vacancy diffusion in oxides are very difficult to measure because of the high energies needed.As a consequence, oxygen diffusion processes can be masked by other mechanisms with lower energy requirements.In fact, experimental results evidence a large dispersion [21][22][23].Energy barriers for vacancy diffusion, E dif a , can be estimated from the corresponding saddle point in the potential energy surface.The saddle point location can be performed through a constrained search applying symmetry considerations.For -Al 2 O 3 , several possible paths exist [see inset (B) in Fig. 3] [24] and a systematic search has been carried out.To validate the present approach we use KCl, isostructural to MgO, for which accurate measurements exist because anion vacancy diffusion occurs at lower temperatures.Our value is E dif a 1:2 eV, close to the experimental value 1:0 0:2 eV [25].For MgO, E dif a 4:2 eV, whereas values as different as 2.2, 3.4, and 5.6 eV have been proposed from experiments [21,22].Recent measurements lead Popov et al. [21] to argue that values lower than 3.4 eVaccount for competing processes.For -Al 2 O 3 225502-2 at 1673 K the experimental value E dif a 4:5 eV [23] is slightly larger than the lowest of the calculated ones, 3.7 eV.For ZnO, E dif a 2:1 eV, indicating a larger mobility of O atoms in this covalent system compared to the highly ionic oxides due in part to the smaller covalent radius of O compared to the O 2ÿ ionic one.
To rationalize the trend in the diffusion barrier for oxygen vacancies in highly ionic oxides, the calculated E dif a is plotted as a function of the Madelung potential at the saddle point,V Mad;TS , obtained from the Ewald summation assuming formal charges for all ions in the crystal (see Fig. 3).In that way we introduce in a compact way the dependence with the ionic radius and the charge suggested in Ref. [9].This approach is not suitable for ZnO where other mechanisms related to the directionality of the bond apply.A symmetric charge distribution based on the density at the transition state [shown in inset (A) of Fig. 3] is taken for the charges in the initial and final sites of the migrating anion (ÿ1e) while a charge of ÿ2e is kept at the saddle point representing the moving O anion.The correlation between the Madelung potential and the calculated E dif a suggests that, for ionic crystals, it is possible to predict E dif a values using the structure of the perfect crystal and an ansatz for the transition state of the migrating oxygen.If several paths are possible, i.e., -Al 2 O 3 , the most favorable diffusion path is given by the structure with the lowest Madelung potential.
Vacancy diffusion can generate dimers and larger aggregates.Experimental information on the existence of dimers is very scarce, although they have been proposed to be present in some MgO preparations [26].The relative stability of the possible divacancies was analyzed by deriving pair potentials for the vacancy-vacancy interaction: where E VÿV is the average formation energy for the divacancy structure, E 0 is the isolated vacancy energy formation, and n d is the number of neighboring vacancies at a distance d from the reference vacancy.The u VÿV pair potentials are reported in Fig. 4 as a function of the distance between the vacancies in the ideal crystal and numbered as I-VII (I stands for the shorter distance).These values have been obtained employing a 2% defect concentration in the supercells.For all three ionic oxides, short distance aggregates are energetically favorable with respect to isolated vacancies.Stabilization is based on the formation of the bonding and antibonding states which reduces the mutual repulsion between the electron pairs while it maximizes the electrostatic stabilization with the neighboring cations.MgO type I shows a larger attraction than CaO due to the shorter distance between the atoms which favors the overlap between the electronic states of the vacancies.In the case of -Al 2 O 3 , divacancies II and IV are more stable than those with shorter distances, I and III, respectively.This is because types II and IV experience a larger Madelung potential and, therefore, a larger stabilization.Again, the Madelung potential at the site is calculated by considering a ÿ1e charge situated at the positions of the  For ZnO, vacancy aggregates are repulsive in nature.The growth of larger metallic nuclides will not occur due to both a sizable diffusion barrier and the thermodynamically unstable dimer state.However, in Fig. 4, minimum repulsion structures for ZnO are found at 5 and 7 A. The reason for this long-range stabilization lies in the changes in geometry for nearest neighbors.For ZnO, type I dimers show a conflict in relaxations where one of the Zn atoms can only relax inwards by 3% compared to 9% in the isolated vacancy (see inset in Fig. 4).Stress diminishes the ability of this Zn center towards bonding and results in a more repulsive interaction than the situation with isolated O vacancies.For other structures, like type III configuration with vacancies align on the (0001) direction, the relaxation spheres around the vacancies do not interact at first neighbors and therefore repulsion is smaller than in type I.This direction (0001) can be regarded as a crystallographic shear (CS) plane where O vacancies can be accumulated without a significant energy cost.Preferential directions, with minimum repulsion, have also been observed for rutile (110) surfaces through STM [10], and the pair potentials obtained there show a similar behavior to those reported in Fig. 4 for ZnO.Similar relaxation processes can provide a reason for the stability of the partial oxides known as Magneli phases for Mo, V, and Ti oxides.
In summary, the highly ionic or more covalent character of the oxides defines the properties of the oxygen vacancies.For highly ionic oxide energy formation, migration and aggregation can be described by analyzing the charge distribution and the corresponding Madelung energies or potentials only.For the ZnO, more covalent material, the large lattice relaxation upon O-vacancy formation and the presence of directed bonds leads to the nearest neighbor dimer instability and to a rather low diffusion barrier.Local structural relaxations induced by the presence of oxygen vacancies in this material strongly indicate why crystallographic shear planes could appear for significant concentrations of oxygen vacancies.

FIG. 4 (
FIG. 4 (color online).Vacancy-vacancy pair potentials plotted with respect to the distance between the O atoms in the regular lattice.For -Al 2 O 3 the Madelung potential is shown for divacancies I-IV, in V, calculated considering a ÿ2e charge at the center of the divacancy and ÿ1e at the vacancy positions.Insets (A) and (B) show the geometric rearrangements (shown as arrows) induced by the presence of the divacancy in ZnO.Red arrows indicate that the center is involved in two different relaxation spheres.

FIG. 3 (
FIG. 3 (color online).Energy barriers for O-vacancy diffusion with respect to the Madelung potential at the saddle point for ionic compounds.SrO and KCl are also included (Cl ÿ vacancy diffusion is plotted for KCl).The Madelung potential on the diffusion saddle point configuration is calculated employing a symmetric distribution of the charge as shown in inset (A) for the electronic density of MgO.Charges employed are ÿ1e; ÿ2e; ÿ1e for the initial, migrating oxygen, and final points.Inset (B) depicts the paths considered for diffusion on -Al 2 O 3 , ordered by increasing O-O distance.Path one shows the lower energy barrier.

TABLE I
e-e the electron repulsion of the electron pair left in the cavity.Notice that from the 2n electrons formed in the ionization process, 2n ÿ 1 are used to form n ÿ 1 anions and two non interacting electrons (2e ÿ ) are left which form an interacting electron pair (e 2 ) in the final system.