Theoretical study of bulk and surface oxygen and aluminum vacancies in aAl 2 O 3

The formation energy, geometry, and electronic structure of isolated oxygen and aluminum vacancies in bulk and on the~0001! surface of corundum ( a-Al2O3) have been investigated by means of periodic calculations in the framework of density functional theory within the generalized gradient approximation and large supercells. The energy cost to form an oxygen vacancy in the bulk is estimated to be of the order of 10 eV, whereas that corresponding to the formation of Al vacancies is found to be at least a 30% larger. The relaxation of the material is rather small for both defects. The removal of an oxygen atom in bulk a-Al2O3 is accompanied by the appearance of an impurity level in the gap, which is a strong indication of electron localization. This has been further confirmed by integration of the density of states in the energy interval corresponding to the impurity level and by several other theoretical analyses. For the a-Al2O3(0001) surface, the formation of oxygen and aluminum vacancies exhibits many similarities with the bulk; the energy cost to form Al vacancies is much larger than for O vacancies and, in both cases, it is accompanied by rather small atomic displacements. However, there are also significant differences between bulk and surface oxygen and aluminum vacancies. Thus the formation energy of one of these point defects in the surface is rather smaller as expected and, more importantly, the degree of electronic delocalization is also larger.


I. INTRODUCTION
2][3][4] In particular, metal oxide surfaces are widely used in catalysis as inert supports for the active component.However, it has been recently proposed that the support has a crucial role in the properties of the catalyst.In fact, oxides are able to stabilize metal particles, the size and morphology of which may change with the oxide nature. 5Thus oxide surfaces are responsible for the dispersion grade and for the electronic and chemisorption properties of the active particles in a catalyst. 2,5However, despite their importance, there are still fundamental questions concerning the atomic structure of both clean and metal-covered metal oxide surfaces. 6][11] Depending on the preparation and treatment, a number of different point defects are known to be present in oxides.In particular, the presence of oxygen vacancies in MgO and other alkaline-earth oxide single crystals has been established long ago and is one of the point defects most extensively studied.  In t case of MgO, the oxygen vacancies are formed either by thermochemical reduction or by neutron irradiation: 34 these point defects are generally described as F centers, although different kinds of these centers exist depending on whether neutral atomic oxygen ͑F cen-ter͒, O Ϫ (F ϩ center͒, or O 2Ϫ (F 2ϩ center͒ is removed from the solid.In each case there are formally two, one, or zero electrons trapped in the cavity left by the removal of the oxygen species.The presence of F centers in MgO is char-acterized by a well-defined spectroscopic feature at ϳ5.0 eV which is interpreted as the lowest spectroscopic transition of the two electrons trapped in the cavity by the external Madelung field.This peak has been resolved into two components due to F ϩ centers at 4.96 eV and to neutral F centers at 5.03 eV, 12,15 and is indeed used to titrate the defect concentration and their diffusion through the solid. 35The optical peak is also accompanied by a luminescence emission at a lower energy ͑see Ref. 36 and references therein͒.Much less is known about the surface oxygen vacancies in MgO.This is because of the lower concentration of these point defects at the surface, which renders their detection is more difficult.Nevertheless, surface oxygen vacancies in MgO have been characterized and some tentatively assigned.17,18,20,37 Recently, some of these assignments have been confirmed through comparison with ab initio calculations on the optical properties of these centers.[38][39][40] Oxygen vacancies are not the only point defects observed in MgO; under neutron irradiation, other bands appear, some of which have been attributed to the presence of aggregates of F centers, 12,41 and again this has been confirmed by accurate ab initio cluster model calculations.42 The occurrence of oxygen vacancies is not exclusive of simple oxides: they are present in more complex materials such as SiO 2 , Li 2 O, ZrO 2 , TiO 2 , and Al 2 O 3 .31,36 However, because of their structural or chemical complexity, information about the properties of oxygen vacancies in these materials is scarce.Over these materials, Al 2 O 3 , or corundum, has been of special interest because of its use in supported model catalysis [43][44][45] and the ␣-Al 2 O 3 (0001) surface may be considered as a paradigm of a complex surface and, consequently, has been the subject of many experimental and theoretical studies.46 -53 Corundum is one of the many possible polymorphic forms of aluminas, which have a broad range of applications in the automobile industry, catalytic reactions, and as coating in several technologies.Corundum itself is the catalytic support for the ethylene epoxidation reaction.54 In the past two decades a considerable amount of research has been devoted to the optical properties of materials based on ␣-Al 2 O 3 .These are mainly due to the presence of radiationinduced defects and, as in the case of MgO, oxygen vacancies are one of the most abundant point defects.There is also evidence for another type of defect which is originated by the removal of cations from the oxide.[55][56][57][58][59][60] Nevertheless, oxygen vacancies are thought to be the most stable defects with characteristic optical absorption spectra showing an absorption band at 6.05 eV, which was found to be formed by two peaks at 5.91 and 6.22 eV and, again as in the case of MgO, are attributed to F ϩ and F centers in ␣-Al 2 O 3 .56 The number of electrons present at these point defects and their spatial extent affects the electronic properties of adsorbed species and consequently they may be responsible for the catalytic activity of the support.For example, the unexpected catalytic behavior of small palladium and gold particles deposited on MgO for CO oxidation has been partially ascribed to the presence of this kind of defects at the oxide surface.[61][62][63] However, while the localized character of the F-centertrapped electrons in MgO is well established, little is known about the electronic structure of these Frenkel defects in ␣-Al 2 O 3 .Obviously, this is a consequence of the much more complicated structure of aluminum oxides compared to those simple cubic oxides.Stashans, Kotomin, and Calais 64 have used semiempirical methods and an Al 26 O 39 stoichiometric cluster model approach to study the optical properties of F ϩ and F centers in corundum crystals.From a limited geometry optimization procedure, these authors conclude that the four nearest-neighbor Al atoms relax outwards whereas the two oxygen atoms nearest to the cavity relax inwards-i.e., decreasing the distance to the center of the vacancy.These semiempirical calculations also suggest that in the electronic ground state the one-electron functions describing the states of the trapped electrons are well localized within the vacancy cavity.One important point is that the analysis of the atomic population indicates a substantial degree of covalence.This is, however, in contradiction with ab initio Hartree-Fock studies which indicate an almost complete ionic character of corundum.65,66 The different description of the ground-state electronic structure of this material can be attributed to two different effects: on the one hand, to the neglect of the Madelung potential in the model and, on the other hand, to the use of empirical parameters.More recently, Xu et al. 67 used a 120-atom cell to study the relaxation and energy formation of a single O vacancy in ␣-Al 2 O 3 within a particular implementation of the local density approximation ͑LDA͒ of density functional theory ͑DFT͒ based on the use of orthogonalized linear combination of atomic orbitals.In this case, the 4 Al and 16 O atoms nearest to the oxygen vacancy were fully relaxed, although the energy minimization procedure is carried out using a numerical grid to approach the potential energy surface.Xu et al. claim 67 that the optimized atoms are all relaxed outwards up to 16% of their initial positions in the bulk.These authors report that the energy required to form the oxygen vacancy point defect is 38.36 eV for the unrelaxed material and drops up to 5.83 eV when relaxation is included.These results are quite surprising.First, one would naively expect that even with two electrons trapped in the vacancy cavity the removal of an oxygen anion will lead to a decrease in the Pauli repulsion with the neighboring atoms, which will tend to relax towards the cavity. Secon, the oxygen vacancy energy formation for the unrelaxed material is enormous and that corresponding to the relaxed geometry is by far too small.Xu et al. 67 justify this value by comparing to the oxygen vacancy of MgO, which is predicted to be 4.7 eV.This is about half of the value obtained by any standard ab initio electronic structure approach, either DFT or Hartree-Fock.This cast reasonably doubts on the reliability of this approach even if it is able to properly predict the parameters of the crystal structure of bulk corundum.Finally, the LDA calculations of Xu et al. 67 predict the F-center electrons of corundum to be largely delocalized. Hoever, there is compelling evidence that the LDA has a strong tendency to delocalize the electron density, resulting in an unphysical description contrary to well-established experiments.68 From the previous discussion it appears that the electronic structure of the oxygen vacancies or F centers of ␣-Al 2 O 3 are far from being understood, even from a qualitative point of view.Therefore, despite the fact that it has been previously found that concerning the relaxation of the ␣-Al 2 O 3 (0001) surface, the LDA approach yields similar results to those obtained by employing other more sophisticated exchange-correlation functionals: 51 it is clear that in the view of the above comments a more detailed study is needed.This is precisely the aim of the present work.To this end, a broad range of theoretical methods are used to study the energy formation, structural relaxation, and electronic structure of bulk and surface oxygen vacancies in ␣-Al 2 O 3 ; aluminum vacancies have also been considered for comparison purposes.Most of the calculations are based on different implementations of DFT, but Hartree-Fock calculations have also been carried mainly to establish the degree of the localization of the electrons trapped in the cavity, a property which can be difficult to measure within DFT because of the well-known tendency of the LDA and generalized gradient approximation ͑GGA͒ to overestimated electronic delocalization.68,69 A periodic approach was used to model the formation of Al and O vacancies in bulk ␣-Al 2 O 3 and in the ␣-Al 2 O 3 (0001) surface; large supercells were used to minimize the interaction among vacancies.However, a cluster model approach has also been used to better analyze the electron density around an oxygen vacancy.The computational details and theoretical methods used to carry out the analysis of electron density are described in Sec.II: further details about the material models are given in Sec.III.The relaxation of the nondefective surface is discussed in Sec.IV A; results concerning geometrical structure and formation energy of aluminum and oxygen bulk and surface ␣-Al 2 O 3 vacancies are reported in Sec.IV B; the analysis of electron density at oxygen vacancy sites, carried out by means of various topological analyses, is also discussed in Sec. IV C. Finaly, the conclusions of the present work are presented in Sec.V.

II. COMPUTATIONAL DETAILS AND METHODS OF ANALYSIS
Periodic DFT calculations were carried out using two different computational packages and two different exchangecorrelation functionals.In both cases, the Bloch functions are expanded in a plane-wave basis set, but with a different cutoff energy because of the different treatment of the effect of core electrons.The first set of DFT calculations use the Perdew-Burke-Erzenhof ͑PBE͒ implementation of the generalized gradient approximation, the Kohn-Sham Bloch functions were expanded up to an energy cutoff of 952 eV, a Goedecker-type pseudopotential was used to mimic the effect of the inner electrons, and, finally, a Monkhorst-pack set of 4ϫ4ϫ4 k points was used for integration in reciprocal space.These DFT calculations were carried out using a parallel version of the CPMD code 70 running in an IBM SP4 machine.A second set of GGA DFT was carried out, but now within the exchange-correlation functional of Perdew and Wang ͑PW91͒.In this case, the core electrons were described by the projected augmented-wave ͑PAW͒ method of Blo ¨chl. 71This is essentially an all-electron frozen-core method that combines the accuracy of all-electron methods with the computational simplicity of the pseudopotential approach, especially in the implementation of Kresse and Joubert. 72The use of this approach permits to obtain accurate results with an energy cutoff of 415.0 eV, significantly smaller than the one used in the pseudopotential PBE calculations.][75] The qualitative and quantitative description of the charge distribution and localization of the pair of electrons left upon the removal of a neutral oxygen atom has been also explored by taking advantage of a topological approach that combines the theory of atoms in molecules 76 ͑AIM͒ and the electron localization function ͑ELF͒ initially introduced by Becke and Edgecombe 77 and recently redefined by Silvi. 78In the AIM theory the surfaces of zero electron density gradient define atomic basins and integration of electron density in a given atomic basin provides an accurate estimate of the number of electrons that this atom has in the molecule.Among the numerous applications of AIM we quote those used to estimate the degree of ionicity of simple oxides and corundum 66 and more recently to rigorously characterize oxygen vacancies in MgO as pseudoatoms. 79The original formulation of the ELF is based in the ratio between the curvature of the electron pair density for electrons with identical spin in the system of interest and that of the homogeneous electron gas.Silvi has recently proposed a new and equivalent definition which is more physically appealing since it is based simply in three simple ratios: that of the parallel spin-pair concentration to the antiparallel spin-pair concentration, to the total spin-pair concentration, or to the single spin-pair concentration. 78In any case, the analytical expression of the ELF restricts its values between 0 and 1; large ELF values indicate that at the current point in real space the electrons are more localized than in the uniform electron gas of identical density.The topology of the ELF permits us to define different types of basin and in particular monosynaptic valence basins which represent electron lone pairs: numerical integration within this basins provides a reliable estimate of the electronic population involved in this particular basin; for more details, see also Ref. 79.
AIM and ELF calculations have been carried employing the densities obtained from periodic and cluster model representations of the ␣-Al 2 O 3 bulk and surface oxygen vacancies as obtained at different levels of theory and considering always the optimized structures obtained with the PW91 functional.The topological analysis have been carried out using the TOPMOD program 80 using the electron density obtained from a single-point calculation on a suitable cluster model ͑see next section͒ using the GAUSSIAN98 suite of computer codes. 81

III. MATERIAL MODELS
In this work, periodic models have been used to represent bulk and surface ␣-Al 2 O 3 without or with either aluminum or oxygen vacancies.In each case periodic calculations were carried out with the material modeled by a suitable supercell constructed in the hexagonal lattice by repetition of the cell in three dimensions.The size of the supercell used to model the oxide was varied during the calculations in order to ensure that convergence was fully achieved.These supercells were obtained by multiplication of the lattice vectors in the basal plane, a and b, by a factor of 2 or 3 resulting in a 2 ϫ2ϫ1 or 3ϫ3ϫ1 cell.The total number of atoms in these supercells ranges from 80 up to 270 atoms, depending on the size and type, bulk or surface, of the unit cell.In each case, the oxygen vacancy was introduced by removing one atom from the supercell.An important point here is that with the size of the present supercells, the minimum distance between repeated vacancies in the nearest-neighbor cell is larger than 9.0 Å; this ensures negligible interaction effects between adjacent vacancies.Notice that these are electrically neutral vacancies.
For each type of vacancy considered in the present work, different optimization strategies were used.For bulk ␣-Al 2 O 3 , and in an initial step, the four Al atoms closer to the removed oxygen atom and its 12 nearest-neighbor oxygen atoms were fully relaxed; see Fig. 1.Looking at the ideal bulk structure, these four aluminum atoms are oriented following tetrahedral directions around the vacancy, with two atoms at ϳ1.86 Å and the other two at ϳ1.97 Å.In subsequent steps, more next-neighbor atoms were allowed to relax and in some cases a full geometry optimization of the totality of atoms in the supercell was carried out.The energy of defect formation was calculated as the energy difference of the perfect bulk supercell and that of the relaxed supercell with the vacancy plus the energy that of an isolated oxygen atom computed in a 10ϫ10ϫ10 Å 3 box and using the spinpolarized version of DFT for two unpaired electrons in oxygen.For convenience, a slightly different box was used for the PBE calculation, but as we will show below, this has no effect on the calculated results.Obviously, using the same box leads to exactly the same results.Surface calculations were performed in a similar way, but adding a vacuum width of 12 Å in the direction perpendicular to the surface.The vacuum cuts the bulk in the ͑0001͒ plane terminating with an aluminum layer.The repeated slabs of finite thickness thus formed contain several aluminum and oxygen layers in an a-b-c packing of ͑Al-O-Al͒ layers, which is repeated 5 and 6 times in the PBE and PW91 calculations, respectively.The energy required to the formation of a surface or subsurface defect was computed using both the relaxed Al-terminated surface, without and with the point defect, and the atomic energy of oxygen as described above.
In the case of oxygen vacancies, cluster models were also used, mainly to carry out the topological analysis.3][84] In the present work, the cluster used to model the oxygen vacancy in bulk Al 2 O 3 contains 26 cations and 12 anions in the quantum region, but the cations were treated at different levels depending on the distance to the vacancy.Hence the four Al cations nearest to the oxygen vacancy were treated using an all-electron basis set while the remaining 22 cations were treated as total ion potentials ͑TIP's͒. 85These 26 cations and 12 anions form the cluster inner region which is further embedded in 3557 point charges ͑PC's͒ with formal values of Ϫ2 and ϩ3 to account for the long-range Madelung potential effects. 86,87The geometrical form of this cluster is represented by a sphere centered in the oxygen vacancy.The cluster used to model the vacancy on the Al-terminated ␣-Al 2 O 3 surface is the oxygen-centered Al 24 O 22 cluster used in a previous work, 84 but in this case the central oxygen atom has been removed to create the F center.These clusters were used to perform single-point calculations, using the pertinent geometry, unrelaxed or relaxed depending on the case, aimed at obtaining the density to be used in subsequent topological analysis or to obtain the isotropic hyperfine constants of the charged F ϩ point defects.In the view of the tendency of the LDA and, to some extent, the GGA to excessively delocalize spin 68 and charge density, 69 the topological analysis was also carried using the Hartree-Fock ͑HF͒ method and the hybrid B3LYP exchange-correlation functional.The latter employs the gradient-corrected Becke's three parameters hybrid exchange functional 88 in combination with the correlation functional of Lee, Yang, and Parr. 89This method offers the best results in terms of bond strengths in molecules and solids, and in general provides a very accurate description of several properties at the molecular level 90 and provides a reasonable description of the band gap in ionic solids. 91,92ll calculations have been carried out without including spin polarization, except to compute the atomic energies.In the case of neutral oxygen vacancies this does not represent a problem because there is always an even number of electrons in the unit cell.In the case of Al vacancies the removal of a neutral Al atom implies the formation of three holes on the valence oxygen 2p band.In this work it is assumed these three holes are equally distributed along the oxygen 2p band.In chemical terms this is equivalent to assume that all oxygen atoms are slightly reduced.However, one must advert that in the case of MgO neutral cation vacancies lead to paramagnetic defects in the bulk 93,94 and in the surface. 95In any case, this possibility was not considered here because the main interest lies in the analysis of the electronic structure of oxygen vacancies, Al vacancies being considered mainly to compare the formation energies.To this end non-spinpolarized calculations seem to be accurate enough.A detailed study of the electronic structure of Al vacancies is out of the scope of the present work.

A. Relaxation of the perfect ␣-Al 2 O 3 "0001… surface
In order to check the consistency of the present computational approach we first discuss the structure of the bulk and of the perfect Al-terminated ␣-alumina surface and compare to experiment and to recent accurate all-electron 50,51 and pseudopotential 52 periodic calculations, respectively.For the bulk corundum the unit cell crystallographic parameters obtained with either PW91 ͑VASP͒ or PBE ͑CPMD͒ are within 1% of experiment.For the surface, a common feature of all these studies is the prediction of a very large degree of relaxation, especially for the outermost Al atomic layer.However, the surface relaxation involves many atomic layers and its direction is not straightforward.1][52] The degree of relaxation predicted by the present GGA/PW91 and GGA/PBE levels of theory for this surface are summarized in Table I.Overall, the present results are in agreement with the above-mentioned previous studies.However, there are two points which merit a further discussion.First, it is important to see that relaxation goes deeper than previously imagined.Only the work of Verdozzi et al. 48previously considered a thick enough slab although within the LDA approach.The present GGA results, either PW91 or PBE, are in very good agreement with the LDA results and with the GGA results of Lodziana and Norskov; 52 the agreement is even better if the results are corrected as usual to remove the artificial z component of the surface dipole.The second point concerns the importance of fully optimizing all coordinates of various surface layers as previously pointed out by Wander et al. 50Finally, it is also worth pointing out that, as far as the surface structure is concerned, the influence of the exchange-correlation functional used in the calculations is rather small.The good agreement between the present GGA results and the previous LDA description, both involving full optimization of up to 11 atomic layers, plus the coincidence with previous all-electron HF and B3LYP studies pointing out the importance in the pseudorotation relaxation, strongly suggest that concerning the structure of the Al-terminated ␣-Al 2 O 3 (0001) surface this is likely to be the final word.To conclude this subsection we mention that the GGA surface energy of this Al-terminated unrelaxed surface is 3.52 J m Ϫ2 .This rather large value is considerably reduced upon relaxation, the final value being 1.55 J m Ϫ2 ͑1.57J m Ϫ2 after including dipole correction͒, in excellent agreement with the GGA results reported by Lodziana and Nørskov. 52

B. Bulk and surface O and Al vacancies in ␣-Al 2 O 3
The energy formation of an isolated oxygen vacancy, either in bulk or surface ␣-Al 2 O 3 , and the corresponding geo-metrical changes induced by the presence of these point defects have been studied following an hierarchical optimization scheme.For the bulk F-center geometry optimization, the atoms around the vacancy are distributed by shells and each shell is labeled according to its coordination sphere with respect to the O vacancy.Hence, in the bulk, the first shell involves 4 Al first neighbors, the second shell 12 O second neighbors, and the third one involves 30 Al atoms.Then the hierarchical optimization procedure involves the full relaxation of atoms in the first plus second or in the first plus second plus third shells, respectively.In the case of a surface oxygen vacancy the hierarchical geometry optimization is defined in a different way.In the first step, the outermost five atomic layers are fully optimized and the 6 atomic layers underneath allowed to vertically relax.In a second step, the 11 outermost atomic layers were fully relaxed.For the bulk Al vacancy the first shell involves 6 oxygen atoms: hence, the first step has been precisely the full optimization of the coordinates of these six atoms; the second optimization step involves all atoms in a sphere of 3.3 Å of radius, which implies a full optimization of 16 atoms; the third step considers the complete third coordination shell ͑58 atoms͒ and, finally, we consider the full optimization of all atoms in a sphere of 5.8 Å of radius which involves 89 atoms in the supercell.
Ignoring geometric relaxation, the removal of a single oxygen atom from bulk alumina is 10.34 and 9.25 eV computed by the GGA/PW91 or the GGA/PBE approaches, respectively, and independently on whether a 2ϫ2ϫ1 or a 3 ϫ3ϫ1 supercell is used.The 10% difference in the energy formation predicted by the two GGA methods provides a measure of the intrinsic uncertainty of these approaches.Most of this uncertainty comes from the different way to estimate the energy of atomic oxygen.In the VASP code, all energies are computed with respect to the atoms in the ͑sometimes unrealistic͒ lowest-spin state; this means a closed shell for all atoms with an even number of electrons.For the oxygen atom, the calculation of a spin-polarized solution leads to an energy lower by 1.61 ͑with the present cutoff energy͒.However, in the PBE calculation using the CPMD code does not refer all energies to atoms and the calculated energy difference between the high-and low-spin states is 2.99 eV.A similar value ͑2.84 eV͒ is obtained using the GAUSSIAN98 code within an extended basis set of Gaussiantype orbitals and the PW91/GGA method.The difference in the estimate of low to high spin between the two methods is clearly due to the particular way the atomic references are constructed in VASP; see also Refs.96 and 97 for a related discussion involving atomic transition metal atoms.Notice that the high-spin state is close to the triplet-state solution for atomic oxygen, but the low-spin state is necessarily a mixture of the two possible spin states ( 1 D and 1 S) corresponding to the 1s 2 2s 2 2p 2 atomic configuration.Therefore, the above discussion permits one to explain the difference between the two GGA approaches.This is confirmed in an straightforward way by comparing the formation energy of the oxygen vacancy predicted by the two GGA approaches for the unrelaxed material, but using the non-spin-polarized solution for oxygen; these are found to be 11.95 and 12.24 eV, respectively.The difference between the two computational approaches drops now to 0.29 eV or less than 3%.These differences are small enough, hence, to avoid redundant results and only the PW91 results will be reported for the oxygen surface and for the Al bulk and surface vacancies.
Allowing relaxation of the first and second shells ͑Fig.1͒, the F-center energy formation is only moderately reduced, 10.28 eV at the GGA/PW91 level.Including relaxation up to the third shell does only reduce the formation energy by 0.02 eV.Therefore, one can consider these results as converged with respect to the deformation in the material induced by the O vacancy.The extent of geometric relaxation of the atoms surrounding the oxygen vacancy, in percentage, is reported in Table II.It is shown that, at the GGA level, the maximum atomic displacements are less than 5%, differences between PW91 and PBE values being less than 1%.In fact, only the two Al atoms nearest to the center of the cavity ͑those with initial Al-O distances of ϳ1.86 Å͒ are significantly displaced from their initial bulk positions and precisely by ϳ5%, while the other 14 ions exhibit much smaller displacements, less than 1%.This behavior contrasts with that reported in the work of Xu et al. 67 where much larger displacements are predicted, 16% and 8% for the first and second shells, respectively.The formation energy reported by these authors, 5.83 eV, is also very different from the values reported above.Therefore, it is likely that the LDA approach used by these authors cannot be trusted, at least for the kinds of systems studied in the present work.The present geometrical structure for the bulk oxygen vacancy is also somewhat in disagreement with the results of Stashans et al. 64 who predict an outward relaxation of the four Al cations surrounding the vacancy.These authors predict that two Al cations relax by ϩ2.8% and the two other cations by ϩ1.8%, whereas the present accurate results predict roughly Ϫ5% and ϩ1%.In any case these are rather small relaxations, especially when compared to the results of Xu et al. 67 TABLE II.Calculated geometry around a single oxygen vacancy in bulk ␣-Al 2 O 3 given as % of relaxation with respect to the bulk.Only the relaxation of the 16 ions nearest to the vacancy ͑4 Al and 12 O: see Fig. 1͒ is given.The number of atoms that have been fully relaxed is also indicated.The remaining ions in the supercell were kept fixed.The number of atoms in each supercell corresponds to that including the oxygen atoms which is later removed.A negative sign indicates relaxation inwards.Next, we turn our attention towards the surface oxygen vacancies of the Al-terminated ␣-Al 2 O 3 (0001) surface for which no previous information exists.This ␣-Al 2 O 3 (0001) surface may be regarded as a repetition of Al-O-Al layers and therefore one can distinguish two different types of surface oxygen atoms, those in the second ͑O-2͒ and those in the fifth ͑O-5͒ atomic layers: both cases have been considered.However, no low-coordinated oxygen vacancies have been studied.The predicted PW91 formation energy of these oxygen vacancies is of 8.91 and 9.53 eV for the oxygen atom in the second and fifth layers, respectively ͑Table III͒.These values are significantly lower than the one predicted for the bulk oxygen vacancy, especially for the oxygen in the second layer.These values include the effect of geometric relaxation in response to the removal of oxygen.Notice that while for the bulk oxygen vacancy the effect of geometric relaxation in the formation energy is rather small ͑0.06 eV͒, in the case of the surface oxygen vacancies the contribution of geometry relaxation is almost one order of magnitude larger ͑0.44 and 0.37 for O-2 and O-5, respectively͒.This larger energy contribution is consistent with a larger degree of relaxation of the ions surrounding the vacancy ͑Table IV͒.Notice, however, that there are some differences with respect to the bulk vacancy.For O-2 there are three nearest-neighbor Al atoms; one moves inwards by 13% and two outwards, but by a different extent ͑6% and 1%, respectively͒.In the case of O-5, there are four nearest-neighbor Al atoms and all relax outwards, at variance with the bulk vacancy, although in any case the atomic displacements are rather small.This result shows that, contrarily to simpler oxides such as MgO, the fifth atomic layer is still far from exhibiting bulk properties and indicates that, for certain purposes, the use of thin slabs or small clusters models of ␣-Al 2 O 3 (0001) should be avoided or used with extreme care.

Method
To end this subsection we consider the energy formation of Al vacancies in the bulk and in the Al-terminated surface.For the unrelaxed ␣-alumina, the PW91 cost to remove a neutral Al atom is of 16.83 eV, allowing the 6 nearest anions to relax drops this energy to 15.91 eV.Further relaxation up to 16, 58, or 89 atoms around the Al vacancy has only a minor effect since the energy formation decreases to 15.63, 15.51, and 15.49 eV, respectively.As in the case of the oxygen vacancy, the larger relaxation effect occurs for the nearest-neighbor ions, which in this case relax roughly by 7%-8% ͑Table V͒.However, for the bulk neutral Al vacancy, the second nearest-neighbor relaxation is also significant, ϳ4%.In the case of the surface Al vacancy, the PW91 energy formation for the unrelaxed substrate, calculated with respect to the corresponding relaxed clean surface, is 14.95 eV.Full optimization of 11 atomic layers ͑Table VI͒ decreases this energy to 14.64 or 14.78 eV if the dipole correction is taken into account.From these values it is clear that the thermal production of vacancies will always favor oxygen vacancies.On the other hand, the removal of an Al atom implies also removing 3 extra electrons per unit cell, slightly reducing the ionicity of the oxygen atoms ͑in this case removing one Al atom from a cell containing 108 Al and 162 O atoms decreases the formal charge of oxygen atoms from Ϫ2.00 to 1.98 a.u.͒.Given the large energy formation of this point defect, its electronic structure will not be further commented.

C. Electronic structure of O vacancies in bulk and surface
For MgO there is evidence that removing a neutral oxygen atom, either from the bulk or the surface, leaves two electrons which are well localized and trapped in the cavity.Finocchi et al. 98 have studied the MgO͑100͒ surface oxygen vacancies by means of LDA periodic calculations.They found that, at low defect concentration, a series of electronic levels is created in the gap of the MgO surface.From a solid-state physics point of view this may be considered the fingerprint of localized electrons.The electronic localization in MgO oxygen vacancies has also been qualitatively found by inspection of simple isodensity contour plots 30 and, more recently, rigorously characterized from sophisticated topological analysis of the density and of the electronic localization function. 79or Al 2 O 3 the only available information is the one from the semiempirical study of Stashans et al. 64 who from Mulliken population analysis found that for the bulk there are almost two electrons trapped in the cavity.No information seems to exist for the surface oxygen vacancies.In order to investigate in more detail the electronic structure of the O vacancy in bulk and surface ␣-alumina we rely on the density of states ͑DOS͒.For the nondefective material the bulk PW91 band gap is of 5.9 eV.This is slightly smaller than the 7.2 eV value obtained using an identical approach, but where the core electrons are described by ultrasoft pesudopotentials instead of the PAW method used here.On the other hand, the surface band gap obtained with the two approaches is on the 4.5-5.0eV interval. 52,84Notice that, in any case, the PW91calculated band gap is rather smaller than the experimental value of 8.75 eV, 99 although since GGA approaches tend to underestimate the band gap value, 91 this may be considered a reasonable value.In any case the important point here is that for the undefective material there is a well-defined gap in the energy levels.However, when an oxygen atom is removed from the 270-atom supercell a clear state appears in the gap that is better illustrated in the DOS reported in Fig. 2.This is a first indication of the localized character of the electrons trapped in the vacancy.The integration of the DOS in the energy region corresponding to the impurity level is 2.00.This is a second strong evidence of localization.However, this electron pair has a finite probability of being located TABLE V. Calculated PW91 geometry around a single Al vacancy in bulk ␣-Al 2 O 3 given as % of relaxation with respect to the bulk.Only the relaxation of the 16 ions nearest in a sphere of 3.3 Å radius is given.The number of atoms that have been fully relaxed is also indicated.The remaining ions in the supercell were kept fixed.The number of atoms in each supercell corresponds to that including the Al atoms which is later removed.A negative sign indicates relaxation inwards.

Method PW91
Atoms near the surrounding Al atoms.This is illustrated by the difference between the integral of the total DOS and of the sum of the DOS projected on the atoms which is of 1.20e.This value means that although there are two electrons trapped in the vacancy, around 40% of the corresponding electron density tends to delocalize in the nearest-neighbor atoms.Surprisingly enough, this prediction is in full agreement with the semiempirical results of Stashans et al. 64 The presence of a trapped electron pair is further supported from the analysis of the ELF, which exhibits a welldefined localization basin as shown in Figs. 3 and 4. Unfortunately, integration of the ELF inside this basin is not possible for the periodic calculation.Therefore, we rely on the corresponding analysis carried out on the electronic density obtained from the cluster models discussed in Sec.III.The use of a cluster model has some advantages, although it has also important limitations.Among the advantages one can quote the possibility to check different theoretical methods ͑HF, B3LYP, or PW91͒ and the possibility to study F ϩ centers, those which are formed when O Ϫ instead of neutral O is removed.To establish the adequacy of the cluster model approach we first discuss the formation energy of the bulk vacancy for the unrelaxed structure; this permits a direct comparison with the periodic calculations described in Sec.IV B. At the PW91 level the O vacancy energy formation predicted by the cluster approach is 9.90 eV: compared to the 10.34 eV value in the supercell calculation, this represents a ϳ4% effect only.It shows that the energy formation predicted by a cluster model has a large uncertainty, but it does also show that the cluster model is able to provide a reasonable, qualitative, description.The energy formation predicted by HF and B3LYP is 8.50 and 10.00, respectively, indicating that the three approaches qualitatively agree.The ELF analysis of the cluster density always predict a localization basin, in agreement with the periodic calculations.However, the integration of the electron density predicted by the two DFT approaches inside the volume of this basin appears to be strongly basis set dependent.This is quite a strange behavior because for the undefective cluster model the integration of the density inside the atomic basins is quite similar, Ϫ1.79e for the oxygen central atom at the HF level compared to Ϫ1.70e for PW91 and Ϫ1.73e for B3LYP, and rather independent of the basis set used for the atoms in the cluster model.Therefore, we have decided to report only the HF results for the integration of the density inside the ELF basin.This appears to be 1.28e, which is consistent with all previous results.Next, we consider the charged oxygen vacancy or F ϩ center.The cluster model calculations predict that, independently of the level of theory, the formation of an F ϩ center is ϳ2 eV higher than the cost to produce a neutral vacancy.In spite of this rather large formation energy F ϩ centers have been experimentally observed in ␣-Al 2 O 3 , although in samples irradiated with high-energy ͑14 MeV͒ neutrons. 100To conclude the description of bulk F and F ϩ centers we describe the calculated isotropic hyperfine constants A iso related to the F ϩ centers.A iso is a measure of the electron density at the nucleus and may serve to see whether the HF, B3LYP, or PW91 descriptions are consistent.These values have been measured and theoretically computed for Mg atoms surrounding an F ϩ center of MgO, 101 but we are not aware of similar studies for corundum.Therefore, we just report the A iso on the Al atoms surrounding the O vacancy.The HF values are ͑58.0,16.1 G͒ for the two different Al atoms, whereas the B3LYP and PW91 are ͑56.0,18.3G͒ and ͑52.1,18.6G͒, respectively.The close similarity between the values computed using different methods indicates that the electron density around the vacancy predicted by the three methods is similar and reinforces the above argument for the use of the HF density to carry out the integration on the ELF basin.Notice that the A iso is larger for the two Al nearest the center of the oxygen vacancy, as expected.Here it is important to point out that the above values for A iso have been computed using the unrelaxed geometry, using a 6-311G** basis 102 in the Al atoms surrounding the vacancy and without including floating basis functions at the center of the vacancy.The latter approach has been tested in MgO where there are almost six Mg cations surrounding the F center, whereas in the case of Al 2 O 3 there are four Al, but placed at different distances.Therefore, in this case the use of floating basis sets is likely to introduce artifacts and has been avoided.The A iso values can be compared directly to electron paramagnetic resonance ͑EPR͒ experiments.We have not been able to find evidence of EPR measurements in ␣-Al 2 O 3 .However, Kurtz et al. 103 have studied the F ϩ centers in Na, K, and Li ␤-Al 2 O 3 .In this form of alumina, oxygen atoms exhibit a local environment rather similar to ␣-Al 2 O 3 .Hence it is expected that F ϩ centers in ␣-Al 2 O 3 will at least have an order of magnitude of the hyperfine constants in ␤-Al 2 O 3 .The hyperfine coupling constants reported by Kurtz et al. 103 exhibit a slight degree of anisotropy, but interestingly enough the values of F ϩ centers in Li, Na, and K ␤-Al 2 O 3 are all in the 37-47 G interval, very close to the values predicted in the present work for the F ϩ centers of ␣-Al 2 O 3 .Comparison to other materials also permits to us obtain important information about the electronic structure of this point defect.To this end it is better to use the spin density at the nucleus.For MgO, Pacchioni et al. 101 report a value of 0.04 a.u.for the six Mg cations surrounding an F ϩ center.For ␣-Al 2 O 3 the present PW91 values are 0.21 and 0.17 for the two distinct Al atoms nearest to the vacancy ͑0.22 and 0.17 at the B3LYP level, 0.26 and 0.18 at the HF one͒.This is a clear indication that in the case of ␣-Al 2 O 3 the single electron in the oxygen vacancy is substantially more delocalized.However, one must advert that the basis set in the Al atoms is large and contains diffuse functions; hence, their average spatial extent is rather large.In any case, the use of A iso , which is an observable quantity ͑or of the spin densities͒, permits quite direct evidence that at least a part of the electron density of the single electron is localized in the cavity.
From the preceding discussion one can conclude that the electronic structure of an isolated neutral oxygen vacancy on ␣-alumina is very similar to those encountered in MgO, both materials having two electrons trapped in the cavity.The question now is whether surface oxygen vacancies of ␣-Al 2 O 3 have also a similar character to those of MgO.
However, to answer this question is not so simple.First, in the slab description of the ␣-Al 2 O 3 (0001) surface one must be aware that some spurious states do appear in the gap which are due to the bottom unrelaxed surface.These states have of course been removed in the calculation of the surface band gap and are the responsible for not being able to predict only a gap in the 4.5-5.0eV interval rather than a more precise value.When an oxygen vacancy in the second layer ͑O-2͒ is introduced in the supercell, a new state appears as in the case of the bulk.However, this new state is only visible when the DOS includes only the states with contributions from the outermost five atomic layers ͑Fig.5͒.The impurity state due to the oxygen vacancy is very similar to that corresponding to the bulk surface and, hence, one expects a similar electronic structure with two electrons trapped and well localized in the cavity.However, the presence of the spurious surface states does also preclude the integration of the DOS in the energy range corresponding to the oxygen vacancy state.Nevertheless, the ELF analysis predicts a well-defined basin as in the case of the bulk, although at present there is no way to integrate the electron density in the ELF basis predicted by the periodic calculations.
The cluster model approach description for the O-2 vacancy is qualitatively similar to that above described for the bulk point defect.The PW91 energy formation computed for a relaxed surface model of the clean surface and without taking into account the geometrical rearrangement in response to the presence of the defect is 6.88 eV, again ϳ2 eV lower than the 9.35 eV value obtained in the periodic calculations ͑Table III͒.The ELF analysis carried out in the cluster model also predicts an electronic localization basin, and in- tegration of the HF density within the volume of this basis is 0.86e, smaller than in the bulk, but still indicating a rather strong localization of the electron pair in the oxygen vacancy cavity.The more delocalized character of the surface F center as compared to the bulk is consistent with the calculated A iso constants for the surface F ϩ center of ␣-Al 2 O 3 (0001); the calculated PW91 values are 86.2,17.5, and 9.0 G for the three distinct Al atoms nearest to the cavity.The corresponding spin density values-0.28,0.18, and 0.06-indicate a significant degree of delocalization in the nearest Al atom.Notice again that no basis functions are included at the vacancy center and that the Al basis set contains diffuse functions.

V. CONCLUSIONS
The formation energy, geometry, and electronic structure of isolated oxygen and aluminum vacancies in bulk and on the ͑0001͒ surface of ␣-Al 2 O 3 have been investigated by means of periodic calculations in the framework of density functional theory within the generalized gradient approximation.Large supercells containing up to 270 atoms for the bulk and 120 for the slab surface model have been used to minimize the interaction of the vacancies.Cluster models have also been used, but for interpretation purposes only.
The energy cost to form an oxygen vacancy in the bulk is estimated to be of the order of 10 eV.This is very close to the energy required to form the same point defect in MgO, 9.55 eV, as obtained in the present work using a similar supercell.The energy formation of vacancies is found to be at least a 30% larger.The relaxation of the material is rather small for both defects.This is contrary to the large relaxation predicted by Xu et al. 67 based on a particular implementation of the LDA.The present results cast doubts on the reliability of such approach.Another important point concerns the displacements of the four Al nearest to the cavity in the bulk, two inwards and two outwards.This is also contrary to the prediction of Stathans et al. 64 who based on semiempirical calculations predicted an outwards displacement of these four Al atoms, but by a small extent as in the present firstprinciples calculations.The removal of an oxygen atom in bulk ␣-Al 2 O 3 leads to the appearance of an impurity level in the gap, a strong indication of electron localization.This has been further confirmed by integration of the density of states in the energy interval corresponding to the impurity level and by several other theoretical analysis including the study of the ELF function in the periodic supercell and in a cluster model.The formation of F ϩ centers is more energy costly: A iso values have been predicted.The analysis of the spin density reveals some important differences with respect to MgO, although one concludes that the single electron is essentially localized in the cavity.Here we would like to men-tion that the present results for the bulk neutral oxygen vacancies are in very good agreement with those reported by Matsunaga et al. 104 while the present paper was being reviewed.
For ␣-Al 2 O 3 (0001) the formation of oxygen and aluminum vacancies exhibits many similarities with the bulk.Hence the energy cost to form Al vacancies is much larger than for O vacancies and, in both cases, it is accompanied by rather small atomic displacements.However, there are also significant differences between bulk and surface oxygen and aluminum vacancies.Thus the formation energy of one of these point defects in the surface is rather smaller as expected and, more importantly, the degree of electronic delocalization is also larger.
To conclude, let us extract some consequences from the present study and dare to make some predictions.][40] Finally, we notice that the presence of oxygen vacancies on the ␣-Al 2 O 3 (0001) surface may have consequences on the spin state of adsorbed particles on this support as recently predicted for Ni atoms on different classes of oxygen vacancies in MgO͑001͒. 108 FIG. 1. ͑Color online͒ The local environment of the oxygen vacancy (V) in bulk Al 2 O 3 .Only the 4 Al nearest neighbors and 12 O ions are shown.The four aluminum cations occupy 2/3 of the anion octahedral interstices.

FIG. 2 .
FIG. 2. Density of states of the 270-atom supercell used to describe a single oxygen vacancy in bulk ␣-Al 2 O 3 .The impurity state is signaled with an arrow.

FIG. 5 .
FIG. 5. Sum of the projected density of states corresponding to the atoms on the five outermost atomic layers of the 18-layer slab supercell used to describe a single oxygen vacancy in the ␣-Al 2 O 3 (0001) surface.The impurity state is signaled with an arrow.

TABLE I .
Relaxation of atomic interlayer spacing ͑in % with respect to the unrelaxed geometry͒ of the top outermost layers of the Al-terminated ␣-Al 2 O 3 (0001) surface from different sources.Z-opt indicates the number of atomic layers vertically relaxed and Full-opt the number of atomic layers for which a full optimization of all atomic coordinates has been taken into account.
a Reference 51. b FP-LAPW calculations reported in Ref. 53. c Reference 48.d Reference 52. e Including dipole correction.

TABLE IV .
Calculated PW91 geometry of the 4 Al and 16 O atoms around a single oxygen vacancy in the Al-terminated ␣-Al 2 O 3 (0001) surface given as % of relaxation with respect to the relaxed clean surface; d i is the distance from the removed oxygen to a given atoms and d f the same value after removal and further relaxation.The 11 topmost layers of the slab model have been fully relaxed and the dipole correction taken into account.A negative sign indicates relaxation inwards.

TABLE III .
Calculated energy formation of surface oxygen vacancies in ␣-Al 2 O 3 ͑0001͒.The oxygen vacancy in the second or fifth layer is denoted by layer 2 and layer 5, respectively.

TABLE VI .
Calculated PW91 geometry around a single Al vacancy in the Al-terminated ␣-Al 2 O 3 (0001) surface given as % of relaxation with respect to the relaxed clean surface.The eleven topmost layers of the slab model have been fully relaxed and the dipole correction taken into account.A negative sign indicates relaxation inwards.