Concentration of Vacancies at Metal Oxide Surfaces: Case Study of MgO (100)

We investigate effects of doping on formation energy and concentration of oxygen vacancies at a metal oxide surface, using MgO (100) as an example. Our approach employs density-functional theory, where the performance of the exchange-correlation functional is carefully analyzed, and the functional is chosen according to a fundamental condition on DFT ionization energies. The approach is further validated by CCSD(T) calculations for embedded clusters. We demonstrate that the concentration of oxygen vacancies at a doped oxide surface is largely determined by formation of a macroscopically extended space charge region.

Metal oxides are key materials for many technological applications. For example, MgO is used as a catalyst for methane oxidation, TiO 2 plays an important role as a photocatalyst, and RuO 2 catalyzes the oxidation of carbon monoxide. It is generally accepted that intrinsic point defects, in particular oxygen vacancies (also called F or color centers), play a decisive role in catalysis at oxide surfaces [1][2][3][4], but significant controversy exists regarding their formation energy, concentration, and charge state. In this paper we study these issues for MgO bulk and the MgO (100) surface in contact with an O 2 gas phase at realistic temperature and pressure. Furthermore, we consider that realistic metal oxides are typically doped, either intentionally or unintentionally. Although the experimental band gap of MgO is 7.8 eV [5], realistic samples are typically neither clear transparent nor insulating. Defects such as intrinsic point defects, impurities, and defect complexes can give rise to electron or hole conductivity [1,6,7]. In this paper, we neglect defect complexes (e.g. dopant plus vacancy). Thus, for our study the role of dopants is to create a Fermi level, i.e. a reservoir for electrons and holes. This is termed "the global effect of doping". We focus our discussion on p-type MgO, since it exhibits intriguing catalytic properties [1,8,9]. Still our theoretical model covers also n-type material, where the concentration of F centers is very low. Our main finding is that for p-type samples, surface O vacancies get doubly positively charged which lowers their formation energy and results in a significant defect concentration. In fact, the free energy of formation in thermodynamic equilibrium is negative under typical catalytic temperatures and pressures. We show that the limiting factor to formation of surface oxygen vacancies is the formation of a space charge region (Fig. 1). Although well-known for doped semiconductor surfaces, the effect of space charge and band bending on the concentration of surface defects has not been discussed so far.
The values of calculated defect energy levels and total energies are sensitive to the employed treatment of exchange and correlation (xc) of the many-electron system. Therefore, special attention is given below to this aspect: (see text) of F centers at MgO (100) for T = 1 000 K and normal pressure of oxygen, as a function of Fermi energy, ǫF, between valence-band maximum, VBM, and conduction band minimum, CBm. Realistic dopant concentration ND=10 18 cm −3 and surface charge density σ = 2.6 · 10 12 e cm 2 (solid lines) and the dilute limit σ → 0 (dashed lines) are shown. b) In p-type MgO (ǫF=VBM) under realistic conditions, band bending, due to formation of a space charge region, limits the formation of surface F 2+ s centers.
Our approach is to determine the best xc functional of the HSE family by the condition that the ionization potentials obtained with the optimized (opt-) HSE functional agree with the results of a G 0 W 0 calculation. In exact DFT such condition must be fulfilled exactly [10]. However, due to the limited flexibility of the HSE family of functionals, which range from the PBE [11] generalized gradient approximation via the HSE06 [12] hybrid functional to the PBE0 [13] hybrid functional, our approach is obviously not exact, but it is "the best compromise". In terms of the HSE [14] exchange and range-separation parameters (α, ω), the three mentioned functionals correspond to (0, arbitrary), (0.25, 0.11 bohr −1 ), (0.25, 0). Our approach is validated for neutral, embedded Mg x O y clusters using the CCSD(T) method.
In this work, we use the all-electron FHI-aims code [22] for the periodic structures (bulk and surfaces) and some embedded cluster calculations. FHI-aims employs atom-centered numeric basis functions and various xc functionals as well as the GW approach. The basis set and numerical grids are of high quality as defined by the tight [22] settings. For all periodic surface models full atomic relaxation is calculated using PBE at the PBE bulk lattice parameter (4.258Å). HSE calculations are performed at these geometries, since using HSE geometries for two smallest unit cells considered results in negligible changes in the calculated formation energies. Vibrational energy ∆F q vib (T ) and vdW contributions to the formation energies were analyzed as well [23], but found to be insignificant for this study.
Furthermore, we employ the TURBOMOLE program package [24] for various embedded clusters using different xc functionals and CCSD(T). Triple-zeta valence plus polarization basis sets [5s3p2d1f ] / [5s4p3d] are used [25]. For the CCSD(T) computations we correlate also electrons in the Mg 2s and 2p shells, using core-valence correlation consistent basis sets, cc-pCVXZ (X = D, T, Q). On the O 2− ions we use the aug-cc-pVXZ basis sets [26]. In both CCSD(T) and DFT calculations, the basis set superposition error was evaluated following the Boys-Bernardi counterpoise correction [27]. MgO clusters are embedded in a periodic point charge array using the periodic electrostatic embedded cluster model [28] in TUR-BOMOLE and a converged finite set of point charges in FHI-aims. (See supplemental information (SI) at [URL will be inserted by publisher] for more details.) To minimize non-physical polarization of peripheral oxygen anions by the embedding point charges, pseudopotentials were added to the first shell of embedding Mg 2+ cations (all-electron Hay&Wadt effective core potentials [29] in TURBOMOLE, and Troullier-Martinstype norm-conserving non-local pseudopotentials [30,31] in FHI-aims). The PBE lattice constant has been used for the embedded clusters. Apart from the outermost frozen shell of atoms, full relaxation is allowed for in the cluster calculations, except for the CCSD(T) and GW calculations and respective DFT values.
In realistic samples defects may get charged due to electron transfer between the dopant-induced Fermi level and the defect states [15]. The neutral oxygen vacancy in MgO bulk and at the MgO (100) surface has an energy level deep in the band gap. This state has s-like symmetry at the defect site, and is fully occupied by two electrons. Thus, a singly and even a doubly charged vacancy is possible. In FHI-aims this situation can be modeled in two ways: either by adding a constant charge density to the density entering the Hartree term, or by slightly modifying the nuclear charges of the atoms in the unit cell. Either concept enables us to describe a charged defect while the supercell is kept neutral. The "constant density approach" is the standard treatment in other periodic codes. For surfaces, where much of the supercell corresponds to vacuum, this approach is obviously unphysical, although it can be partially remedied by a posteriori correction schemes [16]. The other treatment corresponds to the virtual-crystal approximation (VCA) [17,18] of a crystal with dopants. We change the nuclear number of all Mg atoms in the supercell by ∆Z Mg = −q/N Mg , where q is the charge of the oxygen vacancy (+1 or +2), and N Mg is the number of Mg atoms in the supercell. This means that 1 or 2 electrons are transferred to the VBM, which is the Fermi level in the virtual crystal. Once known for one particular Fermi level, the O vacancy formation energy can be trivially calculated for an arbitrary Fermi level (see Eq. 1).
When removing atoms from the bulk or from the surface of a material, we need to consider also a reservoir to which the atoms are brought [19,20]. We assume a gas phase of O 2 molecules which is characterized by an oxygen chemical potential, µ O (T, p) [20]. For an isolated oxygen vacancy, the Gibbs free energy of formation is: Here, E q vac and E host are DFT total energies of defected and undefected systems, respectively, ∆F q vib is the change in vibrational Helmholtz free energy of the crystal upon defect formation, q is the defect charge, and ǫ F is the Fermi energy. The oxygen chemical potential is where ∆µ O contains the vibrational and other T -and pdependent terms [20]. We use the experimental binding energy without zero point energies E bind O2 = 5.22 eV [21] to reduce artifacts originating from the generalized gradient approximations for the binding of the O 2 molecule, but the calculated total energy of the free atom, E O , is calculated with the corresponding electronic-structure approach. ∆µ O = 0 defines the oxygen-rich limit.
First, we address formation energies for isolated vacancies in the bulk, G bulk,q f . We extrapolate our DFT formation energies to the dilute limit using Taylor expansion in terms of reciprocal supercell lattice constant. In agreement with related work [32], we find that the charge-transition levels (2+/+) and (+/0), as well as formation energies G bulk,q f , are almost independent on the xc functional within the HSE family, when referenced to the vacuum level. However, given the more realistic situation of p-type material, where ǫ F is at VBM, G bulk,q f does depend strongly on the choice of HSE parameters (α, ω) for q = 0. The formation energy of the neutral defect is insensitive to the functional. We find that the main error in charged defect formation energies is the error in the VBM position with respect to vacuum (also pointed out in [33]). For fixed ω the formation energies depend practically linearly on the exchange parameter α, which can be traced back to a linear dependence of VBM with respect to vacuum on α.
Next, we identify the optimal xc functional to describe the formation energies of F centers in MgO. The ion-ization potential at a fixed defect geometry for a given functional HSE(α, ω) is where both E q vac and E q+1 vac are extrapolated to the dilute limit. For ǫ F =VBM, I q→q+1 ∆SCF depends on (α, ω). In exact DFT the Kohn-Sham highest occupied orbital (HOMO) does not change with occupation and agrees with the ionization energy. A more practical request is that the HOMO, calculated by G 0 W 0 on top of the HSE electronic ground state, should agree with the HSE ionization energy identifying what we call the optimized HSE functional, opt-HSE that correctly describes the charge excitation of the defect. We determine the opt-HSE functional for fixed ω = 0.11 bohr −1 . The ionization energy I 0→+ for ǫ F at VBM at F 0 geometry is calculated for an embedded Mg 6 O 9 cluster model using FHI-aims. The Fermi level ǫ F is obtained as VBM = E +1 host − E host using HSE, and from the HOMO of the host system in the corresponding G 0 W 0 @HSE calculations. The ionization potential shows a near-linear dependence on the exchange parameter α (Fig. 2) for both ∆SCF and GW method. The intersection of the two linear fits is at α=0.27, very close to HSE06 with parameter set (α=0. 25, ω = 0.11 bohr −1 ). We therefore use HSE06 as our opt-HSE functional that correctly describes the charge excitation of the defect. The difference in formation energies with α=0.25 instead of α=0.27 is negligible for F 0 , less than 0.1 eV for F + , and less than 0.2 eV for F 2+ . FIG. 2. Ionization potential at F 0 geometry calculated for an Mg6O9 embedded cluster by ∆SCF with HSE xc functionals (black symbols) and from the HOMO of a G0W0@HSE calculation (blue symbols). The screening parameter is ω = 0.11 bohr −1 . Solid lines show linear fits to the ionization energies as a function of exchange parameter α.
We perform a validation for the F 0 formation energy using an unrelaxed Mg 6 O 9 cluster model and the CCSD(T) method. This results in a correction ∆CCSD(T) of the DFT formation energies of −0.09 eV for PBE, 0.07 eV for PBE0, and −0.28 eV for B3LYP. Adding these corrections to the DFT formation energies (∆µ O = 0) for a converged cluster size Mg 16 O 19 yields DFT+∆CCSD(T) results of 6.85, 6.88 and 6.89 eV, respectively. These numbers are in good agreement with the HSE06 F 0 formation energies 7.04 eV and 7.05 eV obtained from the same converged embedded cluster and periodic calculations, respectively, using FHI-aims.
Thus, HSE06 is the opt-HSE functional in accordance with GW as well as coupled-cluster results. Our results show that the experimental value for the bulk F 0 center formation enthalpy in MgO of 9.29 eV with respect to the O 2 molecule [34,35] is a significant overestimate. A likely reason is that thermodynamic equilibrium has not been reached in this experiment.
We are now on solid grounds to provide an accurate estimate of G surf,q f for isolated oxygen vacancies at the surface using our periodic model and the HSE06 xc functional. F 0 s , F + s , and F 2+ s formation energies in the dilute limit are 6.34 eV, 2.76 eV, and 0.55 eV, respectively, for ∆µ O = 0 and ǫ F at VBM. For more realistic conditions, the formation energies are lower, as shown in Fig. 1a). Formation energies of neutral O vacancies depend weakly on their concentration (up to approx. 3% for bulk and 12% for surface defects in MgO). On the contrary, due to the slow decay of Coulomb interaction with distance, the formation energy of charged defects will strongly depend on their concentration, as well as the distribution of the compensating charge. Thus, concentration of dopants N D and their distribution (doping profile) should have a strong global effect on the defect formation energies. The equilibrium concentrations η q of oxygen vacancies in three different charge states (q = 0 − 2) are determined by the minimum of the total free energy G of the system with interacting defects: is an effective formation energy of a vacancy in charge state q in the presence of other vacancies. The configurational entropy s conf accounts for energetically degenerate arrangements of the defects (see SI).
Surface defects are charged by accommodating charge carriers from the bulk. This results in depletion of the charge carriers and creation of a space charge layer in the subsurface region. The resulting electrostatic potential causes band bending and prevents more charges from the bulk to reach the surface, increasing the energy cost per defect. As a result, there are two leading electrostatic contributions to the formation energy of charged defects: (i) attraction to the compensating charge, and (ii) band bending. The first contribution reduces the formation energy, while the second contribution increases it. The thickness of the space charge layer, z SC , depends on the doping profile, and may be limited by the thickness of the material. We consider the case of uniformly distributed dopants and unconstrained z SC , but the discussion can be straightforwardly generalized to the more complex situations. To stay focussed on electrostatic effects, we also do not consider a possible (T ,p) dependence of the bulk charge carrier density N bulk e/h , i.e. N bulk e/h ≡ N D is a constant external parameter.
The dependence of G surf,q f on the surface charge density σ is calculated as follows. First, we calculate formation energies G surf,q f (σ, z SC ) at a fixed z SC , equal to the slab thickness d, using VCA (see SI). The calculated formation energies include both electrostatic effects mentioned above. The actual z SC is determined by N D as follows: where e is the absolute value of the electron charge. The formation energy as a function of σ and z SC is is the classic expression for the energy of the space charge region formation at a semiconductor surface. The temperature dependence of z SC and E SC (σ, z SC ) at fixed σ is neglected. The meaning of the last two terms in Eq. 8 is to replace the band bending contribution to the formation energy calculated for z SC = d with the one obtained for the actual z SC from Eq. 7. The remaining dependence on σ after subtracting qE SC (σ, d) from G surf,q f (σ, d) is due to the electrostatic attraction between the defect and the compensating charge.
We can now calculate the equilibrium concentration of O vacancies at a p-doped MgO (100) surface, using Eq. 5. G eff,q f (σ) is calculated from Eq. 6 with σ = eη 1 + 2eη 2 . The concentrations of F 0 s and F + s are found to be negligible at realistic T , p O2 , and N D . The F 2+ s concentration η 2 and corresponding z SC as functions of N D are shown in Fig. 3 for different temperatures and p O2 = 1 atm. Although the F 2+ s Gibbs free energy of formation at σ → 0 is small or even negative at elevated temperatures, the equilibrium defect concentration does not exceed ∼ 0.5% at N D ≤ 10 18 cm −3 . Thus, space charge layer formation can be a mechanism by which wide-band-gap semiconductor surfaces remain stable at high temperatures. The band bending profile for ǫ F = VBM at T = 1 000 K for p O2 = 1 atm and N D =10 18 cm −3 is shown in Fig. 1b). Under these conditions the bulk bands bend downwards by 0.6 eV, and the (2+/+) charge transition level is lowered from 2.2 eV to 1.7 eV above the Fermi level. At T = 1 000 K and p O2 = 1 atm, the contribution of the electrostatic attraction term is small for small N D , but becomes comparable to the formation energy in the dilute limit for N D > 10 18 cm −3 .
We have presented a methodology for calculating charged defect formation energies and concentrations at surfaces, taking into account electrostatic effects due to charge transfer between surface and bulk. Doped material has been simulated using the VCA, and an optimal DFT functional has been identified by validation with coupled-cluster and GW methods. Our analysis shows that the concentration of F 2+ s centers at the (100) terrace of p-type MgO can be as high as 1% at realistic conditions, while relative F + s and F 0 s concentrations are negligible. We demonstrate that the formation of charged vacancies creates a localized, although macroscopically extended, space charge region. This decreases charge transition levels with respect to Fermi level at the surface, raising the formation energy by up to 1 eV and, therefore, limiting the defect concentration. We conclude that electrostatic effects can largely control oxygen vacancy formation at the surface of metal oxides. Experimental information on doping profiles may provide new insights on catalytic activity of doped oxide surfaces.