Morphological Features and Band Bending at Nonpolar Surfaces of ZnO

: We employ hybrid density functional calculations to analyze the structure and stability of the (101 ̅ 0) and (112 ̅ 0) ZnO surfaces, con ﬁ rming the relative stability of the two surfaces. We then examine morphological features, including steps, dimer vacancies, and grooves, at the main nonpolar ZnO surface using density functional methods. Calculations explain why steps are common on the (101 ̅ 0) surface even at room temperature, as seen in experiment. The surface structure established has been used to obtain the de ﬁ nitive ionization potential and electron a ﬃ nity of ZnO in good agreement with experiment. The band bending across the surface is analyzed by the decomposition of the density of states for each atomic layer. The upward surface band bending at the (101 ̅ 0) surface a ﬀ ects mostly the valence band by 0.32 eV, which results in the surface band gap closing by 0.31 eV; at the (112 ̅ 0) surface, the valence band remains ﬂ at and the conduction band bends up by 0.18 eV opening the surface band gap by 0.12 eV.


■ INTRODUCTION
−4 The chemical and physical properties related to the surface structures of ZnO are of fundamental interest and are also key to the material's applications.Accurate characterization of surface structure and properties is therefore essential.In the bulk, the wurtzite structure (zincite, Figure 1) is the most stable polymorph of ZnO over a wide range of temperature and pressure. 3It has four principal low-index surfaces: two side faces that are nonpolar, (101̅ 0) and (112̅ 0); and two opposite polar, (0001)-Zn and (0001̅ )-O.The two nonpolar surfaces are composed of equal numbers of cations and anions in each layer, whereas the polar surfaces have monolayers of cations and anions alternating along the c-axis.−8 Despite extensive computational and experimental study, there remain many uncertainties in both the structure and electronic properties of these surfaces.We therefore report a detailed theoretical study of the two main nonpolar wurtzite ZnO surfaces: (101̅ 0) and (112̅ 0).We analyze the atomic structure of the clean surfaces, the stability of both morphological features (steps and grooves) and dimer vacancies, and the effect on the ionization potential and surface band bending.Our analysis leads to clear and coherent models for these two key surfaces of this widely studied material.
■ STRUCTURE AND ELECTRONIC PROPERTIES OF NON-POLAR SURFACES Until recently, experimental techniques had difficulty in producing high quality large single crystals with well-defined crystalline surfaces, as ZnO-cleaved crystals suffer from stress and external forces that may affect the crystallinity and flatness of the surfaces; as a result, experimental reports concerning zinc oxide surface structures showed significant variations as discussed below.Preparation of clean surfaces is an essential requirement for a number of characterization techniques and many applications.Films grown epitaxally have proved to be highly crystalline 6,7,9−14 and are now widely used in surface studies.In the study of the nonpolar surfaces of ZnO, there is a further complication in that ZnO crystals usually grow along the polar hexagonal direction with opposite polar surfaces exposed at the top and bottom of the film.−21 Choosing an appropriate substrate (such as Al 2 O 3 ) 6,10,12 −14,22 that helps the growth of ZnO crystals along the nonpolar directions is necessary to expose surfaces of interest, which would still be strained due to a large lattice mismatch.With a large mismatch, a very large strain energy may build up in the epilayer, thus creating a series of different defects.5][6][7]9,10,[12][13][14]22,23 Experimental characterization reports on nonpolar surfaces still show the surfaces to have a high density of defect sites (vacancies), steps, and more complex morphological features. 24,25 An atomistic study by Whitmore et al. 26 showed that the creation of vacancies, steps, and nonflat ZnO surfaces is in fact energetically inexpensive. Therefor it is expected that under stress, ZnO surfaces show a certain degree of roughness as a result of surface preparation, including cleaving and polishing.
In gaining an understanding of the surface structural properties, there are two major issues to address.The first is the character of the atomic relaxation at clean surfaces compared to the bulk and the second concerns the electronic structure of such surfaces.Structurally (Figure 2), the (101̅ 0) ZnO surface has been investigated extensively using both theoretical and experimental techniques; however, the nature of surface relaxations remains controversial.Computational studies have predicted an uppermost zinc relaxation toward the bulk by 0.15−0.57Å, 26−35 and experimental reports range from 0.06 to 0.45 Å. 25,36−38 Different structures have been suggested including (a) almost bulk-like terrace structures with no pronounced atomic relaxation, 25,39 (b) terraces showing strong inward O relaxation with the cation lying at the surface, 25,40 (c) terraces with strong Zn relaxations toward the bulk, which also results in a shortening of the corresponding topmost Zn−O bonds, [26][27][28][29][31][32][33][34][35][36][37][38]41 and (d) highly defective surfaces with steps and vacancies.25,26,42 The latter has been supported by experimental work, 25,42 which has shown the presence of steps together with only partial occupation (due to the presence of vacancies) of the first two surface atomic planes (layers), with "occupancies" of 0.77 ± 0.02 and 0.90 ± 0.04 (the value from 0 to 1 giving the probability of site occupation), in the first and second layers, respectively. Whitmore et al. 26 concluded that the energy cost for creating these steps and broadening the terraces was low.Thus, both experiment and theory suggest that steps and terraces are common features of the surface.
The (112̅ 0) nonpolar surface has been so far characterized in much less detail, with pertinent structural results regarding this surface being as follows.Low energy electron diffraction (LEED) measurements 43 on annealed (112̅ 0) surfaces indicate that its atomic structure is bulk-terminated within the accuracy of the measurements. 15According to an early ab initio study 44 the (112̅ 0) surface remains close to a bulk-terminated position, in agreement with both the LEED measurements 43 and the tight-binding (TB) model calculations by Ivanov and Pollmann. 45In the former computational study, however, only three degrees of freedom per surface layer were relaxed, which, as noted by Meyer and Marx, 32 is only a first approximation.The density functional theory (DFT) study by Meyer and Marx 32 has found instead that the atomic relaxations on the two nonpolar surfaces are similar.Finally, as with the (101̅ 0) surface, Dulub et al. 24 report a high density of small terraces running along the ⟨0001⟩ direction and long grooves (ca.250 Å wide and 50 Å deep) along the ⟨11̅ 00⟩ directions as seen from scanning tunneling microscope (STM) images.From LEED and low-energy ion scattering (LEIS) analysis, it was concluded that the (112̅ 0) is the roughest of all the four main low-index ZnO surfaces.Hence, the surface still remains significantly undercharacterized.
As mentioned earlier, the second major issue concerns the electronic structure of the surface.With the ZnO atomic structure well defined, its electronic properties can be predicted using ab initio methods.A correct positioning of the band edges of ZnO is essential to calculate a great range of physicochemical properties such as electron affinity (A), ionization potential (I), work function (Φ), etc., which are crucial in the design of electronic devices.−49 Band gap, Fermi level, and correct positioning of the valence band maximum (VBM) and conduction band minimum (CBM) are fundamental electronic properties which control TCO behavior.Until now, different theoretical approaches have been used to calculate the bulk ionization potential and band alignment; however, each has implicit difficulties. 50Sokol et al. 51 calculated a bulk ionization potential The Journal of Physical Chemistry C Article (I b ) value of 7.71 eV using a hybrid QM/MM approach, which is in close agreement with the experimental value reported by Swank 52 (7.82 eV).A recent method developed by Logsdail et al. 50to calculate I b at the plane wave DFT levelshowed agreement with experimental data for a range of different rocksalt ionic oxides using the PBEsol0 functional.This method includes the surface polarization effects using simple polarizable shell-based interatomic potentials.At the surface, near to the vacuum, there is a shift in the electronic energies, or the band structure (surface band bending) due to the change of the Madelung potential, bandwidth, and surface polarization effects.Therefore, the ionization potential will differ between the bulk and the surface.In a recent study, Hinuma et al. 53 have calculated the surface ionization potential (I s ) and A for the (101̅ 0) and (112̅ 0) ZnO surfaces using GWΓ 1 @HSE (GW approximation with vertex corrections in the screened Coulomb interaction using Heyd−Scuseria−Ernzerhof hybrid functional): viz.an I s value of 8.15 eV (exp.8.00 eV) 54 and 8.17 eV (exp.7.82 eV) 52 and an A value of 4.28 eV (exp.4.60 eV) 54 and 4.30 eV (exp.4.38 eV) 52 for (101̅ 0) and (112̅ 0), respectively.
We now continue with a description of the computational methods used for this work.In presenting our results, we first discuss the structure and stability of the clean (101̅ 0) and (112̅ 0) surfaces, followed by the stability and structure of steps and dimer vacancies on the (101̅ 0) surface.We then introduce our calculations on the bulk ZnO ionization potential and the effect on it of the relaxed structure and the presence of dimer vacancies, grooves, and steps, and finally, we analyze the surface band bending on these two clean nonpolar surfaces.

■ THEORETICAL METHODS
The periodic DFT code VASP (Vienna Ab initio Simulation Package) 55,56 was employed in this study.VASP uses a plane wave basis set to describe the valence electronic states.
Exchange and correlation energy was treated with two different generalized gradient approximations (GGA): Perdew−Burke− Ernzerhof (PBE) 57 and the PBE functional revised specifically for solids: PBEsol. 58PBE0 and PBEsol0 hybrid exchange− correlation (xc) functionals were used with a 25% of the exact exchange from the Hartree−Fock (HF) theory.Interactions between the cores (Zn:[Ar] and O:[He]) and the valence electrons were described using the projector-augmented wave (PAW) 59,60 method.
In our study of defective surfaces we have also used a complementary atomistic approach.The interatomic potential (IP) code GULP (General Utility Lattice Program) 61,62 was employed to study surface morphological features, including steps, dimer vacancies, and grooves, at the (101̅ 0) ZnO surface, with the Born, polarizable shell model potentials developed for ZnO by Whitmore, Sokol, and Catlow, which show excellent agreement with a range of experimental data (see Table 1 in ref 26).Moreover, DFT theory was used to verify the structure and stability of the defective ZnO surfaces.The detailed description of the defective ZnO surface models derived from both the DFT and IP techniques is given in the Results section.
Bulk ZnO.For the structural optimizations, we checked convergence of the total energy with respect to k-mesh sampling and plane wave energy cutoff; the total energy was converged to 1 meV.For GGA functionals, good convergence was achieved with a cutoff of 700 eV and a k-mesh of 11 × 11 × 9 was used for both bulk relaxations and calculation of the density states (DOS).The iterative relaxation of the ions was not stopped until the forces on the ions were all less than 0.01 eV Å −1 .For hybrid functionals, the total energy criteria were kept as in GGA functionals.A computationally less demanding cutoff of 500 eV and a k-mesh of 9 × 9 × 7 was found to be sufficient to converge for bulk relaxations and DOS.The structures were deemed to converge when the force on every ion was less than 0.01 eV Å −1 .In general, all the calculations show good agreement with the experimental lattice parameters (Table 1), although for both a and c, a better agreement was found with the PBEsol and PBEsol0 functionals.The hybrid functionals also improve  The formation energy for zinc oxide and the cohesive energies with respect to Zn metal and O 2 are presented in Table 2. Analyzing the ZnO formation energy, the hybrid functionals give better results than GGA, for the PBEsol0 functional; the difference with respect to the experimental value is ca.11%.Comparing the cohesive energy of Zn metal , the PBEsol0 functional also gives better results with ca. 5% of deviation.The experimental bond energy of O 2 is −5.12 eV, whereas that for the PBE0 functional is −5.188 eV (∼1.3% of difference) and −5.584 eV for PBEsol0 (∼9%).The GGA functionals poorly reproduce the bond energy of oxygen.
Clean (101̅ 0) and (112̅ 0) ZnO Surface Models.For the clean ZnO surfaces, cell parameters and atoms in the middle layer were kept fixed, whereas all other ions were allowed to relax.The surface energy (E surf ) was converged to 1 mJ/m 2 with respect to the thickness of the slab and k-mesh sampling.Convergence was fully achieved for a cell with 15 double layers (60 atoms) and 15 layers (60 atoms) for the (101̅ 0) and (112̅ 0) surfaces, respectively.The slabs were separated by a vacuum gap of 15 Å.A k-point sampling of 7 × 7 × 1 was found to be sufficient.The cutoff energy was kept as in the bulk.

■ RESULTS AND DISCUSSION
Clean (101̅ 0) ZnO Surface.The atomic displacements of the first two layers of (101̅ 0) surface atoms are shown in Table 3 and illustrated in Figure 2. Atoms in deeper layers remain almost in bulk positions.In general, the same pattern was observed whether the GGA or hybrid functional approaches were used.From the first layer, Zn ions show strong relaxation inward and a displacement parallel to the surface (y direction); O atoms remain almost in bulk positions,with just small relaxations away from the surface.In the second layer, both ions relax toward the surface: the O relaxation is very small, whereas that of Zn is more substantial.
Our results are in agreement with the LEED analysis by Duke et al., 37 which predicted Zn inward relaxation of −0.45 ± 0.1 Å and likewise a movement of the top-layer oxygen by 0.1 ± 0.2 Å toward the surface; moreover, the same structure was seen in the high-resolution transmission electron microscopy (HRTEM) images by Ding and Wang. 38An angle-resolved photoemission spectroscopy (ARPES) study 36 also showed this strong Zn relaxation by −0.40 Å.In addition, there is good agreement with previous theoretical studies: 26,28,31,32,35 and there is a general consensus that uppermost zinc atoms relax by (or greater than) −0.21 Å toward the bulk and around 0.16 Å in the y direction, toward the O; second-layer zinc atoms relax by (or greater than) 0.132 Å away from the bulk.We also noted small relaxations of first and second layer oxygen atoms away from the surface.
Clean (112̅ 0) ZnO Surface.In general, as with the (101̅ 0) surface, we found that GGA functionals show larger relaxations in comparison with hybrid GGA (Table 4 and illustrated in Figure 3).We observed larger Zn relaxations in the (101̅ 0) surface along the z axis; nevertheless, lateral displacements are larger in the two directions (x and y) for (112̅ 0).The distortion observed for surface ions in this study can contribute to the surface roughness seen in LEED and LEIS analyses. 24ur calculations predict strong top-layer zinc relaxations along all three crystallographic directions; oxygen ions remain, as in (101̅ 0), almost in their bulk positions, resulting in an anion terminated surface.Ions in deeper layers remain almost in bulk positions.
Our results are in good agreement with previous theoretical work produced using interatomic potential methods by Nyberg et al., 40 and with the ab initio studies of Meyer and Marx 32 and Marana et al. 31 As with the (101̅ 0) surface, we calculated strong top-layer Zn movement to the bulk (from −0.190 to −0.229 Å).Zn also shows substantial displacement in both directions parallel to the surface.In the y direction, towards the oxygen atom, relaxations are from −0.199 to −0.233 Å; and in x from −0.095 to −0.124 Å.Table 4 also shows Zn−O distances to permit comparison with the extensive study on ZnO surfaces made by Meyer and Marx. 32Figure 3 shows the atomic structure of the (112̅ 0) surface produced from our theoretical study using GGA and hybrid functionals.
Stability of the Clean (101̅ 0) and (112̅ 0) ZnO Surfaces.We calculated the surface energy for the nonpolar (101̅ 0) and (112̅ 0) surfaces, which is in good agreement with previous theoretical work (Table 5).We found the (101̅ 0) surface to be more stable for all the DFT functionals used.However, the difference in the surface energy between the two surfaces is small (0.04−0.08 J/m 2 ), so we infer that under thermodynamic equilibrium these two surfaces would coexist in almost equal proportions.
It has been seen that the PBE functional tends to underestimate the surface energies; 32 PBEsol partially corrects the PBE deficiency and hybrid GGA funtionals show more accurate results.LDA and hybrid B3LYP functionals tend to overestimate surface energies.The slight difference between the two main nonpolar ZnO surfaces might be due to the more distorted (112̅ 0) surface, as was shown in STM images by Dulub et al. 24 We also calculated the same surface energy for both surfaces using IPs and obtained similar results to our DFT calculations.
Wander and Harrison 44 calculated a surface energy of 2.05 J/ m 2 for the (112̅ 0) surface, which is in disagreement with the rest of the work published.They used a small slab of seven layers and not all parameters were allowed to relax, which could account for the high calculated surface energy.
Marana et al. 31 report a smaller surface energy difference between (101̅ 0) and (112̅ 0) surfaces than that given in ref 32.However, Marana et al. 31 compared the surface energy difference of 0.1 J/m 2 with the cleavage energy difference of 0.2 J/m 2 given in ref 32.Because the cleavage energy represents double the surface energy, the cleavage energy difference of 0.2 J/m 2 mentioned in ref 32 is actually a surface energy difference of 0.1 J/m 2 , which is the same value as calculated by Marana. 31e surface structures discussed represent definitive models for the clean (101̅ 0) and (112̅ 0) surfaces, which are consistent with experiment.
Steps and Vacancies at the (101̅ 0) ZnO Surface.Following the experimental results of Jedrecy 25 and of Parker 42 and the interatomic-potential calculations by Whitmore, 26 we use IP methods to examine models consistent with a fractional surface site occupancy of 0.75 in the first layer, which has been suggested from experiment.As already seen in ref 26, Zn and O vacancies scattered randomly over the surface are less energetically favorable than nearest-neighbor Zn−O dimer vacancies.Therefore, we concentrate on the problem of the location of dimer vacancies; henceforth, the term "vacancy" will denote a dimer vacancy.
We built a one-sided 2D periodic surface model using a tworegion approach, which has been widely employed in modeling surface structures with interatomic potential based methods. 61,62This approach allows relaxation of the ions in the region next to the vacuum, whereas the second substrate region is held fixed representing the bulk crystal.The (101̅ 0) ZnO surface calculations converged using five layers (20 atoms, ≈13 Å thick) in both region one and region two.To simulate 75% occupation of the topmost surface layer, we constructed a 4 × 4 supercell and removed 4 of the 16 Zn−O dimers in the first layer.We used an in-house developed python code to build all possible reconstructions; all different configurations were fully relaxed using the GULP code employing interatomic potentials.We analyzed the first five lowest energy structures.Key features of the results are as follows as (Figure 4): •  1.00 1.00 LDA 32 1.15 1.25 PBE 32 0.80 0.85 LDA 28 1.19 1.23 PW91 35 1.04 1.06 B3LYP 31 1.30 1.40 IP 40 1.10 1.20 IP 26 1.00 IP 28b 1.20 a Using potentials reported by Whitmore. 26b Using Binks 66 potentials.(Zn) and negative (O) ions, which may explain the high instability of structure z.Next, we calculated the vacancy formation energy as a function of vacancy concentration.Comparison of the energies among the different surface structures is possible; however, it only tells us which structure is more stable, but does not provide any information about the stability of the vacancy.The vacancy energy formation per ZnO dimer was calculated by where E surface def is the energy of the defective surface, n is the number of dimer vacancies ( 4), E dimer is the energy of a ZnO unit in the bulk, and E surface nondef is the energy of the nondefective surface.The energy cost per vacancy values in Figure 4 are in a very good agreement with previous calculations. 26he low energy of formation obtained in the IP calculations for the linear defect is suggestive of the ease of step formation, a suggestion that is supported by experiment, 25,42 where it has been shown that it is almost impossible to create a (101̅ 0) ZnO surface free from steps.We decided to make use of IP and ab initio methods to calculate the energy needed to create a step along [010] and [001] directions.
We find that vacancies aligning along the [001] direction show high stability.We use GGA/PBEsol calculations to investigate the energetic cost of forming a step along the [010] direction (vacancies aligning along the [001] direction) at the DFT level, as such steps were observed in experiment using grazing incidence X-ray diffraction (GIXD) techniques, 25  ).With these models, a vacancy will represent a step along the [010] and [001] direction, respectively.Half of the first and second layer Zn−O dimers were removed from all the structures.Dimers were removed in such a way that dimer vacancies stay together, increasing the size of the step as we increased the size of the supercell (Figure 5).For the larger systems, using 16 unit cells results a lateral separation of ca. 26 Å and ca.44 Å between periodic images of steps along the [010] and [001] direction, respectively.Extrapolating the curve in Figure 6 to 0 (n = ∞), we find that, in the limit of infinite separation, E step[010] = 0.029 eV/Å (cf.0.027 eV/Å using IP).The calculated E step[010] suggests strongly that steps would be seen even at room temperature (kT at room temperature is equivalent to 0.025 eV).
The energy cost of a step along the [001] direction shown in Figure 6 does not seem to converge and could not be extrapolated to the limit of infinite step separation.However, we find that, where the size of the slab is (1 × 16), the step energy is about 9 times greater than that of creating a step along the [010] direction.The high energy cost per step along the [001] direction is as a result of revealing charged atomic rows on each side of the step, creating a strong electric field, which is highly unstable without any major atomic reconstruction, as seen in polar (0001) and (0001̅ ) ZnO surfaces.Steps along the [001] direction show strong displacements with respect to bulk positions and will not be discussed further due to their high instability.
Atomic Relaxation Close to a Step Along [010] Direction.We now summarize the most notable structural features of the (101̅ 0) ZnO surface with a step in the [010] direction (Figure 5).In general, we observed larger relaxations for Zn than for O atoms.For the first layer, a strong average Zn relaxation of 0.312 Å toward the bulk and 0.155 Å along the [001̅ ] direction is calculated (similar to the relaxation on a clean (101̅ 0) surface, summarized in Table 3).Moreover, we find a relaxation of 0.116 Å outward from the step only for the Zn ions that are on the edge of the step (denoted as Zn [1*] , where the number refers to the layer and the " * " indicates step edge positions), whereas smaller relaxations are seen in this direction for the rest of Zn atoms.The maximum displacement of 0.063 Å (with an average of 0.036 Å) is observed for O [1*] away from the bulk; the relaxations are less pronounced in the other two directions.For the remaining O atoms there are no significant relaxations.
In the second layer, Zn [2*] atoms showed the strongest relaxations of 0.292, 0.115, and 0.105 Å toward the step, in the [001̅ ] direction and inward, respectively.The rest of the Zn [2]  ions behave in a completely different way with an average strong relaxation of 0.143 Å outward and smaller relaxations along the other two directions.Again, the O [2] ions show smaller relaxations.
The atoms in the third layer are divided into two sets: atoms with reduced coordination number (e.g., Zn [3−] , where the minus sign indicates the reduction in the number of bonds) and fully coordinated atoms (Zn [3] ).As expected, Zn [3−] atoms show the larger relaxations, whereas the O [3] atoms remain close to their bulk positions.The Zn [3-] ions showed an average relaxation of 0.170 and 0.320 Å in the [001̅ ] direction and toward the bulk, respectively (similar to the relaxations in Zn [1] ).Smaller relaxations were observed for Zn [3] ions.
The ions in the fourth layer are also divided into two sets: ions below under-coordinated atoms (e.g., Zn [4 §] ) and ions below fully coordinated atoms (e.g., Zn [4] ).In general, only Zn [4 §] atoms relax strongly outward (0.140 Å) in a manner similar to that for the Zn [2] ions, with the only difference being that the Zn [4* §] showed a smaller relaxation (0.045 Å) in the same direction.
Ionization Potential and Band Alignment.The bulk ionization potential (I b ) is the energy required to remove an electron from the system.Ideally, we would like to determine this as a "bulk" property, which is independent of the surface termination.When I b is known, the effect of a particular surface, i.e., the surface band bending, or offset can be determined, yielding the surface ionization potential, I s , which is in turn related to a wide range of surface physical/chemical properties.We calculated the ZnO bulk ionization potential (with a The Journal of Physical Chemistry C Article recently developed method) 50 and then studied the effect of different surface features (including point defects in the form of dimer vacancies and a line of four vacancies along the [001] direction, grooves and steps) on the ionization potential.Table 6.ZnO Ionization Potential, Where I is the Ionization Potential, "b" Refers to the Bulk, "s" to the Surface, and "D" to the New Method Used (Taking into Account Surface Polarization Effects) 50    .44eV for the bulk band gap; 80,81 for the (101̅ 0) and (112̅ 0) a I s value of 8.00 eV 54,67 and 7.82 eV, 52 respectively."(D)" 50 has the same meaning as in Table 6.The positioning of the CBM bands for the relaxed (101̅ 0) and (112̅ 0) surfaces made by adding the band gap calculated for each relaxed surface presented in Table 6 to the VBM value.The Journal of Physical Chemistry C

Article
The particular surface morphology will affect the positioning of the bands.As was seen by Jacobi et al. 54 and Uhlrich et al., 67 ideal surfaces result in higher band bending values.Due to the surface treatment given by Swank et al. 52 in (112̅ 0) and the extrapolation to time zero after ion bombardment and annealing by Jacobi et al. 54 in (101̅ 0), the values reported there closely correspond to the bulk values.
Bulk and Surface Ionization Potential.Table 6 shows the calculated bulk (I b ) and surface (I s ) ionization potentials.I b was obtained by: (i) calculating the difference in energy for the O 1s electrons in the bulk and in the middle of an unrelaxed slab; (ii) shifting all electronic levels in the bulk by this difference to align them with the vacuum level obtained in the surface calculation; (iii) calculating the distance between the shifted bulk VBM and the vacuum level.For I s , we followed the same procedure but using the energy of the O 1s electrons in the middle of a relaxed surface instead of those in an unrelaxed slab.The method proposed by Logsdail et al. 50improves I b,D (where "D" makes reference to the multipolar shift) by ca. 1 eV compared to the widely used "band alignment" technique. 68As noted, this method takes into account surface polarization effects using polarizable-shell based IPs.There is not a significant difference between PBE and PBEsol functionals or between PBE0 and PBEsol0 (Figure 7).Moreover, the calculated I b,D values using hybrid functionals show a good agreement with experiment (7.82 eV) 52 (Figure 7), and calculations by Sokol 51 using a QM/MM approach (7.71 eV).Furthermore, this new method 50 is shown to be practically surface independent.
Table 6 shows that the surface effect on I s is stronger when hybrid functionals, rather than GGA functionals, are used.Uhlrich et al. 67 observed an I s value of 8.1 ± 0.1 eV for dry annealed (101̅ 0) crystals, whereas Klein et al. 46 reported a value of 7.7 eV.Hinuma et al. 53 (using GWΓ 1 @HSE) calculated values of 8.15 and 8.17 eV for (101̅ 0) and (112̅ 0) surfaces, and Stevanovićet al. 69 (with DFT and GW calculations) calculated 7.53 and 7.60 eV.However, the calculated numbers by Stevanovićet al. might have converged to the incorrect value as suggested by the Klimešet al. 70 in their study of the energy convergence in the GW approximation using the PAW method.
The different values of I s reported in the theoretical work depends on the method used.Hinuma et al. 53 aligned the CBM by adding the calculated VBM to the experimental bulk band gaps; however, VBM and CBM bend differently, as noted for ZnO in Figure 9 and in CdO. 71The creation of steps and grooves at the (101̅ 0) ZnO surface has a very small effect on the ionization potential (a decrease of ca.0.04 eV), which might be an explanation of the stability of such features as suggested by this and previous work. 25,26However, creating a 25% dimer vacancy at the same surface has a bigger impact on I s (a decrease of ca.0.13 eV).
Next, we have calculated the band gap for the surfaces using the hybrid functionals (shown in Figures 7 and 8).For the (101̅ 0) surface, the band gap is smaller (ca.2.9 eV) than in bulk (ca.3.13 eV); whereas, for the (112̅ 0) surface, the band gap is slightly larger (ca.3.2 eV).The latter finding may imply that using 15 layers thickness is not enough to represent the band gap correctly.We attribute this overestimation to a quantum confinement effect in a relatively thin slab surface model, which is particularly strong for delocalized conduction states.This rationale was confirmed using a computationally less expensive GGA/PBEsol functional: for (112̅ 0) slabs with 29 layers the band gap decreases from 0.89 to 0.77 eV to be compared with the ideal bulk value of 0.70 eV.The rate of the band gap decrease in these calculations, however, is significantly smaller than that in our hybrid functional calculations and, we note, the hybrid functional Gaussian function based study reported in ref 31, has not found any difference between their (112̅ 0) surface and bulk band gap values.
The smaller band gap at the (101̅ 0) surface is attributed to the upward bending of the VBM, as seen in Figures 8 and 9.With the atomic structure and surface energy converged, the slightly larger band gap for the (112̅ 0) surface is caused by the quantum confinement effects: whereas the (101̅ 0) surface is composed of 15 double layers with a pair of ZnO in each single layer and a slab thickness ca.40 Å, the (112̅ 0) surface has 15 layers with two pairs of ZnO on each layer and only half of the (101̅ 0) thickness.Therefore, for the (112̅ 0) surface there are the same number of ions as for (101̅ 0), but the slab has only half the thickness, resulting in a stronger quantum confinement effect on the (112̅ 0) surface.
Work Function.Previous work has given strong emphasis to the determination of the work function (Φ) of ZnO.This property is defined as the energy needed to take an electron from the Fermi level, which like I s , is very sensitive to the sample history, surface preparation procedure, the method of measurement and the facet involved. 72,73For example, for the (101̅ 0) face, different work function values have been reported by experiment and theory: ca.4.6 eV, 52,54,67,74−76 78 Typically, in intrinsic ZnO, the Fermi level lies just below the CBM.Therefore, the CBM values reported in Figure 7 are comparable to the work function of ZnO.
Density of States (DOS).In general, as suggested in Figure 8, cleaving a nonpolar ZnO surface along the (101̅ 0) and (112̅ 0) planes results in a local rise of the VBM and CBM positions, respectively.Figure 9 represents the shift of the band edges at the nonpolar surfaces across the slab.For these calculations, we separated the total DOS into projected DOS per layer.The VBM of the middle layer was set to the corresponding negative I s value (as this layer represents the bulk).The total DOS for the bulk are shown for comparison and only the hybrid calculations are shown because GGA functionals tend to underestimate severely the band gap.No significant differences were observed when using PBE0 and PBEsol0 functionals.For the (101̅ 0) termination we see a split in the O 2s band at the surface due to the lower coordination of the surface atoms.There is also a decrease in the intensity of the Zn d band.For the valence band of the "surface" layer, there is a shift to the VBM of the highest peak, indicating a relative destabilization of the majority of the valence electrons near the surface, which places their energies close to the VBM.Another important feature is the peak that is seen in the CBM (indicated with a dash line in Figure 9); at the surface, this peak overlaps completely with the adjacent peak.
For the (112̅ 0) surface, different behavior is observed.There is no splitting in the O 2s band at the surface.However, there is a displacement of the top of this band.The Zn d band is pushed away from the VBM and the O 2p band is more concentrated close to the VBMas in (101̅ 0).The band gap shrinks at the surface for both (101̅ 0) and (112̅ 0) terminations due to displacement of the top of the VBM; the position of the bottom of the CBM is preserved.The first peak in the CBM overlaps completely with the adjacent peak as in the (101̅ 0) surface.

■ SUMMARY AND CONCLUSIONS
Our hybrid DFT calculations for both clean nonpolar surfaces of ZnO confirm earlier GGA reports of the strong inward cationic relaxation in the topmost atomic layers of the material, although with smaller amplitudes, accompanied by pronounced lateral displacements of cations and anions especially on the (112̅ 0) surface.
The calculated surface energies indicate a higher stability of the (101̅ 0) surface over the (112̅ 0) surface.However, the difference in the surface energy between the two surfaces is small, which implies the importance of thermal vibrational contributions to the free energy, which could determine the crystal morphology under thermodynamic equilibrium, or of the kinetic crystal growth factors.However, the (101̅ 0) surface is seen predominantly in experiment.
Our calculations provide the first computational rationale for the extensive stepping observed on the nonpolar surfaces of ZnO at the ab initio level of theory.The energy cost of creating a step along the [010] direction, E step[010] , on the (101̅ 0) was found to be 0.029 eV/Å (cf.0.027 eV/Å using interatomic potentials).Thus, we expect that steps would be seen even at room temperature (k B T ≈ 0.025 eV).

Article
The surface structures obtained in this work have been employed to determine the bulk and surface ionization potentials of ZnO along with the electronic band bending.The calculated bulk ionization potential values using the method described in ref 50 show an improvement of about 1 eV with respect to the widely used "slab alignment" method.For hybrid functionals, I b is calculated as ≈7.6 eV, compared to the experimental value of 7.82 eV.Surface features such as steps and grooves are shown not to have a strong effect on the ionization potential (a decrease of ca.0.04 eV), whereas a 25% dimer vacancy formation at the surface would decrease the ionization potential by 0.13 eV.
The surface electronic properties are shown to converge much slower with the slab model thickness.Using hybrid functionals, the band gap at the (101̅ 0) surface is still smaller (ca.2.9 eV) than in the bulk (ca.3.13 eV).For the (112̅ 0) surface, a slab with 15 layers was not enough to represent the band gap correctly, which could be attributed to quantum confinement.It is expected that with a thicker slab, the band gap of the surface will become closer to the bulk band gap, as confirmed by our GGA calculations.
To characterize the band bending, we have deconvoluted the density of states into atomic layer contributions.The two nonpolar surfaces are seen to behave markedly differently with a local rise of the VBM for (101̅ 0), which remains nearly flat for (112̅ 0).In contrast, the CBM rises for the (112̅ 0) surface.Therefore, band gap closing will be seen on the first surface by 0.31 eV and band gap opening on the second by 0.12 eV.No significant differences were seen between PBE0 and PBEsol0 functionals.
An extension of this work to polar surfaces of ZnO will be reported in the near future.

Figure 2 .
Figure 2. Schematic representation of the relaxed and unrelaxed (101̅ 0) ZnO surface.Black lines show the bulk position structure.The stick representation is the relaxed (101̅ 0) structure.Only atoms in the first double layer were represented on the top view; darker colors, in (b), were used to represent ions in the first layer.Layers 1, 2, 3, and 4 are represented as L1, L2, L3, and L4, respectively.

Figure 3 .
Figure 3. Schematic representation of the relaxed and unrelaxed (112̅ 0) ZnO surface.Black lines show the bulk position structure.The ball-and-stick representation is the relaxed (112̅ 0) structure.Darker colors were used to represent ions in the first topmost layer.Δ values represent the distance between the selected ions along the specified direction.
The lowest energy structure is a line of four vacancies along the [001] direction, as was seen in ref 26 (structure 1).• The maximum number of vacancies in the same row along [001] is energetically preferable.• Connected vacancies have lower energies (see 2−4 structures) than isolated (see structure 5).• Zigzag patterns (e.g., structure 4) and diagonal connections among the vacancies are observed in the lowest energy configurations.• Two or more vacancies in the same row along the [010] direction are not present in the first ten lowest energy structures.

Figure 4 .Article•
Figure 4. Top view of the lowest five (1−5) and highest (z) energy (101̅ 0) ZnO structures with 75% first layer occupancy.Energy cost per vacancy is represented below each structure.Red, gray, and blue balls represent oxygen, zinc, and vacancies, respectively.Darker balls are first layer atoms.The black box represents the 4 × 4 supercell.

Figure 5 .
Figure 5. Schematic representation of the relaxed and unrelaxed (101̅ 0) ZnO surface showing a step along the (a, b) [010] direction and along the (c, d) [010] direction.(a) Side view along the [001] direction of the (101̅ 0) ZnO surface, (b) top view along the hexagonal axis of the (101̅ 0) ZnO surface, (c) side view along the [010] direction of the (101̅ 0) ZnO surface, and (d) top view along the hexagonal axis of the (101̅ 0) ZnO surface.Black lines show the bulk position structure.Ball and stick representation is the relaxed structure.Only atoms that are in the green region were shown in the top view; darker colors, in (b) and (d), were used to represent uppermost ions.

Figure 6 .
Figure 6.Energy of step formation on the (101̅ 0) ZnO surface, where n is the size of the supercell along the [010] or [001] direction.The dashed line is a linear fit for the step along the [010] direction.

a
Experimental bulk ionization potential, I b = 7.82 eV.52 b Ionization potential for the surface features was calculated only at GGA level (PBEsol) as the size of the supercells makes hybrids very computationally expensive.

Figure 7 .
Figure 7. ZnO band alignment based on the ionization potential.The horizontal lines represent the experimental reported values: 7.82 eV for I b 52and 3.44 eV for the bulk band gap;80,81 for the (101̅ 0) and (112̅ 0) a I s value of 8.00 eV54,67 and 7.82 eV,52 respectively."(D)" 50 has the same meaning as in Table6.The positioning of the CBM bands for the relaxed (101̅ 0) and (112̅ 0) surfaces made by adding the band gap calculated for each relaxed surface presented in Table6to the VBM value.

Figure 8 .
Figure 8. ZnO band bending of the CBM and VBM.Bulk values for CBM and VBM were taken from experiment.E vac , −eV CBM , −eV VBM , E F , χ, and Φ represent the vacuum, band bending at the CBM, band bending at the VBM, the Fermi level, the electron affinity, and the work function, respectively.Each diamond/square represents a layer.Values were taken from Figure 9.

Figure 9 .
Figure 9. (a) Representation of the total density of states for wurtzite ZnO structure using the PBEsol0 functional, for comparison see Figure 2 in ref 82.Surface band bending across the ZnO surfaces: (b, c) (101̅ 0) and (d, e) (112̅ 0).Surfaces structures from (b, d) PBE0 and (c, e) PBEsol0 functionals were taken to produce the surface band bending.Layers 1 and 8 are labeled as "Surface" and "Bulk" layers, respectively.The vertical red line represents the top of the valence band.The valence band (VB) and conduction band (CB) were amplified by a factor of 3 and 25 to make the changes more visible.

Table 1 .
Bulk Properties of the Wurtzite Structure of Zinc Oxide Experimental values were taken from the lowest available temperature, neutron single crystal diffractometry at 20 K. 64 b Width of the oxgyen 2p band.c Energy difference between Zn 3d and O 2p bands. a

Table 2 .
Formation and Cohesive Energies a a The values reported in this table correspond to a cutoff energy of 700 eV.Experimental enthalpies of formation were used as comparison.

Table 3 .
Atomic Relaxations of the First Two Layers of the Nonpolar (101 ̅ 0) Surface a 33The nomenclature used in this table is shown in Figure2.Subindexes represent the layer and the direction of the relaxation.All relaxations and distances are given in Å. b Zn−O angle.cTakenfrom Wander et al.33d Calculated from distances.e Grazing incident X-ray diffraction.

Table 4 .
Atomic Relaxations of the First Two Layers of the Nonpolar (112 ̅ 0) Surface a of parameter u (within 0.26%).Higher cutoff energies (900 eV) do not improve the results.Band gap energies for both PBE0 (3.142 eV) and PBEsol0 (3.128 eV) underestimate experimental data (3.44 eV) by ca.10%.Hybrid fuctionals also show better agreement with the experimental width of the oxygen 2p band and the energy difference between Zn 3d and O 2p bands.
a The nomenclature used in this table is shown in Figure3."Bulk" rows represent values for the unrelaxed surface.All relaxations and distances are given in Å. reproduction

Table 5 .
Surface Energy, E surf (J/m 2 ), of the Nonpolar ZnO Surfaces with the Different Functionals and in Comparison with Previous Calculations (Energies in eV) ca. 4.3 eV.19,69,77−79Kuo et al. reported a work function value of 3.74, 3.95, and 4.21 eV for as-deposited ZnO films, after Ar sputter cleaning and after exposure to oxygen plasma.