Efficient many-body calculations of 2D materials using exact limits for the screened potential: Band gaps of MoS 2 , hBN, and phosphorene
Filip A. Rasmussen, Per S. Schmidt, Kirsten T. Winther, Kristian S. Thygesen
EEfficient many-body calculations of 2D materials using exact limits for the screenedpotential: Band gaps of MoS , hBN, and phosphorene Filip A. Rasmussen, Per S. Schmidt, Kirsten T. Winther, and Kristian S. Thygesen
Center for Atomic-scale Materials Design (CAMD),Department of Physics and Center for Nanostructured Graphene (CNG),Technical University of Denmark, DK-2800 Kongens Lyngby, Denmark (Dated: June 21, 2016)Calculating the quasiparticle (QP) band structure of two-dimensional (2D) materials within theGW self-energy approximation has proven to be a rather demanding computational task. Themain reason is the strong q -dependence of the 2D dielectric function around q = that calls for amuch denser sampling of the Brillouin zone than is necessary for similar 3D solids. Here we use ananalytical expression for the small q -limit of the 2D response function to perform the BZ integral overthe critical region around q = . This drastically reduces the requirements on the q -point meshand implies a significant computational speed-up. For example, in the case of monolayer MoS ,convergence of the G W band gap to within ∼ . × q -points rather thanthe 36 ×
36 mesh required with discrete BZ sampling techniques. We perform a critical assessment ofthe band gap of the three prototypical 2D semiconductors MoS , hBN, and phosphorene includingthe effect of self-consistency at the GW and GW level. The method is implemented in the opensource GPAW code. PACS numbers: 71.15.Dx, 71.15.Qe, 73.22.-f
I. INTRODUCTION
The past few years have witnessed an explosion in re-search on atomically thin two-dimensional (2D) mate-rials. Of particular interest are the 2D semiconductorsincluding the family of transition metal dichalcogenides,which have been found to exhibit a number of uniqueopto-electronic properties.
For understanding and pre-dicting these properties the electronic band structure ofthe material is of fundamental importance. The GWmethod, introduced by Hedin in 1965 and first ap-plied to real solids in an ab-initio framework by Hybert-sen and Louie and Godby, Sham, and Schl¨uter, hasbecome the “gold standard” for calculating quasi-particle(QP) band structures. Over the years its performancehas been thoroughly established for bulk materials and more recently also for molecules. In comparison,critical assessments of the accuracy and numerical con-vergence of GW calculations for 2D materials are ratherscarce.
Nevertheless these studies have shown that(i) it is essential to use a truncated Coulomb interac-tion to avoid long range screening between periodicallyrepeated layers which reduces the QP band gap, and (ii)when a truncated Coulomb interaction is used, the con-vergence of the GW calculation with respect to the num-ber of k -points becomes much slower than is the case forsimilar bulk systems.The slow k -point convergence of the GW band struc-ture is directly related to the nature of electronic screen-ing in 2D which is qualitatively different from the wellknown 3D case. Specifically, while the dielectric func-tion, ε ( q ), of bulk semiconductors is approximately con-stant for small wave vectors, the dielectric function ofa 2D semiconductor varies sharply as q → . As aconsequence, the number of q -points required to obtain a proper sampling of the screened interaction W ( q ) overthe Brillouin zone (BZ) is much higher for the 2D ma-terial than what would be anticipated from the 3D case.For example, the band gap of bulk MoS is converged towithin ∼ . k -point grid of 12 × requires agrid of 36 ×
36 when standard BZ sampling schemes areapplied.Importantly, supercell calculations only describe thecharacteristic features of screening in 2D materials cor-rectly when the unphysical screening between the peri-odically repeated layers is removed, e.g. using a trun-cated Coulomb interaction. Without this, the dielec-tric function does not approach unity for q → . with 10 ˚A vacuum ) as a re-sult of over screening. Since in this case, the screeningis really intermediate between 3D and 2D, the GW gapshows faster convergence with k -points than is observedusing a truncated Coulomb interaction. This is presum-ably the reason why most earlier GW calculations for 2Dsemi-conductors have been performed with k -point gridsranging from 6 × ×
15 which are much too coarsefor describing the truly isolated 2D material.Here we show that the slow k -point convergence of theGW self-energy in 2D materials can be alleviated by per-forming the BZ integral of W ( q ) analytically in the criti-cal region around q = where ε ( q ) varies most strongly.The analytical expression for W ( q ) is obtained from alowest order expansion in q of the head, χ ( q ), andwings, χ ( q ), of the non-interacting density responsefunction. This simple scheme reduces the required num-ber of q -points by an order of magnitude without loss ofaccuracy. a r X i v : . [ c ond - m a t . m t r l - s c i ] J un II. THE GW SELF-ENERGY
We split the GW self-energy into the exchange and cor-relation part, respectively. The former does not presentparticular problems in 2D materials and is calculated us-ing a Wigner-Seitz truncated Coulomb interaction as de-scribed elsewhere. In a plane wave expansion the corre-lation part of the self-energy, evaluated for an electronicstate | n k (cid:105) , takes the general form (cid:104) n k | Σ c ( ω ) | n k (cid:105) = 1(2 π ) (cid:90) BZ d q (cid:88) GG (cid:48) i π (cid:90) ∞−∞ dω (cid:48) W GG (cid:48) ( q , ω (cid:48) ) × (cid:88) m [ ρ m, k + q n, k ( G )][ ρ m, k + q n k ( G (cid:48) )] ∗ ω + ω (cid:48) − (cid:15) m k + q − iη sgn( (cid:15) m k + q − µ ) , (1)where m runs over all electronic bands, (cid:15) m k + q are singleparticle energies, µ is the chemical potential and η is abroadening parameter. The pair densities are defined as ρ m k + q n k ( G ) = (cid:104) n k | e i ( G + q ) · r | m k + q (cid:105) , and W GG (cid:48) ( q , ω ) isthe correlation part of the dynamical screened potentialgiven by W GG (cid:48) ( q , ω ) = v G ( q ) (cid:2) ε − GG (cid:48) ( q , ω ) − δ GG (cid:48) (cid:3) , (2)where v G ( q ) = 4 π/ | G + q | is the Coulomb interaction.In most implementations the BZ integral is evaluatednumerically with a standard quadrature method usinga regular q -point grid matching the k -point grid of theground state DFT calculation. Since the screened po-tential, Eq. (2), diverges for G = 0 , q = (for bulkmaterials) this point is generally handled separately, sothe integral may be written (cid:90) BZ d q S ( q , ω ) → Ω N q (cid:88) q n (cid:54) = S ( q n , ω ) + (cid:90) Ω d q S ( q , ω ) , (3)where S ( q , ω ) denotes the entire integrand, Ω is the vol-ume of the BZ, N q is the total number of q -points in thegrid and Ω denotes a small region around q = . Forbulk systems Ω is normally defined as a sphere centeredat q = . For a 2D material, the BZ integral in Eq. (1)reduces to a 2D integral with in-plane sampling of q , andΩ represents an area. We now focus on how to calculatethe contribution to the integral around the special point q = in the 3D versus the 2D case.Within the random phase approximation (RPA) thedielectric matrix is given by ε GG (cid:48) ( q , ω ) = δ GG (cid:48) − v G ( q ) χ GG (cid:48) ( q , ω ) . (4)For a solid with a finite band gap it can be shownthat the head of the non-interacting response function χ ( q , ω ) ∝ q for small q . Since v ( q ) = 4 π/q itfollows that in 3D the head of the dielectric function ε ( q , ω ) converges to a finite value > q → ε ( q , ω ) in a relatively large region around q = . This means that in the BZ integration in Eq. (1) aroundthe singular point G = G (cid:48) = q = all factors, except4 π/q , can be assumed to be constant and the integralcan be performed analytically over a sphere centred at q = or numerically on a fine real space grid. For GW calculations of 2D materials performed withperiodic boundary conditions in the out-of-plane direc-tion, the direct use of Eq. (1) leads to significant over-screening due to the interaction between the repeatedimages. One way of dealing with this is to subtract theartificial interlayer contribution calculated from a classi-cal electrostatic model. A more direct way of avoidingthe spurious interactions is to truncate the Coulomb po-tential in the direction perpendicular to the layers. Thusin Eqs. (4) and (2), v G ( q ) should be replaced by v G ( q (cid:107) ) = 4 π | q (cid:107) + G | (cid:104) − e −| q (cid:107) + G (cid:107) | L/ cos( | G z | L/ (cid:105) , (5)where only in-plane q are considered. L is the length ofthe unit cell in the non-periodic direction and the trun-cation distance is set to L/
2, which simplifies the expres-sion. In the long wavelength limit, G = , q (cid:107) → ,the truncated interaction becomes v ( q (cid:107) ) ≈ πLq (cid:107) . Wesee that the q = divergence in the truncated Coulombpotential is reduced by a power of q compared to thatof the full Coulomb interaction. As will be shown in thefollowing this has important consequences for the form ofthe screened interaction. However, before presenting theform of the screened interaction of a 2D semiconductorevaluated using the full expression for the response func-tion and truncated Coulomb interaction, it is instructiveto consider a simplified description of the 2D material.Let us consider a strict 2D model of a homogeneousand isotropic semiconductor. In the small q limit, thedensity response function takes the form χ ( q ) = α q where α is the 2D polarizability. Using that the 2DFourier transform of 1 /r equals 2 π/q , the leading orderof the dielectric function becomes ε D ( q ) ≈ πα q. (6)Some examples of macroscopic dielectric functions fora representative set of 2D semiconductors are shown inFigure 1 (see Ref. 21 for a precise definition of thisquantity). The linear form (6) is clearly observed inthe small q regime. Importantly, if we use the samestrategy for evaluating the BZ integral in Eq. (1) as in3D, i.e. assuming (cid:15) − ( q ) to be a slowly varying func-tion and evaluating it on the discrete q -point grid, weobtain zero contribution for the q = 0 term because1 /(cid:15) D − q = 0, see Eq. (2). On the otherhand, it is clear that the screened interaction takes theform W D ( q ) = − π α / (1 + 2 πα q ) for small q . Inparticular, W D ( q ) takes a finite value for q = 0 whichis qualitatively different from the 3D case where W ( q )diverges for q → .In Appendix A we show, following an analysis similarto that of Ref. 32 adapted to the case of a truncated . . . . . . . . . q [˚A − ]0246810 ε D ( q , ω = ) h-BNPhosphorene H-MoS ( M )T-HfO ( M )T-NiS ( M )h-BN ( M )Phosphorene ( X )Phosphorene ( Y ) FIG. 1. (Color online) Static macroscopic dielectric functionsof a representative set of 2D semiconductors as a function of q along the Γ → M direction for the hexagonal structures andalong the path from Γ to X or Y in the case of phosphorene. Coulomb interaction, that for a general non-isotropic 2Dmaterial, the small q (cid:107) limit of the head of the screenedpotential takes the form W ( q (cid:107) ) = − (cid:18) π (1 − e −| q (cid:107) | L/ ) | q (cid:107) | (cid:19) × ˆ q (cid:107) · A ˆ q (cid:107) π (1 − e −| q (cid:107) | L/ ) ˆ q (cid:107) · A ˆ q (cid:107) , (7)where ˆ q (cid:107) = q (cid:107) / | q (cid:107) | and A is a second rank tensor whichalso depends on the frequency. We see that we have W ( q (cid:107) = ) = − (2 πL ) ˆ q (cid:107) · A ˆ q (cid:107) . The expression ˆ q (cid:107) · A ˆ q (cid:107) is closely related to the slope of the dielectric function andthe 2D polarizability but includes local field effects. Inaddition to Eq. (7) there are similar expressions for thewings and body of the screened interaction, see Eq.(A25)to Eq.(A28). These expressions must be integrated overthe Ω -region, that we now define as the primitive cellin the 2D BZ that surrounds the q (cid:107) = point. Theexpression is simplified to one that can be integrated an-alytically as shown in Appendix A.The full expression for W in Eq. (7) is therefore eval-uated numerically on a discrete sub-grid, constructed asa Monkhorst-Pack grid within Ω , and the simplified ex-pression in Eq. (A31) is only used for q (cid:107) = on thesub-grid. The limit of the integral is now given by theradius r Ω , defined as πr = Ω /N q , where N q is thenumber of grid points in the sub-grid. This approach en-sures a smooth evaluation of W , that converges fast with N q . It is found to be necessary to have N q ≈
100 when q (cid:107) = is evaluated using Eq. (A31) for both iso- andanisotropic materials where as N q ≈ is needed if theanalytical correction at q (cid:107) = is omitted. FIG. 2. (Color online) The head of the static componentof the screened potential (subtracted the bare interaction) ofmonolayer a) h-BN b) MoS and c) phosphorene as a functionof q along the Γ → M direction or Γ → X and Γ → Y in the case of phosphorene. The crosses are the numericalvalues obtained on a fine q -point grid while the circles ortriangles represent the values obtained on a coarse q -pointgrid. The bars represent a simple numerical approximationto the BZ integral of W ( q ) performed on the coarse q -pointgrid. The value of the screened potential for q = is set tothe analytical result Eq. (7). The full curve represents theanalytical small q approximation, Eq. (7), and the hatchedarea shows its contribution to the integral. A. Results
To investigate how this method performs we have car-ried out test calculations for the three monolayers h-BN,MoS and phosphorene, which have quite different dielec-tric functions as seen on Figure 1. h-BN is a large gapdielectric with low screening ability leading to a smallslope of the dielectric function at q = , while MoS hasa larger dielectric function and quite steep slope at q = .Phosphorene has a dielectric function similar to MoS insize and steepness but is anisotropic with slopes varyingby ∼
40% between the two high symmetry directions,Γ → X and Γ → Y .All the calculations were performed using the GPAWelectronic structure code. The structures used inthe present calculations are relaxed with DFT using thePBE exchange-correlation (xc) functional. The result-ing lattice constant for h-BN is 2 .
504 ˚A, the in-plane lat-tice constant for MoS is 3 .
184 ˚A with a S-S distanceof 3 .
127 ˚A. For phosphorene the in-plane unit cell is4 .
630 ˚A by 3 .
306 ˚A, the in-plane P-P-P angle is 95 . ◦ and the layer thickness is 2 .
110 ˚A. A convergence testwith respect to the amount of vacuum between repeat-ing periodic images was carried out and 15 ˚A of vacuumwas necessary for h-BN and Phosphorene where as only10 ˚A was needed for MoS . The PBE eigenvalues andwavefunctions were calculated with a plane wave basiscut-off energy of 600 eV and used as input in the GWcalculations. For the initial investigation of the q -pointconvergence, the dielectric function and the correlationself-energy were calculated using a cutoff of 50 eV. Thiscutoff is insufficient to ensure properly converged quasi-particle energies, but it is adequate to describe the trendsrelated to the improved q -point sampling relevant for thisstudy. The following fully converged calculations werecarried out using a 1 /N pw extrapolation to the completebasis set limit using cutoff energies of up to 200 eV. In Figure 2 we compare the analytical small q ex-pression, Eq. (7), for the head of the screened potential W ( q ) with the numerical values obtained using a fineand coarse q -point sampling. In all the cases the q = value has been set to the analytical value. It is evidentthat the screened potential falls off quickly and thus fora coarse q -point sampling the q = contribution to theintegral is by far the largest and should therefore not beneglected. Similarly, using only the exact value in q = could also pose a problem as the contribution will begrossly overestimated due to the convex nature of poten-tial. We note that the analytical expression follows thenumerical results quite closely and is even accurate faraway from the Γ-point – for MoS we have an almostperfect agreement for the points shown. Thus using theanalytical limit within the region around q = is rea-sonable. We notice that the anisotropy of phosphorenemakes W ( q ) ill-defined at q = (different limit valuesdepending on the direction of q ). For larger q the dielec-tric anisotropy becomes negligible. However, because ofthe relatively large weight of the q = contribution tothe BZ integral, the anisotropy should be taken into ac-count for accurate GW calculations.We note that a similar approach to the treatment ofthe q = term of the screened potential was suggestedin Ref. 20. That particular method was based on fittingto an empirical expression for ε ( q ) calculated from thevalue at a small but finite q . The method outlined hereis different in that the analytical expression for W ( q ) is obtained from a lowest order expansion of the head, χ ( q ), and wings, χ ( q ), of the non-interacting den-sity response function and thus can be obtained with-out fitting or using empirical parameters. This also en-sures that the effect of in-plane dielectric anisotropy isexplicitly included. (a) h-BN × × × × × × /N k Q P g a p ( e V ) (b) 2H-MoS × × × × × × /N k N k (Number of k -points) (c) Phosphorene × × × × × × /N k No q =0 term Analytical q → limit FIG. 3. (Color online) The G W quasi-particle band gapof monolayer (a) 2H-MoS (b) h-BN and (c) phosphorene,calculated using two different treatments of the q = termin Eq. (1). The dashed (green) line shows the contributionobtained when the head and wing elements of the q = termare neglected corresponding to the standard treatment usedfor 3D systems. The solid (blue) line shows the contributionobtained when using the analytical results, Eq. (7), to performthe integral over the q = element. The insets shows theresults for the largest k -point grids on a reversed linear scalein 1 /N k . Notice the zero point is at the right side of the x -axis. In Fig. 3 we show the minimum QP band gap of mono-layer h-BN, MoS and phosphorene as a function of 1 /N k where N k is the total number of k -points in the BZ sam-pling (the q point grid for the GW integration is the sameas the k -point grid used in DFT). We compare the resultsobtained using two methods: i) neglecting the q = con-tribution to head and wings of the screened potential andii) evaluating Eq. (7) as described. It is clear that methodi) in all cases underestimates the correlation self-energydue to the underestimation of the screening; In order toget the band gap converged to within ∼ . k -point sampling of minimum 36 ×
36 forh-BN, 36 ×
36 for MoS and 22 ×
30 for phosphorene. Wealso note that for large k -point grids the band gaps us-ing this method converge approximately as 1 /N k as themissing contribution is almost proportional to the areaof the q = region. Clearly, the latter approach variessignificantly less with the k -point grid and in fact thegap is converged to within 0 . k -pointgrid in the order of 6 × ∼ . ×
12 grid (in the worst case). We have performed testcalculations for other 2D semiconductors and obtainedsimilar conclusions although the number of k -points re-quired to reach convergence within 0 . q = term neglected) is some-what system dependent; materials with efficient screen-ing, e.g. MoS and NiS , require larger k -point gridsthan materials with poor screening, e.g. h-BN and HfO (see Fig. 1). N k (Number of k-points) . . . . . . . . . Q P g a p ( e V ) × × × × × × Vacuum = 10 ˚AVacuum = 15 ˚AVacuum = 20 ˚A No q = term Analytical q → limit FIG. 4. (Color online) The band gap of monolayer 2H-MoS calculated with different amounts of vacuum betweenrepeated layers. The solid and dashed lines are with andwithout the q = 0 correction, respectively. As the vacuum isincreased, the weight of the correction is decreased and it isnecessary to use denser in-plane k -point sampling to achieveconvergence. To obtain converged band gaps it is necessary to use aunit cell with enough vacuum between repeated layers toavoid an artificial interaction. This is true even when atruncated Coulomb interaction is used as the finite vac-uum affects wave functions and energies, in particular forhigher lying unbound states. As the amount of vacuumis increased, the Brillouin zone shrinks and the analyti-cal correction around q = 0, applied only for G = 0, hassmaller weight. This means a slower convergence withrespect to in-plane k -points. This is shown in Fig. 4 for MoS , where it is clear that the correction is less effec-tive for larger vacuum. The calculations converge towardthe same value indicating that for MoS
10 ˚A of vacuumis sufficient. The most efficient procedure to obtain con-verged band gaps is therefore to first converge the amountof vacuum at a low k -point sampling without applyingthe correction and then afterwards converge the k -pointsampling with the correction at the given vacuum.In Table I we report the converged values for the quasi-particle band gaps. For h-BN the band gap is indirectbetween the K - and Γ-point, for MoS and phosphoreneit is direct at the K - and Γ-point respectively. For thesecalculations we used 18 × k -points for h-BN, 18 × k -points for MoS and 10 ×
14 for phosphorene with theanalytical integration of W ( q ) around q = . Accord-ing to Fig. 3 this is sufficient to ensure convergence towithin 0 .
05 eV. We note that spin-orbit interactions arenot included in the reported values. Inclusion of spin-orbit interactions split the valence band of MoS at the K point by 0 .
15 eV thereby lowering the QP gap by around0 .
07 eV.
Spin-orbit interactions have no effect for h-BN and phosphorene.For MoS the converged G W @PBE band gap of2 .
54 eV agrees well with our previously reported valueof 2 .
48 eV (with spin-orbit coupling) obtained using aWigner-Seitz truncated Coulomb interaction and 30 × k -points. Other reported gaps range from 2 .
40 eV to2 .
82 eV.
However, these calculations were performedi) without the use of a truncated Coulomb interactionand including 15-25 ˚A vacuum, ii) employing relativelysmall k -point grids of 6x6 to 18x18, and iii) using differ-ent in-plane lattice constants varying between 3.15 and3.19 ˚A. These different settings can affect the band gapby as much as 0.5 eV, and therefore we refrain fromproviding detailed comparison of our result to these ear-lier calculations. An overview of previous GW results forMoS can be found in Ref. 35.In Ref. 22 a G W @LDA band gap for MoS of 2 .
70 eVis reported using a truncated Coulomb interaction and acalculation of the screened potential at q = based onthe method in Ref. 20. In that study, the lattice con-stant of MoS was 3 .
15 ˚A. With this lattice constant weobtain a gap of 2 .
64 eV, which is in fair agreement withRef. 22. Our result is very close to the experimental valueof 2 . Per-forming partially self-consistent GW the band gap in-creases to 2 .
65 eV (2 .
58 eV including spin-orbit coupling).For h-BN, we obtained a G W band gap of 7 .
06 eVwhich increases to 7 .
49 eV with GW . In Ref. 49 theG W band gap was calculated to be 7 .
40 eV. Insteadof a truncated Coulomb interaction the band gap wasextrapolated to infinite vacuum. The treatment of the q = 0 term is not mentioned nor is the size of the k -point grid. Despite the difference at the G W level,they report a similar increase of the band gap of 0 . calculation.For phosphorene we calculate a G W band gap of2 .
03 eV which agrees well with the previously reported
Transition DFT-PBE G W @PBE GW @PBEh-BN K → Γ 4.64 7.06 7.49 K → K K → K → Γ 0.90 2.03 2.29TABLE I. Band gaps in eV calculated with DFT-PBE,G W @PBE and GW @PBE using the PBE-relaxed struc-tures. The GW calculations were performed using analyticintegration of W ( q ) around q = without including spin-orbit interactions. 10 ˚A of vacuum was used for MoS and15 ˚A for h-BN and phosphorene. The following k -point gridswere used; h-BN: 18 ×
18, 2H-MoS2: 18 ×
18 and phosphorene:10 × value of 2 . using the method of Ref. 20. The bandgap increases to 2 .
29 eV with GW . III. CONCLUSION
In conclusion, we have discussed the connection be-tween the form of the q -dependent dielectric function ofa 2D semiconductor and the slow k -point convergence ofthe GW band structure. We have derived an analyticalexpression for the q → limit of the screened potential ofa semiconductor when a 2D truncation of the Coulombpotential is used. The method accounts for dielectricanisotropy and does not rely on any additional param-eters or fitting. Using this expression we have shownthat convergence of the GW self-energy with respect tothe size of the k -point grid is drastically improved. Forthe specific case of monolayer MoS , we found that theuse of the analytical form alone reduces the k -point gridrequired to achieve convergence of the GW self-energycontribution to the band gap to within ∼ . ×
36 to 12 ×
12 – a reduction in the numberof k -points by an order of magnitude. This method maytherefore enable future large-scale GW calculations for2D materials without compromising accuracy. ACKNOWLEDGMENTS
We acknowledge support from the Danish Council forIndependent Research’s Sapere Aude Program, GrantNo. 11-1051390. The Center for NanostructuredGraphene is sponsored by the Danish National ResearchFoundation, Project DNRF58.
Appendix A: Calculation of the q → In the following we derive the analytical form of thescreened potential, Eq. (2), for 2D materials in the limit q (cid:107) → . We largely follow the approach of Ref. 32 where the same limit for bulk systems was considered. As ex-plained in the main text we use a truncated Coulombinteraction of the form v ( r (cid:107) , z ) = θ ( R − | z | ) (cid:112) | r (cid:107) | + z . (A1)Using this potential we effectively turn off interactionbetween electrons on different 2D layers of the supercellcalculation. We choose R to be half the height of theunitcell, R = L/
2, so that an electron in the center ofthe layer will not interact with electrons located in theneighboring unitcell. This means that the 2D truncatedcoulomb interaction of Eq. (16) in reduces to v G ( q (cid:107) ) = 4 π | q (cid:107) + G | (cid:104) − e −| q (cid:107) + G (cid:107) | L/ cos( | G z | L/ (cid:105) , (A2)where only in-plane q is considered. We note that in thelimit L → ∞ it takes the usual 3D form, v G ( q ) = π | q + G | .In the long wavelength limit it has the asymptotic behav-ior v ( q z = 0 , q (cid:107) → ) = 2 πL | q (cid:107) | , (A3)diverging slower than the full Coulomb potential withprofound consequences for the properties of 2D materials.In the long wavelength limit q → the non-interactingdensity response function or irreducible polarizability hasthe following behavior χ (cid:48) ( q → ) = q · P q = | q | ˆ q · P ˆ q (A4) χ G0 ( q → ) = q · p G = | q | ˆ q · p G , , (A5) χ ( q → ) = q · s G = | q | ˆ q · s G , (A6)where P is a second rank tensor, p G and s G are propervectors and ˆ q = q / | q | . The density response function,and therefore also P and p G , has a frequency dependencewhich here and through the rest of this section has beenleft out to simplify the notation. Within the randomphase approximation the dielectric function is given by(schematically) ε = 1 − vχ . (A7)Due to technical reasons it is easier to work with asimilar symmetrized version given in Fourier space by˜ ε GG (cid:48) ( q ) = δ GG (cid:48) − (cid:112) v G ( q ) χ GG (cid:48) ( q ) (cid:112) v G (cid:48) ( q ) . (A8)Inserting the Coulomb potential, Eq. (A2), and the ex-pressions for the non-interacting response function Eqs.(A4)-(A6), the head and wings of the symmetrized di-electric function are˜ ε ( q (cid:107) → ) = 1 − v ( q (cid:107) ) | q (cid:107) | ˆ q (cid:107) · P ˆ q (cid:107) = 1 − π (1 − e −| q (cid:107) | L/ ) ˆ q (cid:107) · P ˆ q (cid:107) (A9)˜ ε G0 ( q (cid:107) → ) = − (cid:113) v G ( ) (cid:113) v ( q (cid:107) ) ˆ q (cid:107) · p G = − (cid:113) v G ( ) (cid:113) π (1 − e −| q (cid:107) | L/ ) ˆ q (cid:107) · p G (A10)˜ ε ( q (cid:107) → ) = − (cid:113) v G ( ) (cid:113) v ( q (cid:107) ) ˆ q (cid:107) · s G = − (cid:113) v G ( ) (cid:113) π (1 − e −| q (cid:107) | L/ ) ˆ q (cid:107) · s G . (A11) To determine the inverse dielectric function we writethe dielectric function as a block matrix in the G , G (cid:48) components with head, wings and body of the form˜ ε = (cid:18) H w (cid:124) v B (cid:19) (A12)The inverse is then given by˜ ε − = (cid:18) ( H − w (cid:124) B − v ) − − ( H − w (cid:124) B − v ) − w (cid:124) B − − B − v ( H − w (cid:124) B − v ) − B − + B − v ( H − w (cid:124) B − v ) − w (cid:124) B − (cid:19) (A13)From this we see that ˜ ε − = ˜ ε − (cid:88) G , G (cid:48) (cid:54) = ˜ ε B − GG (cid:48) ˜ ε G (cid:48) − (A14)˜ ε − G0 = − ˜ ε − (cid:88) G (cid:48) (cid:54) = B − GG (cid:48) ˜ ε G (cid:48) (A15)˜ ε − = − ˜ ε − (cid:88) G (cid:48) (cid:54) = ˜ ε (cid:48) B − G (cid:48) G (A16)˜ ε − GG (cid:48) = B − GG (cid:48) + ˜ ε − (cid:88) G (cid:48)(cid:48) (cid:54) = B − GG (cid:48)(cid:48) ˜ ε G (cid:48)(cid:48) (cid:88) G (cid:48)(cid:48) (cid:54) = ˜ ε (cid:48)(cid:48) B − G (cid:48)(cid:48) G (cid:48) (A17)Introducing the vectors a G , b G and the tensor A givenby a G = − (cid:88) G (cid:48) (cid:54) = B − GG (cid:48) (cid:113) v G (cid:48) ( ) p G (cid:48) (A18) b G = − (cid:88) G (cid:48) (cid:54) = (cid:113) v G (cid:48) ( ) s G (cid:48) B − G (cid:48) G (A19) A = − P − (cid:88) G (cid:54) = (cid:113) v G ( q (cid:107) ) s G ⊗ a G , (A20)where ⊗ denotes the tensor product, the long wavelengthlimit of the inverse dielectric function is seen to be given by ˜ ε − ( q (cid:107) → ) = 11 + 4 π (1 − e −| q (cid:107) | L/ ) ˆ q (cid:107) · A ˆ q (cid:107) (A21)˜ ε − G0 ( q (cid:107) → ) = − (cid:113) π (1 − e −| q (cid:107) | L/ ) ˆ q (cid:107) · a G π (1 − e −| q (cid:107) | L/ ) ˆ q (cid:107) · A ˆ q (cid:107) (A22)˜ ε − ( q (cid:107) → ) = − (cid:113) π (1 − e −| q (cid:107) | L/ ) ˆ q (cid:107) · b G π (1 − e −| q (cid:107) | L/ ) ˆ q (cid:107) · A ˆ q (cid:107) (A23)˜ ε − GG (cid:48) ( q (cid:107) → ) = B − GG (cid:48) + 4 π (1 − e −| q (cid:107) | L/ )( ˆ q (cid:107) · a G )( ˆ q (cid:107) · b G (cid:48) )1 + 4 π (1 − e −| q (cid:107) | L/ ) ˆ q (cid:107) · A ˆ q (cid:107) (A24)Inserting these expression in the equation for thescreened potential, Eq. (2), we see that the head andwings are given by W ( q (cid:107) → ) = v ( q (cid:107) ) (cid:2) ˜ ε − ( q (cid:107) ) − (cid:3) = − (cid:18) π (1 − e −| q (cid:107) | L/ ) | q (cid:107) | (cid:19) ˆ q (cid:107) · A ˆ q (cid:107) π (1 − e −| q (cid:107) | L/ ) ˆ q (cid:107) · A ˆ q (cid:107) (A25) W G0 ( q (cid:107) → ) = (cid:113) v G ( )˜ ε − G0 ( q (cid:107) ) (cid:113) v ( q (cid:107) )= − π (1 − e −| q (cid:107) | L/ ) | q (cid:107) | (cid:113) v G ( ) ˆ q (cid:107) · a G π (1 − e −| q (cid:107) | L/ ) ˆ q (cid:107) · A ˆ q (cid:107) (A26) W ( q (cid:107) → ) = (cid:113) v G ( )˜ ε − G0 ( q (cid:107) ) (cid:113) v ( q (cid:107) )= − π (1 − e −| q (cid:107) | L/ ) | q (cid:107) | (cid:113) v G ( ) ˆ q (cid:107) · b G π (1 − e −| q (cid:107) | L/ ) ˆ q (cid:107) · A ˆ q (cid:107) . (A27)and the body also gets a correction and becomes W GG (cid:48) ( q (cid:107) → ) = (cid:113) v G ( ) v G (cid:48) ( ) (cid:2) ε − GG (cid:48) ( q (cid:107) ) − δ GG (cid:48) (cid:3) = (cid:113) v G ( ) v G (cid:48) ( ) (cid:34) B − GG (cid:48) − δ GG (cid:48) + 4 π (1 − e −| q (cid:107) | L/ )( ˆ q (cid:107) · a G )( ˆ q (cid:107) · b G (cid:48) )1 + 4 π (1 − e −| q (cid:107) | L/ ) ˆ q (cid:107) · A ˆ q (cid:107) (cid:35) . (A28)Taking the limit q (cid:107) → we see that the head of thescreened potential is W ( q (cid:107) ) → ) = (2 πL ) ˆ q (cid:107) · A ˆ q (cid:107) which is a finite value.Defining the dimensionless quantity x = q (cid:107) L/ A = π (cid:82) π ˆ x ( φ ) · A ˆ x ( φ ) dφ it ispossible to rewrite the head of the screened potential as˜ w ( x ) = W (2 x /L )(2 πL ) A = − (cid:18) − e −| x | | x | (cid:19)
11 + 4 πA (1 − e −| x | ) . (A29)It is evident that the polar integral of Eq. (A29), over asmall circle with radius r Ω , cannot be evaluated analyt- ically: (cid:90) Ω ˜ w ( x ) d x = − π (cid:90) r Ω0 (cid:18) − e − x x (cid:19)
11 + 4 πA (1 − e − x ) x d x . (A30)It is however noticed that the function ˜ y ( x ) = πA ) x agrees very well with the integrand for small x . It hasthe same first order Taylor expansion and it is integrable.This yields (cid:90) Ω ˜ w ( x ) d x ≈ − π (cid:90) r Ω0 x πA ) x d x = − π (4 πAr Ω + r Ω − ln(4 πAr Ω + r Ω + 1)(4 πA + 1) . (A31)Since the expression in Eq. (A31) only holds for small x , it is generally not valid in the entire Ω region. Also,the expression is not valid for non-isotropic materials,where it is not justified to take the rotational average ofthe A tensor. K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz,Phys. Rev. Lett. , 136805 (2010). A. Splendiani, L. Sun, Y. Zhang, T. Li, J. Kim, C.-Y.Chim, G. Galli, and F. Wang, Nano Lett. , 1271 (2010). H. Zeng, G.-B. Liu, J. Dai, Y. Yan, B. Zhu, R. He, L. Xie,S. Xu, X. Chen, W. Yao, and X. Cui, Sci Rep (2013), 10.1038/srep01608. C. Zhang, A. Johnson, C.-L. Hsu, L.-J. Li, and C.-K. Shih,Nano Lett. , 2443 (2014). Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman,and M. S. Strano, Nature Nanotechnology , 699 (2012). L. Britnell, R. M. Ribeiro, A. Eckmann, R. Jalil, B. D.
Belle, A. Mishchenko, Y.-J. Kim, R. V. Gorbachev,T. Georgiou, S. V. Morozov, A. N. Grigorenko, A. K. Geim,C. Casiraghi, A. H. C. Neto, and K. S. Novoselov, Science , 1311 (2013). M. Bernardi, M. Palummo, and J. C. Grossman, NanoLett. , 3664 (2013). W. G. Aulbur, L. J¨onsson, and J. W. Wilkins, in
SolidState Physics , Vol. Volume 54, edited by Henry Ehrenreichand Frans Spaepen (Academic Press, 1999) pp. 1–218. F. Aryasetiawan and O. Gunnarsson, Rep. Prog. Phys. ,237 (1998). L. Hedin, Phys. Rev. , A796 (1965). M. S. Hybertsen and S. G. Louie, Phys. Rev. Lett. ,1418 (1985). R. W. Godby, M. Schl¨uter, and L. J. Sham, Phys. Rev.Lett. , 2415 (1986). M. Shishkin and G. Kresse, Phys. Rev. B , 035101(2006). T. Kotani and M. van Schilfgaarde, Solid State Communi-cations , 461 (2002). A. Marini, C. Hogan, M. Gr¨uning, and D. Varsano, Com-puter Physics Communications , 1392 (2009). C. Rostgaard, K. W. Jacobsen, and K. S. Thygesen, Phys.Rev. B , 085103 (2010). F. Caruso, P. Rinke, X. Ren, A. Rubio, and M. Scheffler,Phys. Rev. B , 075105 (2013). F. Bruneval, Phys. Rev. Lett. , 176403 (2009). X. Blase, C. Attaccalite, and V. Olevano, Phys. Rev. B , 115103 (2011). S. Ismail-Beigi, Phys. Rev. B , 233103 (2006). F. Huser, T. Olsen, and K. S. Thygesen, Phys. Rev. B ,245309 (2013). D. Y. Qiu, F. H. da Jornada, and S. G. Louie, Phys. Rev.Lett. , 119901 (E) (2015). F. A. Rasmussen and K. S. Thygesen, J. Phys. Chem. C , 13169 (2015). L. Keldysh, Jetp Lett. , 658 (1979),wOS:A1979JN00600015. P. Cudazzo, I. V. Tokatly, and A. Rubio, Phys. Rev. B , 085406 (2011). R. Sundararaman and T. A. Arias, Phys. Rev. B ,165122 (2013). M. S. Hybertsen and S. G. Louie, Phys. Rev. B , 5390(1986). R. M. Pick, M. H. Cohen, and R. M. Martin, Phys. Rev.B , 910 (1970). J. Deslippe, G. Samsonidze, D. A. Strubbe, M. Jain, M. L.Cohen, and S. G. Louie, Computer Physics Communica-tions , 1269 (2012). C. Freysoldt, P. Eggert, P. Rinke, A. Schindlmayr, andM. Scheffler, Phys. Rev. B , 235428 (2008). C. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, and A. Rubio, Phys. Rev. B , 205119 (2006). C. Freysoldt, P. Eggert, P. Rinke, A. Schindlmayr, R. W.Godby, and M. Scheffler, Computer Physics Communica-tions , 1 (2007). J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen, Phys.Rev. B , 035109 (2005). J. Enkovaara, C. Rostgaard, J. J. Mortensen, J. Chen,M. Du(cid:32)lak, L. Ferrighi, J. Gavnholt, C. Glinsvad,V. Haikola, H. A. Hansen, H. H. Kristoffersen, M. Kuisma,A. H. Larsen, L. Lehtovaara, M. Ljungberg, O. Lopez-Acevedo, P. G. Moses, J. Ojanen, T. Olsen, V. Petzold,N. A. Romero, J. Stausholm-Møller, M. Strange, G. A.Tritsaris, M. Vanin, M. Walter, B. Hammer, H. H¨akkinen,G. K. H. Madsen, R. M. Nieminen, J. K. Nørskov,M. Puska, T. T. Rantala, J. Schiøtz, K. S. Thygesen, andK. W. Jacobsen, Journal of Physics: Condensed Matter , 253202 (2010). F. H¨user, T. Olsen, and K. S. Thygesen, Phys. Rev. B ,235132 (2013). J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996). M. L. Tiago, S. Ismail-Beigi, and S. G. Louie, Phys. Rev.B , 125212 (2004). B.-C. Shih, Y. Xue, P. Zhang, M. L. Cohen, and S. G.Louie, Phys. Rev. Lett. , 146401 (2010). J. Klimeˇs, M. Kaltak, and G. Kresse, Phys. Rev. B ,075125 (2014). J. Yan, J. J. Mortensen, K. W. Jacobsen, and K. S. Thyge-sen, Phys. Rev. B , 245122 (2011). Z. Y. Zhu, Y. C. Cheng, and U. Schwingenschl¨ogl, Phys.Rev. B , 153402 (2011). A. Ramasubramaniam, Phys. Rev. B , 115409 (2012). T. Cheiwchanchamnangij and W. R. L. Lambrecht, Phys.Rev. B , 205302 (2012). H.-P. Komsa and A. V. Krasheninnikov, Phys. Rev. B ,241201 (2012). A. Molina-S´anchez, D. Sangalli, K. Hummer, A. Marini,and L. Wirtz, Phys. Rev. B , 045412 (2013). H. Shi, H. Pan, Y.-W. Zhang, and B. I. Yakobson, Phys.Rev. B , 155304 (2013). H. J. Conley, B. Wang, J. I. Ziegler, R. F. Haglund, S. T.Pantelides, and K. I. Bolotin, Nano Lett. , 3626 (2013). A. R. Klots, A. K. M. Newaz, B. Wang, D. Prasai,H. Krzyzanowska, J. Lin, D. Caudel, N. J. Ghimire, J. Yan,B. L. Ivanov, K. A. Velizhanin, A. Burger, D. G. Mandrus,N. H. Tolk, S. T. Pantelides, and K. I. Bolotin, Sci. Rep. (2014), 10.1038/srep06608. N. Berseneva, A. Gulans, A. V. Krasheninnikov, and R. M.Nieminen, Phys. Rev. B , 035404 (2013). V. Tran, R. Soklaski, Y. Liang, and L. Yang, Phys. Rev.B , 235319 (2014). A. Baldereschi and E. Tosatti, Phys. Rev. B17