Finite-size correction for slab supercell calculations of materials with spontaneous polarization

The repeated slab approach has become a de facto standard to accurately describe surface properties of materials by density functional theory calculations with periodic boundary conditions. For materials exhibiting spontaneous polarization, we show that the conventional scheme of passivation with pseudo hydrogen is unable to realize a charge-neutral surface. The presence of a net surface charge induces via Gauss’s law a macroscopic electric field through the slab and results in poor size convergence with respect to the thickness of the slab. We propose a modified passivation method that accounts for the effect of spontaneous polarization, describes the correct bulk limits and boosts convergence with respect to slab thickness. The robustness, reliability, and superior convergence of energetics and electronic structure achieved by the proposed method are demonstrated using the example of polar ZnO surfaces.


INTRODUCTION
Density functional theory (DFT) has evolved into the work-horse approach for electronic structure calculations allowing to study of almost all aspects of materials and their properties. A common assumption underlying most DFT codes is periodic boundary conditions (PBC). These conditions naturally reproduce the translational symmetry of crystals to describe an infinite bulk crystal. Thus only the unit cell has to be represented, allowing for highly efficient calculations. PBC can also be easily extended to study structures where translational symmetry is broken. For surfaces, where translational symmetry applies only along the lateral dimensions, a repeated finite-slab configuration is used. Similar approaches are used for point and extended defects 1 . Common to all these approaches is that one or more dimensions that are infinite are approximated by finite cell dimensions. Thus, a prerequisite to obtaining accurate results is to ensure that the size of this dimension(s) is sufficiently large by performing careful convergence checks. Otherwise, finite-size errors may become substantial and obscure the results.
A major reason why PBC have become so popular and successful is that finite-size errors quickly decay even for modest system sizes. For example, accurate surface calculations can be performed using slabs consisting of only a few (typically 5-10) atomic layers 2,3 . Finite-size convergence, however, dramatically slows down when long-range interactions are present. Such Coulombic interactions are present in various technologically relevant systems. Prominent examples are charged defects at surfaces [4][5][6] in the context of catalysis or Fermi level pinning at semiconductor devices, and electrified surfaces in strong electric fields [7][8][9] or in electrochemical environments [10][11][12] . Indeed, slow convergence for such systems has been reported 13 . The origin of the slow convergence can be easily understood when considering that a slab with two oppositely charged surfaces represents electrostatically a sheet capacitor. Thus, a constant electric field throughout the entire slab arises which does not decay when going away from the surface toward the bulk. This asymptotic behavior of surfaces is fundamentally different from that for charged point defects where the potential decays inversely proportional with the distance to the defect.
To overcome these extremely slow convergence rates for charged surface calculations, two main strategies are followed. The first strategy is to derive and compute the asymptotic behavior of the target quantity, usually the total energy E tot , as a function of slab thickness. The converged limit can then be simply obtained by the finite-size energy plus the change of the asymptotic function from that size to infinity. This strategy has been originally developed for and highly successfully applied to charged point defects 1,[14][15][16] . It has later been extended to surfaces 4,6,17 . The second strategy is to employ boundary conditions such that the correct asymptotic limit is already reproduced in the finite-size supercells. For surfaces, this strategy has been successfully applied to remove long-range electric fields that arise when studying unpassivated polar surfaces 9,18,19 .
In the present study, we will follow the second strategy and derive a concept that efficiently overcomes the extremely slow 1/d convergence (with d the slab thickness) when computing surface properties, such as surface states or surface energies, for systems exhibiting spontaneous polarization. This type of polarization is present in many technologically relevant semiconductors and insulators. A prominent example are the group-III-nitrides, which are the technological backbone of energy-efficient lighting 20,21 . Also, many oxides discussed in the context of topological insulators 22 , opto-and micro-electronic devices 23,24 , catalysis [25][26][27] , and other applications exhibit spontaneous polarization. Thus, schemes that boost size convergence for these technologically important and scientifically exciting material systems are highly desirable.
While the discussions and the proposed methodology are generic and can be applied to any material exhibiting spontaneous polarization we will focus on ZnO as a prototype system, because of the wealth of previous studies that exist for ZnO surfaces 2,27-29 .
The paper is organized as follows. Section "Convergence behavior in the conventional scheme" shows calculations performed using the conventional scheme to passivate partially filled surface dangling bonds (db). This scheme is shown to have deficiencies when applied to materials exhibiting spontaneous polarization. In Sections "Origin of convergence failure", "Fundamentals of the proposed scheme", and "Implementation and parameterization of the proposed scheme", we analyze the origin of these deficiencies, propose a scheme to overcome them, and give practical advice to its application. Finally, in Section "Applications: surface phase diagram" and in Section "Applications: electronic structure" the proposed scheme is applied to study surface properties of ZnO, such as surface phase diagrams and surface band structures.

RESULTS AND DISCUSSION
Convergence behavior in the conventional scheme The common approach in DFT is to model surfaces using a repeated slab approach. Within this approach, the translational symmetry along the direction perpendicular to the surface is broken by the introduction of a vacuum region within the supercell. To decouple consecutive slabs, as well as the two surfaces bounding the slab, the vacuum thickness has to be carefully controlled to guarantee size convergence. If the thickness of the substrate or epitaxial film to which the surface is attached is thicker than the DFT slab, which is the case for typical technological applications and devices, also convergence with respect to the slab thickness needs to be carefully controlled. We note that when describing thin nanostructures consisting of only a few atomic layers or 2D materials the slab thickness of the real structure and the DFT slab can be exactly matched. In this case, the buildup of electric fields through the slab and charge transfer between front and back surface are real, i.e., slab size convergence and related corrections should not be applied.
Since surface properties are commonly given for the case of zero external electric fields, we employ the dipole correction 30 to remove the spurious external electric field that arises when the front and back surfaces of the slab is structurally inequivalent and have different work functions. We note that the dipole correction can be also used to model a slab embedded in an external electric field. In particular, its generalized form 9 is now being used to study surfaces even in extremely high electric fields 31 .
Creating a slab for a wurtzite (WZ) material, like ZnO, by cleaving the bulk perpendicular to its polar direction (the c-axis) results in two nonequivalent surfaces. These are the ZnO(0001) surface, which is terminated by Zn atoms, and the ZnOð0001Þ surface, which is terminated by O atoms. Because all partially filled Zn dbs are located on one side of the slab, while all partially filled O dbs are located on the other side of the slab, the charge is transferred from the Zn to the O side of the slab. To prevent this spurious charge transfer we first apply the conventional scheme and passivate the partially filled O db on the ZnO(0001) surface of the slab by pseudo hydrogen atoms 18 of valence 1/2 (denoted in the following as psH 0.5 ).
For many semiconducting and insulating materials, this scheme shows a highly efficient size convergence 18 . However, as we demonstrate in the following, materials exhibiting spontaneous polarization exhibit a poor size convergence. We consider the ZnO set-up described above, in which the O-terminated side of the slab is passivated by psH 0.5 . For the WZ ZnO(0001) surface of the slab we focus on two structures that satisfy electron counting (EC), i.e., structures in which all low-lying db states (usually states of the anion) are filled and all high-lying db states (usually cation states) are empty 32 . Structures that satisfy EC represent local minima on many surfaces and can be constructed simply by counting the number of electrons in the db states. Here, we use the (2 × 2)-O ad and the (2 × 2)-V Zn structures, in which an O atom is adsorbed on, respectively a Zn atom is removed from, a (2 × 2) surface unit cell of the ZnO(0001) surface. For these (and all other cases reported in this work) we perform DFT calculations using the Vienna Ab initio Simulations Package (VASP) 33,34 with projector augmented waves (PAW) 35 . Both the generalized gradient approximation due to Perdew, Burke, and Ernzerhof (PBE) 36,37 and the hybrid functional due to Heyd, Scuseria, and Ernzerhof (HSE06) 38,39 are used for exchange-correlation functionals (for further details on the set-up see Section "Methods").
An important physical quantity in surface science is the energy difference between two surfaces. This difference provides direct information about which surface is thermodynamically more stable and is thus the key quantity to construct surface phase diagrams. The energy difference between the two EC structures of interest is defined as (1) where α indicates the bulk phase (e.g., WZ) and E ZnO tot accounts for the fact that to convert the O ad into a V Zn surface an O and a Zn atom, i.e., a ZnO formula unit, have to be removed. The convergence of ΔE α with respect to slab thickness is shown in Fig. 1a. An accuracy of 1 meV/uc (meV/surface unit cell) is typically needed to identify the thermodynamically most stable surfaces. The two green points, obtained for the thinnest considered slabs consisting of 6 bilayers (BL), i.e., six Zn and six O atoms in a (1 × 1) unit cell, respectively seven BL, give the impression that the energy difference is already converged to within 1 meV/uc. However, increasing the slab thickness (blue points) up to 14 BL (the thickest slab we considered) reveals a qualitatively very different convergence behavior. In fact, a linear fit (dashed blue line) indicates that a slab thickness d slab of 61.2 Å (corresponding to more than 23 BL) is required to converge the energy to within 1 meV/uc. The resulting surface energy difference ΔE WZ is vastly different (by more than 10 meV/uc) from the one obtained for the conventionally sized 6 BL (7 BL) calculations (green line).
These results clearly demonstrate a poor size convergence which directly translates to large system sizes that are computationally expensive for DFT. Particularly in the context of surface studies, where lateral dimensions of reconstructions can be very large, this would lead to calculations that become practically intractable. For example, the study of surface reconstructions on ZnO(0001) surfaces requires (at least) a ( ffiffiffiffiffi 48 p ffiffiffiffiffi 48 p ) surface unit cell 40 , resulting in a DFT calculation with more than 2000 atoms to achieve a 1 meV/uc accuracy in the surface energy.
The kink in the energy dependence observed in Fig. 1a for WZ ZnO is directly related to the change in convergence behavior discussed above. Using the averaged electrostatic potential we can directly evaluate the electric field as a function of position in the slab; the technical details are described in Section "Implementation and parameterization of the proposed scheme". Plotting the resulting electric field values at different slab thicknesses for a slab with the (2 × 2)-O ad surface reconstruction (Fig. 1b) reveals a qualitatively similar convergence behavior: there is again a kink in the curve. Extrapolation of the values (blue dashed line) shows that the field approaches zero as the slab thickness approaches infinity.
The field through the slab can also be visualized in the layerresolved density of states (LDOS). We show electronic structures and electrostatic potentials calculated with HSE06 since the effects are better visible due to the improved description of the bandgap. The LDOS together with the total DOS and the averaged electrostatic potential for a ZnO(0001) (2 × 2)-O ad slab (Fig. 2a) are shown in Fig. 2b and c for two different slab thicknesses (9 BL and 13 BL). The presence of the electric field tilts the band edges. As a consequence, the effective overall bandgap of the system decreases significantly (Fig. 2b) or even vanishes (Fig. 2c), as seen in the total DOS. Beyond a critical thickness, the valence-band maximum (VBM) on one side of the slab and the conduction-band minimum (CBM) on the other side of the slab overlap (Fig. 2c). This overlap, which leads to a disappearance of the bandgap, is referred to as breakdown: electrons from the VB at the right surface are transferred into the CB at the left surface.
The quantitative changes observed as a function of slab thickness indicate that within the conventional calculation scheme, extremely large slab thicknesses are required to obtain converged results. The occurrence of breakdown in such thick slabs also makes it impossible to analyze the electronic band structure directly. Even more problematic is the fact that different convergence regimes occur, which may give the incorrect impression that results for small slab thicknesses are already converged, but they actually extrapolate to the wrong limit. In the next section, we will shed light on the origin of this behavior by deriving asymptotic bulk limits, which a slab set-up has to emulate. Based on those insights we will propose a scheme to achieve the correct limits.

Origin of convergence failure
In this section, we derive the asymptotic limit of electronic structure for infinitely thick slabs (semi-infinite systems). We focus on two representative and well-studied polar surface orientations (Fig. 3): the (111) orientation for zinc blende (ZB) (no spontaneous polarization) and the (0001) orientation for the WZ structure (with spontaneous polarization). These close-packed surface structures show an identical atomic structure in the top surface layer; the difference is that ZB ZnO does not have a spontaneous polarization. Indeed, the symmetry of the ZB structure is such that all nearest-neighbor bonds are equivalent. At the bulktruncated surface, where the cut is performed such that the number of cuts (broken) bonds is minimized, each atom has one db pointing away from the surface.
Defining Z c and Z a as the valencies of a cation, respectively anion atom, we note that in the case of the ZB structure the O dbs located on one side of the slab are partially filled with Z a /4 (i.e., 3/ 2) electrons, while Zn dbs on the other side of the slab is filled  The computed electrostatic potential energy profile and the layer-resolved density of states (LDOS) for slabs with this atomic arrangement. The slab thickness d slab has been chosen to show b no breakdown (9 BL) and c breakdown (13 BL). The planar-and the macroscopic-averaged electrostatic potential energies are shown as blue and red solid lines, respectively. The LDOS along the z axis (axis perpendicular to the slab surfaces) is resolved for each Zn-O BL and shown as a gray region. The energy of the highest occupied state is set to 0 eV.
with Z c /4 (i.e., 1/2) electron, as shown in Fig. 3a. Because cation dbs have higher orbital energy compared to anion dbs, in a finitesize slab a charge transfer from the Zn to the O side of the slab occurs, denoted as q t in Fig. 3b. This results in two oppositely charged surfaces that form a sheet capacitor. As a consequence of the charged capacitor, a macroscopic electric field arises, as shown in the band structure in the lower part of Fig. 3b. The schematic picture is confirmed by an actual DFT calculation (cf. Supplementary Fig. 1). After charge transfer both ZB ZnO surfaces fulfill EC 32 , i.e., all cation/anion dbs states are completely empty/filled. However, the slab surfaces are no longer charge-neutral.
The macroscopic electric field induced by the charge transfer results in canted valence and conduction bands (Figs. 3b and 4a). Fig. 3 Schematic picture. Schematics illustrating how to construct and passivate an asymmetric semiconductor slab within a DFT code in order to mimic a semi-infinite bulk system. Panels a-d describe the case of a material not exhibiting spontaneous polarization, while panels e-i describe the case of a material exhibiting spontaneous polarization. Corresponding band diagrams are shown right below the schematic of the slab: for a finite slab formed by truncating the bulk before charge transfer from cation dbs to anion dbs (a, e) and after charge transfer (b, f), for a semi-infinite bulk structure (c, g) and for a slab where the left side is passivated by the conventional (d, h) and the proposed (i) scheme. Surface atoms are shown as balls containing a chemical element (O or Zn). The green/cyan colors within a db at the surface indicate filling, while parts remaining white/gray indicate missing electrons (i.e., empty dbs). The corresponding amount of electrons in a db is written in the green box in units of electrons (e − ). Passivation of the left side of the slab is denoted by psH. In each schematic for the structure of a slab the core charge (white rectangle with Z tot ) and the valence electrons (blue rectangle with N e ) are shown as separate entities to visualize the charge displacement. q t , q p and q s stand for transferred charge, polarization charge and surface charge, expressed in units of elementary charge (e). The valence of the pseudohydrogen atom (Z psH ) is shown as a "blow up". "EC" stands for electron counting.  1) five BL slab to which no passivation is applied and breakdown has occurred, i.e., the charge has been transferred from one side of the slab to the other. b WZ ZnO(0001) (2 × 2)-O ad 13 BL slab for which the O dbs on the left surface has been passivated using the proposed passivation scheme. The planar-averaged and the macroscopic electrostatic potential energies are shown as blue and red solid lines, respectively. The LDOS along the z-axis (axis perpendicular to the slab surfaces) is resolved for each Zn-O BL and shown in gray. The energy of the highest occupied state is set to 0 eV.
As a consequence, the effective bandgap is dramatically reduced and may even disappear. This sets an upper limit on the charge transfer q t : once the occupied valence band on the left side moves above the empty conduction band on the right side (which formally corresponds to the formation of a negative band gap) a back transfer of the electronic charge from the left to the right surface occurs. This back transfer of charge reduces the magnitude of the electric field. In the limit of an infinitely thick slab, the full charge has been transferred back. As a consequence, the field disappears and the Zn/O dbs on the surface become again partially occupied by 1 2 e À / 3 2 e À (Fig. 3c). Therefore, following the conventional strategy to passivate the O surface of a ZB ZnO(111) slab by providing the missing 1/2 electron in an O db by a pseudo hydrogen atom 18 with valence 1/2 (psH 0.5 ) correctly reproduces the asymptotic bulk limit (Fig. 3d): there is no electric field in the slab and both surfaces are charge-neutral. The absence of a field that depends on the slab thickness is the fundamental origin of the fast convergence behavior for ZB shown in Fig. 1a (orange line).
We now discuss changes that occur when spontaneous polarization is present using the prototype example of WZ ZnO. The induced changes are summarized in Fig. 3e-i. Spontaneous polarization in WZ originates from a shift in the electron cloud, relative to the ZB structure, leading to a negative charge q p on the surface as shown in Fig. 3e. This charge corresponds to the difference in formal polarization P f between the WZ and the ZB structure, i.e.,

Fundamentals of the proposed scheme
Having identified the polarization charge q p as the origin of the failure of the conventional passivation scheme we will now derive (i) the correct limit for infinitely thick slabs and (ii) extend the conventional passivation scheme to systems exhibiting spontaneous polarization. To derive the correct asymptotic limit we use the same arguments as in the polarization-free case: in this limit, any macroscopic field must disappear, i.e., the slab surfaces must be charge-neutral. For cases with spontaneous polarization, this means that the charge displacement q p due to spontaneous polarization must be compensated for so that the surface attains a net charge of 0.
In order to achieve a field-free structure (Fig. 3g), we need to remove the charge q p from the (0001) surface. This can be accomplished by transferring this charge to the left surface [the Oterminated (0001) surface], to which end we modify the pseudo hydrogen atom (denoted in the following as psH * ) such that the psH-O bond seen in Fig. 3h becomes an acceptor bond 41 . Note that the modified pseudohydrogen atom psH * remains a chargeneutral object, i.e., valence and nuclear charge are equal, and as such does not induce any uncompensated free carriers and/or a compensating background charge in the supercell, which remains neutral. Practical advice on how to adjust the charge of the pseudohydrogen is given in the Supplementary Notes Section B. Before charge transfer, such a bond is not completely occupied and thus can accept a fraction of an electron. For the specific case of a bond between psH * (charge 1 2 þ q p e ) and oxygen (which contributes 6/4 = 1.5 electrons to each of its four bonds) the occupation is two plus the polarization charge, i.e., Since the energy level of this bond is below the VBM it will be filled to a total of two electrons, giving rise to a net negative charge q p . By construction, this charge is the exact opposite to the polarization charge at this surface. The two charges thus exactly cancel and make the surface charge-neutral (Fig. 3i). In contrast to the conventional compensation scheme, the valence of the pseudo hydrogen is no longer simply a function of valence and bond coordination of the surface atom. Instead, the polarization charge, which is an inherent bulk property, also enters (cf. inset in Fig. 3i).
The success of the proposed scheme in eliminating the macroscopic electric field is seen in Fig. 1 and Fig. 4b. Figure 1 shows a dramatically improved size convergence: the energy difference is converged to within 1 meV/uc already for the thinnest considered slab of six BLs (Fig. 1a). Indeed, as shown in Fig. 1b there is no electric field left. The values for the energy difference between the two considered reconstructions (Fig. 1a) obtained by extrapolation to an infinitely thick slab (blue line) and by applying the proposed passivation scheme (red line) agree to within 1.5 meV/uc. The proposed scheme also overcomes the artifacts in band-structure and DOS calculations: band edges are no longer canted but perfectly flat (Fig. 4b) revealing again the absence of an electric field. Furthermore, the surface state originating from the partially filled db of the adsorbed oxygen atom on the Zn-terminated surface, which was difficult to discern in the slab passivated with psH 0.5 (Fig. 2c), is now clearly visible as a peak in the DOS at~0 eV in Fig. 4b. A similarly improved convergence behavior is observed for the energy, structural properties, and electronic structure of GaN and AlN when using the proposed scheme, as seen in Supplementary Figs. 1-4.
An interesting implication of introducing valencies that are not a simple integer fraction, such as 1/2 or 3/2, is that surface reconstructions that satisfy the EC rule, which requires bonds with integer fillings (either 0 or 2), can only be obtained for large surface unit cells. Specifically, the product q p ⋅ N surf (with N surf the number of surface atoms in the cell) needs to be an even integer. In the example discussed here, q p = −0.020, and therefore a surface unit cell with about 100 surface atoms would be required to achieve EC.
The proposed scheme is readily implemented in DFT calculations: all that is needed is, in the case of passivating an O db, to replace the pseudo hydrogen with valence ( 1 2 by pseudo hydrogen with valence 1 2 þ q p e , i.e., psH Ã 0:5þq p =e ). As a consequence, both surfaces of the WZ ZnO(0001) slab become charge-neutral and the slab remains fieldfree. Practical guidelines for determining the pseudo hydrogen charge will be developed in the following section.

Implementation and parameterization of the proposed scheme
A key input quantity to our proposed passivation scheme is the proper valence of the pseudo hydrogen. This value is obviously related to the spontaneous polarization constant of the respective material and there are different ways by which it can be determined. One option is to perform a Berry phase calculation employing the modern theory of polarization 42,43 . This allows to explicitly calculate the spontaneous polarization constant of a specific material and thus the polarization charge q p . For our example of ZnO the calculated spontaneous polarization constant referenced to ZB ZnO is −0.020 e/uc. Since we are interested in studying the (0001) surface, we want to passivate the ð0001Þ side of the slab. The pseudo hydrogen atom will need to have a valence Z psH Ã ¼ 1=2 À 0:020 ¼ 0:480 to passivate the O db; we will denote it as psH Ã 0:480 . A second approach to obtain the spontaneous polarization constant is to use the electric field we observed in the finite-sized slab calculations (cf. Fig. 2). In this approach, we perform a calculation using a slab, in which the partially filled dbs on each of the two slab surfaces are passivated using conventional psH atoms, i.e., Z psH = Z c /4 on the anion-terminated side and Z psH = Z a /4 on the cation-terminated side. The slab in this calculation has to be thinner than the critical thickness at which breakdown occurs. Furthermore, the atoms are fixed at their bulk positions, which is appropriate since we are interested in bulk polarization. Keeping the atoms fixed and thus excluding atomic relaxations allows us to use the electronic part of the dielectric constant (i.e., the highfrequency dielectric constant) when extracting the polarization charge from the observed electric field. Based on this set-up we can directly calculate the spontaneous polarization constant with respect to the ZB phase from the surface-bound charge obtained by using Gauss' law (q s ¼ E Á ε, where q s , E and ε are the surfacebound charge, electric field along the polar direction, and the dielectric constant). We note that according to the surface theorem of the modern theory of polarization 44 the surface-bound charge obtained in this fashion corresponds to the polarization charge when the latter is determined using the ZB phase as a reference. For the example of ZnO(0001) we measure the electric field in a 7 BL slab as E z = 75.51 mV/Å. Using our calculated high-frequency dielectric constant for ZnO ε = 5.25, we obtain for the surfacebound charge q s = −0.021 e/uc. This value is in very good agreement with the spontaneous polarization constant of WZ with respect to the ZB phase obtained from a Berry phase calculation, −0.020 e/uc. Based on the obtained value, we can now use psH Ã 0:480 for passivating the O db on the ZnO(0001) surface, or psH Ã 1:520 for passivating the Zn db on the ZnO(0001) surface.

Applications: surface phase diagram
In the following, we apply the proposed passivation scheme to study surface properties such as surface phase diagrams and surface band structures, which are important in the context of, e.g., catalysis 2,27 , corrosion 40 , or semiconductor devices 45 . We use again the example of a Zn-terminated (0001) surface of WZ ZnO. Typical calculations for ZnO surfaces use slabs of about 4-8 Zn-O BLs thickness, as well as the conventional passivation scheme 2,28,29,46 . The energy convergence analysis (Fig. 1a) shows that within this conventional scheme slabs of 8 BL produce energy differences that are inaccurate by 10 meV/uc; slabs thicker than 23 BL are needed to converge energy differences to within 1 meV/uc. Therefore, we investigate how the convergence behavior affects commonly computed surface properties, such as surface phase diagrams. To perform this benchmark we consider several surface reconstructions of the ZnO(0001) surface 28,40 . These reconstructions have been selected to be energetically highly stable. They also have a small surface unit cell to enable extensive convergence checks with respect to slab thickness. All these structures obey EC. Specifically, we consider the 0.5 monolayer (ML) OH structure, in which OH groups are adsorbed on 1/2 of the surface atoms of the Zn-terminated surface in a regular pattern, the already discussed (2 × 2)-O ad and (2 × 2)-V Zn structures, and a ( ffiffiffiffiffi 12 p ffiffiffiffiffi 12 p )-n3 triangular reconstruction, in which 6 Zn and 3 O atoms are removed from the surface Zn-O BL in a ( ffiffiffiffiffi 12 p ffiffiffiffiffi 12 p ) surface cell. For each of these structures α we evaluate the excess Gibbs free surface energy ΔG α with respect to a reference structure. Here we chose the (2 × 2)-O ad structure as the reference, i.e.
E α tot is the DFT total energy of the surface phase α, μ i is the chemical potential (i.e., the energy of the reservoir) of element i and N i is the number of atoms of element i exchanged with the reservoir. For the construction of the surface phase diagrams, we use energies obtained from calculations for two different slab thicknesses. Specifically, we employ a 5 BL structure, where no breakdown occurs, as well as a thicker 9 BL structure that shows breakdown. For each structure, both the conventional and the proposed passivation schemes are applied. The resulting four surface phase diagrams are shown in Fig. 5. These diagrams are not intended to identify the globally stable surface reconstructions, since we limit our benchmark to only four surface structures rather than studying all known ones (which would be computationally extremely expensive). Still, the constructed diagrams contain sufficient information to reach conclusions regarding convergence with respect to slab thickness.
Each diagram shows only two stable phases over the thermodynamically allowed range of chemical potentials. Interestingly, these two phases are not always the same. The fact that they are not identical implies that slab convergence has not only a quantitative but even a qualitative impact. As shown in Fig. 5 this qualitative change occurs only when applying the conventional passivation scheme (psH 0.5 used to passivate the O dbs). For the five BL slab the Zn vacancy reconstruction is predicted to be thermodynamically stable in the left bottom corner of the phase diagram (cf. Fig. 5c). However, the 9 BL slab shows that this reconstruction is unstable and a triangular structure becomes energetically preferred (Fig. 5d). This qualitative change is a direct consequence of the slow convergence of the conventional scheme with slab thickness. In contrast, using our proposed scheme we observe the correct phase diagram already for small and thus computationally highly efficient slab thicknesses: both the nature of the low-energy phases and the boundary in surface phase diagrams remain the same when the slab thickness is increased (Fig. 5a, b).
We conclude that the proposed scheme allows obtaining converged results at significantly reduced computational costs. For some of the systems studied here the required CPU time could be reduced by more than an order of magnitude.
Applications: electronic structure Finally, we demonstrate that the proposed passivation scheme not only improves surface energetics, but enables the description of accurate and converged surface band structures. Figure 6 shows the surface band structures calculated for the (2 × 2)-O ad 13 BL slab (the same slab as in Fig. 4b), using the conventional (Fig. 6a) and the proposed (Fig. 6b) passivation schemes.
Using the conventional scheme we obtain a surface band structure that is severely distorted because the conduction band on one side of the slab overlaps with occupied surface states near the valence band on the other side of the slab. We note that this behavior cannot be improved by using thicker slabs: due to the breakdown, the bandgap remains zero even in the asymptotic limit of infinitely thick slabs. The proposed scheme removes this artifact and the energy difference between the highest occupied state (the surface state) and the lowest unoccupied state (the CBM) becomes 2.42 eV. We know this value is reliable because of the absence of an electric field in the slab, as seen in Fig. 4b. The surface band structure in Fig. 6b also allows us to determine the energy difference between the CBM and the VBM, which turns out to be 3.36 eV. This is slightly larger than the bulk bandgap (3.27 eV) due to the expected presence of quantum confinement in a finite slab.
We have carefully analyzed the convergence behavior of slab calculations for materials exhibiting spontaneous polarization. Based on this analysis we show that the conventionally used passivated-slab scheme fails: applying it results in a poor convergence with respect to slab thickness and an inaccurate description of surface states. By analyzing the asymptotic bulk limit, we extend the conventional slab passivation scheme to accurately reproduce this limit. The proposed scheme is shown to rapidly converge with slab thickness and to provide accurate and robust results when computing surface energies, surface phase diagrams, and surface states. The proposed formalism can be easily employed in all existing codes since it only requires replacing the passivating pseudo hydrogen atom with valence Z psH = Z c /4 or Z psH = Z a /4 by one with valence Z psH Ã ¼ Z c =4 þ q p =e or Z psH Ã ¼ Z a =4 À q p =e. The proposed scheme will thus allow performing slab calculations for systems exhibiting spontaneous polarization with the same quality and computational resources as for systems where spontaneous polarization is absent.

DFT calculations
The VASP code 33,34 , employing the PAW approach 35 , is used for all presented DFT calculations. The kinetic energy cutoff employed for the plane-wave basis set is 500 eV. A Γ-centered (8 × 8 × 6) k-point grid is used for WZ bulk and a (8 × 8 × 1) grid is used for the (1 × 1) surface unit cell to perform Brillouin-zone integrations. Equivalently folded k-point meshes are used for larger surface cells. In order to obtain accurate electronic structure results, the electronic relaxation is carried out until the total energy convergence is less than 10 −7 eV per surface unit cell. For the exchange-correlation functional we use the generalized gradient approximation PBE 36,37 for analysis related to energetics and the hybrid functional HSE06 38,39 for analysis of the electronic structure. The reason is the decay behavior of the potentials: When employing (semi-)local exchange-correlation functionals the effective electron potential of a surface decays exponentially toward the vacuum. Exact exchange (EXX) reproduces the exact asymptotic behavior proportional to 1/z, with z being the distance of a point in a vacuum from the surface [47][48][49][50] , which results in a slow decay of the potential extending far into the vacuum and hence a significantly slower supercell-size convergence.
As expected, we observe a slow 1/z convergence characteristic for EXX for small and intermediate vacuum sizes. Since the EXX contribution is Fig. 5 Phase diagram constructed using selected ZnO(0001) surface phases. The DFT energies used for the construction are obtained from calculations for two different slab thicknesses and for the two discussed passivation schemes. Different colors indicate different surface phases: 0.5 ML OH adsorbed on a bulk-terminated surface (blue), a Zn vacancy in a (2 × 2) unit cell (green), and a triangular reconstruction ð ffiffiffiffiffi 12 p ffiffiffiffiffi 12 p Þ À n3 in which 6 Zn atoms and 3 O atoms are removed from the first two surface layers (red). The respective surface structures (top view) are shown as insets with O, Zn, and H atoms shown as red, white, and blue balls. smoothly truncated in HSE06 at a distance of 10 Å one might expect at the largest vacuum sizes that the slow converging EXX is absent. The surface energy data shown in Supplementary Fig. 6b clearly shows the expected slow convergence for supercells with a small to medium vacuum region. For the supercells with the largest vacuum size (>30 Å) the slow convergence seems to disappear. However, the much larger scatter in the HSE06 computed surface energies prevents validation of convergence with the accuracy achieved when using semilocal DFT functionals such as PBE ( Supplementary Fig. 6b). Therefore, we use the HSE06 functional for the electronic structure analysis (where having correct bandgaps is important, and the size of the supercells allows achieving electronic convergence) and the PBE functional for the energetics.
With our chosen convergence parameters we obtain for ZnO bulk lattice parameters a = 3.286 Å, c/a = 1.613 and u = 0.379 using PBE and a = 3.248 Å, c/a = 1.609 and u = 0.380 for HSE06 (with α = 0.36). These values show good agreement with experimental 51,52 and theoretical results 2,28,46 . Based on the calculated lattice parameters, the calculated spontaneous polarization constants for bulk ZnO are − 0.041 C/m 2 (equivalent to −0.0203 e/uc) for PBE and − 0.046 C/m 2 (equivalent to −0.0207 e/uc) for HSE06, respectively, which lie within the range of previously reported values (−0.022 to −0.057 C/m 2 ) 53-56 . Our calculated bulk bandgap for ZnO is E g = 0.73 eV with PBE and E g = 3.27 eV with HSE06 (α = 0.36), in good agreement with other theoretical work [E g = 0.74 eV (PBE) 57 , E g = 3.29 eV (HSE06, α = 0.36) 58 ]. The experimental ZnO bandgap is reported to be E g = 3.37 59 to E g = 3.44 eV 60 .

DATA AVAILABILITY
All data used in this work is available in the Pyiron repository (https://www.mpie.de/ 3880353/pyiron) and can be given access to upon request. Fig. 6 Band structure calculated for the WZ ZnO(0001) (2 × 2)-O ad structure. Band structures are obtained using a the conventional and b the proposed passivation schemes for the O-side of the slab. The energy of the highest occupied state is set to 0 eV, and the CB, VB, and surface states are shown as orange, blue and green solid lines, respectively. The energy differences between the highest occupied state (a surface state) and the CBM, and between the VBM and the CBM are indicated.