Rigorous characterization of oxygen vacancies in ionic oxides

P. Mori-Sánchez, 1,2 J. M. Recio, 1,3 Bernard Silvi, 3,4 C. Sousa, 2 A. Martı́n Penda ́s, Vı́ctor Luaña, and F. Illas Departamento de Quı ́mica Fı́sica y Analı ́tica. Universidad de Oviedo, 33006-Oviedo, Spain Departament de Quı ́mica Fı́sica, Universitat de Barcelona, 08028-Barcelona, Spain Departament de Cie `ncies Experimentals, Universitat Jaume I, 12071-Castello ́n, Spain Laboratoire de Chimie The ́orique. Universite ́ Pierre et Marie Curie, UMR-CNRS 7616-Paris, France ~Received 27 November 2001; revised manuscript received 3 April 2002; published 1 August 2002 !


I. INTRODUCTION
The study of F centers in insulating and semiconducting crystals has received significant experimental and theoretical attention for more than 30 years ͑see for example Refs.1,2, and references therein͒.F centers in metal oxide crystals result from oxygen vacancies that are not naturally associated to metal vacancy counterparts such as in Schottky defects.In the laboratory, defective solids containing an excess of oxygen vacancies can be obtained by chemical dehydroxylation of the surface, irradiation, doping with metal vapors, or Ar ϩ bombing. 2 As a result, new interesting defectinduced optical and electronic properties appear that are mainly controlled by the crystalline structure of the host material.At the surface, oxygen vacancies are associated to reactive centers in adsorption and catalytic processes of technological relevance. 2,3The catalysis experiments are interpreted by assuming that F centers in MgO are localized regions with a well defined chemical behavior. 3Consistent with this image, the electron-paramagnetic-resonance ͑EPR͒ spectra of the F ϩ center in MgO also point towards a substantial localization of the unpaired electron. 4These and other observations suggest that F centers in ionic crystals are localized objects calling for a rigorous chemical characterization and a quantitative measure of their electron localization.Both questions are very hard to be settled experimentally and pose a challenge to theory.
The subject has been the focus of several theoretical investigations, indeed.Supercell and embedded cluster studies [5][6][7] have accurately predicted the energetics and spectral features of F centers in MgO.However, attempts to quantify the electron localization failed to give conclusive results due to the use of Mulliken population analysis as well as to the observed basis dependence of the projected charges. 8Besides, charge density differences and spin density maps 4,9,10 provide a clear but only qualitative picture of the electrons trapped in the hole cavity.This image requires a quantitative description to fully understand the dependence of the center properties upon the local environment and the type of the host lattice.Therefore, we believe that a definitive, physically founded characterization of the nature of F centers in ionic crystals is still in need.
In this work, we investigate this problem following a topological approach that examines the crystal electron density combining the theory of atoms in molecules 11 ͑AIM͒ and the electron localization function ͑ELF͒. 12,13The AIM theory provides a rigorous partition of physical space into open quantum subsystems ͑atoms or functional groups͒ 14 such that physical properties can be determined for a subsystem and all subsystems contribute additively to the crystal properties.The ELF function was designed by Becke and Edgecombe to provide an orbital independent description of the electron localization. 12The ELF is defined as

͑1͒
where D and D ‫ؠ‬ represent, respectively, the curvature of the electron pair density for electrons of identical spins ͑the Fermi hole͒ for the actual system and a homogeneous electron gas with the same density.The analytical form of ELF confines its values between 0 and 1. ELF is close to 1.0 in those regions where the antiparallel spin-pair probability is high and close to zero where it is low.
Using the AIM method, Bader and Platts 15 discussed in 1997 a Li 14 F 13 ϩ finite cluster model of the LiF F center, while Madsen et al. 16 extended in 1999 the analysis using also the ELF to the F center arising from an extra sodium atom in sodalite.Both works suggested the F centers to be very localized quantum subsystems or pseudoatoms.
Here, we show that F and F ϩ centers located in the bulk and on clean and defective surfaces of MgO can be identified with pseudoatomic entities within the host lattice.Our main aim is to provide the physical grounds to the chemical be-havior exhibited by these pseudoatomic subsystems.In particular, we clarify: ͑i͒ whether they are to be described as dangling bonds, lone pairs, or single electrons; and ͑ii͒ whether there is any electron delocalization involving the rest of the lattice or a defined part of it.After these points are solved, we turn to the evolution of the properties of the defect as a subsystem, considering in detail its electronic population, spin multiplicity, and coordination at the vacancy site.

II. THE TOPOLOGICAL APPROACH: BASIC CONCEPTS
The study of the topology of scalar fields constructed from the multielectron wave function, like the electron density (r ជ ) or the ELF (r ជ ), has proven a very fruitful pathway to extracting chemical information from quantum mechanics.The topological jargon is introduced through the associated gradient vector fields ٌ ជ and ٌ ជ .Very succinctly, every gradient line of these vector fields starts, or ends, at a critical point of the scalar field, i.e., a point where the gradient is null.These special points are further classified according to the nature of the Hessian of the scalar field, which in Bader's terminology, 17 involve a pair of integer indices: the rank of the Hessian matrix, and the sum of the signs of its eigenvalues at the critical point.A (3,Ϫ3) critical point, for example, is a local maximum or three-dimensional attractor of the field, and the geometrical locus of all the points whose field lines end up at that site defines the attraction basin of the critical point.Three-dimensional attractors of the electron density are usually located at the nuclear positions, and their basins are referred to as atomic basins, but occasionally (3,Ϫ3) attractors not associated to any nuclei, or nonnuclear maxima ͑NNM͒ appear.Different atomic basins are necessarily separated by two-dimensional surfaces ͑separatrices͒ characterized by local zero flux of the gradient field.Every other critical point of the field is located onto these surfaces, and whenever two atomic basins share a twodimensional portion of a separatrix, a (3,Ϫ1) critical saddle point appears that signals the presence of a chemical bond between both basins.These points are called bond critical points ͑bCP's͒, and two unique gradient lines start at them that end at the two bonded nuclei.The union of these two lines is called the bond path.
The topology of the ELF provides a real space look at the pairing of electrons, and allows an easy identification of all the chemical concepts traditionally emanated from the Lewis model.The three-dimensional attractors of the ELF may occur at nuclear ͑core basins͒ or non-nuclear ͑valence basins͒ positions.Different classes of valence basins are distinguished by means of the synaptic order, the number of core basins sharing separatrices with a given valence basin, provided that all these cores belong to a common localization domain. 18A localization domain is the volume enclosed by one or more closed isosurfaces of (r ជ ).When a localization domain contains more than one attractor within it, it is called reducible.Otherwise it is said to be irreducible.Given this terminology, monosynaptic valence basins, V(X), correspond to electron lone pairs or to groups of electron lone pairs of atom X, disynaptic basins V(X,Y ) to the electrons of two-center bonds between X and Y, trisynaptic basin V(X,Y ,Z) to three-center bonds, and so on.
Since three-dimensional basins fill the space, relevant observables may be additively partitioned into basin contributions by integrating the relevant observable densities, say A (r ជ), over the basin under study, ⍀ i , As a simple application of the above prescription, basin populations, N(⍀ i ), are obtained as and their variances through the following formula: where (r ជ 1 ,r ជ 2 ) is the spinless pair density. 19t has been shown that the population variance of a given basin i can be readily written as a sum of contributions arising from all the other basins ͑covariance͒:

͑5͒
In this expression N ¯(⍀ i )N ¯(⍀ j ) is the number of electron pairs classically expected from the basin population, whereas N ¯(⍀ i ,⍀ j ) is the actual number of pairs obtained by direct integration of the pair function over the ⍀ i and ⍀ j basins.The variance is thus a measure of the quantum-mechanical uncertainty of the population of a given basin as a consequence of electron delocalization, whereas the pair covariance indicates the degree of correlation between the population fluctuations of two given basins.

III. MODELING AND COMPUTATIONAL DETAILS
We have chosen F and F ϩ centers in MgO covering the coordinations and symmetry environments depicted in Fig. 1.They are prototypical examples of oxygen vacancy systems suited to be investigated through a sensible computational modeling.In order to asses the performance of our results, we have completed extensive cluster in the lattice and supercell periodic quantum-mechanical calculations using the GAUSSIAN98 ͑Ref.21͒ and CRYSTAL98 ͑Ref.22͒ codes, respectively.The electron densities so generated have been analyzed with the EXTREME ͑Ref.23͒ and CRITIC ͑Ref.24͒ codes, while the ELF functions have been studied with the TOPMOD ͑Ref.25͒ program.To illustrate in a simple manner our main findings, we concentrate on the most representative cluster-type results.
Our clusters consist of a quantum-mechanical all electron region containing the vacancy and the neighboring Mg and O shells shown in Fig. 1.These atoms are described with flex-ible enough Gaussian-type double-basis sets plus polarization functions (͓11s11p1d/4s4p1d͔ for O and ͓13s13p1d/4s4 p1d͔ for Mg͒.The inclusion of basis functions at the vacancy position has been checked to introduce negligible effects.In particular, we have considered the ͓3s2 p1d͔ basis set of Ref. 7 with optimized exponents for the bulk F ϩ case.Peripheral oxygens are surrounded by a set of Mg 2ϩ effective core potentials in order to avoid the unphysical electron density polarization that otherwise occurs.The electrostatic lattice potential in the cluster region has been accurately represented by the Evjen's method, 26 using an appropriate set of fractional point charges that also provide neutrality to the system.Similar embedding models have been successfully used to determine the optical spectra in bulk and surface F centers. 7 A two-step procedure has been followed.In a first stage, the clusters have been examined at fixed bulk geometries under both the unrestricted Hartree-Fock ͑HF͒ and the nonlocal ͑Ref.27͒ density-functional approximations.Geometrical relaxations of the lattice around the defects were later taken into account at the HF level.As an interesting outcome, we have found that our basic results derived from the properties of the electron density and the ELF function appear to be mostly insensitive to the theoretical level and to the computationally expensive geometry optimizations.Thus, our discussion carried out on cluster HF results obtained at the bulk frozen geometry.

IV. RESULTS AND DISCUSSION
The F centers in MgO turn out to be proper quantum subsystems, or chemical objects.In all cases studied, a local maximum of the electron density ͑or non-nuclear maximum, NNM͒ appears near the vacancy center.This NNM is necessarily surrounded by a zero flux surface of the electron density gradient field that separates it from the rest of the lattice ions.The NNM shares an important number of properties with the host oxide it replaces, including basin size, net charge and topology, behaving as a coreless or valence-only anion ͑see Fig. 2͒.These are remarkable features if we notice the bulk electron density values at the O 2Ϫ nucleus and at the NNM site: 291 versus 0.02-0.03e/bohr 3 .Despite this disparity, the NNM basin is isomorphic to that of the pure host oxide, and both basins show bCP's that link them to their nearest (Mg 2ϩ ions), and next nearest neighbors (O 2Ϫ ions).Moreover, the electron density at the NNM-Mg 2ϩ bCP (0.021 e/bohr 3 for a bulk F center͒ is strikingly near the value found at its O 2Ϫ -Mg 2ϩ counterpart (0.034 e/bohr 3 ).The nature of the bonds is also unaffected upon replacing the oxides by the pseudoanions.In both cases the interactions are of the closed-shell type.This is a significant departure from the shared-shell bonding nature of NNM found for example in metallic systems, but agrees with that previously reported in ionic systems. 15,16,28s differences among different centers are regarded, the NNM basins turn out to be slightly smaller than the O 2Ϫ basin either in the bulk or at the surface sites (S), as a detailed exam of Fig. 2 shows.The distance from the basin attractors to their nearest bCP's, i.e., the topological radii, provide a precise and significant measure of the basin sizes: r(O 2Ϫ )ϭ2.28, r(F)ϭ2.01,r(F S )ϭ1.93, r(F ϩ )ϭ1.91, and r(F S ϩ )ϭ1.85 bohr.As topological radii correlate with atomic electronegativities (), 29 the relative size of the host and defect basins can then be understood by assuming a well defined sequence of electronegativities: ).The patterns of electronic localization near the defects may be easily gained by analyzing the ELF scalar function ͓(r ជ )͔, which approaches 0 and 1 in the limits of weak and strong localization.F centers in MgO emerge as very localized objects, always associated to a localization basin around a punctual attractor.The values at the attractors are close to 1, revealing a true Lewis pair and an unpaired electron in the F and F ϩ cases, respectively.ELF basins are connected through (3,Ϫ1) CP's.The ELF value at the (3,Ϫ1) CP on the surface separating the pseudoatom valence, V(vac), from the magnesium core basin, C(Mg), is very low, whereas that on the V(vac)-V(O) separatrix is higher than the V(O)-C(O) separation.In the bulk case, for instance, these values are 0.031, 0.134, 0.029, respectively.The vacancy and the valence shells of the n nearest oxides thus belong to the same localization domain.They define a combined object with stoichiometry O n (2nϩ2)Ϫ in F centers and O n (2nϩ1)Ϫ in F ϩ defects that we call a superanion.As a result, we find that the ELF delocalization is small between V(O) and V(vac) and negligible between C(Mg) and V(vac).This means, in turn, that the superanions are mainly stabilized by the electrostatic interactions with the neighboring magnesium cations.This superanionic ELF picture clearly distinguishes the O-vac from the Mg-vac interactions, and complements qualitatively the AIM results.
According to our calculations, the nature of the defective centers is governed by the actual formal charge and by the location of the vacancy.We have found clear trends along the Bulk → Surface → Step → Corner ͑BSSC͒ sequence.The electron density inside the pseudoatom becomes flatter and flatter, and the value at the NNM site approaches that at the neighboring bCP's.As a consequence, the pseudoatom gets softer and small geometrical and vacancy-lattice polarization effects appear along this series.These are enhanced in the charged F ϩ systems.A careful comparison of the ٌ ជ vector field maps plotted in Fig. 2 reveals that the nearest-neighbor atoms expand towards the vacancy as the coordination index decreases, the effect being larger for the F ϩ than for the F center.Accordingly, the NNM-Mg bCP's, and to a less extent the NNM-O bCP's, approach the position of the NNM, and FIG. 2. Flux lines of the electron-density gradient vector field ٌ ជ around the bulk ͑top͒ and bare surface ͑bottom͒ defects.From left to right we plot the perfect crystal, the F, and the F ϩ centers.The O 2Ϫ vacancy site lies at the center of all the plots, oxide ions occur along the main diagonals, and magnesium ions along the ͓100͔ directions.Thick lines represent the bond paths and interatomic separatrices.
the pseudoatom becomes less spherical and more polarized towards the oxide ions.It is to be noticed that these results are preserved if lattice relaxations around the defects are allowed.
The ELF localization domains enclosed by the ϭ0.8 isosurfaces are shown in Fig. 3 for the neutral F defects.Only the vacancy domain changes substantially through the BSSC sequence.In the bulk, it is a region of high and nearly constant electron localization.But as the coordination index decreases and space free of electrons becomes increasingly accessible, the superanion and V(vac) basins expand out of the crystal surface, and the defect turns more delocalized.As a direct effect of the increasing delocalization in space of the V(vac)-V(O) localization domain, the values at the O-vac (3,Ϫ1) CP increase, while those at the Mg-vac CP decrease along the BSSC sequence.These trends are smoothed in the single electron F ϩ cases.
A quantitative measure of the degree of localization is provided by the electron populations shown in Table I, obtained by integrating the electron densities over the NNM and ELF defect basins.A similar magnitude providing information about spin localization is found by integrating the spin density ͗S z ͘.The pseudoatom population is smaller than its formal value, and decreases along the BSSC sequence in agreement with the progressive delocalization of the defec-tive electrons among their neighbors.From 1.4 to 0.9 electrons dwell within the F NNM's, and from 0.7 to 0.6 in the F ϩ case.The step and corner geometries yield rather close populations for the F centers, which are about 30% smaller than those calculated for the bare surface and the bulk defects.The integrated spin density in the paramagnetic F ϩ centers confirms that the unpaired electron is mostly located in the pseudoatom basin even in the low coordination sites.It is interesting to observe that the ELF V(vac) populations are much closer to the perfectly localized picture than the NNM ones, even though both the electron density and ELF basins follow the same trends.The greater ELF populations are due to nonnegligible lattice contributions to the V(vac) population.This image gives additional support to the superanion concept.It is remarkable that the V(O)-V(vac) covariances increase along the BSSC sequence, thus enhancing the delocalization within the superanion valence shell.
One of our most interesting findings concerns the symmetry constrained position of the NNM and ELF attractors for the surface defects ͑see Table II͒.The NNM sinks substantially from its nominal position into the crystal, the displacement being larger as the coordination index decreases.The ELF attractor also moves from its nominal position, but this time in the opposite direction.This divergent behavior illustrates how the image provided by the electron density and the ELF function complement each other and increase the interpretative capabilities of the topological approach.Two different forces acting in opposite directions, the electrostatic potential and the Pauli repulsion, mould this behavior.Whereas the NNM is driven inside the bulk material by the dominant Madelung attraction to the nearest Mg 2ϩ basins, the ELF attractor moves outward where the Pauli repulsion from neighboring shells is minimized.
Interestingly enough, these results allow us to suggest some physical ideas concerning the surface reactivity of F centers.Firstly, these defects should be unfavorable positions for physisorption, that mostly relies on the balance of attrac-   tive electrostatic and Pauli repulsion forces between the surface and the adsorbate.The different location of the NNM and ELF attractors splits the optimum distance at which the adsorbate should rest over the surface into a classical ͑elec-trostatic͒ and quantum ͑Pauli͒ pair.The adsorption energy becomes smaller as compared to the case where both distances coincide.Secondly, chemisorption of Lewis acids should be favored on F centers, as the V(vac) basin has a strong dangling bond character and should, therefore, tend to form dative bonds with electron deficient reactants.Analogously, F ϩ centers might chemisorb radicals.In this case the vacancy basin merges with the Lewis's lone pair basin of the substrate within which most of the spin density is localized.The total number of basins of the surface plus adsorbate is lowered by one, and the total Pauli repulsion decreases, thus explaining the observed formation of superbase sites associated to the defects. 3Finally, both F and F ϩ centers should be favored protonation sites.Inserting a proton close to the attractor of the V(vac) increases the attractive potential energy and leaves the number of basins unchanged.Hence, the least topological change principle is fulfilled. 31The protonation of a F ϩ center should yield the formation of a weakly bounded hydrogen available for further chemical reaction, and a sec-ond proton insertion in a F center could lead to the formation of an adsorbed hydrogen molecule. 31,32e conclude that MgO F centers are clear quantum subsystems, definitely characterized by a high degree of electronic localization.When the anion vacancy is formed, a coreless pseudoanion, which behaves as an activated host site appears, and its valence basin merges with those of the nearest oxides to form a polyatomic superanion.The properties of these objects are mainly dependent on the actual coordination index of the defect.The topological approach provides the essential analytical tools to understand the nature of the F centers and to predict their chemical role in the bulk and surface reactivity of the material.

FIG. 1 .
FIG. 1. Quantum-mechanical regions of our embedded cluster models.From top to bottom, left to right: bulk (O h symmetry͒, corner (C 3v ), bare surface (C 4v ), and step (C s ).Small dark balls represent the vacancy position, white balls the Mg ions, and big dark spheres the oxide ions.

FIG. 3 .
FIG. 3. Localization domains of the F vacancy in the bulk ͑top left͒, bare surface ͑top right͒, step ͑bottom left͒, and corner ͑bottom right͒ in MgO.The limiting isosurfaces correspond to (r ជ )ϭ0.8.V(vac) basins are clearly identified as the biggest light domains, whereas C(Mg) ͑smallest domains͒ and V(O) basins appeared located in their corresponding crystallographic positions.

TABLE II .
Distances ͑in bohr͒ from the actual position of the NNM and ELF attractors to the nominal vacancy position.Symmetry precludes displacements in bulk defects.