Doping dependence of the surface phase stability of polar O-terminated (000) ZnO

In this section, we briefly summarize how we apply ab initio thermodynamics to determine the stability of the various surface structures when in contact with an atomic reservoir, i.e. a surrounding gas. For the sake of brevity, we refer to the supplementary material, which is available online at stacks.iop.org/NJP/19/083012/mmedia, or [41] for further details.

When a surface is in thermodynamic equilibrium with an external reservoir at temperature T and pressure p, the most stable structure and composition of the surface is the one that minimizes the Gibbs free surface energy $\gamma (T,p)$. As we are interested in the relative surface stability, the thermodynamic quantity of interest is the change in surface free energy ${\rm{\Delta }}\gamma$ with respect to a given reference structure. We assume that Zn is always in thermal equilibrium with the surface, and restrict the problem to the two independent gas phase reservoirs H₂ and O₂. The corresponding chemical potentials of the reservoir gases ${\rm{\Delta }}{\mu }_{{\rm{H}}}$ and ${\rm{\Delta }}{\mu }_{{\rm{O}}}$ are referenced to the total energies of the isolated H₂ (${E}_{\mathrm{gas}}^{{{\rm{H}}}_{2}}$) and O₂ (${E}_{\mathrm{gas}}^{{{\rm{O}}}_{2}}$) molecules:

$$\begin{eqnarray}&&{\rm{\Delta }}{\mu }_{{\rm{H}}}={\mu }_{{\rm{H}}}-\displaystyle \frac{1}{2}{E}_{\mathrm{gas}}^{{{\rm{H}}}_{2}},\,\,\,\,\,\,{\rm{\Delta }}{\mu }_{{\rm{O}}}={\mu }_{{\rm{O}}}-\displaystyle \frac{1}{2}{E}_{\mathrm{gas}}^{{{\rm{O}}}_{2}}.\end{eqnarray} \tag{ 1 }$$

More than two external reservoirs, i.e. gas components, can be considered in principle (see e.g. [42]). We refrain from this here so as not to unnecessarily complicate the discussion of the main focus of our manuscript—the dependence on the bulk doping concentration.

The change of surface free energy per unit area for a certain surface structure with respect to the chemical potentials of the reservoir gases ${\rm{\Delta }}{\mu }_{{\rm{H}}}$ and ${\rm{\Delta }}{\mu }_{{\rm{O}}}$ is then given as

$$\begin{eqnarray}&&{\rm{\Delta }}\gamma =\displaystyle \frac{1}{A}\left({E}_{\mathrm{slab}}^{}-{E}_{\mathrm{slab}}^{\mathrm{ref}}+{n}_{{\rm{H}}}^{\mathrm{rem}}\left[{\rm{\Delta }}{\mu }_{{\rm{H}}}+\displaystyle \frac{1}{2}{E}_{\mathrm{gas}}^{{{\rm{H}}}_{2}}\right]+{n}_{{\rm{O}}}^{\mathrm{vac}}\left[{\rm{\Delta }}{\mu }_{{\rm{O}}}+\displaystyle \frac{1}{2}{E}_{\mathrm{gas}}^{{{\rm{O}}}_{2}}\right]\right)\end{eqnarray} \tag{ 2 }$$

where ${E}_{\mathrm{slab}}^{}$ is the total energy of the surface slab with ${n}_{{\rm{H}}}^{\mathrm{rem}}$ removed hydrogens and ${n}_{{\rm{O}}}^{\mathrm{vac}}$ oxygen vacancies present at the surface; ${E}_{\mathrm{slab}}^{\mathrm{ref}}$ is the energy of the reference surface. Taking into account all the considered surface structures, the most stable one for a specific set of chemical potentials $({\rm{\Delta }}{\mu }_{{\rm{H}}},{\rm{\Delta }}{\mu }_{{\rm{O}}})$ is the one with the lowest ${\rm{\Delta }}\gamma$. With this knowledge, we can build two-dimensional phase diagrams with the axes being the two chemical potentials. Since the energy differences between the surfaces can be small, it is crucial to correctly calculate the total energies of the investigated surface structures. Therefore, it is essential to account for the energy contributions associated with the charge transfer or band bending, which we can now do with the CREST approach.

The CREST approach was introduced by Sinai et al in [34]. In this section, we review the main idea and explain our extensions, which allow us to use a smaller slab to facilitate the calculation of the total energies. For ease of understanding, we explain the application of CREST with the example of the ZnO (000) surface. At this surface, the oxygen atoms are only coordinated three-fold, rather than four-fold as in the bulk. Hence, they provide an extra half an electron each, or one extra electron for every two surface O atoms. Under typical ultra-high vacuum (UHV) conditions, the surface reacts with hydrogen, which adsorbs on every other surface O atom to form the $(2\times 1)$-H surface structure shown in figure 1. This structure satisfies the electron counting rule and exhibits flat band conditions, as shown in figure 2(a). When individual hydrogen atoms are removed, they expose a dangling bond with an associated surface state in the ZnO band gap. Due to the intrinsic n-type doping of ZnO, the Fermi level is close to the conduction band minimum. Electrons from the bulk can therefore lower their energy by filling this surface state. The migration of electrons to the surface results in a net surface charge and leaves behind a depletion of charge carriers in the region adjacent to it. The depletion region is also referred to as a space charge region, and depending on the doping concentration, can reach macroscopically far into the crystal. Electrostatically, this leads to the formation of band bending, which shifts the energy levels near the surface relative to their bulk positions until they are in resonance with the bulk Fermi level, as shown in figure 2(b).

Zoom In Zoom Out Reset image size

With the CREST approach, band bending can be included in the first-principles calculations without having to treat the entire space charge region with first principles [34]. The concept of CREST is to only treat the near-surface region explicitly (e.g. by means of DFT), and use a classical electrostatic approach for the rest of the system. This is conceptually shown by the shaded region in figure 3(a). In the slab representing the surface, the free charge carriers are modeled using the virtual crystal approximation (VCA) [43] (see section 3 for details on how the VCA is employed in this work). The potential drop in the classical, electrostatic region that arises from ionized dopants can be adequately captured by homogeneously distributed point charges. However, for technical and numerical reasons, the homogeneously distributed charge is collapsed and is implemented as a two-dimensional charge sheet on the back side of the semiconducting slab, as shown in figure 3(b). Thus, the parabolic potential originating from the remaining ions in the space charge region is replaced by a linear drop of the potential between the bottom side of the slab and the charge sheet. Compensating charge-carriers are introduced in the DFT calculations of the surface slab, rendering the whole system effectively charge neutral.

Zoom In Zoom Out Reset image size

Sinai et al determined the position and charge of the charge sheet by fitting the plane-averaged electrostatic potential within the slab and extrapolating into the vacuum region below. However, when the slab thickness is reduced (e.g. to facilitate the use of computationally expensive hybrid functionals) such a fit might become unreliable. We therefore present an alternative, simplified approach to define the properties of the charge sheet, which is crucial for correctly accounting for band bending effects. Although we will explain our version of CREST using the example of n-type doped ZnO, we emphasize that our approach is completely general.

The major parameter of our CREST method is the work function of the 'bottom' side of the slab. In the slab approach, one of the essential convergence parameters is the slab thickness. Adsorption or desorption on one side of the slab should not affect the properties on the other side of the same slab. It is thus sensible to define the vacuum level on the 'bottom' side of the slab as our energy zero. Furthermore, it is clear that the 'bottom work function' ${\rm{\Phi }}{{\prime} }_{{\rm{L}}}$ (i.e. the difference between the Fermi energy and the vacuum level) must not be affected by any change to the upper side of the slab. This fact can be exploited to determine the electrostatic potential drop over the space charge region and calculate the charge and position of the charge sheet. We can therefore define a target bottom work function ${{\rm{\Phi }}}_{\mathrm{target}}$ for an ideal system that does not exhibit band bending (as shown in figure 2(a)). This target will guide us in the modified CREST scheme.

In an ideal system (e.g. the ZnO $(2\times 1)$-H surface before hydrogen desorption), where doping is modeled using the VCA, the Fermi-level is well above the mid-gap and no band bending occurs. When a hydrogen atom is now desorbed, a new acceptor-like state forms in the band gap. This state may be capable of holding one full electron. Since the state is below the Fermi energy, all free charge carriers will flow into this state. However, unless the slab is very thick, it will typically contain less than one full electron. As a result, the state becomes fractionally occupied, while the conduction band in the finite slab is completely depleted of electrons. Furthermore, the Fermi energy of the entire system (which is calculated from the electron distribution) is now in resonance with the gap state, which results in a change in the work function of the lower side of the slab ${\rm{\Phi }}{{\prime} }_{{\rm{L}}}$. This change indicates a deviation from the target value ${{\rm{\Phi }}}_{\mathrm{target}}$ and requires action.

The difference between ${\rm{\Phi }}{{\prime} }_{{\rm{L}}}$ and ${{\rm{\Phi }}}_{\mathrm{target}}$ is equal to the electrostatic potential drop associated with band bending outside the quantum-mechanical region, i.e. outside the slab (see figure 3(a)). We refer to this part of the band bending as the tail ${\rm{\Delta }}{\varphi }_{\mathrm{tail}}^{\mathrm{BB}}$,

$$\begin{eqnarray}&&{\rm{\Delta }}{\varphi }_{\mathrm{tail}}^{\mathrm{BB}}={\rm{\Phi }}{{\prime} }_{{\rm{L}}}-{{\rm{\Phi }}}_{\mathrm{target}}.\end{eqnarray} \tag{ 3 }$$

${\rm{\Delta }}{\varphi }_{\mathrm{tail}}^{\mathrm{BB}}$ will later be used to determine the position and total charge of the charge sheet that mimics band bending effects. This is where our extended CREST deviates from the procedure of Sinai et al [34]. Instead of fitting the electrostatic potential at the boundary of the slab we use the fact that we know the amount of band bending outside the slab and apply a self-consistent scheme to adjust the properties of the charge sheet to reproduce this amount. A flow chart covering all steps within our modified CREST approach is shown in figure 4. In the next paragraphs, we will derive the mathematical relations to define the correct properties of the charge sheet and explain the self-consistent scheme in more detail. However, for the sake of understanding, we first explain what effect the inclusion of the charge sheet has on the bottom side work function of the slab.

Zoom In Zoom Out Reset image size

By including a charge sheet in the DFT calculation, we will change the bottom side work function ${{\rm{\Phi }}}_{{\rm{L}},n}$, because the electric field between the charge sheet and the bottom side of the slab results in a potential drop, which affects the bottom side work function. Our objective in CREST is now to adjust the properties of the charge sheet iteratively (the subscript n denotes the iteration) such that ${{\rm{\Phi }}}_{{\rm{L}},n}={{\rm{\Phi }}}_{\mathrm{target}}$. If this is the case, the potential drop due to band bending would be correctly reproduced. We can interpret this procedure as finding the right amount of charge that has to be transferred to the surface, such that the distance of the conduction band to the globally defined Fermi level is correctly described. The slab states are now offset with respect to the lower vacuum level by the correct band bending, and the right amount of mobile electrons is included in the calculation.

To explain the mathematical details of our modified CREST method determining the charge and position of the charge sheet, it is useful to briefly recap the classical electrostatic equations associated with band bending in extended space charge layers. Let us consider an infinite surface parallel to the x-y plane with the semiconductor extending in the z direction. z_d denotes the extent of the space charge region in the bulk starting from the surface at z = 0. In the depletion approximation, we assume all donors to be ionized in the space charge region, which gives a space charge density of $\rho ={{eN}}_{{\rm{D}}}$, with N_D being the doping concentration. The resulting electrostatic potential associated with the electric field in the space charge region is determined by the one-dimensional Poisson equation,

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}}^{2}\varphi (z)}{{\rm{d}}{z}^{2}}=-\displaystyle \frac{\rho (z)}{{\epsilon }_{r}{\epsilon }_{0}},\end{eqnarray} \tag{ 4 }$$

where $\varphi (z)$ is the electrostatic potential in the direction of the surface normal. The vacuum permittivity is given by ${\epsilon }_{0}$, and ${\epsilon }_{r}$ is the relative permittivity of the material. With the boundary condition that the electric field vanishes at the edges of the space charge region, i.e. at z = 0 and at z_d, the solution of the Poisson equation for the electrostatic potential yields

$$\begin{eqnarray}&&{\varphi }_{{\rm{S}}}(z)=-\displaystyle \frac{{{eN}}_{{\rm{D}}}}{2{\epsilon }_{r}{\epsilon }_{0}}{(z-{z}_{d})}^{2}.\end{eqnarray} \tag{ 5 }$$

The bands bend parabolically, with a total potential drop over the space charge region of

$$\begin{eqnarray}&&{\rm{\Delta }}{\varphi }_{{\rm{S}}}={\varphi }_{{\rm{S}}}(0)-{\varphi }_{{\rm{S}}}({z}_{d})=-\displaystyle \frac{{Q}_{{\rm{S}}}{z}_{d}}{2{\epsilon }_{r}{\epsilon }_{0}A};\end{eqnarray} \tag{ 6 }$$

here we replaced the charge density $\rho ={{eN}}_{{\rm{D}}}$ by the total amount of charge of the ionized donors ${Q}_{{\rm{S}}}$ divided by the volume of the space charge region ${V}_{\mathrm{SCR}}=A\cdot {z}_{d}$, where A is the base area of the slab. These equations describe the situation for homogeneously distributed dopants. In the next paragraph, we derive the equations for a system with a charge sheet as used in the CREST method.

For technical reasons, in the CREST method it is beneficial to collapse the homogeneous doping region into a single charge sheet, positioned within the space charge region. This charge sheet, combined with the quantum mechanical region, can be viewed as an auxiliary system. The properties of this auxiliary system will be defined in such a way that the original situation of the homogeneous doping region is described as accurately as possible. One difference between the original and the auxiliary system is the difference in relative permittivities. We position the charge sheet in the vacuum region on the back side of the slab calculation, and not in a medium with permittivity ${\epsilon }_{r}$. The slope of the potential therefore changes abruptly at the interface between the back side of the slab and the vacuum region due to the change in the relative permittivity from dielectric to vacuum. Within the slab, the relative permittivity of the material ${\epsilon }_{r}$ is given by DFT. It would also be desirable to set ${\epsilon }_{r}$ in the vacuum region between the lower slab surface and the charge sheet, which, however, has not been possible within most common DFT codes so far. Therefore, we correct for the influence of the relative permittivity in our equations, which determine the charge and distance of the sheet by considering the ${\epsilon }_{r}$ of the material in the equations of the homogeneous charge density, and setting it to the vacuum value in the following equations of the auxiliary system. With this, we introduce ${\epsilon }_{r}$ as an input parameter for our CREST (see figure 4 for reference) and ensure that the potential drop and the transferred charge is correctly reproduced.

Mathematically, the electrostatic environment of the auxiliary system can by described again by the solution of the Poisson equation. The constant charge density is now replaced by a delta function at the position d of the charge sheet

$$\begin{eqnarray}&&\rho (z)=\displaystyle \frac{{Q}_{{\rm{C}}}}{A}\cdot \delta (z-d).\end{eqnarray} \tag{ 7 }$$

The two parameters that determine our auxiliary system are the surface charge density of our charge sheet ${Q}_{{\rm{C}}}/A$ and the distance of the sheet from the back side of the slab d. The solution of the Poisson equation returns a linear drop in the potential between the back side of the slab at z = 0 and the position of the charge sheet at z = d

$$\begin{eqnarray}&&{\varphi }_{{\rm{C}}}(z)=\displaystyle \frac{{Q}_{{\rm{C}}}}{{\epsilon }_{0}}\displaystyle \frac{z-d}{A}.\end{eqnarray} \tag{ 8 }$$

To ensure that our auxiliary system describes the intended system correctly, it has to fulfill two requirements to fix the two adjustable parameters of the charge sheet: first, the potential difference due to band bending ${\rm{\Delta }}{\varphi }^{\mathrm{BB}}={\rm{\Delta }}{\varphi }_{S}$ in the true system is reproduced by the auxiliary system

$$\begin{eqnarray}&&{\rm{\Delta }}{\varphi }_{{\rm{S}}}\mathop{=}\limits^{!}{\rm{\Delta }}{\varphi }_{{\rm{C}}}={\varphi }_{{\rm{C}}}(0)-{\varphi }_{{\rm{C}}}(d).\end{eqnarray} \tag{ 9 }$$

Second, the electric field acting on the slab is correctly described. The second criterion is fulfilled, if the charge sheet contains the same amount of charge ${Q}_{{\rm{C}}}$ as the total charge of the volume of the space charge region ${Q}_{{\rm{C}}}\mathop{=}\limits^{!}{Q}_{{\rm{S}}}=A\cdot {z}_{d}\cdot {{eN}}_{{\rm{D}}}$. Equation (9) for condition one then determines the second free parameter, which is the distance d of the charge sheet. Inserting the potential drop for the homogeneous charge distribution (6), and for the auxiliary CREST system (8) into equation (9) yields $d=\tfrac{{z}_{d}}{2{\epsilon }_{r}}$. This ensures that the potential difference associated with the band bending is correctly reproduced.

Finally, to obtain the correct total energy in CREST we need to correct for the incorrect potential shape that was introduced by collapsing the space charge region into a charge sheet. We subtract the energy of the effective plate capacitor ${U}_{\mathrm{sheet}}$ introduced by the charge sheet and add the contribution for the parabolic shape ${U}_{\mathrm{SCR}}$ of the desired homogeneous doping distribution. This gives a total energy correction of

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{\mathrm{tot}}=-{U}_{\mathrm{sheet}}+{U}_{\mathrm{SCR}}=\left(-\displaystyle \frac{1}{2}+\displaystyle \frac{1}{3}\right){\rm{\Delta }}{\varphi }^{\mathrm{BB}}\cdot {Q}_{{\rm{C}}}=-\displaystyle \frac{1}{6}{\rm{\Delta }}{\varphi }^{\mathrm{BB}}\cdot {Q}_{{\rm{C}}}.\end{eqnarray} \tag{ 10 }$$

This correction only depends on the potential drop due to band bending ${\rm{\Delta }}{\varphi }^{\mathrm{BB}}$ and the charge within the space charge region, which is equal to the charge of the introduced sheet Q_C. The total energy of the slab system, which enters equation (2) for the phase diagram calculations, can now be written as ${E}_{\mathrm{slab}}={E}_{\mathrm{slab}}^{\mathrm{DFT}}+{\rm{\Delta }}{E}_{\mathrm{tot}}$. Additional information about the implementation of our modified CREST approach for FHI-aims can be found in the supplementary material.

The CREST method explained in the previous section allows us to study the influence of the doping concentration on the stability of semiconductor surfaces. As a first step, we explore the defect-free O-terminated (000) ZnO surface and the influence of the doping concentration on the hydrogen termination. The $(2\times 1)$-H surface with a half-monolayer of adsorbed hydrogen, shown in figure 1, completely compensates the charge from the dangling bonds on the surface and therefore has a very stable surface structure. This means the partially occupied surface bands are filled by the adsorption of hydrogen on every second oxygen atom. Alternatively, free charge carriers introduced by dopants fill and thus passivate the surface gap states. This alternative passivation introduces an energy penalty, which results from the associated build-up of band bending. The question is whether the bare, bulk- truncated O-terminated (000) ZnO surface could ever be stabilized by such free charge carriers. To investigate this question we gradually remove hydrogen from the $(2\times 1)$-H surface for different doping concentrations. We performed our calculations for a 4 × 4 supercell, which implies that the hydrogen coverage is reduced by $\approx 6$% for every hydrogen atom that is removed.

Naively, one may expect the hydrogen coverage to reduce significantly for high doping concentrations, in accordance with the results reported by Moll et al [33] (who did not include band bending). However, figure 5, in which the coverage of the lowest energy structure is plotted as a function of the hydrogen chemical potential and doping concentration, shows that this is not necessarily the case. Converting the hydrogen chemical potential into the temperature and partial pressure of the hydrogen gas phase shows that for natively doped ZnO (i.e. approximately ${10}^{17}\,$ cm^–3), the unperturbed $(2\times 1)$-H surface structure is still thermodynamically preferred for most realistic experimental conditions. Only at much higher doping concentrations or at very high temperatures do we find a deviation from the ideal structure. For UHV conditions and at 600 K, which is the temperature typically used for sample annealing, reduction to a surface coverage of 44% occurs for a doping concentration of $\approx 5\times {10}^{19}\,\,\,{\mathrm{cm}}^{-3}$. At lower temperatures, the transition would require an even higher number of free charge carriers. The stability of the $(2\times 1)$-H surface against hydrogen removal can be explained by the cost of band bending. To compensate for the desorption of hydrogen from the surface requires a certain amount of charge to be transferred from the bulk to the surface state that is generated in the desorption process. However, charge transfer from the bulk to the surface states is limited by the electrostatic field originating from band bending. This energy cost stabilizes surface structures close to 50% coverage under realistic conditions and illustrates the importance of including band bending in the surface calculations. At these conditions a hydrogen reduction of $6 \%$ is expected at a doping concentration of approximately ${10}^{20}$ cm⁻³. For this case, the generated surface state becomes filled by about 0.7 electrons. For lower doping concentrations, this critical value cannot be reached because band bending limits the charge transfer. However, one needs to consider that because of the limited size of the unit cell, we can only model hydrogen reductions in increments of 6%. With this in mind, we can argue that diminishing H-coverage of $\lt 6$% is to be expected already at lower doping concentrations. So far our findings cannot be compared to experimental results, because to our knowledge there is no experimental study that systematically investigates ZnO surfaces depending on the doping concentration.

Zoom In Zoom Out Reset image size

Oxygen vacancies at or near the O-terminated ZnO (000) surface have been discussed in the literature [12, 51], but from ab inito thermodynamics they are not considered to be stable in thermal equilibrium under realistic conditions. Still, they have been observed in experiments and might form in reality as a result of physical treatment during surface preparation procedures, which include sputtering and annealing steps. Therefore, we take a closer look at the stability of the ZnO surface including oxygen vacancies with the CREST approach. Introducing one oxygen vacancy in our 4 × 4 unit cell corresponds to an oxygen vacancy concentration of $1/16\approx 6$%. Figure 6 shows the doping dependency of the hydrogen coverage in this case.

Zoom In Zoom Out Reset image size

One can see that there are now two competing surface structures in UHV conditions. For conditions that are richer in hydrogen, the (2 × 1)-H structure is still the most favorable. However, at smaller chemical potentials the most stable surface structure for low doping concentrations is now the one with two fewer hydrogens in the unit cell. This can easily be understood, since a structure with an oxygen vacancy and two hydrogens missing again fulfills the electron counting rule by completely filling the surface band. Thus no charge transfer from the bulk is required to saturate a dangling bond at the surface. Again, an increase in the doping concentration leads to a decrease in the hydrogen coverage. However, the desorption of hydrogen sets in at lower doping concentrations than for the defect free surface. While for the defect-free surface, hydrogen desorption is only observed for doping concentrations above ${10}^{20}$ cm⁻³, oxygen vacancies provide additional charge carriers and lower this critical value to approximately ${10}^{19}$ cm⁻³. This is because the oxygen vacancy pins the Fermi level at a higher energy level, a concept we recently explored for self-assembled organic monolayers at ZnO surfaces [58], resulting in less band bending and a spatially smaller space charge region. The schematic band diagram explaining the position of the generated surface states and Fermi level pinning can be found in figure 7. Still, as in the defect-free case, band bending limits the charge transfer to the surface and so hinders a further decrease of its H concentration.

Zoom In Zoom Out Reset image size

Looking at the two phase diagrams with and without an oxygen vacancy, the question arises as to when the surface with an oxygen vacancy concentration of 1/16 becomes stable? To answer this question, we combine the results of the previous two sections and assume that the surface is simultaneously in equilibrium with two independent H₂ and O₂ gas reservoirs. The change in surface free energy now depends on the two chemical potentials of the surrounding gases and the doping concentration adds a third dimension. Such a three-dimensional plot would be hard to follow and therefore we representatively plot the phase diagrams for the two doping concentrations ${N}_{{\rm{D}}}={10}^{18}$ cm⁻³ and ${N}_{{\rm{D}}}={10}^{20}$ cm⁻³ as a function of the chemical potentials ${\rm{\Delta }}{\mu }_{{\rm{H}}}$ and ${\rm{\Delta }}{\mu }_{{\rm{O}}}$. Additionally, we added the corresponding partial pressures at a temperature of 600 K to the opposing axes of the hydrogen and oxygen chemical potentials, respectively.

Figure 8 shows that for natively doped ZnO, the $(2\times 1)$-H surface structure is still the thermodynamically preferred surface structure under UHV conditions at 600 K. This means that we expect surfaces with oxygen vacancies and reduced hydrogen coverage to only become stable in extremely oxygen poor conditions or at very high temperatures. The doping concentration does have an influence on the hydrogen coverage once an oxygen vacancy is present, but does not distinctively influence the stability of the vacancy itself. So, how can we justify the fact that it is still relevant to consider surface oxygen vacancies? First, an oxygen vacancy with a fifty percent coverage of hydrogen results in an empty gap state at the surface, and therefore might be the cause of Fermi level pinning, which is seen experimentally; and second, vacuum chambers are not a thermodynamic equilibrium. In the majority of surface preparation techniques, the ZnO polar surface is sputtered at least once with Ar${}^{+}$, for example. This sputtering could create oxygen vacancies that remain as metastable structures, even though they are, in principle, thermodynamically unstable.

Zoom In Zoom Out Reset image size