Quasi-degenerate states and their dynamics in oxygen deficient reducible metal oxides

Transition and rare-earth metal oxides have attracted much attention because of their versatile chemical and physical properties, spanning the whole range from conductors to semiconductors and wide gap insulators.^1,2 This results in a large manifold of exciting catalytic, electrocatalytic, and even photo-catalytic activity in important reactions, as first illustrated in the water oxidation by Fujishima and Honda.³ Oxygen defects are crucial in defining these properties, which has been highlighted in many reports.^4–7 While entropic reasons dictate that all solid crystals must contain a certain proportion of defects at non-zero temperature,^2,8 this is particularly acute for oxides because they are typically synthesized via heating at high temperatures.^1,9,10 This procedure generates oxygen depletion, which provokes severe changes in the material properties and can be monitored, for instance, in the reduction of the work function.¹¹ 

Consequently, the understanding of the electronic structure of defective oxides is key, not only to optimizing and fine-tuning their performance but also aiding the design of novel materials with tailored properties.^2,4 Yet, many computational approaches are based on describing pristine (non-defective) structures and, thus, most of the rich chemical complexity is simplified below the levels of proper predictions.^12,13 Importantly, most of these simulations on non-defective materials lie at the core of search engines of large databases and machine learning algorithms.^14–20 In order to properly account for chemical and physical attributes, better theoretical tools to assess structural defects have to be developed and implemented in those high-throughput approaches.^21–24 In the present perspective, taking some established examples from the literature, we aim at defining new rules to include defects in such studies.

Oxides can be analyzed following the standard classification in terms of reducibility,²⁵ as depicted in Fig. 1. Here, the bandgap acts as a general descriptor for the metal–oxygen bond and reducibility. Highly ionic, non-reducible binary oxides, as those derived from alkaline earth metals or aluminum, present simple electronic structures with wide bandgaps and are, thus, effective insulators. The fully ionic model adequately predicts their properties, such as vacancy formation energies and diffusion barriers, which tend to be high in energy for those materials.^26,27 For these oxides, once the oxygen atoms are removed from the lattice, the two excess electrons are stored in the cavity left by the O atom.⁴ Such bound electrons are found, for instance, in defective MgO and give rise to additional optical absorption lines, from which their name color centers or F centers (from German Farbe = color) was derived.^2,8

In covalent oxides, such as SiO₂, the additional electrons left upon oxygen removal remain as dangling bonds of sp³ character at the Si center or the defect is quenched by forming homoatomic bonds, e.g., Si–Si units, as shown in Fig. 1.^28–31 

In reducible oxides, the electronic structure is more intricate and oxygen vacancies tend to be electronically active.^5,25 Due to the significant amount of covalency, cationic and anionic contributions cannot be easily disentangled, and a clear one-to-one mapping of formal charges and computed ones is not possible.³² Many reducible oxides behave as semiconductors, i.e., the valence and conduction bands are separated by a small energy gap (below 4 eV). In these materials, the vacancy formation energies are smaller than in non-reducible metal oxides as the electrons left by oxygen removal end up reducing the metal centers and populating otherwise forbidden electronic states. This change in electron occupation typically causes a mismatch in the ionic radii and induces a local lattice deformation, also known as a polaron.^33,34 Depending on the strength of the electron–phonon coupling and the energy difference between the vacancy-induced electronic states and the conduction band, electron transport may occur via polaron hopping or band-like conduction.^13,35 The oxygen vacancy can therefore act as an effective n-type dopant, which is beneficial for many (photo-)electrochemical applications.⁶ Not only electronic but also ionic conduction may take place via oxygen vacancy migration.³⁶ 

One of the best known reducible metal oxides is TiO₂ in its three polymorphs, rutile, anatase, and brookite.^3,37,38 Upon oxygen removal, two Ti centers are reduced to Ti³⁺, which have electronic states close to the conduction band,^39–41 forming mobile polarons, as was observed during ab initio molecular dynamics simulations.⁴² The extent of the formed polarons and the interplay between trapping and detrapping of the excess charge is largely dependent on the exposed surface facet.^39,41,43 Furthermore, TiO₂ and many other reducible oxides can easily lose more than one oxygen and store the deficiencies along one particular direction (glide planes), known as Magnéli phases.^2,44 The clustering of oxygen vacancies leads to broad distributions of electronic states in the bandgap and has a crucial impact on the applications as catalysts,⁴⁵ photoelectrodes, or charge transport layers.^46–48 

However, the situations described above are far too simple for the characterization of vacancies in many other defective metal oxides. Removing oxygen may lead to different geometric as well as electronic configurations with similar energies, i.e., quasi-degenerate states. In other words, one single vacancy state is not enough to account for the complexity of many oxides at finite temperatures, where a dynamic equilibrium between multiple configurations is present. In this perspective, we showcase four relevant metal oxides, namely, CeO₂,^50–52 In₂O₃,^53,54 MoO₃,⁵⁵ and BiVO₄,^56,57 which exhibit quasi-degeneracy. The understanding of their complex electronic structures is essential not only to improve their catalytic efficiencies but also to expand that knowledge to other reducible metal oxides. With this perspective, we aim at giving a guideline to improve current high-throughput approaches, which shall screen for possibly advantageous defective oxides and include the possibility of quasi-degenerate structures.

Density functional theory (DFT) is the workhorse of electronic structure simulations due to its favorable scaling and high (close to chemical) accuracy, providing practical estimates for the ground state energy and electron density of many-electron systems.^58–60 Its implementations differ in the description of the so-called exchange–correlation functional that is analytically unknown. Most commonly applied for metal oxides are functionals based on the Generalized Gradient Approximation (GGA), including PW91,⁶¹ PBE,⁶² and their revised versions revPBE,⁶³ RPBE,⁶⁴ and PBEsol.⁶⁵ The latter two typically yield more accurate adsorption energies or lattice parameters, respectively.⁶⁶ However, GGA often fails in accurately describing the electronic properties of reducible oxides due to its inherent self-interaction error,⁶⁷ in which one electron interacts with itself in an unphysical manner. This leads (i) to overdelocalization of the electronic density, overestimating the degree of covalency in the system, and (ii) to a systematic underestimation of the bandgap. The latter is particularly severe for small bandgap semiconductors since they are often incorrectly predicted to be metals. The error becomes relevant for elements with strongly correlated electrons, such as transition metals and rare-earth elements with d- and f-electrons.

Here, the sum runs over the occupation values $ρl$ of the atomic orbitals with the angular quantum number l. DFT+U is able to accurately describe reducible transition and rare-earth metal oxides; yet, the choice of U is not trivial and largely depends on the system and atomic occupancies and, therefore, on the charge associated with the cations.^75,76 Its value can be derived in a stand-alone fashion using a linear response approach to balance out the derivative of the atomic chemical potentials.^75,77 In many practical cases though, U is instead fitted to the bandgap⁷⁸ or redox energies⁷⁹ from experiments or higher levels of theory (e.g., hybrid functionals), either alone or coupled to the thermodynamics of some processes.^76,80,81

One of the drawbacks in assessing oxygen vacancy formation using the Hubbard correction is the U dependence on the local electronic environment.⁷⁶ In principle, the U parameters for the metal sites surrounding the vacancy should differ from those that are unaffected by it. Hybrid functionals circumvent this issue, introducing non-locality in the exchange-correlation energy by incorporating a certain amount of exact exchange $EXHF$ from Hartree–Fock theory.⁸² Compared to DFT+U, this is a rather computationally expensive procedure (1–2 orders of magnitude larger); however, the implementation of these approaches in periodic codes has improved over the last years and are nowadays becoming routine (except for dynamic properties). Out of those, commonly used for solids are the full-range PBE0⁸² and range-separated HSE (HSE03/HSE06)⁸³ functionals, the latter replacing long-range exchange by DFT exchange, which allows for smoother convergence. Like most hybrid functionals, they employ a fraction of 25% $EXHF$⁠.^82,83 This has shown to be too much for several transition metal containing oxides as it leads to overestimated band splittings, including exceedingly large bandgaps, and to a nonphysical stabilization of high spin states.^20,84–87 Defect levels^88,89 and polaron hopping barriers,⁸⁶ and therefore, the localization of excess charge, are crucially dependent on the amount of $EXHF$ introduced. The optimum fraction $a(EXHF)$ of exact exchange can be obtained through empirical fitting to experimentally measured constants, such as bandgaps⁸⁵ or polaron hopping barriers,⁸⁶ but it can also be calculated self-consistently^90,91 or approximated by $a(EXHF)≈1ϵ∞$ for many oxide semiconductors.^88,89 Typically, 10%–15% of $EXHF$ was found to describe well the electronic structures of transition metal oxides,^85–87,92 as well as rare-earth metal oxides.^91,93

To achieve a more accurate description of defect states in reducible oxide semiconductors, it is necessary to go beyond density functional theory. By definition, the electronic structure of excited states, i.e., the conduction band properties, cannot be assessed by ground state DFT. Many-body perturbation theory, such as the GW method, can determine the quasiparticle band structures and excited state properties with high precision,^94,95 but its implementation is yet limited due to its large computational cost and required proficiency. Instead of fully self-consistent GW, typically non-self-consistent G₀W₀ calculations are applied on top of DFT or hybrid functional results. Those may, however, not necessarily improve the accuracy of defect levels, when compared to hybrid functional calculations.⁹⁶ 

Bulk CeO₂ presents a simple fluorite structure, as shown in Fig. 2(a). Its experimental bandgap is 3.2 eV,^97,98 and vacancy formation was computed to require 1.85 eV with respect to gas-phase molecular oxygen.⁵² Upon oxygen removal, two separate Ce atoms are reduced, populating the 4f subshell in between the valence and the conduction band. Consequently, their ionic radii increase. The material can be employed either as a support, for instance, in a three-way car-exhaust catalyst,⁹⁹ or as a catalyst itself, for example, in high temperature HCl oxidation to chlorine¹⁰⁰ and other chlorinated structures,¹⁰¹ or the production of carbamates.¹⁰² 

Many of the appealing catalytic properties of ceria stem from its facile reducibility and innate oxygen transport ability at the surface. Similar to TiO₂,^39,41,103 the surface orientation determines the most likely position of the oxygen vacancy on ceria. In the case of CeO₂(111), the most exposed surface on polycristalline⁵² and rod-shaped¹⁰⁴ nanoparticles, its oxygen vacancies migrate to the subsurface positions, while the electrons reside at the surface cerium centers to avoid nearest neighbor interactions.^105–108 There they affect intrinsic properties including the basicity of the surface and the redox potential, as well as mesoscopic parameters, such as its wettability.¹⁰⁹ Some of these properties can be fine-tuned via doping.⁵² 

Nanocubes, by contrast, preferentially exhibit the polar (100) termination, shown in Fig. 2. Cuts along these directions generate a type 3 surface¹¹⁰ according to the Tasker polarity rules.¹¹¹ They require a substantial reconstruction to become stable so that the top surface layer expels half of its oxygen atoms. Since this process strives toward overall surface neutrality and a minimal dipole moment, no cation reduction occurs as a result. As such, the surface does not present any defect state in the bandgap, which can be seen in Fig. 4. The oxygen depletion of the purely O-terminated surface was originally considered to yield a regular pattern, where rows of top-surface oxygen atoms run parallel and are interspersed with the vacancies that act like trenches.

The high oxygen vacancy mobility in bulk ceria extends to the (100) surface vacancies. The low coordination number along with newly liberated space allows for an increase in the thermal fluctuations of the low-lying soft modes of the oxygen atoms. With the number of particles (O atoms) reduced to half the number of boxes (potential O-positions) and low kinetic barriers, other potential O-configurations become accessible on the reconstructed (100) surface. Ab initio molecular dynamics support this claim, demonstrating that the diffusion starts at temperatures below 500 K, as shown in Ref. 112.

The configuration phase space of these patterns does not have one single low-energy outlier. Rather, all the structures appear with a very small energy spacing, resulting in a densely packed configuration space. For instance, the second lowest configuration in energy presents a zigzag O-pattern and has a surface energy of only 0.01 J/m² higher than the row pattern. Because its multiplicity (the number of identical zigzag structures) is larger than one, it contains a new, configurational entropy term. The surface, thus, has many degrees of freedom in terms of oxygen configurations, and the contribution from their entropies renders the oxygen-termination more stable.¹¹²  In situ HRTEM experiments confirmed the presence of variable terminations.¹¹³ 

Finally, the coexistence of different local geometries also provides a wider array of anchor points for complex reactants to adsorb onto the surface, radically changing the reactivity of this surface. This is the case when stabilizing single metal adatoms.¹¹⁴ The geometric isolation from the other metals on the surface, together with the strong coupling to the electronic structure of the oxide material, provides a new mechanism for the metal/support interactions.

In₂O₃ is an example of intrinsic structural complexity, which is mainly used in optoelectronic devices in combination with Sn (ITO) to improve its electrical conductivity.^53,116,117 It can be a transparent, conductive oxide meaning that it has a high electrical conductivity and that it is nearly transparent in the visible region.⁵³ In₂O₃ is also a selective catalyst for reactions such as propane dehydrogenation to propene,¹¹⁸ acetylene semi-hydrogenation to ethylene,¹¹⁵ and, most importantly, CO₂ hydrogenation to methanol.^54,119,120 It is a more cost-efficient catalyst than the noble metals as the commodity price of indium is 324 USD/kg.¹²¹ 

The bulk solid can appear in cubic (bixbyite),¹²² rhombohedral (corundum),¹²³ and amorphous¹²⁴ forms. The bixbyite, shown in Fig. 2, is formed by a Jahn–Teller deformation of the more regular corundum structure. Yet, it can also be understood as a reconstruction of the much simpler face-centered cubic fluorite structure (CaF₂) by removing one-quarter of the oxygen atoms in a (2 × 2 × 2) supercell.¹¹⁷ The In atoms in this distorted structure appear in five non-equivalent crystallographic sites, each one surrounded by six oxygen atoms.¹²⁵ 

Its most stable surface, (111),^117,126 was studied using the pure PBE functional, after having verified that the U_eff contribution is negligible.¹¹⁵ It is highly anisotropic and can easily lose its oxygen atoms.^127,128 A close look to the surface reveals that it contains four chemically non-equivalent oxygen sites, as shown in Fig. 2.¹¹⁵ In H₂-rich atmospheres, these oxygen atoms can desorb as water, with energies ranging from exothermic (−0.8 eV) to mildly endothermic (+0.4 eV). This reaction is fully reversible as the vacancies heal in the presence of water even at trace concentrations.^119,129

Amidst the anisotropy, the three topmost positions are energetically equivalent. Vacancy diffusion to the neighboring sites shown in Fig. 2 is, therefore, thermoneutral.¹²⁸ All of these vacancies remain confined to the surface¹³⁰ because their formation in the bulk is highly endothermic. Upon formation, the surface is able to accumulate the residual electrons.¹²⁷ They partially localize onto the low-lying s- and p-states of the three exposed In atoms surrounding the vacancy,^115,119 although most of it remains delocalized,¹³⁰ as illustrated in Fig. 4. It has been shown that variations of the synthesis protocol can lead to surface hydroxide formation.^131,132 Doping strategies can modify the number of vacancies accumulated on this surface. For instance, the incorporation of Pd atoms via coprecipitation into the lattice stabilizes the formation of a cluster of vacancies.¹³³ 

α-MoO₃, which is shown in Fig. 2, is used for several chemical and physical applications⁵⁵ as a catalyst and photosensor. Oxygen vacancy formation in this compound is essential for maintaining the catalytic activity of the oxide in important chemical transformations, such as the Formox process (selective methanol oxidation to formaldehyde)^134–137 or hydrodeoxygenation (HDO) reactions.^138–141 These defects are in fact unavoidable, even at high oxygen pressures.¹⁴² The formation of vacancies also influences its physical applications since its optical properties are correlated with the vacancy concentration.¹⁴³ 

The alpha phase of molybdenum trioxide is a semiconductor (3.0–3.2 eV for the experimental bandgap^144–146), and it is formed by bilayers of distorted octahedra (along the xz plane) packed through van der Waals interactions (along the y axis). This is the cleave plane that leads to the lowest energy surface α-MoO₃(010). The structure contains three symmetrically nonequivalent O atoms in its lattice of which the terminal one is the easiest to lose.^81,147

A PBE+U study (with U_eff = 3.5 eV optimized against the hybrid HSE06 functional) found that the vacancy can yield two different electronic structures. The ground state configuration of $VOt$ (2.61 eV for the formation of 1/2 O₂ at a coverage of $θVOt$ = 1/9 ML) corresponds to the reduction of two Mo⁶⁺ centers to Mo⁵⁺ (2 × Mo⁵⁺), the undercoordinated Mo (Mo_cus) and a neighboring Mo center along the [001] direction. The DOS analysis shows that each Mo center has one electron located in a d-orbital, being lower in energy than the conduction band, as shown in Fig. 4. In the second configuration (2.84 eV for the formation of 1/2 O₂ at $θVOt$ = 1/9 ML), Mo_cus is doubly reduced to form a Mo⁴⁺ center. The DOS shows that the electron pair is divided over two different d orbitals of Mo_cus and have the same spin (Fig. 4). The two electronic states were fully characterized showing clear differences in bond distances, Bader charges, XPS displacement, and vibrational frequencies, each usable as fingerprints. Both states and their vacancy formation energies were confirmed by hybrid HSE06 functional, which we used as a benchmark.

The ground state, thus, corresponds to a bipolaron with the two electrons localized at two Mo centers. The metastable state, on the other hand, is defined as a single dipolaron, where both electrons are stored at the same center. The coexistence of these states can be rationalized by the combination of redox and electrostatic properties of the oxide. By the addition of one and two hydrogen atoms to the surface, we found that for Mo⁵⁺, the energy to oxidize it (to Mo⁶⁺) and to reduce it (to Mo⁴⁺) is equal but with opposite sign. Therefore, the redox process to the conversion between the 2 × Mo⁵⁺ and Mo⁴⁺ states is thermoneutral. This is also confirmed by the DOS analysis in Fig. 4 since the vacancy-induced gap electronic states have the same energy. Only the electrostatic interactions between both cations and the surrounding oxygen atoms are responsible for the different stability of both states.⁸¹ 

The predicted electronic states are in agreement with the experimental characterization of α-MoO_3−x surfaces,^11,148,149 where three different reduced Mo centers (six-coordinated Mo⁵⁺, five-coordinated Mo⁵⁺, and five coordinated Mo⁴⁺) were observed. At low reduction levels, Mo⁵⁺ centers predominantly exist, whereas only at high reduction levels, also Mo⁴⁺ centers were observed.^11,138,149 We have studied the variation of vacancy formation energies with $VOt$ coverage for both electronic states. The energy needed to form 2 × Mo⁵⁺ increases linearly with coverage, whereas it remains constant to form Mo⁴⁺, until a coverage of 1/2 ML, where both energies converge. For $θVOt$ = 1/9, the Maxwell–Boltzmann distribution (at 300 °C) predicts 22.0% of Mo⁵⁺ centers and only 0.1% of Mo⁴⁺ centers on the catalyst surface. However, for $θVOt$ = 1/2 ML, more similar distributions of 19.2% Mo⁵⁺ and 15.4% Mo⁴⁺ centers are predicted.

Experimentally, it was found that polarons are highly mobile in α-MoO₃.^150,151 We found that both states are interconnected by electron hopping with an energy barrier of 0.3 eV, which can easily be overcome at typical reaction temperatures. The Born–Oppenheimer molecular dynamics were performed to illustrate the dynamic behavior of electronic states, as shown in Fig. 3. Simulations at 573 K and $θVOt$ = 1/4 ML resulted in 95% and 5% population of 2 × Mo⁵⁺ and Mo⁴⁺ states, respectively.

These two states have an immediate impact when studying the chemical properties of molybdenum. The difference in adsorption energies of methanol and formaldehyde, the reactant and product of the industrial Formox process,¹³⁶ on a defective center of either configuration largely influences the catalytic activity. The adsorption of methanol through its lone pair on the undercoordinated center is independent of the oxidation state. However, formaldehyde also interacts with a neighboring oxygen atom, which results in an adsorption energy of 0.23 eV higher on Mo⁴⁺ centers. Therefore, to analyze and understand the chemistry of these materials, at least those two states are needed to present a comprehensive analysis of the reaction pathways.

A level of complexity is added in multinary oxides, for instance, in oxide perovskites, such as SrTiO₃,^152–154 or complex light absorbers, such as BiVO₄, which we discuss as a case study herein and which is shown in Fig. 2. During the last decade, BiVO₄ has emerged as a promising candidate for photo-electrochemical water oxidation catalyst due to its suitable bandgap and valence edge positions, as well as superior hole mobilities and charge separation efficiencies.^{56,57,155,156} Nonetheless, it suffers from poor electron transport, which is ascribed to small polaron formation self-trapping the electrons.^57,86,157 Oxygen vacancies were found to increase the catalytic activity, which often is attributed to an increase in n-type conductivity typical for many metal oxides.^6,158–161 Being effective n-type dopants, they need to form electronic states at the conduction band minimum or close to it, hence becoming shallow donors. They can, however, also interact with the intrinsic polaronic levels and lead to enhanced polaronic transport, although they form deep states.³⁵ Oxygen vacancies may also have a detrimental effect on the photo-electrochemical behavior since they can act as charge-carrier traps and/or electron–hole recombination centers.^46,162–164 The exact role of oxygen vacancies in BiVO₄ is still under debate.^35,165,166 To a large extent, this is due to artifacts of the employed methods, i.e., differences in the DFT functionals, which largely impact the computed electronic structure of the vacancy.

Oxygen vacancies in the bulk and at the (001) surface have been screened with different density functionals including GGA, GGA+U, and hybrid level of theory.⁹² For the hybrid study, we used the modified PBE0-10 functional that is PBE0 with only 10% of the exact exchange introduced; as for GGA(+U), we used PBEsol and, respectively, PBEsol+U with U_eff = 2.5 eV that has yielded good results in previous studies.^92,165,167 We assessed the electronic structure of the defective oxide, including a detailed comparison of the influence of each functional.

Importantly, we found distinct geometric arrangements of oxygen vacancies in BiVO₄:^92,168 localized vacancies, $VOloc$⁠, with undercoordinated $Vcus4+O3$ and $Bicus2+O7$ units and split vacancies, $VOsplit$⁠, in which the vacancy is quenched by forming aV₂O₇ dimer with an adjacent VO₄ unit. This is due to the flexibility of the lattice (and the ability of V to generate oligomers with oxygen), which allows for facile rotation of the VO₄ tetrahedra and, therefore, breaking and forming of V–O–V bonds between the metallic planes along the [001] direction. Such bridge formation can be viewed as VO₄ orthovanadates condensing to produce polyoxovanadates, which has been widely observed in solution,¹⁶⁹ as well as in solid V₂O₅, where it was found to reduce the cost of $VO$ formation.¹⁷⁰ 

In the bulk, the different types of vacancies have similar formation energies with small functional-dependent differences, whereas $VOsplit$ is largely stabilized (≳1 eV) in the subsurface of the thermodynamically most stable (001) surface. Moreover, the transition from $VOloc$ to $VOsplit$ proceeds without significant barrier when O is removed from the (001) surface.⁹² In the subsurface, the V–O–V formation is so favorable because (i) surface relaxation decreases the distance between the upper two metallic layers and, thus, the distance between the bridge forming V atoms and (ii) the VO₄ tetrahedra that terminate the (001) surface move and rotate more easily than fully coordinated bulk VO₄ units, regarding solely for geometrical reasons. The (001) subsurface, therefore, presents an energy sink for accumulation of oxygen vacancies in the form of V–O–V bridges, which is also supported by recent experiments detecting reduced metallic surface sites by LEED and photoemission spectroscopy.^171,172

Formation of those different structures has tremendous consequences for the electronic properties and, therefore, for the photo-electrocatalytic activity. The resulting electronic configurations are shown in Fig. 4 on the right. In the localized vacancy, $VOloc$⁠, the excess electronic charge is trapped as an electron pair at the vacancy site and at the surrounding Bi and V atoms, therefore reducing Bi³⁺ to Bi²⁺ and V⁵⁺ to V⁴⁺. This gives a doubly occupied electronic state of around 0.9–1.2 eV below the conduction band, which is composed of a Bi 2p and V 3d hybrid orbital. Interestingly, a similar electronic configuration was obtained for all GGA, GGA+U, and hybrid functionals; however, with GGA, the trap state is less deep, lying only 0.6–0.7 eV below the conduction band.

When forming a V–O–V bridge in $VOsplit$⁠, the surplus electrons are no longer bound as a pair but split and delocalized onto different vanadium centers. The extent of delocalization is highly dependent on the functional and local geometry. Considering the hybrid functional PBE0-10 as highest in accuracy, two isoenergetic electronic configurations were obtained, depending on the relevant V–O distances, i.e., whether the input geometry was taken from a prior PBEsol or PBEsol+U optimization. One electron always remains localized in the 3d orbitals (t₂ symmetry) of the bridge forming V with a slightly longer V–O bond, which is reduced to V⁴⁺, forming a deep intragap state at ≥1.2 eV below the conduction band (with hybrid and GGA+U methods; ∼0.7 eV with GGA). The second electron is either delocalized onto different V centers within the same metallic plane or the adjacent planes along the [001] direction when calculated with GGA or hybrid with GGA geometry, or localized onto the other bridging V (V⁴⁺) with GGA+U or hybrid with GGA+U geometry (see Fig. 4). Other electronic configurations, where the second electron is localized on a V atom farther away from the vacancy site and, thus, forms small polarons in the bulk region of BiVO₄, have also been observed.^35,92 The interaction between the two additional electrons in $VOsplit$ is negligible, since both singlet and triplet states with equal energies were found. In the case where the extra electron is two-dimensionally delocalized in the metallic plane (and adjacent planes), an additional electronic state at the bottom of the conduction band (in the bulk) or very close to the conduction band ((001) surface) was found. This would favor band-like conduction via conduction band states. In the case where the second electron localizes and reduces a different V atom, a deeper intragap state was observed and hopping conduction via small polaronic carriers would be expected.^35,92 Apart from the electronic conductivity, also the ionic conductivity is expected to increase through O vacancy migration in BiVO₄ bulk due to the low transition barriers for consecutive forming and breaking of V–O–V bonds.^36,92

Those different geometries and electronic configurations have a large impact not only on the conductive properties but also on the catalytic activity. For example, such an extended delocalized charge state due to oxygen vacancies at the surface has also been found for WO₃ and is predicted to enhance the surface reactivity by promoting the formation of intermediate radical species at the surface.¹⁷³ 

The search for new materials relies on extensive computation of the main properties in different crystalline materials. The realization that databases act as the starting point for these searches prompts the FAIR¹⁷⁴ principles of Findability, Accessibility, Interoperability, and Reusability to ensure a more ergonomic design. The Materials Genome database is a good example of a practical implementation.^14,21,175 Algorithms faced with a large candidate pool usually start by filtering out desired properties based on lower level theory.²³ This allows for a fast comparison of structural and electronic properties of a large set of different materials. Then, more specific filters are applied to this dataset, targeting particular properties, such as bandgaps or valence band levels.^17,20 In doing so, the number of potential materials is reduced, and in a next step, the desired properties are calculated at a higher level of theory for a smaller set of compounds. In this way, the structures can be employed for fast materials screening.^14–16,24 This is a great advance in the search for new materials, but in light of our studies, some subtleties, such as the role of defects and defect states, are only starting to be included.²⁴ 

Our results for the oxides described above that can be consulted in the ioChem-BD database^{168,176–180} show that vacancies can completely change the electronic structure, even for simple binary oxides. In particular, those defects which allow for quasi-degenerate states may be interesting for the desired applications and should therefore be carefully checked. They typically induce cations to show multiple redox pairs, configurations, or dynamicity, among other effects that might be amplified in high-entropy oxides.¹⁸¹ 

All these aspects should be included as filters in the new generation of screening engines in order to reduce the number of false positives, as proposed in Fig. 5: if the pristine material fulfills a given property at a particular level of the screening process, it should not be immediately accepted. Instead, the thermodynamical feasibility to form a vacancy has to be tested (considering a supercell). This defective structure then re-enter the funnel. Then, the possibility of a wide configurational space (many structures) or the formation of new bonds within the solid has to be checked, as it may introduce large differences in the material’s properties. Finally, some metal centers may present multiple oxidation states and therefore lead to several possible electronic configurations, which have to be accounted for. Of course, the location of this new constraint in the list of filters needs to be evaluated with care. Typically, the funnel structure is built so that the larger number of structures are evaluated at a cheaper computational cost. Therefore, the particular location of our new filter will require to take into consideration the computational level needed and the number of configurations to be assessed.

Oxygen vacancies are ubiquitous in oxides, in part because of the thermal treatment employed in their synthesis. Traditionally, the residual electrons left by the O defect were assigned either to the vacancy site, to the metal cations (if those were reducible), or they delocalized in the case of a metallic oxide. In this perspective, by contextualizing a few examples from our group, we demonstrated how this picture is exceedingly simple. For instance, some surfaces are intrinsically defective and do not lead to the formation of electron centers. This is the case for some of the massive surface reconstructions that regulate the surface as with CeO₂ (100). The ultimate consequence of the lack of half of the atoms is that (i) the diffusion barrier of the remaining O atoms substantially decreases and (ii) several surface termination patterns coexist, resulting in a large stabilization due to configurational entropy. The most stable In₂O₃ surface, (111), exhibits a pattern where only some of the oxygen atoms can be removed preferentially, leaving spatially resolved isolated centers. Upon vacancy formation, the excess electrons are partially delocalized between three vacancy-surrounding In cations. In other systems, such as MoO₃, the cation can attain more than one oxidation state easily. The corresponding vacancy defect has two quasi-degenerate energy states that interconvert through polaron dynamics. A final observation regards more complex structures, such as BiVO₄. There, the vacancy can present two different geometries with the formation of a V–O–V bridging unit. As such, it can adapt to different electronic configurations, in which the excess electrons may localize on different V centers or instead delocalize in a two-dimensional extended charge state. The flexibility of the geometry and the redox character of both metals explains the versatility of vacancy structures in these ternary oxide materials.

In summary, the most common defect in oxides, the oxygen vacancy, presents a wide variety of behaviors that affect their electronic structure. In many cases, several near-degenerate electronic structures are needed to fully model the defect state. Even symmetric configurations have been found to break their symmetry under polaron interactions. Therefore, the rational design of new materials via high-throughput screening should also consider structures with one or more vacancies, the possibility of reconstruction, and that some constituents may adopt more than one oxidation state. Thus, any survey of reactivity and physical properties of reducible metal oxides, such as the ones discussed here, will require a wide sampling of the various cationic oxidation states and their locations. This aspect is often still overlooked in contemporary materials screening but will prove to be imperative in order to catch false positives and improve the accuracy of future predictions.

The authors acknowledge the MCIU/AEI/ERDF-EU RTI2018-101394-B-I00, MINECO BES-2016-076361, and AGAUR 2017-SGR-90 projects for financial support, as well as the Barcelona Supercomputing Center (BSC-RES) for providing generous computer resources.