Screening in 2D: GW calculations for surfaces and thin films using the repeated-slab approach

In the context of photoelectron spectroscopy, the $GW$ approach has developed into the method of choice for computing excitation spectra of weakly correlated bulk systems and their surfaces. To employ the established computational schemes that have been developed for three-dimensional crystals, two-dimensional systems are typically treated in the repeated-slab approach. In this work we critically examine this approach and identify three important aspects for which the treatment of long-range screening in two dimensions differs from the bulk: (1) anisotropy of the macroscopic screening (2) $\mathbf k$-point sampling parallel to the surface (3) periodic repetition and slab-slab interaction. For prototypical semiconductor (silicon) and ionic (NaCl) thin films we quantify the individual contributions of points (1) to (3) and develop robust and efficient correction schemes derived from the classic theory of dielectric screening.

In the context of photoelectron spectroscopy, the GW approach has developed into the method of choice for computing excitation spectra of weakly correlated bulk systems and their surfaces. To employ the established computational schemes that have been developed for three-dimensional crystals, two-dimensional systems are typically treated in the repeated-slab approach. In this work we critically examine this approach and identify three important aspects for which the treatment of long-range screening in two dimensions differs from the bulk: (1) anisotropy of the macroscopic screening (2) k-point sampling parallel to the surface (3) periodic repetition and slab-slab interaction. For prototypical semiconductor (silicon) and ionic (NaCl) thin films we quantify the individual contributions of points (1) to (3) and develop robust and efficient correction schemes derived from the classic theory of dielectric screening.

I. INTRODUCTION
Surface science has developed into a highly active and multidisciplinary area of research that has produced a rich variety of surface-sensitive experimental techniques. In this context theory and calculations have become an invaluable tool for interpreting the often complex and indirect data and for making predictions for structures or properties still inaccessible experimentally.
For many experimentally and technologically relevant questions the electronic structure of a system is of central importance. Experimentally, it can be probed, for example, by photoemission spectroscopy. Theoretically, many-body perturbation theory in the GW approximation, where G is the Green's function and W is the dynamically screened Coulomb interaction, has developed into the method of choice for describing electronic excitations as measured in direct and inverse photoemission in weakly correlated solids and their surfaces. 1,2,3,4,5,6,7 Most GW implementations, notably such that employ pure or augmented plane-waves as basis functions, 8,9,10,11 but also others, 12 rely on three-dimensional periodic boundary conditions, which are appropriate for crystalline bulk systems. To treat systems with a reduced periodicity like surfaces, thin films, nanowires, clusters, or molecules, three-dimensional periodicity is imposed artificially. For this, the system of interest is placed in a three-dimensional unit cell, with an empty (vacuum) region to separate the physical system from its periodic images in the broken-symmetry direction(s); see Fig. 1. This procedure is frequently referred to as "supercell approach," but for systems with a two-dimensional periodicity we prefer the more descriptive term "repeated-slab approach." In this paper we will focus on the two-dimensional case, i.e., the description of surfaces and thin films. Already in the early days of modern GW calculations, Hy- bertsen and Louie applied the then new methodology to simple semiconductor surfaces using the repeated-slab approach. 13 In their paper they raised several fundamental and technical questions. Of importance here is their speculation that "the crucial change in the self-energy operator at the surface may be largely contained in the Green's function," while "the screened interaction may more closely follow the variation in the local density." They further emphasize that their preliminary conclusions "will require further examples . . . as well as a critical evaluation of the slab approach for the surface selfenergy operator". Surprisingly, these issues have only rarely been addressed in later studies and the "critical evaluation of the slab approach" is still missing in the literature. This paper is a contribution to fill this gap, illustrating that the screened interaction does not simply follow the local density but unfortunately is substantially influenced by the repeated-slab geometry.
In principle it is straightforward to investigate how the artificial periodic repetition (in one or more dimensions) affects a given computational method by increasing the vacuum region until the properties of interest exhibit no dependence on this parameter anymore. This is easily achieved in density functional theory (DFT) with the commonly employed local or semilocal functionals, provided that appropriate correction schemes for the low-order electrostatic multipole moments are applied, 14 if necessary. The decoupling becomes more difficult for GW due to the long-range nature of the (screened) Coulomb interaction, as has been explicitly shown for a Na 4 cluster. 15 To solve the decoupling problem, it has been proposed to cut off the Coulomb interaction in the broken-symmetry directions. 15,16,17,18 However, to ensure that no interactions within the physical system are truncated, the vacuum region must be at least as large as the physical system in these schemes, which makes the approach unproportionally expensive (in terms of the computational effort) for thicker slabs. Moreover, the k-point sampling of the Brillouin zone, the computational parameter naturally associated with these longrange effects, becomes a critical convergence parameter along the nontruncated directions. 16,18 As we will show below, even this seemingly technical issue is related to the screening and hence the physical properties of the system. An obvious way to avoid these spurious interactions for systems with broken translational symmetry is to abandon the concept of periodic boundary conditions altogether in the relevant directions and perform the calculations entirely in real space. For semi-infinite jellium surfaces such a GW embedding scheme has been successfully implemented. 19,20 Its extension to realistic surfaces, however, is computationally still too expensive. A real-space implementation for finite systems has also been reported, 21,22 but its applicability to systems with periodicity in one or more directions remains to be shown.
The purpose of this work is to reevaluate the performance of the repeated-slab approach when no Coulombtruncation technique is employed. We address the nature and the magnitude of the effects introduced by the periodic repetition and propose a robust correction scheme that allows to extract the isolated-slab limit already from very small vacuum separations. To achieve this, it proved necessary to explicitly take the anisotropy of the macroscopic screening into account. In addition, we have found that the k-point sampling of the Brillouin zone becomes a more critical parameter for GW repeated-slab calculations than recognized so far. Appropriate samplings sufficient for bulk GW or slab DFT calculations cannot be transferred to surfaces in general.
For weakly correlated bulk systems GW calculations are now routinely performed. Remaining open issues, such as computational efficiency, influence of pseudopotentials or self-consistency, are actively being addressed 3,9,12,21,23,24,25,26,27 but do not affect the conclusions drawn in this paper. For more strongly correlated systems, it becomes necessary to go beyond GW , but these schemes often include the GW self-energy diagrams as lowest order. 28,29 Similarly, the Bethe-Salpeter approach to electron-hole excitations, as probed in optical absorption or electron energy-loss spectroscopy, builds on the GW self-energy. 30 In these schemes, screening plays a similar role as for GW calculations, but they go beyond the scope of the present work.
We note also that the surface electronic structure of a slab may differ from that of a semi-infinite solid. Surface resonances and plasmons, for example, may be affected by confinement effects. The finite thickness of the slab becomes an additional potential source of error in surface calculations, but will not affect free-standing films or heterostructures. Converging the surface electronic structure with respect to slab thickness and/or developing robust correction schemes is a separate problem that would require a detailed study in its own right. We will briefly address the relevant issues in Sec. III.
The remainder of this paper is organized as follows: In Sec. II we present the methodological and physical aspects of the repeated-slab approach. After briefly summarizing the GW space-time method (Sec. II A), we focus first on the anisotropy of macroscopic screening (Sec. II B) and the k-point convergence parallel to the surface (Sec. II C). Then we discuss the influence of the periodic repetition and the associated convergence of the band energies with respect to the vacuum separation between the slabs (Sec. II D). In Sec. III we summarize our main findings and put our conclusions into the context of other computational GW schemes. In the Appendix we present the computational scheme for the classical dielectric models employed.

II. LONG-RANGE SCREENING
A. The GW space-time method All calculations have been performed with the GW space-time method. 9,31,32 Due to its advantageous linear scaling behavior with respect to the k-point sampling of the Brillouin zone it is ideally suited for performing the extensive convergence studies presented below. We have recently extended the code to include the anisotropy in the long-range screening, 33 a crucial point for the repeated-slab approach, as we will demonstrate below. In the following, we will only briefly sketch out the computational scheme and refer the interested reader to the previous papers for further details.
In the space-time method, the Green's function G is constructed in real space (r and r ′ ) and imaginary time (iτ ) from the Kohn-Sham wave functions φ nk and energies ǫ nk (the Fermi level is set as the energy zero).
where Ω denotes the unit-cell volume and the integral over k runs over the first Brillouin zone. Depending on the signs taken, the state summation over n runs over unoccupied (occupied) states for positive (negative) imaginary times. The polarizability is calculated in the random-phase approximation by and is then Fourier transformed to reciprocal space and imaginary frequency (iω). The polarizability P is used to compute the screened interaction W via the symmetrized dielectric matrix (3) From its matrix inverse for each k and iω, the screened interaction is obtained and then Fourier transformed back to real space and imaginary time. The self-energy is then computed as To obtain the self-energy corrections to the Kohn-Sham energies ǫ nk , the matrix elements of the perturbation operator Σ − V xc (where V xc denotes the local exchangecorrelation potential) are computed and Fourier transformed to the imaginary frequency axis. These matrix elements are then analytically continued to the real frequency axis by fitting a multi-pole function on the imaginary frequency axis. 31 Finally, the quasiparticle energies ǫ qp nk are given by the solution of ǫ qp where ϕ nk ||ϕ nk denotes matrix elements with respect to the Kohn-Sham wave functions.

B. Anisotropy of the macroscopic screening
An important point for computing the screened interaction in reciprocal space [Eq. (3)-(4)] is the 1/k 2 singularity of the Coulomb potential as k → 0. For the G = G ′ = 0 element of the symmetrized dielectric matrix [Eq. (3)] it is cancelled by the k 2 behavior of the polarizability. 34 In practice, we enforce this cancellation analytically by performing a Taylor expansion for the polarizability around the Γ point k = 0. However, the expansion depends in general on the spatial direction in which the Γ point is approached. This directional dependence reflects the anisotropy in the macroscopic screening and introduces a nonanalytic, yet finite contribution to the inverse dielectric matrix.
where L(iω) denotes the macroscopic dielectric tensor. The tensor expression in the denominator reduces to a scalar when the macroscopic screening is isotropic. To our knowledge, most GW implementations explicitly exploit this simplification (exceptions being the exact treatment of Ref. 36 and the approximate but accurate 33 tensor treatment of Ref. 37) and therefore implicitly assume that the screening is isotropic. However, this is not the case for repeated-slab systems, as is easily demonstrated.
Assuming that the slabs in Fig. 1 are homogeneous dielectrics of finite thickness s with an isotropic bulk dielectric constant ε b and separated by a vacuum region with thickness v, the dielectric tensor components of the repeated-slab system are given by where c = s+ v denotes the total height of the simulation cell. ε (parallel to the surface) and ε z (perpendicular to it) agree for isotropic systems (s/c = 0 or 1) but deviate considerably for s/c ratios between these limiting cases (cf. Fig. 2). The ratio ε /ε z becomes largest for s/c = 1/2. In practice, the parameters of repeated-slab systems are often close to this ratio of maximum anisotropy.
In the GW calculations the singular elements of W require a special treatment in the subsequent computational steps, which formally involve an integration over the Brillouin zone. Since the singularity is integrable, the result is formally well defined, but cannot be approximated by a finite summation. The solution is to split the screened interaction according to  into a long-range part W lr and a short-range remainder W sr . For W lr , a simple analytic form is chosen that exhibits the same singularity as W , but for which the integral can be computed analytically. 8,10,31,33,36,37,38 The remainder W sr then becomes nonsingular, and its integral can safely be replaced by the sum over a finite k-point grid.
We demonstrate the importance of an anisotropic treatment for W lr for the case of a hydrogen-saturated four-layer Si(100) slab with a vacuum region equivalent to four layers of silicon. The calculated nonzero elements of this repeated-slab system's dielectric tensor are ε xx =5.1, ε yy =5.5, and ε zz =2.2 at the smallest imaginary frequency ω=0.036 hartree, highlighting the large anisotropy predicted by the electrostatic considerations above. In Fig. 3, we show the convergence of the fundamental gap 39 with respect to the k-point sampling N z perpendicular to the surface. It is obvious that isotropic averaging for the screened interaction, i.e., using a scalar rather than a tensorial expression for the singularity of W lr , leads to an unphysical linear increase in the band gap, which will not level off when N z is increased fur- ther. The reason for this linear divergence lies in the inadequate treatment of the singularity in Eq. (7), which is not fully removed. 33 In contrast, the anisotropic treatment converges rapidly.
Only the proper anisotropic treatment allows us to investigate the importance of the k-point sampling in the direction perpendicular to the surface. This convergence has not been discussed for repeated-slab systems before. Using the anisotropic treatment in the GW space-time method and a N × N × N z sampling, we find that the self-energy corrections exhibit a 1/N z behavior, the magnitude of which rapidly decreases with increasing N . 33 Such a behavior might result from the remaining approximations made for the Γ point. In practice, the convergence with respect to N z plays a relevant role only for coarse k samplings. We will address convergence with respect to this parallel sampling in Sec. II C. Test calculations indicate that the qualitative behavior is not affected by the choice of N z . For simplicity, we have therefore used N z = 1 for the calculations reported in the following.

C. k-point sampling parallel to the surface
The k sampling is a further critical aspect in GW calculations for surfaces and thin films, which may not have attracted sufficient attention so far. We demonstrate the importance of this issue for two representative slab systems: a two-layer NaCl(100) slab and a hydrogensaturated four-layer Si(100) slab. In Fig. 4 the gap for a N × N × 1 sampling is plotted as a function of 1/N . In both cases we observe a similar shape of the curves that does not follow a simple power law as observed for N z . We will show in the following that this can be understood in terms of the particular decay behavior of the screened interaction in a slab system and can be modeled using a simple analytic function.
To better understand the approximations introduced by a finite k-point sampling, we will first derive the general connection between the k-point sampling and the range of nonlocality, and then discuss its relevance for the present case. As an introductory remark, we note that in periodic systems the two-point functions F = G, W , or Σ reflect the lattice periodicity where R denotes a vector of the real-space lattice defined by the unit-cell vectors a 1 , a 2 , and a 3 . We can then introduce the representation where r and r ′ are restricted to the unit cell. The corresponding reciprocal-space representation F k (also denoted as mixed-space representation) 31 is obtained from a Fourier transformation We denote the unit vectors of the reciprocal lattice by b i , i ∈ {1, 2, 3}. They are defined by those of the realspace unit cell via a i b j = 2πδ ij . The connection to the R representation provides a real-space picture of the kpoint discretization. For this, we note that a regular, Fig. 5). This supercell coincides with the interaction cell of the space-time method, 9 or in other words the range of nonlocality for the two-point function F = G, W , and Σ. Discretizing k for the function F k is therefore equivalent to imposing a translational symmetry in real space, where S denotes a lattice vector of the interaction-cell lattice (N i a i ). 0000000000000  0000000000000  0000000000000  1111111111111  1111111111111  1111111111111   0000000000000  0000000000000  0000000000000  1111111111111  1111111111111  1111111111111   0000000000000  0000000000000  0000000000000  1111111111111  1111111111111  1111111111111   000000000000000  000000000000000  000000000000000  111111111111111  111111111111111  111111111111111   000000000000000  000000000000000  000000000000000  111111111111111  111111111111111  111111111111111   000000000000000  000000000000000  000000000000000  111111111111111  111111111111111  111111111111111 image charges From this connection we conclude that convergence in the k sampling can only be achieved when the interaction cell associated with the sampling is large enough to encompass the dominant nonlocality in the relevant functions (G, P , W , and Σ). This also holds true for GW schemes that do not explicitly rely on the real-space interaction cell, or which employ a representation other than plane waves for the unit-cell coordinates (r, r ′ ). Let us now consider the behavior of the different two-point functions. In nonmetallic bulk systems, the Green's function decays exponentially with increasing distance of its arguments, typically over a few bond lengths. 40 We see no reason to believe that repeated-slab systems show a qualitatively different behavior in this respect. However, the surface may exhibit a different band gap than the bulk, or it may even be metallic. The decay properties of the Green's function will then change accordingly. 40,41 In particular, a reduction in the gap at the surface implies a slower decay of G at the surface compared to the bulk. The decay properties of G are directly transferred to P [Eq. (2)] and Σ [Eq. (5)]. The screened interaction W , on the other hand, exhibits a very slow 1/r decay. However, this long-range limit W lr is subtracted from the full W in the reciprocal-space singularity treatment discussed previously. For the numerical convergence only the remainder W sr is relevant, i.e., the difference between W and the model function W lr used in the singularity treatment. W sr is dominated by the variation of the electron density and thus by the atomic structure, and should become negligible when the typical length scale of these variations is exceeded. In bulk systems, this corresponds to a few bond lengths. However, the decisive structural variation in a repeated-slab system is the slab itself (at the slab boundary the dielectric constant drops from ε b to 1). The thickness of the slab or the vacuum therefore defines the length scale for the interaction-cell convergence also parallel to the surface.
To understand how the k discretization modifies the screened interaction, we employ the classical theory of dielectric screening. Similar to what was done above for the macroscopic dielectric tensor, we approximate the slab by a homogeneous dielectric with a sharp interface to the vacuum as sketched out in Fig. 6. We then apply the method of image charges (cf. Appendix) to compute the screening function. When a charge is placed at point z ′ , the resulting image charges lie on a line through the original charge perpendicular to the surface (cf. Fig. 6). The effect of the k discretization is simulated by imposing translational symmetry for the image charges in the direction parallel to the surface according to Eq. (14). In other words, we explicitly consider a lattice of perpendicularly aligned image charges [cf. Fig. 6(b)] with the interaction cell's lattice constant of N times the original cell's lattice constant. To simplify the discussion we focus on one characteristic aspect of the screened interaction: the potential V im generated by the image charges at the position of the original charge, i.e., at z = z ′ . We obtain where the sum over S runs over the lattice vectors of the interaction cell N a . d and σ enumerate the image charges that are computed as described in the Appendix. From this the long-range model function 1/(ε r) and the contribution of the original image charges at S = 0 are subtracted. Exploiting the sum rule Eq. (35), we obtain for the error introduced by the parallel discretization In Fig. 7 we show ∆V im (z) for the dielectric model corresponding to the two-layer NaCl(100) slab with a 3× 3 k sampling for both the isolated and the repeatedslab case. The absolute error in the potential is very large. We note that the average shift in the potential is implicitly corrected for in the full GW calculations by setting the W 0 (0, 0) element to zero after the subtraction of the model function. Therefore, only the deviation from this average (dashed line in Fig. 7) contributes to the k-point convergence, which is considerably smaller but far from negligible. We conclude that the convergence of a GW calculation with respect to the k sampling is crucially influenced by the discretization error of the screened interaction.
The variation of ∆V im from its average value may be taken as an indicator for the discretization error in the GW self-energy. We investigated how this changes as a function of slab and vacuum thickness within the dielectric model. When the vacuum region is increased, the variation increases because ∆V im quickly decays in the vacuum region (cf. Fig. 7), thereby pulling down the average over the computational cell. The periodic repetition (small vacuum) thus introduces a highly advantageous compensation effect: While the absolute value of ∆V im is quite large in the slab, the average is dominated by the contribution of the slab, which then cancels out a large part of the total discretization error. The physical origin behind the slab and vacuum dependence of the k convergence lies in the scale dependence of the screening in such an inhomogeneous system. For distances much smaller than the slab thickness, a bulk-like screening is expected inside the slab. However, in cells with large vacuum separations the components of the dielectric tensor may deviate considerably from this bulk limit, as illustrated in Sec. II B. Since the treatment for the long-ranged part of the screened interaction assumes this behavior at all length scales, a portion of this difference in screening enters the numerical treatment of the short-ranged part W sr and affects the k || convergence. If, on the other hand, the vacuum separation is small, the dielectric screening of the repeated-slab system is close to bulklike and the magnitude of the discretization error is reduced correspondingly. We find this anticipated behavior fully confirmed for the slab systems considered in this work ( Figure 10 will show this explicitly for the two-layer NaCl(100) slab).
To estimate the influence of the W discretization on the quasiparticle energies computed in the GW space-time method more precisely, we have developed a fit function to describe the dependence of the band energies on the kpoint sampling N . For this purpose, we assume that the G 0 W 0 corrections show the same functional dependence on N as ∆V im in Eq. (16). Retaining a single term of the summation, we arrive at a three-parameter function: The parameters Q n , D n , and ǫ n (∞) are determined for each state n by fitting to the GW data. We find that this relatively simple form accurately describes the convergence with respect to N for all slab systems studied. Typical examples of the quality of the fit are shown in Fig. 4. We employ the fitting procedure to estimate the remaining error at finite k sampling, or to extrapolate the converged value ǫ n (∞) in cases where the error at finite sampling is unacceptably large.

D. Periodic repetition: slab-slab interaction
Until now, we have focused the discussion on technical and numerical aspects of GW calculations for repeatedslab systems. We now turn to the physical effects associated with the interslab polarization. For this, we computed the band structures for hydrogen-saturated Si(100) slabs with four, six, and eight Si layers and varying amounts of vacuum. The k-point convergence was tested for each system, based on the fitting procedure described above. An 8 × 8 × 1 (7 × 7 × 1) sampling was found to be sufficient for the four-layer (six-and eight-layer) case. The numerical results for the fundamental gap presented in Fig. 8 show a slow convergence with respect to the vacuum size. The changes within the numerically accessible range are noticeable (ranging between 0.05 eV and 0.2 eV per 10 bohr of vacuum). More importantly, we estimate that vacuum regions of 80-100 bohr would be required to achieve absolute convergence within 0.05 eV of the isolated-slab limit (which can be determined with the correction scheme that will be described later). However, only the gap between occupied and unoccupied states is affected. The changes within the occupied or unoccupied part of the spectrum are minor. As we will show, this is a further effect of the long-range screening that affects the occupied and unoccupied bands as a whole: The presence of the neighboring slabs gives rise to an additional contribution in the screened potential, which decreases when the distance between the slabs is increased.
We investigated this within the dielectric model. The model parameters are obtained consistently with the full GW calculation from the dielectric tensor components ε and ε z and the total cell height c. From Eqs. (8) and (9), we find To simulate the effect of the periodic repetition, we take into account a finite number of slabs (≈20) in our model calculation and compare this with the limiting case of an isolated slab, which corresponds to increasing the vacuum size to infinity. The change in the self-energy that results from the long-range screening effects can be estimated from the static Coulomb-hole screened-exchange (COHSEX) approximation 42 : Here W denotes the statically screened interaction and v the bare Coulomb interaction. To isolate the contribution arising from the neighboring slabs we separate the screened interaction into a bulklike part W bulk and an additional contribution W pol from long-range polarization effects. Since Σ is a direct product of W with the Green's function, the self-energy reflects this separation: (23) To simplify this expression further, we split off purely local contributions from W pol according to implicitly defining the purely nonlocal remainder W pol,nl . We now argue that the self-energy Σ pol,nl arising from this separation is negligible for long-range polarization effects. In particular, we have by construction, The corresponding COH contribution therefore vanishes. We now turn to Σ pol,nl SEX . Assuming that W pol is generated by a set of image charges at distance d, W pol,nl (r, r ′ ) ≈ 0 for |r − r ′ | ≪ d. The magnitude of d is determined by the distance to the relevant dielectric interface. On the other hand, the second factor in the SEX term is the Green's function at t → 0, t > 0. Its exponential decay limits the magnitude of |r − r ′ | to a few interatomic distances. 40 Therefore, Σ pol,nl SEX is expected to be small for polarization effects taking place at a length scale of more than a few bohr, and we can safely neglect its contribution. For Σ pol,loc arising from the local part [Eq. (24)], the state summation reduces to a projection onto the occupied states, and the expectation values of Σ pol,loc become ϕ n |Σ pol,loc |ϕ n ± 1 2 ϕ n |W pol (r, r)|ϕ n , where the plus (minus) sign applies to unoccupied (occupied) states. For the dielectric model, we now identify W bulk with 1/(ε(z)|r−r ′ |). W pol (r, r) then reduces to the potential V im (z) induced by a charge at its own position, the S = 0 term in Eq. (15). Equation (27) corresponds to an adiabatic switching on of "image charge" effects for the charged final state, taking into account the sign convention of single-particle energies. We note that Delerue et al. 43 have derived the same final formula from an electrostatic model and used it to estimate self-energy shifts in isolated nanoparticles. In Fig. 9, we compare V im (z) for an isolated slab with that of the repeated slabs. We find for both cases a positive contribution inside the slab, and a negative one outside, in agreement with Delerue. 43 The divergence at the interface is an artifact of the steplike dielectric profile employed. In the repeated-slab case, however, the potential is shifted downward because of the additional polarization in the neighboring slabs. It is precisely this polarization of the neighboring slabs that creates the undesired perturbation in the central slab and is responsible for the observed dependence of the band gap on the vacuum size. The shift in the image potential (indicated by the dashed line in Fig. 9) is essentially constant over the slab and continuous across the slab-vacuum interface. The Σ pol self-energy corrections according to Eq. (27) therefore do not vary significantly between different states, giving rise to a scissor-like change in the gap in agreement with the observations made for the full G 0 W 0 calculation. Moreover, this fortuitously simple situation facilitates a quantitative comparison between the dielectric model and the dependence of the band gap found in the G 0 W 0 calculation. To demonstrate this, we have extracted the magnitude of the finite-vacuum effect from the dielectric model and corrected the numerical G 0 W 0 data by these values for every vacuum thickness. The result has been included in Fig. 8 as the dashed curves. We find that this corrected data no longer depends significantly on the vacuum size. Using this correction scheme, it becomes possible to determine the isolated-slab values for arbitrary vacuum sizes. These a posteriori corrections do not depend sensitively on the microscopic details of the material. We demonstrate this for an ionic material, a two-layer NaCl(100) slab. The dependence of the k convergence behavior on the slab/vacuum ratio (left side of Fig. 10) agrees very well with the predictions from the dielectric model: With increasing vacuum size, the convergence becomes more difficult. For this reason, we found it necessary to extrapolate the k convergence for each vacuum thickness using Eq. (17), because its magnitude becomes comparable to that of the finite-vacuum effect. Using the extrapolated values, we find a very good agreement for the finite-vacuum effect compared to the dielectric model. This is demonstrated by the quasiparticle correction ∆ qp to the band gap shown at the right-hand side of Fig. 10. Including the finite-vacuum correction, ∆ qp coincidentally converges even faster with respect to the vacuum size than the band gap of the underlying DFT-LDA calculation in this case.

III. DISCUSSION AND CONCLUSIONS
In this work we have investigated the performance of the repeated-slab approach in G 0 W 0 calculations. Using the classic theory of dielectric screening we have found that the relevant effects can be understood from simple dielectric slab models. We have demonstrated that long-range polarization effects introduce several important differences in slab systems compared to the cor-responding bulk systems for prototypical semiconductor (Si) and ionic (NaCl) materials. The periodic repetition of the slabs introduces an additional scissorlike change of the one-particle spectrum that can be quantitatively corrected for a posteriori. These corrections are derived from dielectric slab models.
Since long-range effects are described by the small-k behavior in reciprocal space, this limit requires special attention for slab calculations. The anisotropy of the macroscopic screening at k → 0 and the slow convergence with respect to the k grid must be taken into account. If this is not done correctly, it may not be possible to assess the importance of the vacuum separation. However, a careful treatment of these aspects does not pose principal problems and can be achieved with moderate computational effort, in particular for small vacuum sizes. Using the finite-vacuum correction scheme provides a highly efficient technique to extract accurate values for the isolated-slab limit from only a few calculations for small vacuum separations.
We will briefly discuss how far the issues raised in this work for the GW space-time method are relevant also for other computational schemes. The anisotropy of the dielectric tensor is specific to the repeated-slab approach when no Coulomb-truncation techniques are employed. We have demonstrated recently 33 how the anisotropy can be treated to arbitrary precision in plane-wave approaches by expanding the angular dependence in spherical harmonics. In practice, already a maximum angular momentum of 2 is sufficient to capture most of the effect. 33 We expect that a similar strategy is applicable for other basis sets, too.
The slow convergence with respect to the k sampling is rooted in the reciprocal-space computation of the screened interaction. To our knowledge, all GW implementations with periodic boundary conditions rely on this strategy. However, studying the convergence with respect to this parameter is computationally very efficient in the space-time approach since it scales linearly in the number of k points, compared to a typically quadratic scaling in convolution approaches. The discretization error arises from the part that is treated numerically, i.e., the deviation of the screened interaction from the model function used for the 1/k 2 singularity treatment. Alternative integration schemes, such as the improved integration scheme by Pulci et al. 37 or the offset Γ-point method, 10 will probably not improve this significantly, because they aim at better integrating the 1/(k T Lk) model function, not the deviations from it. The k convergence is a critical issue also for isolated slabs when treated by Coulomb-truncation techniques. In fact, it becomes even more dramatic in these cases because the long-range tail then equals the unscreened 1/r Coulomb interaction and the small-vacuum compensation effect is excluded. Indeed, Coulomb-truncation techniques have been found to require drastically enlarged k-point grids in the nontruncated directions. 16,17,18 We note that the anisotropy and the k-point convergence are relevant also for condensed slablike systems, such as multilayers or quantum-well superlattices, but their magnitude depends critically on the dielectric constants of the materials involved.
The dependence of the gap on the vacuum thickness, on the other hand, is a physical property of repeatedslab systems. Since not only the closest neighboring slabs contribute to this effect, its magnitude reduces with increasing slab thickness. 44 It is reproduced with the model calculation and is thus independent of the computational scheme used for the G 0 W 0 calculations. Any implementation that does not find this trend must employ additional approximations that suppress this behavior. Whether such an implementation then provides the isolated-slab limit depends on the actual approximations and cannot be foreseen in general.
For future GW calculations in repeated-slab systems, we recommend the following procedure: 1. Use an anisotropic treatment of the screened Coulomb singularity for systematic k-point convergence studies.
2. Check the k-point sampling parallel to the surface separately for each slab and vacuum thickness. If necessary, extrapolate using a physically justified function.
3. Correct band energies/gaps for the finite-vacuum thickness according to Eq. (27) and the Appendix.Equations (18) and (19) provide model parameters consistent with the actual dielectric tensor.
For surface calculations, it should be kept in mind that the surface electronic structure of a slab differs from that of a semi-infinite solid due to confinement of the electronic states. This aspect is already present at the level of DFT and is directly transferred to the GW band structure. Quantum confinement may for instance split a surface resonance into a series of quantized states. In addition states localized at the two slab surfaces may couple through the slab at insufficient slab thicknesses. GW corrections might become important if they shift the energies of surface states or resonances relative to the bulk, as this alters the decay behavior.
For the surface GW self-energy itself, the role of the slab approximation has, to our knowledge, not been studied. An in-depth analysis of this question goes beyond the scope of the present paper, but we will briefly discuss the relevant aspects for G and W . For semiconductors the exponential spatial decay of the Green's function leads to a rapidly decreasing influence of the second surface with increasing slab thickness. However, if this decrease is fast enough to be unimportant at the commonly employed slab thicknesses of only a few atomic layers cannot be guaranteed a priori. For the screened interaction, on the other hand, the influence of the dielectric discontinuity decays only slowly. Using simple electrostatic models similar to those detailed in Sec. II C (see also Delerue et al. 43 ), we find that the presence of the additional surface at a distance s gives rise to an image-charge-like contribution ≈ ε b −1 2ε b (ε b +1)s . Unlike for the finite-vacuum effect (Sec. II D), however, quantitative corrections for this finite-slab screening effect depend strongly on the spatial extent of the electronic states.
In summary, we find that the repeated-slab approach provides a computationally efficient way to calculate electron and hole energies in the G 0 W 0 approximation for two-dimensional systems like surfaces and thin films. Three manifestations of long-range screening effects in two dimensions -the anisotropy of the macroscopic screening, the k sampling and the periodic repetition of the slabs -have been identified and analyzed in terms of the classic theory of dielectric screening. Their effect on G 0 W 0 calculations in the repeated-slab approach has been demonstrated for thin Si and NaCl films. Robust and efficient correction schemes have been developed. With these, the isolated-slab limit can be easily obtained from repeated-slab calculations.
In order to develop a computational scheme for a multilayer system, a proper book-keeping is crucial to track the reflection and propagation of the various image charges. For this purpose, we divide the system into individual layers. Each layer has the same thickness L and a layerspecific dielectric constant ε(z). 46 For the model calculations in this work, we use L=1 bohr. We will restrict the notation to the z coordinate for the position of the charges and layers. The parallel coordinate ρ becomes only relevant for the computation of the potential; see Fig. 6(a). The origin of the coordinate system is chosen such that the layers are centered around integers L, and the interfaces are at half integers.
For reasons that will become clear below, we denote an image charge that contributes to the potential in layer z and is located at z + σd by q(z, d, σ), where d ≥ 0 is the distance from the layer and σ = ±1. To show that the image charges can be determined iteratively we will now derive the iteration for d → d + L. Consider a charge q(z, d, σ) relevant for the potential in layer z. Due to the interface at z − 1 2 σL two additional image charges appear. The reflected image charge, located at z − σ(d + L), describes the potential in layer z and is given by [cf. (Eq. 30)] q rf (z, d + L, −σ) = ε(z) − ε(z − σ) ε(z) + ε(z − σ) q(z, d, σ) .
The propagated image charge remains at the position z + σd and describes the potential in layer z − σL. Using our book-keeping notation and Eq. (30), it can be written as q pr (z − σ, d + L, σ) = 2ε(z − σ) ε(z) + ε(z − σ) q(z, d, σ) . (32) Obviously, the distance parameter is increased by L for each interface taken into account. The image charges for the distance d + L can thus be computed iteratively from those at distance d and will in general combine a reflected and a propagated contribution q = q rf + q pr . The iterations are started by setting q(z ′ , 0, ±1) = 1 .
We restrict z ′ to integer L, i.e., the original charge is placed at the center of a layer, thereby avoiding the divergence of the image potential at the dielectric discontinuities between adjacent layers. In practice, the iterations are stopped at some d max , which thereby becomes a convergence parameter. In addition, we truncate the number of layers considered explicitly and neglect image charges that fall outside. This truncation becomes a second convergence parameter. The convergence for both parameters was tested by doubling the parameter until the changes became negligible. Summing the Coulomb potential of all the image charges relevant for this layer, one obtains The first term in Eq. (34) describes the potential of the original unit charge for z = z ′ and corresponds to that of a bulk material with the local dielectric constant ε(z).
The sum over image charges q(z, d, σ) then introduces the long-range polarization effects due to the variations in the dielectric constant.
As a final remark, we mention a useful connection to the macroscopic anisotropy. The long-range limit in the direction parallel to the surface 1/(ε ρ) is obtained for ρ ≫ d by neglecting the perpendicular distance d in the denominator of Eq. (34). From this, we obtain the sum rule