Why Are Polar Surfaces of ZnO Stable?

: We probe and rationalize the complex surface chemistry of wurtzite ZnO by employing interatomic potential calculations coupled with a Monte Carlo procedure that sampled over 0.5 million local minima. We analyze the structure and stability of the (0001) and (0001 ̅ ) ZnO surfaces, rationalizing previous patterns found in STM images and explaining the (1 × 1) periodicity reported by LEED analysis. The full range of Zn/O surface occupancies was covered for a (5 × 5) supercell, keeping | m Zn − m O | / N ≈ 0.24 where m and N are the numbers of occupied surface sites and total surface sites, respectively. Our calculations explain why the (5 × 5) reconstructions seen in some experiments and highlight the importance of completely canceling the inherent dipole of the unreconstructed polar surfaces. The experimentally observed rich reconstruction patterns can be traced from the lowest occupancy, showing the thermodynamically most stable con ﬁ gurations of both polar surfaces. Triangular and striped reconstructions are seen, inter alia , on both polar surfaces, and hexagonal patterns also appear on the O terminated surface. Our results explain the main experimental structures observed on these complex surfaces. Moreover, grand canonical simulations of ZnO polar surfaces reveal that disorder is favored and, thus, con ﬁ gurational entropic factors is the the cause of their stability.


■ INTRODUCTION
Modeling or predicting surface structures often relies on intuition and input from experiment. Even if successful, for complex surfaces, this approach is often insufficient to determine their real structure and fails to explain the effect of experimental conditions that are vital in understanding the surface chemistry of oxides. To address this failure, we have implemented unbiased Monte Carlo methods for canonical and grand canonical ensembles and applied them to the intriguing case of the polar surfaces of wurtzite ZnO.
ZnO is an affordable, earth abundant, wide band gap transparent conducting oxide with applications in electronics, optoelectronics, sensors, pharmaceuticals, and catalysis. 1,2 These applications are greatly affected by or even crucially depend on the ubiquitous polar terminations of ZnO crystals. Exposed polar surfaces should not be able to have stable flat faces with bulk periodicity, yet nevertheless experiment often shows that they do. 3,4−12 There is still no satisfactory explanation for this stability, which is vitally important not only from the fundamental point of view but is also key for our understanding the synthesis of methanol, photocatalysis, and hydrogen sensing. 1,13,14 Our results lead to a new understanding of the complex and intriguing surface chemistry of this widely studied oxide. The proposed approach is novel in the field of surface structure prediction and, in fact, in modeling realistic extended systems across the field of the chemistry of materials, where grand canonical Monte Carlo methods have been primarily confined to Boltzmann lattice simulations and related techniques. Our work has allowed us for the first time to rationalize and reconcile a large body of seemingly contradictory experimental findings by linking the preferred structural motifs to the environmental and preparation conditions and to explain a number of observed structures including the controversial "unreconstructed", or pristine ZnO polar surface. We do not only attribute the stability of the polar surfaces to the high degree of disorder but calculate the relevant thermodynamic surface potential (grand, or Landau surface potential) and show that it can become comparable with or even lower than the corresponding surface energy of the competing nonpolar surfaces, which is a first proof of the observed surface stability. These results are key both to fundamental and applied surface chemistry of ZnO, which in turn is one of the most important advanced functional materials. A recent study on CeO 2 by Loṕez' group 15 has demonstrated the importance of configurational entropy in stabilizing polar surfaces in general.
ZnO crystals show four dominant low-index surfaces: (101̅ 0), (112̅ 0), (0001), and (0001̅ ). The latter two are categorized as polar ( Figure 1); 16 as noted they play a key role in the applications of this material 13 and have different chemical properties. 17 Their stability is the most puzzling feature of ZnO chemistry, and such stability for polar surfaces is not observed in other oxides. 18 We report here a computational study of the surface reconstruction mechanism that rationalizes this unique behavior. To explain the many contradicting experimental findings, 1,3−12,19−23 we examine in detail the most stable surface structures that will emerge upon crystal growth (under thermodynamic equilibrium) along the polar directions, which coincide with the c axis of the wurtzite crystal structure.
Despite numerous attempts guided by experimental and theoretical studies, to date it has proven difficult to rationalize the contradictory experimental observations on the polar surfaces of ZnO. Considering first the Zn terminated face, low-energy electron diffraction (LEED) 3−5 and helium atom scattering (HAS) 6 studies have shown (1 × 1) surfaces with no evidence of atomic reconstruction. In contrast, the best fitting of a grazing incidence X-ray diffraction (GIXD) study 11,12 was obtained with only 75% of the topmost Zn ions present. This finding was corroborated by a more recent scanning tunneling microscopy (STM) study, 30 where steps forming triangular shapes of different sizes were seen. Moreover, a high-resolution scanning force microscope (SFM) study by Torbrugge et al. 25 shows a Zn terminated morphology dominated by triangular shaped reconstructions coexisting with a (1 × 3) periodicity. The strong deviation of the observed 1/3 surface charge reduction from the ideal value of 1/4 is explained by the authors by, presumably, partial compensation of the ionic charge of Zn cations by hydrogen adsorbed on top of double Zn rows forming hydrides. On the O face, LEED, 4,7,8 lowenergy ion scattering (LEIS), 9 GIXD, 10−12 and surface X-ray diffraction 10 again showed an unreconstructed (1 × 1) structure, whereas a (1 × 3) (0001̅ ) reconstruction with O vacancies lying in the [100] direction was observed using HAS and LEED. 1,20 The O face reconstruction was supported by an X-ray photoelectron spectroscopy (XPS) quantitative analysis of the O 1s intensities, where a deficiency of 39 ± 10% of the topmost O ions was found. As for the Zn face, the nonideal charge reduction can be attributed to the surface hydroxyl formation supported by DFT calculations, 26 although other morphological features including steps have been discussed. 1 Crucially, recent high-resolution STM 21 also revealed (5 × 5) (0001̅ ) hexagonal patterns.
Theoretical work on polar ZnO surfaces has so far typically followed experiment by reproducing either unreconstructed 10,27 or nonstoichiometric surfaces 19,28−31 (compare with ref 32, discovering nonpolarity of ultrathin, few atomic layers thick, ZnO films). Prominently, Kresse and co-workers simulated formation of triangular patterns on the Zn face, 29,30,33,34 and hexagonal reconstructions on the O face 21,31 with a (5 × 5) periodicity. These calculations provided, however, only a partial solution to the problem and left many experimental findings unexplained, exigently, the reported (1 × 1) periodicity.
Stabilization Mechanisms. The intrinsic instability of the flat ZnO polar surfaces is described from two points of view: the covalent and ionic. In the covalent picture, surface atoms lose a quarter of their neighbors and are left with dangling bonds. Because ZnO is a polar semiconductor, the Zn face would show an excess of electrons in the dangling bonds localized on Zn, which would be formed by states that belong to the Zn 4s and 4p conduction bands. At the O face, there would be dangling bonds localized on O, which in turn has electrons missing (holes) from the O 2p valence band. Both dangling bonds and excess charge carriers should be evident from surface metallization or redox activity, which, however, are not observed. From this description, it is also clear that surface cleavage requires at least half the band gap energy equivalent per surface atom for breaking one bond. In the ionic picture, the instability of the ZnO polar surfaces is a consequence of the charge distribution of the ionic monolayers along the c axis. Assembling such monolayers produces a dipole moment in the repeat unit cell normal to the surface 16 (Figure 1). Crystals with unreconstructed flat polar surfaces would therefore resemble charged capacitors. Beyond a critical slab thickness, the accumulated voltage will result in a dielectric breakdown, causing the material to collapse and conduct, although polar surfaces of very small particles or ultrathin films could be stable. 2 In the ionic picture, which we will use in the following discussion, the surface stabilization is achieved by charge transfer: (i) electronic that results in surface metallization; (ii) ionic that yields surface nonstoichiometry; or (iii) surface coverage with charged guest species, or adsorbates (e.g., H, O, or OH).
The "metallization" mechanism 10,27,29,30,35,36 arises naturally from the covalent picture, whereas in the ionic description it implies electron transfer from the O face to the Zn. The metallization has been suggested before as the main stabilization mechanism; 10 a later scanning tunnelling spectroscopy work, however, has clearly demonstrated the prevalent ionic reconstruction mechanism and no experimental support for the existence of surface conduction in this material. 37 Moreover, photoemission experiments did not find evidence for the proposed surface states, 38 whereas calculations yield higher surface energies for the electronically stabilized surfaces compared to the atomically (ionically) reconstructed. 16,29,39 A number of experimental and theoretical reports have attributed the stability of the nonstoichiometric polar ZnO surfaces to Zn and/or O vacancies. 19,28−31 For example, the best fit to X-ray diffraction data was obtained in ref 40 with 25% Zn 2+ vacancies at the Zn terminated surface, which according to the ionic model is approximately the charge needed to compensate the dipole. For the (0001) ZnO surface, STM images and density functional theory (DFT) calculations have shown that the vacant Zn sites complemented by O vacancies can cluster as triangular shaped indentations, whereas the remaining atoms can also form triangular islands. 29,30,33 Triangular terraces terminated by single-layer steps, ca. 2.4 Å in height, 37 were assumed to be decorated by oxygen. Such reconstructions could be produced by removing the topmost Zn and O atoms from the surface, as confirmed by Monte Carlo (MC) simulations using empirical potentials. 34 In general, larger triangles are preferred over smaller ones, 37 and the energetic competition between triangles of different size is found to be very tight, explaining the observations in experiment of macroscopic roughening of the surface. These reconstructions appear to be electrostatically driven and are stable over a wide range of oxygen and hydrogen chemical potentials. 30 In contrast, at the (0001̅ ) surface, honeycomb-like reconstructions are observed occasionally in small portions of STM images and are supported by ab initio calculations. 21,31 The difference in reconstructions between the two polar surfaces has been attributed to a higher flexibility of Zn atoms to form bonds than O atoms. 31 In addition, metal oxide surfaces are usually in contact with some water and its dissociated species. There is experimental and theoretical evidence that the stability of the ZnO polar surfaces can be achieved by adsorption of charged H, O, and OH species, mainly by the creation of stable hydroxylated surfaces. 31,41,42 However, such surfaces are strongly passivated by the presence of hydroxyl groups, which is not desirable if a high concentration of active sites is needed. Different factors will determine the actual stabilization mechanisms of observed ZnO polar surfaces, resulting in differing structures; they include synthetic conditions or residuals in the ultrahigh-vacuum chamber used in STM. For example, Ostendorf et al. 43 and Zheng et al. 44 have proposed an alternative stabilization method for the (0001) surface, via the formation of facets of high step density. Steps reported by Zheng et al. have a (101̅ 4) surface orientation and are created after annealing at ≈850°C: the use of high annealing temperatures (>700°C) causes roughening.
In summary, as seen from experiment and available computational work, the atomic structure of ZnO polar surfaces is very sensitive to mechanical, chemical, and thermal treatments or conditions. 1,21,25,45−47 A consensus seems to be reached that the ionic mechanism is responsible for stabilization of clean ZnO surfaces, but how it works in atomistic detail and what effect external conditions have on these surfaces is unclear.
In this study, we conduct a first unbiased exhaustive search over the full range of possible surface occupancies and atomic configurations for both polar surfaces of ZnO, which results in a definitive prediction of surface reconstruction patterns and occupation of surface atomic planes as a function of external ZnO chemical potential and temperature. We demonstrate that accurate exclusion of the surface model dipole is an essential prerequisite for a correct account of the surface reconstruction mechanisms and associated surface energy. To address the challenge of the diverse structural and chemical conditions, we have applied Monte Carlo methods (e.g., refs 48 and 49) developed in statistical mechanics, and widely used in materials simulations (e.g., refs 50−52). We note that, previously, surface structures have already been a target of a number of Monte Carlo studies using a canonical ensemble, 34 which, however, have been severely limited in the scope, accuracy, and comparability with realistic systems. A more general approach employing the concept of a grand potential has been used in thermodynamic analysis of (i) phase segregation for nonstoichiometric surface terminations of various materials (for example, see refs 53−58), (ii) nonstoichiometric slabs, 59 (iii) interfaces, 60,61 and (iv) surface reconstructions. 62−66 Following these thermodynamic approaches, a local or global deviation from precise 1:1 stoichiometry at surfaces of ZnO can be trivially incorporated using chemical potentials of O and Zn atoms or ions that are simply added to the total energy of preoptimized surface potential energies. Such an approach is particularly beneficial for modeling ZnO growth that in solutions would involve transient nonstoichiometric (and, possibly, high-energy) configurations, with a general addition (or removal) of Zn x O y at each elementary stepcf. refs 67−73. We note that even though kinetics of crystal growth is not considered in this work, the thermodynamic approach can be used to augment and utilize our results for this purpose in the future. The majority of world production of ZnO, however, is based on metallurgical processes, in which crystal growth occurs from the gas phase and, therefore, the primary interest is in the thermodynamically stable pure stoichiometric surfaces.
In this work, we move beyond the previous MC studies of ZnO and, via an extensive sampling of the configurational space extrapolated to the thermodynamic limit, and obtain Helmholtz surface free energies. Using the concept of the chemical potential and the grand canonical ensemble, we exploit our canonical ensemble MC results to predict macroscopic surface plane occupation and resulting morphologies.

■ METHODS AND COMPUTATIONAL DETAILS
Strategy. In view of the complexity of the chemistry of polar ZnO surfaces, here we investigate the origin of their stability using unbiased global search techniques 74 on the energy landscape defined by highly accurate interatomic potentials. 75 The applicability of our method, however, is not restricted to interatomic potentials; any potential energy surface, including for example DFT or any other flavor of ab initio or semiempirical approaches can be explored, which will only be restricted by the computational costs.
We used interatomic potential (IP) based methods with a polarizable shell model 75 to calculate the Zn/O vacancy ratio needed to minimize the magnitude of dipole moment. First, bulk ZnO was optimized and cleaved along the c axis, creating both the (0001) and the (0001̅ ) surfaces. Second, we calculated the dipole due to a single bilayer composed of three sequential layers of point charges: Zn core, O core, and then O shell. Finally, we calculated the charge compensation needed to reduce the dipole to a minimum value. This compensation is done through an ionic transfer between opposing polar sides ( Figure 1).
Our choice of a simulation cell size is based on a compromise between the optimum Zn/O vacancy ratio, which could be supported, and available computer power. In particular, the dipole moment of a ZnO crystal in the polar direction can be considered as a result of accumulation of dipoles from Zn and O bilayers, within which Zn and O atomic planes are separated by the smallest distance of kc, where c is the hexagonal lattice parameter and k = 1/2 − u is the internal atomic parameter of the wurtzite lattice (with u = 0.3806 in our simulations to be compared with 0.3819 from low-temperatures X-ray diffraction, 24 see Figure 1). Since O atoms in the wurtzite structure do not occupy centrosymmetric positions, they experience spontaneous polarization with shells displaced from the cores (Δ c−s ≈ 0.0006 in fractional coordinates, which correspond to a (δ) from one surface to the other would compensate the intrinsic ZnO dipole. To account for different possible mechanisms of compensation including formation of vacancies, adatoms, electron and hole charge carriers, and charged adsorbates, we allow for the δ and −δ charges to occupy arbitrary planes nearby both surfaces and separated from them by small distances vc and v′c (see Figure 2). In principle, v and/ or v′ could be zero, i.e., part of the surface. The resultant dipole is given by where p is the dipole per bilayer: which comprises the contribution from cation−anion charge (q) separation and intrinsic polarization of cations (d + ) and anions (d − ). The condition of nonpolarity, P = 0, in the limit of N z → ∞ yields the ideal ratio (Θ) of charge transferred between the two opposing surfaces and the charge of an atomic plane: For the case of ZnO, which we simulate using the polarized shell model (d + = 0), the ideal ratio Θ is We note that Θ depends on two parameters: the electronic polarization of O ions and the distances between atomic charged layers in the c direction, which will vary slightly according to the adopted parametrization of the IP or the details of alternative computational approaches. All real ZnO samples can present different stabilization mechanisms for the polar surfaces with different morphologies and adsorbates, but they will necessarily have to satisfy eq 3.
To implement such a value of Θ in a molecular mechanical model, we need to employ a sufficiently large slab model, which will contain N surf unit cells in the surface plane. We choose to transfer m ions from one surface to the nearest bulk-like empty layer on the other side of the slab (forming a vacancy on one side and an adatom on the other) with the ideal Θ (eq 4) constraining the occupancy of the reconstructed surfaces: (5) or ca. 24% (compare refs 21, 31, 76, and 77). We note that alternative mechanisms are possible for surface reconstructions, which can include formation of vacancies on both sides of the slab (cf. ref 34). Our approach, however, has an advantage in conserving the number of atoms in the simulation cell. Even with relatively large simulation cells that comprise upward of several thousands of ions, the ideal value of Θ still cannot be achieved and a small residual dipole across the unit cell will persist. Our approach raises, therefore, a crucial question: what is the structural and energetic influence of a weak residual dipole? To answer this question, we have performed two exemplar simulations, one with such a residual dipole and another in which the dipole has been fully compensated with point charges placed on both sides of the slab (dipole = 0 method). Details of this study are given in the Supporting Information (SI). The method used here can be also used in surface response to external field, for example, due to surface accumulation layers. We have found that using an uncompensated model results in an incorrect energy ranking (compared to the ideal zero dipole surface terminations) of most stable surface configurations and therefore erroneous surface energies. We conclude that an accurate dipole cancellation is essential and, therefore, only the dipole = 0 method will be used henceforth.
Surface Models. We have used an atomistic approachas implemented in the IP code GULP (general utility lattice program). 78,79 The use of IP allowed us to optimize a vast number of structures generated by our in-house global optimization code, the knowledge led master code (KLMC), 80−82 and identify the lowest energy configurations.
The surface structures were built and optimized by KLMC and GULP, respectively. We used the Born, shell model potentials for ZnO developed by Whitmore, Sokol, and Catlow, which show excellent agreement with a range of experimental data (see Table 1 in ref 75), as with other calculations.
The bulk structural parameters used were a = 3.2518 Å, c = 5.1969 Å, and u = 0.3806 as calculated in our recent publication. 83 The dipole compensation was performed by removing Zn and/or O atoms from the surface and spreading compensating charge uniformly across the bottom of the slab. Using a graphical representation of a residual dipole, we obtain the best integral approximations (of |m Zn − m O |) for the solution of eq 5 as a function of the cell sizesee

Chemistry of Materials
Article are allowed to be moved. For example, a model with thirteen V Zn and seven V O will be as viable as that with just six V Zn for dipole cancellation.
We employ a one-sided 2D-periodic surface model using a two-region approach, which has been widely employed in modeling surface structures with IP based methods. 79,84−86 Our surface models are six layers thick in the c direction (250 atoms per unit cell, ca. 15 Å), the top three atomic layers of which were allowed to relax (region 1), whereas the bottom three were held fixed representing the bulk crystal (region 2). The bottommost (sixth) layer had compensating charge spread uniformly ( Figure 4). The lowest energy structures were tested with a slab twice as thick as the one implemented in our calculations: no significant structural changes were observed, and the average difference in surface energy was below −0.1 J/ m 2 .
Thorough testing has shown that even as good a dipole cancellation as that achieved using the (5 × 5) surface supercell is still not sufficient to guarantee correct energy ranking of competing surface structuressee the Supporting Information. Therefore, two extra layers of fixed point charges were created to compensate for the remaining dipole: 100 points were distributed over a rectangular grid placed 50 Å above and below the slab; 50 Å was deemed large enough to prevent spurious interactions with the surface atoms.
Stochastic Sampling and Optimization. We used Monte Carlo routines as implemented in KLMC for the stochastic sampling of ZnO polar surfaces. The top layer was made of a grid with Zn and O bulk lattice positions where the number of V Zn and V O was specified and KLMC created the surface structures by swapping the vacancies among the grid positions. Two restrictions were specified: KLMC restricted anions from swapping with cations, and vacancies were only allowed to occupy bulk lattice positions in the top layer. For each structure created, KLMC called GULP to perform a structure optimization using the BFGS algorithm.
We tested all the Zn/O (and O/Zn) possible ratios in a (5 × 5) supercell, a total of 26 different occupancies for each surface termination; every initial structure was subjected to zero dipole. The stochastic sampling using Monte Carlo routines as implemented in KLMC probed more than 10,000 different reconstructions for each occupancymore than 500,000 different structures in total. The reproducibility of the distribution of local minima and that of the lowest energy structures found in these simulations has been confirmed by analyzing several independent sample data sets. The need to analyze the whole range of Zn/O occupancies arises from the fact that the concentration of Zn and O ions on each surface depends strongly on the sample preparation conditions, including temperature, annealing time, sputtering time, and energy. 37 Surface Energy. The surface energy (E surf ) of a material is the energy per unit area required to create a surface. Given a bulk energy containing the same number of atoms as the slab, E bulk , and a relaxed energy of the cleaved system, E slab , then the surface energy E surf is defined by where A is the surface area of the cleaved system. In our 2D interatomic potential model, E surf is calculated as bulk unrelax (7) where E unrelax is half of the energy of an unrelaxed cleaved slab where the charge of the top and bottom layer was spread uniformly minus E bulk . The problem with the many-electron techniques including DFT in general is their inability to account for the surface energies of two opposing polar surfaces in an unambiguous way. Conventionally, the energy of cleaving is equally distributed between two resultant surfaces. The surface model is a slab exposing both polar surfaces, one of a particular current interest and another assumed to be "passive", i.e., not contributing to the surface relaxation energy. But such an assumption is wrong for any calculation using self-consistent procedures for the charge density and/or wave function to account for the electronic relaxation from an initial guess, or trial solution. In the case of quadrupolar surfaces (two equivalent opposing dipolar surfaces), these energies could be safely assumed to be equal to each other; however, for any true dipolar system, the electronic relaxation on two inequivalent surfaces within the same calculation is impossible to apportion uniquely. This problem is not pertinent to the method of interatomic potentials where electronic relaxation is fully controlled by the user. In this sense, a rigorous validation of the interatomic potentials using polar surfaces with DFT is impossible.
Canonical Ensemble and Helmholtz Free Energy. Next, we estimate the partition function (Z) that will allow us to calculate the surface free energy. Given the energy of a state i, E i , then the partition function of a system for a canonical ensemble can be approximated as where P m is the number of site permutations for each occupancy: where M is the number of lattice sites, 25 in our simulations; and n = m + 6 for 0 < m < 19 and n = m − 19 for 19 < m < 25. N m is the total number of attempts (10,000 in all our cases); k B is the Boltzmann constant; and T is the temperature. Using P m / N m , a normalizing factor, allows us to extrapolate results from our sampling to complete configurational space (constrained by the choice of initial lattice and interstitial positions) and thus obtain an estimate of the configurational entropy. The Helmholtz surface free energy of the canonical ensemble model (γ m c ) of a cleaved crystal is defined as This model has, however, the restriction that the surface free energy for a given occupancy is calculated as if that particular occupancy is repeated over the space indefinitely, a situation that is not completely realistic since we would expect that, in nature, a combination of different occupancies will be seen on polar surfaces of zinc oxide. The change in occupancy can be tackled using a more general grand canonical ensemble as detailed below.

Article
Grand Canonical Ensemble and Grand Potential. In a model of variable occupancy, we assume that the macroscopic surface consists of a patchwork of local patterns with different occupancies. All macroscopic variables can be obtained in a usual manner by averaging over all local patterns (occupancies). We rearrange eq 8 for a partition function for a occupancy m as follows (cf. refs 87−89): where ΔE i (m) = E i (m) − E 0 (m). Now we introduce the chemical potential (μ) as a surface energy cost of adding/removing a ZnO unit (a pair of Zn and O atoms) to the surface model (in our simulations, this is equivalent to moving these atoms from the bottom to the topmost layer or vice versa). A partial partition function (Z m ) can then be expressed as In this model, the full partition function will be given by the sum of all partial contributions: Thus, we obtain the distribution function over different occupancies (w m ): The observed occupancy can be found as an average ⟨m⟩ : The free energy of the grand canonical ensemble, or grand potential (Ω), is obtained using the full partition function: which allows us to determine the surface grand potential (γ m g ) similar to eq 10: Although the number of particles (ZnO units) is not conserved in this ensemble, and the chemical potential is used instead, it is still convenient to relate observables such as the surface grand potential to the observed number of particles, or occupancy. To this end, we invert eq 15. We note that, in relation to the latter, numerical difficulties impeded the procedure.
In order to find the values for μ (⟨m⟩ ), we have implemented a discrete approach that approximates the value of μ by minimizing the difference between the target ⟨m⟩ T and the ⟨m⟩(μ) and exploits the monotonic character of the latter function. Due to a near singular behavior of ⟨m⟩(μ), standard numerical optimization routines could not be applied and an elementary step based method was used (simulating the Zeno's Achilles' race with the tortoise). Starting from an initial guess and by varying the step size with respect to the ⟨m⟩ T − ⟨m⟩(μ) difference, the μ values corresponding to all the occupancies were found.
Next, we present results of a statistical study of the two polar ZnO surfaces across all occupancies for a given super cell size, employing (i) the canonical ensemble, with the surface stability characterized by the Helmholtz free energy, and (ii) the grand canonical ensemble and its grand potential. Moreover, we identify the atomic structures of the lowest energy surface configurations for each occupancy. More technical details are given in the Supporting Information including the density of states and the number of unique structures seen in each occupancy as a function of temperature.
Density Functional Theory. The periodic DFT code VASP (Vienna ab initio simulation package) 90,91 was employed in this study to corroborate the accuracy of our ZnO IP. VASP uses a plane wave basis set to describe the valence electronic states. Exchange and correlation energy was treated with the generalized gradient approximations (GGA) PBE functional revised specifically for solids: PBEsol. 92 Interactions between the cores (Zn:[Ar] and O:[He]) and the valence electrons were described using the projector-augmented wave (PAW) 93,94 method.
A cutoff of 700 eV and a k-mesh of 1 × 1 × 1 were used, and the iterative relaxation of the ions was not stopped until the forces on the ions were all less than 0.02 eV Å −1 . The slabs were separated by a vacuum gap greater than 25 Å. Zn and O fractional coordinates taken from the IP structures and PBEsol lattice parameters 83 were used. As in the IP calculations, the top three atomic layers were allowed to relax, whereas the bottom three were held fixed. The corresponding Zn and O ions belonging to the sixth layer were placed randomly across their lattice sites.

■ RESULTS AND DISCUSSION
To establish the mechanism stabilizing polar ZnO surfaces, we have employed unbiased stochastic sampling techniques on the energy landscape defined by accurate interatomic potentials (IP). 75 A two-dimensional (2D) periodic surface supercell model was used. The size of the supercell is an important variable in surface reconstruction as it constrains the extent of reconstruction patterns and the lowest dipole that can be achieved. Figure 3 shows the smallest possible dipole values for supercell sizes up to 400 unit cells. As noted, the (5 × 5) supercell is used for all our calculations as it gives a good balance between size and the magnitude of the dipole, and is also supported by STM images. 21 Attempts to reduce the dipole, e.g., to 10 −6 e Å, would require a computationally infeasible (1000 × 1000) supercell.
The procedure for canceling the dipole after the ionic transfer is illustrated in Figure 4. Our calculations showed that initial structures with weak and no dipole differ in energetics and structure (further details are given below and in the SI), which suggests that an accurate dipole cancellation is essential. Compensating the residual dipoles with point charges allowed us to reduce their values below 10 −6 e Å in all initial configurations.
Canonical Ensemble and the Helmholtz Free Energy. The calculated surface energies (E surf , Equation 7), using the global minimum configurations, are shown in Figure 5 (blue line, T = 0 K) as a function of occupancy, for both polar As we have statistically sampled a significant number of surface structures, here, we gain an insight into which structures are adopted in nature and how they evolve with temperature, first under the constant area (a 2D analogue of 3D volume) conditions, using the Helmholtz free energy as a measure of surface stability (eq 10). We evaluate the extrapolated total number of accessible configurations in each occupancy, using eq 9 for lattice site permutations, which has a maximum near 50% surface occupancy and a minimum at both maximum and minimum occupancy. We note that the minimum occupancy corresponds to one ZnO layer stripped off the surface and the maximumto one ZnO layer grown on the surface. The midrange occupancy offers most opportunities for Zn and O ions to change their lattice positions. In most cases, the number of configurations, especially near 50% occupancy, is too large for an exhaustive evaluation; and we applied Monte Carlo sampling (or global optimization) techniques with the expectation that the sample evaluated is large enough to represent the thermodynamic limit.

Chemistry of Materials
Considering the tentative lowest surface energy and the Helmholtz surface free energy as a function of occupancy, the global and local minima, see Figure 5, suggest stable and metastable occupancy at zero and higher temperatures. These profiles are similar to those drawn for solid solutions, 97,98 and likewise can be interpreted in a similar manner. A maximum of the free energy of a particular element ratio (R, for instance), or occupancy in our system, suggests a miscibility gap, i.e., regions of the solid have ratios larger and smaller than R (with the average of R over all regions). The influence of temperature on the surface structure was evaluated by the Helmholtz surface free energy (F surf ) ( Figure 5). Assuming our sampling for each occupancy is sufficient, we found that the strongest effect is for the mid-range occupancies, which is expected as there are more local minima (LM) here; moreover, the change in F surf is more pronounced for mid-range concentrations as the temperature increases. The size of the supercell is restricted by our computer resources which is not a constraint for nature; for global minima (GM) dominated free energies we would expect the structural pattern of the global minima to be repeated, whereas for the other occupancies the surface structure will be a montage of all the low-energy local minima. For large and complex energy landscapes, locating the global minima by straight Monte Carlo methods is difficult. 99 Although more sample points were taken across the energy landscape for midrange concentrations, the resulting percentage occupancy is still smaller than achieved for high and low concentrations (points at the edge of the graph). However, as our results appeared to have converged with respect to additional sample points, we assume that as the low-energy local minima are difficult to find, their contributions to F surf will also be small.
Grand Canonical Ensemble and Grand Potential. We recall that in the choice of the canonical ensemble we constrain the surface occupancy, or the numbers of Zn and O ions in the topmost surface layer in the simulation cell (5 × 5). In the absence of a known strictly periodic order, the simulation cell

Chemistry of Materials
Article size severely reduces the freedom of ions to adopt such configurations or structural arrangements, patterns, that would minimize the surface energy. Whereas a brute force increase in the cell size is not feasible with current computer resources, we can simulate the large scale disproportionation of ions between areas of different local occupancy by using the grand canonical ensemble method, described in Methods and Computational Details. Combining patterns from the left and right of a maximum in Figure 5 with the constraint of a fixed average occupancy, ⟨m⟩, will effectively, or implicitly, sample configurations of larger unit cell sizes than that of m = ⟨m⟩ already explicitly modeled. By doing so it may be possible to lower the free energy for that degree of occupancy. To ascertain the best combinations, we need to compute just one general partition function, bringing together all surface structures across the whole range of occupancies, in the presence of the chemical potential, μ. As described in Methods and Computational Details, this partition function is used to calculate the average degree of occupancy for any value of μ (eq 15) and the corresponding grand potential (eq 16).
The availability of additional or lack of Zn and O atoms, which can be related to the value chosen for the chemical potential, will affect both the surface pattern and the degree of occupancy. In Figure 6, we show the calculated relationship between the chemical potential and the average occupancy of the topmost layer. As expected, the degree of occupancy decreases with the chemical potential increase, since it is more favorable for Zn and O to bond to the surface. Looking more carefully at the lowest temperature graph (Figure 6), occupancy for the O surface depends on the sign of μ: a positive value of the chemical potential produces zero occupancy (i.e., with only six Zn adatoms on top of a complete layer); otherwise there is a drive to a maximum occupancy of the topmost layer (⟨m O ⟩ = 25 with six Zn adatoms on top of a new layer). Hence, the sign of μ determines whether our simulations predict that the O surface will dissolve or grow with complete layers covered by a few adatoms of opposite polarity being more stable. The Zn

Chemistry of Materials
Article terminated surface, in contrast, is more interesting, as we discover five possible degrees of occupancy. For a small positive value of μ the occupancy stabilizes at 17, which corresponds to the global minimum of the surface free energy (see Figure 5 at 0 K). The degree of occupancy increases by one when μ becomes negative, a occupancy that remains stable for a significant range of negative μ, after which the average occupancy increases to 23, then 24, and finally, a maximum occupancy for μ < −0.025 J/m 2 . Similar ranges of the chemical potential are predicted for the average occupancies of 17 and 18; this range is twice as wide for 24 and an order of magnitude lower for ⟨m⟩ = 23. An appearance of a plateau in this graph indicates a metastable state (degree of occupancy) and is associated with a slower growth, which will be observed in experiments for this surface termination. Apart from the only change in occupancy that occurs for a positive chemical potential, there is no obvious link between the local minima of the surface free energy as a function of occupancy and the predicted values of average occupancy.
To gain a further insight into observed surface structure and emerging atomic patterns, we consider next the occupancy distribution functions (w m , eq 14) that give us access to probabilities of finding occupancies m Zn and m O of the Zn and O terminated surfaces spanning a range of chemical potential, which covers all possibilities from an empty to a complete topmost surface layer (see Figure 7). For the Zn terminated face, we note that, at very low temperatures, contributions from m Zn = 25, 24, and 17 and m Zn = 0 and 17 are identified for negative and positive values of μ, respectively. As the temperature increases, the occupancy distribution function smooths around these peaks, and for very high temperatures, only contributions from middle occupancy ranges are expected, as seen in Figure 6.  Figure 5.
We note that the structural patterns we have sampled are constrained by the choice of the simulation cell. A number of the peaks in our calculated surface free energy as a function of occupancy will be removed by considering mixtures (if a larger unit cell is used, then a new pattern can be constructed with the chosen degree of occupancy from patterns from smaller unit cells that have larger and smaller degrees of occupancy, although in our calculation of the grand potential we do not worry about boundary effects between patterns). From the insight gained from Figure 7, we can conclude that the peaks (seen in Figure 5) in the surface free energy at m Zn = 12 and 15 are artifacts of the chosen simulation cell size. Likewise, the oscillations of the surface energy between m O = 9 and 11, and the peak at m O = 14 are also artifacts. Returning to Figure 6, the sharp edges of ⟨m⟩(μ) for both Zn and O terminated surfaces are smoothed out as temperature increases. Furthermore, it is noticeable that the first nonzero occupancy for the Zn terminated surface decreases (right most plateau lowers in height), whereas for the O terminated surface a plateau centered at μ = 0 develops at ⟨m⟩ = 12. Although the latter does not correspond to the global minimum (GM) in the surface free energy as reported in Figure 5, it is similar in value to the GM, and once artifacts caused by constraints of the unit cell size used in simulations are removed, then we can easily

Chemistry of Materials
Article imagine ⟨m⟩ = 12 is centrally located in the superbasin that it resides in. For extreme temperatures (see Figure 6 for 10,000 K), all plateaux disappear, and the chemical potential of a greater magnitude is required to generate minimum or maximum occupancy of either surface.
Finally, the surface grand potential (γ g , eq 17) is shown in Figure 8 as a function of the chemical potential. For all chosen temperatures, we observe a monotonic increase with μ. At low to moderate temperatures that are expected to be used in experiment, the surface grand potential for the Zn face plateaus from above ca. 0.007 J/m 2 , whereas, for the O face, it reaches a constant value at about 0 J/m 2 . At very high temperatures, illustrated here with T = 10,000 K, the chosen range of μ is not sufficient to maintain the surface occupancy above 13 and below 7 (see Figure 6) on both faces, which is of course only hypothetical as the temperature is way in excess of the ZnO evaporation point. The steepest rise of the surface grand potential of the Zn face corresponds to the reduction of occupancy from 25 to 17. As the Helmholtz surface free energy passes through the minimum, the slope of the surface grand potential becomes more gentle. The simpler shape of γ g (μ) with only one rise and a plateau can be related to the dominant character of the central super basin of the Helmholtz surface free energy seen in Figure 5.
For the Zn face, γ g stabilizes about 0.05 J/m 2 above the ground state minimum of γ c that corresponds to a layer of ZnO stripped off the surface. For the most negative values of μ and the complete occupancy of the topmost layer of ZnO, γ g lowers to about 0.7 J/m 2 , which is even lower than the surface energy of the nonpolar surfaces of ZnO 83 (which is however taken in the athermal limit with μ = 0). The O face proves to be more stable than the Zn face across the whole range of thermodynamic parameters. Based on our thermodynamic analysis of the behavior of the Zn terminated surface at low temperatures, there is a considerably wide range of the chemical potential on both sides of zero where the occupancy is about m Zn = 17 and 18. For more extreme chemical potentials, the preference is for a fully occupied topmost layer with a few Zn vacancies (negative μ) or unoccupied topmost layer with a few O adatoms (positive μ). In contrast, the O face at low temperatures only shows a fully occupied or unnocupied topmost surface layer with O vacancies or Zn adatoms, respectively. Only on a significant temperature increase, an appreaciable but still narrow range of μ (around zero) starts to support an intermediate m O = 12 occupancy. Our calculations, therefore, predict that slight changes to the synthesis conditions would lead to different structures as reported by Torbrugge et al., 25 Beinik et al., 46 Lauritsen et al., 21 Woll, 1 Batyrev and van den Heuvel, 47 and Valtiner et al. 45 Moreover, even for the same synthetic conditions, more than one structure or motif can become thermodynamically accessible (see Supporting Information) and, therefore, as mentioned, the morphological structure of a surface may vary dramatically in the same sample.
Next, we consider prevalent patterns of the surface atomic structure that can be predicted for both faces based on the presented thermodynamic analysis.
Atomic Structure. Figure 9 shows the global minima structures predicted for the Zn face and Table S1 shows the average atomic relaxation for the first two atomic double layers for both terminations. Three patterns are identified: stripes, triangles (with O edges as seen in experiment), and randomly sited vacancies. These findings show good agreement with scanning force microscopy work by Torbrugge et al., 25 where most of the surface area is covered by triangles of different sizes but also holes and rows of missing atoms running along the [100] direction. The (1 × 3) reconstructions reported in ref 25 are essentially similar to what is presented in Figure 9 (m Zn = 14), where the model is constructed by rows of missing Zn and O atoms aligned along the [100] direction. However, the symmetry reported in this study is (5 × 5), which accounts for a complete removal of the dipole moment, whereas in the (1 × 3) model by Torbrugge et al., H adsorption contributes to the stabilization. 100 We note that samples in ref 25 that show stripy surface features were annealed at 1150 Kthe need for annealing is rationalized by our calculations that predict their first occurrence at temperatures above 500 K, and clear persistence at 1000 K (see Figure 7). At lower and higher temperatures, we observe more disordered structures with less pronounced patterning. These features are characteristic for μ > 0, which corresponds to sample preparation and/or utilization under vacuum conditions (low partial pressures of both Zn and O) also in agreement with the report by Torbrugge et al. On  there is no vacancy either in the first or second layer. The first layer presents smooth lateral displacements with six oxygens on top in interstitial sites. Triangular reconstructions of different sizes are observed among the occupancies. We noted as well that most of our reconstructions are surrounded by oxygen edges.

Chemistry of Materials
The O face predicted global minima structures are shown in Figure 10 and include small islands, hexagons and triangles. In contrast to the Zn termination, there is no ionic exchange between layers. For a range of 0 < m O < 6, the lowest energy structures show small islands spread all over the surface with Zn interstitials. For m O = 0, 2, hexagonal patterns appear formed with a combination of zinc and oxygen interstitials; the smaller hexagon is regular with a side length of 5.65 Å, whereas the bigger hexagon has two sides of 11.3 Å and four of 8.59 Å. At m O = 6, there is the first well-defined triangular pattern, which increases smoothly until m O = 11 and 12, where well-defined big triangles connected by oxygen interstitials are formed. Although triangular and hexagonal shapes appeared for some occupancies, we do not observe a preference for them: those reconstructions are a mechanism of stabilization but not the only one. The preferred configuration, as seen in experiment, will depend on synthesis conditions and surface occupancies.
DFT calculations were performed for six different coverages (m Zn = 14, 17, and 22 and m O = 0, 9, and 11) to validate the quality of the interatomic potentials used in this study. Reconstruction patterns have been preserved, and only minor atomic relaxations were observed below 0.083 Å/atom. Relaxations in the surface plane were minor (below 0.037 Å/ atom). The strongest relaxations were observed along the surface normal direction (the hexagonal c axis), with a maximum of 0.068 Å/atom. Similar behavior has been reported in our study of ZnO nonpolar surfaces. 83

■ CONCLUSIONS
To summarize, we report the first unbiased computational study of the stability of the polar ZnO surfaces. Accurate cancellation of the surface slab dipole is key to gain reliable surface models of differing occupancies. The lower polarizability of the Zn face over the O face results in the higher surface energy (γ c ). The O face has a smaller γ c range compared to the Zn face (by 24%), implying a higher probability of disorder, which rationalizes the observed (1 × 1) periodicity 96 and rates of growth. Surface free energy calculations support this conclusion and point to stronger ordering trends for the Zn face. Structurally, the two ZnO polar surfaces behave

Chemistry of Materials
Article differently. The Zn face undergoes significant reconstruction with some Zn ions hopping from the second to the first layer. On the O face, the movements are weaker: the ionic transfer between layers only happens when the top Zn layer is full (m O > 19), and every O ion fills a vacancy. For both surfaces triangular reconstructions are common. The methods we have developed allow complex surface structures of oxides to be successfully predicted.
Considering the surface energies predicted with the method of the grand potential, we observe dramatic stabilization due to the availability of many microscopic states of similar energy, which can be traced back to the chosen large size of the surface simulation cell. Our choice was influenced by the necessity of having to remove the inherent dipole; i.e., the bulk structure of wurtzite ZnO determines the ideal ratio of charge transfer and therefore optimum surface plane occupation that is incompatible with smaller surface cells. We found that with the optimum surface occupancies the large size surface cell is incompatible with unique high-symmetry structural motifs, and instead there is realized several competing motifs. Ideally, we would have liked to have used an even larger surface cell size, e.g., (19 × 19), where the resulting dipole is substantially smaller than that achieved in this study within the (5 × 5) . There will be more possible patterns, or microscopic states, in a larger cell, and we assume that again there will be many low-energy configurations composed of several competing structural motifs. From the grand potential, however, we have considered larger patterns that may be possible in much larger surface cell sizesalbeit ignoring the cost of mismatching across boundarieswhich confirmed that both surfaces are highly disordered and disorder is the source of the stability of polar surfaces.
Nevertheless, certain reconstruction patterns observed experimentally under differing conditions are clearly recognized in our global minima optimizations, including the hexagonal patterns on the O terminated surface in agreement with the work of Lauritsen et al.; 21 and on the Zn terminated surfaces triangular patterns, both disordered as in the work of Dulub et al. 30 and with striped reconstructions reported by Torbrugge et al. 25 The appearance of partially ordered atomic arrangements can be associated with particular values of the chemical potential and associated surface energy distributions (partition functions). Further analysis of such phenomena is necessary. This is the first grand canonical ensemble study that allows us to investigate stable surface structures that emerge on crystal growth under thermodynamic equilibrium, but in this work we constrain our model considering only one layer of ZnO. A more general study focusing on crystal morphology, which goes beyond this approximation and considers the effect of interaction between adjacent areas of varying occupancy, is in progress and will be reported in the near future.