Understanding the role of oxygen-vacancy defects in Cu2O(111) from first-principle calculations

The presence of defects, such as copper and oxygen vacancies, in cuprous oxide films determines their characteristic carrier conductivity and consequently their application as semiconducting systems. There are still open questions on the induced electronic re-distribution, including the formation of polarons. Indeed, to accurately reproduce the structural and electronic properties at the cuprous oxide surface, very large slab models and theoretical approaches that go beyond the standard generalized gradient corrected density functional theory are needed. In this work we investigate oxygen vacancies formed in proximity of a reconstructed Cu2O(111) surface, where the outermost unsaturated copper atoms are removed, thus forming non-stoichiometric surface layers with copper vacancies. We address simultaneously surface and bulk properties by modelling a thick and symmetric slab, to find that hybrid exchange-correlation functionals are needed to describe the oxygen vacancy in this system. Our simulations show that the formation of oxygen vacancies is favoured in the sub-surface layer. Moreover, the oxygen vacancy leads to a splitting and left-shift of the shallow hole states in the gap, which are associated with the deficiency of copper at the surface. These findings suggest that surface electronic structure and reactivity are sensitive to the presence of oxygen vacancies, also when the latter are formed deeper within the film.


Introduction
Due to its availability, non-toxicity, and low cost, cuprous oxide (Cu 2 O) is a well investigated semiconductor material. Thanks to its direct band gap of 2.17 eV [1], it finds applications for photovoltaic conversion in solar cells, photoelectron chemical water splitting [2,3], and other catalytic processes. The high chemical activity of Cu 2 O arises from the Cu + oxidation state, which can be either oxidised to Cu 2+ or reduced to Cu 0 , thereby facilitating a wide range of redox reactions. Further, Cu-based transparent conducting oxides are widely employed in nano-electronics and heterogenous catalysis [4].
The various applications are largely reliant on defects which are often responsible for the electronic performance. In particular, point defects tend to result in the emergence of localized states in the semiconductor band gap, and thereby play a key role on properties like the photoelectronic conductivity [5]. Both experiment [6] and theory [7][8][9] indicate that the defects play a central role in Cu 2 O, however there is no general consensus on the underlying mechanisms on how they work. The experimental evidence suggests that the defect chemistry of Cu 2 O is polaronic in nature [10], with contributions from copper vacancies (V Cu ), oxygen vacancies (V O ) as well as other ionic and electronic defects. In the bulk material, polarons lead to relatively deep hole trap levels [11], with Fermi level pinning in the band gap. This is responsible for the remarkably low open circuit voltages of 0.0-0.4 eV, resulting in poor performance when applied in solar cells [10]. On the other hand, research comparing a series of Cu 2 O facets has found that Cu 2 O(111) presents higher catalytic activity and stability among the low-index facets of cuprous oxide, therefore, the characterization of its surfaces and interfaces has attracted extensive interest both theoretically (ab initio methods) [12][13][14] and experimentally (XPS, LEED TEM) [15][16][17].
According to deep level transient spectroscopy, the acceptor levels caused by V Cu are not localized, but rather form a band in the 0.45-0.55 eV range above the valence band maximum (VBM) [10,18]. Density functional theory (DFT) calculations with the generalized gradient approximation (GGA) [8,9,19], instead, predict V Cu delocalized hole states across the Fermi level, and a quasi-metallic character, rather than a polaronic conducting mechanism for p-type conductivity. The deficiencies of GGA-based methods in properly modelling the polaronic nature of p-type conductivity has been already recognized by Scanlon et al [10], who proposed to apply higher level of theories, such as hybrid functionals to achieve more reliable results. When using a screened hybrid-DFT (Heyd-Scuseria-Ernzerhof (HSE) hybrid functional [20]), distinct single particle levels associated with V Cu defects appear within the band gap, consistent with experimentally known small polaron hopping mechanism. In addition, Isseroff and Carter [21] used DFT-HSE with more numerically stringent parameters and k-point sampling. In doing so, they demonstrated that DFT-HSE reproduces the correct electronic structure of the Cu 2 O cation vacancies and the nature of the hole traps observed in experiments as well.
The presence of V O in transition metal oxides has also an important influence on their semiconducting character. The charge density and atomic structure rearrangement around the vacant site affect the electronic conductivity, thus changing structural, electrical and optical properties as well as the photocatalytic performance [22]. Previous studies find that around each V O site in Cu 2 O the four equivalent copper atoms have unsaturated dangling bonds and act as a source of electron density, thus contributing to n-type character [23]. However, this interpretation has not reached a general consensus [8,19,24]. First-principles calculations have been performed to assess their role as electron-donating defects. Raebiger et al [19] apply DFT-GGA augmented with the Hubbard model (+U) and assert that the oxygen vacancies are unable to act as donors, because no ionization levels within the band gap has been revealed. They also estimate a low defect formation energy under Cu-rich/O-poor conditions, and conclude that V O defects are stable only in the charge neutral state. Also based on GGA-DFT calculations, Soon et al [8] report that the presence of V O introduces an occupied singlet state deep in the valence band and an empty triplet state inside the conduction band, thus resulting in a charge neutral defect, with no ability to compensate produced holes. On the other hand, the hybrid-DFT study of V O defects by Scanlon and Watson [24] found no indication that the oxygen vacancies act as donors in Cu 2 O, because no defect state is located in the band gap. Photoluminescence experiments, however, indicate n-type conductivity in Cu 2 O, with oxygen vacancies acting as donors. Furthermore, a combined experimental and theoretical work by Nandy et al [25] based on GGA-DFT shows that oxygen vacancies can enhance the electron donating ability in Cu 2 O, because of the unshared Cu-3d state which stretch toward the conduction band, forming an impurity energy (donor) state. They also discuss the work function dependence on the V O concentration, which may lead to the Fermi level pinning, and consequent enhancing of the electron tunnelling proficiency of the surfaces. In summary, the available theoretical studies seem to suggest that un-doped n-type conductivity in Cu 2 O is not driven by simple V O defects, but the presence of V O can increase the catalytic ability by decreasing the work function.
Given the sensitivity of predictions for Cu 2 O to the level of theory applied from DFT and the unresolved questions that remain regarding vacancy defects at the surface, a thorough analysis of V O and V Cu in Cu 2 O is desirable. Current experimental findings will be explained by an insightful knowledge of the function of V O , which may be seen theoretically by considering how it operates in semiconductor devices or catalytic reactions situated on the surface or interface of Cu 2 O. It is crucial for further advancements in abundant Cu 2 O-based materials, particularly for the development of p-n heterojunctions for Cu 2 O-based photovoltaic devices [24] and the adsorption and dissociation of gases (O 2 [26], NO [27], CS 2 [28] etc) On photocatalysis, Khasanah et al [29] pointed out that the dominating effect of V O in Cu 2 O remains to be determined, because it exerts opposite effects by facilitating the photocharge recombination and reducing the high reaction overpotential at the interface. Further, Kumar et al [30] demonstrate that due to the suitable band gap, Cu 2 O are widely employed as the photocathode and the photoanode in photoelectrochemical water splitting, and enhance its productivity by tuning nanostructures, such as by introducing defects or dopants. In lithium-ion battery investigation, Cu 2 O-based materials are commonly used as electrodes where V O defects can improve the reaction kinetics by reducing the energy barrier of ion intercalation and accelerating the charge transfer [31]. Meanwhile V O s are native defects in transition metal oxides and their presence has a critical effect on the physicochemical properties of the oxide. The study of defects in Cu 2 O serves as a benchmark for research on other transition metal oxides and unleash their potential in more fields, with better understanding the effects of the presence and diffusion of V O s [32].
In this work, we modelled a large Cu 2 O slab consisting of 8 atomic layers and 736 atoms in the unit cell. This allowed us to model bulk-like and surface properties simultaneously. In order to address the contradictions and uncertainties in V O defects, reliable descriptions of the electronic characteristics of the defects are acquired and accomplished through hybrid DFT. A reconstructed Cu 2 O(111) surface was built by removing all the outermost copper atoms, and forming a non-stoichiometric surface with copper vacancies. We find that at hybrid functional level, the semiconducting features with acceptor levels in the gap are correctly reproduced. Moreover, the oxygen vacancy induces a splitting and left-shift of the in-gap states, which suggest that surface electronic structure and reactivity are sensitive to the presence of oxygen vacancies.

Computational details
All simulations were performed using the Quickstep module of the CP2K program package [33]. The Kohn-Sham DFT [34,35] within the hybrid Gaussian and plane waves framework was applied for the electronic structure calculations, where the interactions with the atomic cores is described through Goedecker-Teter-Hutter (GTH) pseudopotentials [36] and the molecular orbitals of the valence electrons are expanded in Gaussian type orbitals. We used DZVP-MOLOPT-GTH primary basis sets and auxiliary uncontracted basis sets to accelerate the calculation of the electron repulsion integrals for hybrid functionals simulations [37]. For O atoms, cFIT3 auxiliary basis sets were used, and cFIT9 were used for Cu atoms. The plane waves cutoff has been set at 600 Ry. All calculations were periodic in three dimensions, only at Gamma point and spin-polarized (unrestricted Kohn-Sham). We used and compared three different flavours of exchange and correlation schemes: A standard GGA functional, Perdew-Burke-Ernzerhof (PBE) funcational [38], PBE+U [39,40], and the hybrid functional HSE06 [20]. We used a 3 × 3 × 3 supercell (162 atoms For the optimization of geometries the long-range dispersion interactions D3 [41] was used. For the self-consistent field optimization we applied the orbital transformation method with conjugate gradient, converging the energy down to 5 × 10 −7 Hartree. All proposed structures have been fully relaxed at the PBE+U+D3 functional level of theory, with the Broyden-Fletcher-Goldfarb-Shanno [42] scheme and a force threshold of 1 × 10 −3 Hartree/Bohr. At each optimized geometry we refined the electronic structure by an additional calculation using the hybrid functional HSE06. More computational details are provided in the supplementary information (SI) section S.1. The reported illustrations of structures and electronic distributions have been produced with the visualization tools VESTA [43] and VMD [44][45][46][47][48][49][50][51].

Cu 2 O Bulk with V Cu
A first assessment of the computational approach is provided by reproducing some properties of the bulk system, i.e. the electronic structure of pristine Cu 2 O and the formation energy of one V Cu as well as the induced electronic rearrangement. We obtain a bulk energy gap of 0.56 eV using PBE+U (U = 2 eV) and 1.74 eV using HSE06. The band gap of HSE06 is closer to the experimental value (2.17 eV) [1]. Removing one Cu atom from the bulk (figure 1(b)) and relaxing the structure at the PBE+U level, we observe some contraction of the Cu-O and Cu-Cu distances involving the first neighbours of the vacancy site, i.e. the Cu-O distance shortens to 1.84 Å from 1.86 Å, and the Cu-Cu distance shortens to 2.94 Å from 3.03 Å. The PBE+U projected density of states (PDOS), shown in the SI, is not significantly modified by the presence of the vacancy, since no defect states appear within the gap (figure S.5).
By computing the electronic structure of the optimized geometry at the HSE06 level of theory, one unoccupied (hole) spin down state appears within the energy gap, at 0.57 eV above the VBM ( figure 1(b)). From the spatial distribution of this hole state (figure 1(c)), we observe that it is localized mainly on one Cu atom that was coordinated to the same O as the vacant Cu atom, and partially also on two Cu atoms in the vicinity of the defect site. It can be seen that the hole has a Cu 3d character with a small O 2p contribution. These results are in agreement with the previous hybrid functional calculations by Isseroff and Carter [21], who found a single particle state 0.46 eV above the VBM, and Scanlon et al [10], who also discuss a single-particle state 0.52 eV above VBM. The small numerical difference can be attributed to the use of different system size and calculation parameters.
By comparing the band gap values calculated at various theoretical levels-where the hybrid functional predicts a more acceptable band gap-and the electronic characteristics of the defective structure, we were able to demonstrate that HSE06 can produce plausible explanations for defects in related systems that are consistent with both experimental and other theoretical results [10,18,21].   unsaturated oxygen atoms (O u ), a second layer of copper, and a third layer of saturated oxygen atoms (O s ) atoms. The Cu atoms between the two oxygen layers can be classified into two groups, according to their bonding environment, i.e. saturated (Cu s ) and unsaturated (Cu u ) copper atoms. Previous studies [12,23,25,[52][53][54] have considered two possible models for the Cu 2 O(111) surface, the stoichiometric (non-polar) surface and the non-stoichiometric (polar) surface, where the Cu u atoms are removed. Thermodynamic analysis shows that the Cu u poor surface becomes increasingly more stable at higher oxygen chemical potential [53,54]. In agreement with this, DFT predicts a lower surface energy for the reconstructed surface obtained by removing all Cu u ions [16]. Hence, our simulations take as reference a model slab of eight trilayers, where from both surfaces all Cu u atoms have been removed, thus introducing surface Cu vacancies. The lateral dimension is 24.15 × 24.15 Å 2 , corresponding to a 4×4 super cell, where the bulk lattice constant is 4.28 Å, comparing to 4.27 Å in experiment [1]. The resulting structure is illustrated in figure 2 and we refer to it as Cu 2 O(111)-V Cu . We take the Cu 2 O(111)-V Cu model as the reference structure. The resulting slab is no longer stoichiometric, owing to the vacant copper positions. Upon relaxation, both surface layers undergo some structural rearrangements, while the middle trilayers maintain the equilibrium bulk structure, which confirm that the slab thickness is appropriate to describe the surface properties. Panel (a) of figure 2 shows the top layer, where each oxygen atom is coordinated with three copper atoms, forming a pattern of interconnected triangles, where the oxygen atoms sit in the centre. Upon relaxation, the Cu-O and Cu-Cu distances around the unsaturated oxygen become slightly shorter with respect to equilibrium bulk (i.e. 1.83 vs 1.86 Å and 2.88 vs 3.03 Å, respectively). Correspondingly, the Cu-Cu distances around the saturated oxygen are slightly elongated (i.e. 3.19 Å). The structural changes are similar to those in Cu 2 O bulk with a single V Cu .

Cu 2 O(111)-V Cu
The structural relaxation at the surface is combined with the electronic rearrangement. As discussed in previous works [10,11,19,21] and the last subsection, PBE and PBE+U methods fail in predicting the correct band structure and defect energy levels. It should be noted that the PBE exchange-correlation functional is a continuous extrapolation of the exchange-correlated energy of homogeneous electron gas, and in this approximation, self-interaction error (SIE) is inherent in PBE. As a result, the description for insulator and semiconductor systems exhibit flaws, and especially for highly correlated systems involving dand f-electron, as Cu 2 O, the band gaps and band structures are grossly underestimated. Therefore we refine the characterization of the electronic structure of the PBE+U optimized geometries by HSE06 single point energy calculations. In doing so, our work also confirms that different electronic properties are predicted with PBE+U and HSE06 levels of theory for the defective Cu 2 O(111) surfaces. Specifically, HSE06 predicts more localized states than PBE+U. The Mulliken population analysis shows that the Cu atoms at the surface become slightly more positive with respect to bulk (+0.32 vs +0.26 with our settings). Accordingly, the oxygen atoms within the surface trilayer acquire some additional negative charge (see figure 2). Figures 3(a) and (b) display the PDOS computed at the PBE+U and HSE06 levels of theory (see also SI figure S.6). In both cases we observe the appearance of one feature within the gap. These are the 16 acceptor states (empty) in each of the spin channel, corresponding to the 16 V Cu at the surface. These electronic (hole) states are located at the surface and subsurface atoms and not in the bulk region of the slab, as shown by the spatial distributions depicted by isosurfaces in figures 3(c) and (d). PBE+U places the hole states just above the valence band, between the VBM and VBM+0.5 eV, such that the defect band crosses the Fermi level, resulting in a notionally metallic distribution. The spatial distribution ( figure 3(c)) involves all Cu atoms around each V Cu vacancy (six Cu on the surface, and three Cu on the subsurface). On the other hand, HSE06 predicts defect states with mainly Cu 3d character, and only small contributions from O 2p states. The spatial distribution also suggests major contributions on the surface and subsurface copper atoms, though, not all Cu atoms around V Cu are involved. Only five out of nine neighbouring Cu atoms contribute to the defect state, in disagreement with PBE+U. Moreover, according to HSE06, the defect states are acceptor states in the middle of the energy gap of the bulk, at about 0.8 eV from VBM. This is consistent with the experimental findings, revealing states in the range 0.12-0.70 eV above Fermi level, also associated to copper vacancies [10,18,55,56]. These results assess the reliability of our model for the characterization of the Cu 2 O(111) surface and confirm that using the very large simulation cell and HSE06, we can correctly reproduce the properties of defective Cu 2 O for both bulk and surface Cu vacancies.

Oxygen-vacancy defect
Single neutral oxygen vacancies have been introduced by removing one oxygen at three different sites: two of them within the surface trilayer, and one within the bulk region in the middle layer. Since placing one vacancy at the surface generates an asymmetric slab structure, which could lead to some artifact, we also considered the symmetric case, where two almost equivalent vacancies are placed at the two opposite surfaces. The neutral oxygen vacancy formation energy is evaluated as: where the E slab-O is the total energy of the defective slab, E O2 is the energy of free-standing O 2 molecule, and E slab is the energy of the Cu 2 O(111)-V Cu slab. N VO is the number of defects, i.e. one for the single V O and two for the symmetric slab with two surface vacancies. As illustrated in figure 4, we use the following notation: V u O is the vacancy formed at the site of an unsaturated oxygen (surface layer), V s O refers to the vacancy within the surface sublayer, i.e. at the place of one O s , and V b O is the vacancy in the bulk region.

Atomic structures
The structural rearrangements occurring upon full relaxation of the system with PBE+U+D3 are limited to the local region around the vacancy sites. In particular, the nearest neighbour Cu atoms that become unsaturated once the O atom is removed, move towards the vacant site (see figure 4 and SI   while the Cu-O distances to the next neighbour oxygen atoms are elongated. This results in a four-membered Cu-cluster with unshared 3d states and partially reducing the electro-positivity of the cations according to Mulliken charges, which decrease by roughly 0.10 on each member of this Cu-cluster.
Farther from the site of the vacancy no significant structural or charge rearrangement is observed. table 1 lists the formation energies calculated for the three distinct V O defects. The HSE06 energies are computed on the PBE+U+D3 optimized structure. V u O and V b O , turn out to be significantly less stable than the vacancy introduced in the sublayer. In the case of V u O this observation can be rationalized considering the much reduced relaxation around the vacancy site located at the very surface. In this region the structure is already strongly modified due to removal of the unsutured Cu atoms, therefore there is not enough flexibility to compensate for the missing oxygen. On the other hand, upon the formation of the vacancy within the bulk the energy drop along the relaxation is also large, 3.8 eV, due to the contraction of the first neighbour shell. However, there is a higher initial energy cost to removing an oxygen atom from the bulk-like region than from the subsurface and therefore, even after structural relaxation, the subsurface oxygen vacancy is the most thermodynamically favoured site.

V O formation energy
We also testes the case of an symmetric model when the vacancy is created close to the surface, by placing a second vacancy on the opposite site of the slab at an almost equivalent site. The new structures have been relaxed with PBE+U+D3, and the electronic properties have been refined at the HSE06 level. The resulting geometric and electronic properties, as well as the formation energies turn out to differ only marginally from those obtained with the asymmetric slab.

Electronic structures
The HSE06 DOS projected onto the topmost trilayer of the Cu 2 O(111)-V Cu slab containing one of the described vacancy defects are reported in figure 5. We observe that the presence of V O in the structure causes detectable modifications in the region of the VBM and in the in-gap states, associated to the surface Cu vacancies. The vacancy placed at the unsaturated oxygen site results in a slight shift in energy of the VBM and the in-gap states, probably owing to the reduced possibility of rearranging the electronic density among the neighbours. A larger effect is instead observed in the case of V s O , i.e. when the V O is located at the sub-surface site, which results in a reduction of about 0.1 eV of the gap between the VBM and the in-gap states. In particular we observe that some additional Cu d-states get occupied at the VBM. The spatial distribution of these states shows that they are localized on some Cu atoms of the sub-surface layer, just below an unsaturated oxygen site, as shown in figure 5. Hence, the excess negative charge resulting from the removal of one oxygen partially occupies the empty states resulting from the formation of the surface Cu vacancies, thus forming a polaron.

Conclusions
In this work we characterized Cu and O vacancies formed at Cu 2 O(111) by means of DFT calculations at the HSE06 level of theory applied to an unusually large slab model system. It has been confirmed that PBE and PBE+U give conflicting predictions for electronic properties of the defects even for this comparably large Cu 2 O(111) slab. Similar findings have been reported previously, showing that due to the SIE PBE fails to accurately characterize positive defects and predicts a semimetallic material instead of a semiconductor. HSE06, instead predicts localized in-gap defect states, which are associated to the removal of the unsaturated Cu atoms. The formation of one oxygen vacancy, both in the bulk or at the surface, induces a local relaxation around the defect corresponding to contraction of the Cu-atom coordination shell. No additional defect states are revealed within the gap. However, when the vacancy is placed in the sub-surface layer by removing a surface saturated oxygen, which is also the most favourable site for the formation of this defect, we could reveal the localization of a polaron. The polaron is displaced with respect to the vacant site, but still within the subsurface layer, and it has a Cu-d character. We conclude that HSE06 applied to the proposed slab model predict polaronic properties associated with the combined presence of V O and V Cu defects in Cu 2 O. The results suggest that the neutral V O defect induces the hole states splitting, which could enhance the electron mobility and boost the electronic activity at the surface, i.e. V O in p-type Cu 2 O (along with V Cu ) semiconductor can improve electron donating ability and its reactivity. In principle, the polaron can enhance conductivity as it can hop from Cu to Cu site, while exploiting the presence of the in-gap empty states and vibrational dynamics at finite temperature. More interestingly, it can potentially have an effect on the catalytic activity of the material. However, these are just speculation based on the behaviour of polarons in other oxides (see TiO 2 [57,58]), since we have neither carried out simulations of the dynamics of the polaron nor on the effects on adsorbates. These could be topics of future work on this system.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files). All calculation files will be available on the Materials Cloud (materialscloud.org) soon.