How dielectric screening in two-dimensional crystals affects the convergence of excited-state calculations: Monolayer MoS 2
HHow dielectric screening in two-dimensional crystals affects the convergence ofexcited-state calculations: Monolayer MoS Falco H¨user ∗ and Thomas Olsen Center for Atomic-scale Materials Design (CAMD), Department of PhysicsTechnical University of Denmark, 2800 Kgs. Lyngby, Denmark
Kristian S. Thygesen † Center for Atomic-scale Materials Design (CAMD), Department of PhysicsTechnical University of Denmark, 2800 Kgs. Lyngby, Denmark andCenter for Nanostructured Graphene (CNG)Technical University of Denmark,2800 Kgs. Lyngby, Denmark (Dated: January 14, 2014)We present first-principles many-body calculations of the dielectric constant, quasiparticle bandstructure, and optical absorption spectrum of monolayer MoS using a supercell approach. As theseparation between the periodically repeated layers is increased, the dielectric function of the layerdevelops a strong q dependence around q = 0. This implies that denser k -point grids are requiredto converge the band gap and exciton binding energies when large supercells are used. In the limitof infinite layer separation, here obtained using a truncated Coulomb interaction, a 45 × k -pointgrid is needed to converge the G W band gap and exciton energy to within 0 . k -point sampling and the interactions between the repeated layershave opposite effects on the band gap and exciton energy, leading to a fortuitous error cancellationin the previously published results. PACS numbers: 71.20.Nr, 71.35.-y, 73.22.-f, 78.20.Bh, 78.60.LcKeywords: MoS2
I. INTRODUCTION
Atomically thin two-dimensional (2D) materials suchas graphene, hexagonal boron nitride, and transitionmetal dichalcogenides (TMDC) possess unique electronicand optical properties including high intrinsic carriermobilities, tunable band gaps, and strong light-matter interactions.
These features, combined withthe possibility of engineering their electronic propertiesfurther via strain, alloying or stacking, make the 2D ma-terials ideal as building blocks for new opto-electronicstructures and devices with minimal sizes and perfor-mances surpassing present technologies.After the intense focus on graphene, the TMDCs arenow attracting increasing interest. This stems mainlyfrom the greater variation in their electronic propertiesincluding both semiconducting and metallic behavior. Sofar, the most intensively studied single-layer TMDC isthe semiconductor MoS . Nanostructured forms of MoS have previously been explored as potential catalysts fordesulferization of crude oil and more recently for (photo)-electrochemical hydrogen evolution. Bulk MoS iscomposed of two-dimensional sheets held together byweak van der Waals forces and individual sheets can beisolated by exfoliation techniques similar to those used toproduce graphene. Single layers of MoS therefore com-prise highly interesting two-dimensional systems with afinite band gap and have recently been proposed for nano-electronics applications. The optical properties of MoS have been thoroughlystudied experimentally. The absorption spectrumshows two distinct low energy peaks at 1 .
88 eV and2 .
06 eV, which are denoted by A and B, respectively, and derive from direct transitions between a split valenceband and the conduction band at the K point of the Bril-louin zone. Their Rydberg satellites, Zeeman splitting,and dependence on crystal thickness have been investi-gated in detail. Recently, the quantum yield of lumi-nescence from MoS was shown to increase dramaticallywhen the sample thickness was changed from a few layersto a monolayer indicating a transition to a direct bandgap in the single layer.In the past couple of years a number of theoreticalstudies of the electronic band structure and optical exci-tations in monolayer MoS have been published. These studies are based on many-body perturbationtheory in the GW approximation (mainly the non-selfconsistent G W approach) for the band structureand the Bethe-Salpeter equation (BSE) with a staticallyscreened electron-hole interaction for the optical excita-tions. As is standard practice, the calculations have beenperformed on a supercell geometry where the MoS lay-ers have been separated by 10 −
20 ˚A vacuum and theBrillouin Zone (BZ) sampled on grids ranging from 6 × ×
15. With these parameters G W band gaps inthe range 2 . − . W -BSE exciton bindingenergies of 0 . − . and indirect band gaps have been foundat the G W level, while only direct gaps have been ob- a r X i v : . [ c ond - m a t . m t r l - s c i ] J a n tained with self-consistent GW and GW . Whencomparing these values, it should be kept in mind thatboth size and nature of the band gap of MoS dependssensitively on the in-plane lattice parameter, a . One of the most fundamental quantities describing theelectronic structure of a material is the dielectric func-tion. The dielectric properties of atomically thin 2D ma-terials are quite different from their 3D counterparts. For example plasmons in 2D metals have acoustic dis-persion relations ( ω p ( q ) → q → obtained using thesupercell approach lie in the range 4 . − . These values have been used to ra-tionalize the exciton binding energy in MoS using thesimple Mott-Wannier model.In this paper, we present an in-depth study of the di-electric function, quasiparticle band structure and ex-citonic states in monolayer MoS . We focus on sepa-rating the spurious interlayer screening from the intrin-sic intralayer screening in supercell geometries, and theconsequences that the physics of screening in 2D hasfor the convergence of many-body excited state calcu-lations. The 3D macroscopic dielectric constant, as usedfor solids, converges to 1 for all q vectors in the limitof infinite separation of the layers and is thus meaning-less for a 2D material. We use an alternative approachto calculate the dielectric constant by averaging the totalfield over the material rather than the supercell. This 2Ddielectric constant shows strong q -dependence for smallwave vectors and becomes exactly 1 for q = 0. This prop-erty has important consequences for the k -point conver-gence of many-body calculations.In general, the use of a truncated Coulomb interac-tion is essential to avoid interlayer screening which de-cays slowly with the layer separation, L . The interlayerscreening yields too large dielectric constant for wave vec-tors q < /L . As a consequence, the G W band gapsand exciton energies are 0 . q -dependence of the dielectric constantfor small q implies that a k -point grid of at least 45 × . k -point grids below 15 ×
15 the band gapis at least 0 . L → ∞ . Thusthe effect of interlayer screening and too coarse k -pointgrids partially cancel out leading to reasonable values forthe band gap and exction binding energy with undercon-verged parameters as applied in previous studies.The paper is organized as follows. In Sec. II we presentG W band structures and study the convergence of thegap with respect to interlayer separation and k -pointsampling. In Sec. III we show calculations for the 2Ddielectric constant and explain the origin of the slow k -point convergence of the band gap. In Sec. IV we presentmany-body calculations of the lowest excitons and anal-yse their convergence with layer separation and k -point sampling. Our conclusions are given in Sec. V. II. QUASIPARTICLE BAND STRUCTURE
In this section we demonstrate that GW band struc-tures for monolayer MoS converge extremely slow withrespect to the interlayer separation. In order to obtainwell converged results (within 0 . k -point grid of around 45 ×
45. Previously reported calcu-lations with the full Coulomb interaction have employedonly separation between 10 and 20 ˚A and used from 6 × × k points. The resulting band structures are,however, somewhat saved by a fortunate error cancella-tion between the two effects. A. Computational details
All our calculations have been performed with the pro-jector augmented wave method code GPAW.
TheKohn-Sham wave functions and energies of monolayerMoS were calculated in the local density approxima-tion (LDA) using a plane wave basis with cut-off energy400 eV. The 4 s and 4 p semicore electrons of Mo wereexplicitly included in all calculations. Unless otherwisestated the calculations have been performed for the ex-perimental lattice constant of 3 .
16 ˚A. One-shot G W calculations were performed using the LDA wave func-tions and eigenvalues to obtain the G W @LDA quasi-particle energies. A plane wave cut off of 50 eV and200 bands were used for the dielectric function, screenedinteraction and GW self-energy. Convergence with re-spect to these parameters has been checked very carefully.With these values band gaps were found to be convergedwithin around 10 meV. The plasmon pole approximationfor the dielectric function was found to yield QP energieswithin 0 . sheets, we have applied a truncated Coulomb in-teraction of the form v c ( r ) = (1 /r ) θ ( R c − z ), followingRefs. 30 and 31. For details on the implementation ofthe GW method in the GPAW code we refer to Ref. 32.We note that we have used a numerical averaging of thehead of the screened potential W ( q ) for all wavevectors q in the Brillouin zone (similar to Eq. 9 in Sec. III C).This was found to be crucial in all calculations with thefull 1 /r Coulomb interaction.
B. Results
The band structure calculated using 45 × k pointsand the truncated Coulomb interaction is shown in Fig.1. At the LDA level, we find a direct band gap at the Kpoint of 1 .
77 eV while the smallest indirect gap of 1 .
83 eV
FIG. 1: (Color online). Bandstructure of monolayer MoS calculated with LDA and G W @LDA using 45 × k pointsand a truncated Coulomb interaction to avoid interaction be-tween periodically repeated layers. The valence band topshave been aligned.FIG. 2: Definition of the interlayer separation, L . occurs from Γ to a point along the Γ-K direction. Incontrast, G W predicts an indirect gap of 2 .
58 eV andat direct gap at K of 2 .
77 eV.In Fig. 3 we show the convergence of both the directand the indirect band gap with respect to the k -point gridfor a fixed interlayer separation of 23 ˚A (see Fig. 2 forthe definition of L ). It is clear that a very dense k -pointgrid is needed in order to obtain well-converged resultswith the truncated Coulomb interaction. For 45 × k points, band gaps are converged within less than 0 . × k points withthe bare Coulomb interaction. However, the calculatedvalues are too low. The slow convergence with respect to k points when the truncation is used will be discussed in FIG. 3: (Color online). Direct (full symbols) and indirect(open symbols) G W band gaps as function of the number of k points in one of the in-plane directions for a layer separationof L = 23 ˚A.FIG. 4: (Color online). Direct G W band gap plotted as afunction of interlayer distance for different k -point samplingswith the full 1 /r interaction. Dotted lines serve as a guide forthe eye to extrapolate for L → ∞ . They were obtained byfitting all values for L >
30 ˚A, including the results with theCoulomb truncation, to a quadratic function. Dashed hori-zontal lines indicate the calculated values with the truncatedCoulomb interaction. detail in Sec. III C.We see that results do not converge independently withrespect to the number of k points and the interlayer sep-aration. In Fig. 4, we plot the L -dependence of thedirect band gap for different k -point samplings with thebare interaction. The k -point dependence becomes muchstronger for large L . For L → ∞ , the values are expectedto converge to the results calculated with the truncation(indicated by dotted lines). They seem to exhibit a linear1 /L behaviour only for L >
50 ˚A or very dense k -pointsamplings. Fig. 5 shows all results and interpolated val-ues in a contourplot as a function of 1 /L and numberof k points. The effects of using more k points and in- FIG. 5: (Color online). Contour plot of the direct G W band gap as a function of the inverse interlayer distance andnumber of k points in one of the in-plane directions with thefull 1 /r interaction. Contour lines are separated by 0 . creasing L are of different sign and partially cancel eachother. This is the reason, why different choices of thetwo parameters yield the same results. Especially, theband gaps calculated with 9 × k points and L = 23 ˚Aand 15 × k points and L = 43 ˚A are the same as with45 × k points and infinite L . This seems, however,conincidental and we do not expect it to be the case forother systems.We note that all calculations have been performed witha single k point in the direction perpendicular to thelayer. This is, however, insufficient for small interlayerdistances. For L = 13 ˚A, we find an increase of the bandgap of around 0 . . k points areused, for example. For L >
20 ˚A or use of the truncation,this effect is negligible.
C. Comparison with previous work
In table I we show our converged results obtained withthe truncated Coulomb interaction and 45 × k pointstogether with previous G W results from the literature.For each reference we show the values used for the lat-tice constant, the interlayer separation and the k -pointsampling. It can be seen that all the previous calcula-tions have used small layer separations and no trunca-tion method. As pointed out in the preceding discussion,this gives a fast k -point convergence. A properly con-verged calculation, however, requires larger separationsand thereby more k points. But as a consequence of acancellation of errors, a calculation with 19 ˚A layer sepa-ration and 12 × k points yields almost the same bandgaps as our converged result (within 0 .
15 eV). We arethus led to conclude that the reasonable agreement be- tween our results and previous ones are to a large extentfortuitous.Furthermore, the effect of strain can have a large im-pact on the MoS band gap. As demonstrated in Ref.4, using 12 × k points and 19 ˚A layer separation, theG W band gap for the experimental lattice constant of3 .
160 ˚A is indirect. With a lattice constant of 3 .
190 ˚A,corresponding to 1% strain, the gap changes to direct.The lowering of the direct band gap becomes even morepronounced for larger lattice constants. As can be seenfrom the table our converged results predict the sametrend, in particular the decrease of the direct gap as func-tion of strain, with our values for the direct gap beinggenerally 0 . .
255 ˚A, thesmallest indirect transition occurs from the Γ point atthe valence band to the K point at the conduction band.This is also in agreement with Ref. 4. In the partiallyself-consistent GW calculations of Ref. 26, the oppo-site trend was found, namely a transition from a directto indirect band gap for ∼ k points in thein-plane directions were used.In Ref. 22, the band gap was determined by extrap-olating from L = 20 ˚A to infinite layer separation, un-der the assumption that the gap scales linearly with theinverse distance between the layers. The obtained val-ues for the direct and indirect band gaps are ∼ . ∼ . .
18 ˚A and the same (under-converged) k -point grid of 12 ×
12 as in Ref. 22.From our studies, we conclude that the G W @LDAband gap of monolayer MoS is indirect with a value of2 . . .
16 ˚A is used. The question of howwell the one-shot G W @LDA approximation describesthe true electronic structure of this system remains open.Partially self-consistent GW and fully self-consistentGW calculations have been shown to consistently yielddirect band gaps of 2 . − .
80 eV.
III. STATIC SCREENING
In this section we present a detailed investigation ofthe (static) dielectric properties of monolayer MoS . Thisserves a dual purpose. First, it illustrates the origin of theslow convergence of the GW results presented in the pre-vious section (and the BSE results presented in the nextsection). Secondly, it shows that the usual definition ofthe macroscopic dielectric constant of a periodic solid isnot meaningful when applied to a 2D system representedin a periodic supercell. We discuss the difference betweenscreening in 2D and 3D which becomes particularly pro-nounced in the q → TABLE I: Calculated G W band gaps obtained in present work and compared with previous results from the literature. Allour calculations have been performed using a truncated Coulomb interaction. E gap (eV)Ref. starting point a (˚A) number of k -points layer separation (˚A) direct indirectThis work LDA 3.16 45 × × v c ) 2.77 2.58This work LDA 3.19 45 × × v c ) 2.65 2.57This work LDA 3.255 45 × × v c ) 2.41 2.51Ref. 25 LDA 3.15 18 × × ∼ . × × ∼ × × ∼ × × × × × × /L extrapolation 2.97 3.26Ref. 22 PBE 3.18 12 × × ∼ ∼ × × × × ∼ × × A. 3D macroscopic dielectric constant
The microscopic dielectric function determines the re-lation between a weak external potential and the totalpotential in the material, V tot ( r ) = (cid:90) d r (cid:48) (cid:15) − ( r , r (cid:48) ) V ext ( r (cid:48) ) . (1)For a periodic system the dielectric function can be con-veniently expressed in plane waves (cid:15) − ( r , r (cid:48) ) = (cid:88) GG (cid:48) (cid:88) q e i ( G + q ) r (cid:15) − GG (cid:48) ( q ) e − i ( G (cid:48) + q ) r (cid:48) , (2)where G is a reciprocal lattice vector, q belongs to the1. BZ. Within the random phase approximation (RPA)we have (cid:15) GG (cid:48) ( q , ω ) = δ GG (cid:48) − v c ( q + G ) χ GG (cid:48) ( q , ω ) , (3)where χ is the non-interacting density response func-tion. Here, v c can be the Fourier representation of eitherthe full or the truncated Coulomb interaction. For thecalculations in this section we have used a 50 eV cut-off for the reciprocal lattice vectors to account for localfield effects. The non-interacting response function, χ ,was constructed from local density approximation (LDA)wave functions and energies including states up to 50 eVabove the Fermi level. All calculations were performedwith the projector augmented wave method code GPAW.Details on the implementation of the dielectric functionin the GPAW code can be found in Ref. 34.It follows from Eq. (2) that the total potential result-ing from a plane wave external potential V e i q · r has theform V tot ( r ) = (cid:101) V q ( r ) e i q · r (4) where (cid:101) V q ( r ) is a lattice periodic function. We thus definethe macroscopic dielectric constant as1 (cid:15) M ( q ) ≡ (cid:104) (cid:101) V q (cid:105) Ω V = (cid:15) − ( q ) , (5)where (cid:104) . . . (cid:105) Ω denotes a spatial average over a unit cell.Note that in general (cid:15) M ( q , ω ) (cid:54) = (cid:15) ( q , ω ) because of localfield effects. To explicitly demonstrate that Eq. (5) does not ap-ply to low-dimensional materials, we have calculated themacroscopic dielectric constant as a function of the layerseparation, L . The results are shown in Fig. 6 for differ-ent values of the in-plane momentum transfer q . We alsoshow the dielectric constant corresponding to polariza-tion orthogonal to the layer. Clearly the macroscopic di-electric constant approaches unity for all q -vectors in thelimit of large interlayer separation. This occurs becausethe total field is averaged over an increasingly larger vac-uum region.Previously reported values for the macroscopic di-electric constant of monolayer MoS lie in the range4 − In these calculations the MoS layers wereseparated by 10 −
20 ˚A vacuum. As can be seen from (cid:15) (cid:107) ( q = 0) in Fig. 6 this is consistent with our results.However, it should also be clear that numbers depend onthe distance between layers and in fact are not meaning-ful. B. 2D macroscopic dielectric constant
For a 2D material, the average of the total potential inthe definition of the macroscopic dielectric constant mustbe confined to the region of the material. Since Eq. (4)still holds for a 2D material when q is confined to theplane of the material, we average the in-plane coordi-nates ( r (cid:107) ) over the unit cell area A and the out-of-plane FIG. 6: (Color online). The 3D static macroscopic dielectricconstant 1 /(cid:15) − ( q ) of monolayer MoS as a function of theinterlayer separation, L . (cid:15) (cid:107) is the dielectric constant withpolarization parallel to the monolayer and (cid:15) ⊥ is the dielectricconstant for polarization orthogonal to the layer. coordinate ( z ) from z − d/ z + d/ z de-notes the center of the material and d its width. The 2Dmacroscopic dielectric constant then becomes1 (cid:15) DM ( q (cid:107) ) ≡ (cid:104) (cid:101) V q (cid:105) A,d V = 2 d (cid:88) G ⊥ e iG ⊥ z sin( G ⊥ d/ G ⊥ (cid:15) − G0 ( q (cid:107) ) , (6)where the sum is over all G with G (cid:107) = . In this workwe have taken d = 6 .
15 ˚A corresponding to the interlayerseparation in bulk MoS . We shall return to the problemof chosing d below.The results for the static dielectric constant evaluatedfrom Eq. (6) using the bare Coulomb interaction is shownin Fig. 7 for four different layer separations. The re-sult for L = d = 6 .
15 ˚A coincides with the 3D dielectricconstant of bulk MoS given by Eq. (5). The result ob-tained with the truncated Coulomb interaction is shownin black; it represents the case of infinite layer separation.Before discussing the results, it is instructive to considerthe potential arising from a 2D charge density fluctuationof the form, n ( r ) = n e i q (cid:107) · r (cid:107) δ ( z ) , (7)The corresponding potential follows from Poisson’sequation φ ( r ) = n q (cid:107) e − i q (cid:107) · r (cid:107) e − q (cid:107) | z | . (8)It follows that the potential perpendicular to the layerfalls off exponentially over a characteristic distance of1 /q (cid:107) . This explains why in general (cid:15) DM ( q (cid:107) ) coincideswith the isolated layer result for q (cid:107) (cid:38) /L .The variation of (cid:15) DM when the parameter d is changedby ±
10% is indicated by the shaded region in Fig. 7.
FIG. 7: (Color online). Static macroscopic dielectric constantfor a single layer of MoS calculated along the Γ-K line. Thecalculations are performed using Eq. (6) with the microscopicdielectric constant, (cid:15) − GG (cid:48) ( q ), evaluated from Eq. (3) with ei-ther the bare Coulomb interaction (dotted and dashed lines)or truncated Coulomb interaction (full black line). The greyarea represents the result obtained when the averaging regionperpendicular to the layer, d , is varied by ± .
15 ˚A and thus co-incide with the dielectric constant of bulk MoS . The curveshave been interpolated from a 32 × q -point mesh. To the left of the maximum, (cid:15) DM ( q (cid:107) ) is insensitive to d since the induced potential is more or less constantover the averaging region. To the right of the maximum,the variation in (cid:15) DM ( q (cid:107) ) follows the ±
10% variation in d . This is because for these wave vectors the inducedpotential has essentially vanished at the borders of theaveraging region. In general, increasing (decreasing) d will decrease (increase) (cid:15) DM ( q (cid:107) ) in the large wave vectorregion.Another characteristic feature of the potential in Eq.(8) is the 1 /q (cid:107) scaling which should be compared withthe 1 /q form of the Coulomb potential from a 3D chargeoscillation. Since the non-interacting response function, χ ( q ), scales as ∼ q for q → (cid:15) DM (0) = 1, whilethis is in general not the case in 3D. In our calcula-tions, the effect of interlayer interactions is eliminatedby using a truncated Coulomb interaction of the form v c ( r ) = (1 /r ) θ ( R c − z ). For q z = 0 and in the limitof small q (cid:107) , the Fourier representation of the truncatedCoulomb interaction becomes v D ( q ) = πR c | q | , i.e. itscales as 1 /q as the potential from the 2D charge den-sity wave ensuring the correct limit (cid:15) DM (0) = 1.Finally, we note that previous studies have em-ployed a strict 2D model for the dielectric function inthe small q limit of the form (cid:15) ( q (cid:107) ) = 1 + αq (cid:107) . This formis convenient as it leads to an analytical expression forthe screened interaction in 2D. Our definition differs bybeing a 3D (or quasi 2D) quantity valid for general q (cid:107) . C. Screened interaction
In Fig. 8 we show (cid:15) − as a function of q (cid:107) evaluatedwith and without the truncated Coulomb interaction.For small q , the two curves differ significantly due tothe long range nature of the induced potential (8). Atlarge q ( ∼ K/ , since they givethe average potential over the supercell and not over theMoS layer. In particular their value will be highly de-pendent on the size of the unit cell (in the limit of infinitelayer separation both will equal 1 for all q ). Nevertheless,this quantity is a crucial ingredient of both the GW self-energy and the BSE kernel as it provides the screeningof the divergent term of the Coulomb interaction.For q = 0 the Coulomb kernel diverges and we approx-imate W ( q = 0) by the integral W ( q = 0) = 1Ω Γ (cid:90) Ω Γ d q v c ( q ) (cid:15) − ( q ) ≈ Γ (cid:15) − ( q = 0) (cid:90) Ω Γ d q v c ( q ) , (9)where Ω BZ is the Brillouin zone volume and Ω Γ is a smallvolume containing q = 0. In isotropic systems (cid:15) − ( q ) isusually constant in the vicinity of q = 0 and the approxi-mation works well. However, when (cid:15) − is evaluated withthe truncated Coulomb interaction, (cid:15) − acquires muchmore structure for small q as can be seen from Fig. 8.Thus, for coarse k -point samplings we will underestimatethe Γ-point screening since we simply use (cid:15) − ( q = 0) = 1.The linear behavior of the screened interaction forsmall q suggests that a better approximation for W ( q =0) would be W ( q = 0) = 1Ω Γ (cid:90) Ω Γ d q v c ( q ) (cid:104) q · ∇ q (cid:15) − ( q ) (cid:12)(cid:12)(cid:12) q =0 (cid:105) . (10)Since the dielectric matrix in RPA is (cid:15) GG (cid:48) ( q ) = 1 − v c ( q ) χ GG (cid:48) ( q ), we can derive an analytic expression forthe first order Taylor expansion in q , and its inverse.These quantities can be evaluated with vanishing addi-tional cost, but we will leave the assessment of this ap-proximation to future work. IV. OPTICAL ABSORPTION SPECTRUM
In this section we present many-body calculations ofthe optical absorption spectrum of monolayer MoS bysolving the Bethe-Salpeter Equation (BSE) under thestandard assumption of static screening of the electron-hole interaction. As for the GW band gap, we find thatthe use of a truncated Coulomb interaction is essential FIG. 8: (Color online). The 3D static inverse dielectric con-stant (cid:15) − ( q ) of monolayer MoS calculated in the RPA fordifferent values of in-plane momentum transfer q along the Γ-K direction. The separation between layers is L = 20 ˚A. Notethat neither of the quantities can be interpreted as the macro-scopic dielectric constant of the monolayer (this quantity isthe black curve in Fig. 7). to avoid interlayer screening and obtain well convergedexciton binding energies. Furthermore, the very strongq-dependence of the 2D static dielectric function around q = 0, leads to very slow k -point convergence for theexciton binding energy.In order to obtain an accurate absorption spectrum in-cluding excitonic effects we calculate the response func-tion from the Bethe-Salpeter Equation (BSE). Using thestandard assumption of a static dielectric screening of theelectron-hole (e-h) interaction, the BSE can be recastas an effective two-particle Hamiltonian, which is diag-onalized on a basis of electron-hole pairs. In this way theexcitonic eigenstates can be expressed as a linear combi-nation of single-particle transitions | λ (cid:105) = (cid:88) vck A λvck | vck (cid:105) , (11)where v , c , and k denote valence band, conduction bandand Brillouin zone wave vector, respectively. The absorp-tion spectrum is proportional to the imaginary part ofthe macroscopic dielectric function, which in the Tamm-Dancoff approximation can be written (cid:15) ( ω ) =2 π lim q → v c ( q ) (cid:88) λ δ ( ω − E λ ) × (cid:12)(cid:12)(cid:12)(cid:12) (cid:88) vc k A λvck (cid:104) v k − q | e − i q · r | c k (cid:105) (cid:12)(cid:12)(cid:12)(cid:12) , (12)where E λ are the eigenvalues associated with | λ (cid:105) .In all calculations we have included a single valenceband and a single conduction band in the BSE Hamilto-nian. We have tested that the first excitonic peak is com-pletely unaffected if we instead include 6 valence bandsand 4 conduction bands. This is also expected since thehighest (lowest) valence (conduction) band is well iso-lated from the remaining bands at K where the excitonis centered, see Fig. 1. For the screening we have in-cluded 65 bands in the evaluation of the response func-tion, which is sufficient for converged results. Increasingthe number of bands to 300 affects the position of thefirst exciton by less than 10 meV. The plane wave cutofffor the response function (local field effects) was set to50 eV and we checked that the excitonic binding energychanged by less than 10 meV when increasing the cutoffto 200 eV. The dependence on k -point sampling and in-terlayer separation will be examined below. Details onthe implementation of the BSE method in the GPAWcode can be found in Ref. 5. A. Convergence tests
In the lower panel of Fig. 9, we show the excitonbinding energy as a function of interlayer separation cal-culated for different k -point samplings using the bareCoulomb interaction and the truncated Coulomb inter-action. With the bare Coulomb interaction, the obtainedresults are far from convergence, even for L = 50 ˚A. Thedependencies on the layer separation and number of k points is very similar as for the quasiparticle gap dis-cussed in Sec. II B, even on a quantitative level. There-fore, the optical gap, which is given by the differenceof the QP gap and the exciton binding energy, is al-most indenpedent of L and whether or not the truncationmethod is used, as shown in the upper panel. This is con-sistent with the observations in Ref. 22.The convergence of the binding energy with respect tothe k -point sampling is plotted in Fig. 10 for an inter-layer separation of 20 ˚A. The truncated Coulomb ker-nel gives a much slower convergence with respect to thenumber of k -points than the bare Coulomb interaction.However, it should be clear from Fig. 9 that the bind-ing energy obtained with the bare Coulomb interactionconverges to a value which is highly dependent on theinterlayer separation. In Ref. 25, convergence was foundwith 18 × k points, but for a layer separation of only24 ˚A. The obtained exciton binding energy was around0 . k -point convergence observed when usingthe truncated Coulomb interaction is related to the q -dependence of the screening in two-dimensional systems.As demonstrated by Eq. (9) and Fig. 8 (blue curve),a too low k -point sampling leads to an underestimationof the screening in the vicinity of q = 0 and thus anoverestimation of the exciton binding energy. B. Results
From the convergence tests described above we con-clude that the BSE calculations are (nearly) converged if
FIG. 9: (Color online). Optical gap and binding energy of thelowest exciton in monolayer MoS as a function of interlayerseparation calculated from the BSE and the G W quasipar-ticle gap. Results with the full 1 /r Coulomb interaction (fulllines) and the truncated interaction (dashed lines) are shown.Dotted lines give an estimation for extrapolation to infinite L .FIG. 10: (Color online). Binding energy of the lowest excitonin monolayer MoS as a function of k -point sampling for asupercell with a layer separation of L = 20 ˚A. we use a truncated Coulomb interaction and a 45 × k -point sampling. With these setting we have calculatedthe BSE spectrum on top of a G W quasiparticle bandstructure obtained with the same parameters. The re-sulting absorption spectrum is shown in Fig. 11. We alsoshow an RPA calculation, i.e. neglecting electron-holeinteractions in the BSE, performed on top of the sameG W band structure for comparison. With electron-hole interaction included, we obtain an exciton bindingenergy of ∼ . exhibits a spin-orbit split peak around1 . Since we have not included spin-orbit couplingin our calculations, the spectrum Fig. 11 only showsa single peak at low energies. However, it has previ-
FIG. 11: (Color online). Absorption spectrum of single layerMoS calculated with the RPA and BSE using the G W quasiparticle band structure. The calculation has been per-formed with a truncated Coulomb interaction to avoid inter-actions between repeated layers and with a 45 × k -pointgrid. ously been shown that the spin-orbit coupling doesnot have a large effect on the exciton binding energy andonly results in a split excitonic peak. The main peak inthe BSE@G W spectrum is situated at 2 . . W or underestima-tion of the exciton binding energy in BSE with a staticelectron-hole interaction.From the above discussion it should be clear that itis extremely challenging to converge the exciton bindingenergy with respect to interlayer separation and k points.In general, the optical gap is much easier to converge withrespect to interlayer separation than either the quasipar-ticle gap or the exciton binding energy. Nevertheless,for many physical applications it is of importance to ob-tain accurate values for both the quasiparticle gap andthe exciton binding energy separately. In Ref. 22 the ex-citon binding energy was obtained by 1 /L extrapolationof the quasiparticle gap calculated in a range of interlayerseparations between 10 and 20 ˚A and assuming the samedependence for the exciton binding energy. Our resultsindicate that one should be cautious with such extrapola-tions. This is because the screening at different q -pointshas a very different dependence on interlayer separation,which results in different convergence behavior at differ-ent k -point samplings (see Fig. 9 full lines). The extrap-olation procedure may therefore not give reliable results,since higher k -point samplings are required at larger in-terlayer separation. We are aware that the convergenceissues may depend a lot on the implementation of theBSE method. However, we have previously performedthe same calculations with YAMBO code, which pro- duced very similar convergence behavior for quasiparti-cle gap and exciton binding energy (also using truncatedCoulomb cutoff and 45 × k -point sampling). V. CONCLUSIONS
We have presented a careful investigation of the di-electric properties, band gap and excitonic states ina two-dimensional semiconductor exemplified by mono-layer MoS . We have demonstrated that the ”standard”macroscopic dielectric constant used for solids is not ap-plicable (meaningless) to supercells describing the 2Dmaterial as an infinite array of parallel sheets, and there-fore replaced it by a 2D version in which the induced fieldis averaged over the extent of the material rather thanover the entire supercell. We showed that the effect ofinterlayer screening leads to underestimation of the bandgap and exciton binding energy by up to more than 0 . <
30 ˚A. The reason for this is thatinterlayer screening affects (cid:15) ( q ) for q < /L where L isthe distance between layers in the supercell. Since it isthe small q limit of (cid:15) ( q ) that is most important for thescreened interaction W ( q ) = (cid:15) − ( q ) /q , the effect cannotbe neglected. Here we have circumvented the problemby using a truncated Coulomb interaction that explicitlycuts off the interaction between neighboring layers.The properly defined 2D dielectric function (cid:15) DM ( q ) hasa very sharp wave vector dependence for small q and sat-isfies (cid:15) DM (0) = 1 in general. This has the consequencethat quasiparticle- and optical excitations obtained fromthe GW and Bethe-Salpeter Equation, respectively, re-quire much denser k -point grids than experience from3D systems would suggest. For MoS we find that a pre-cision of 0 . k -point grids of at least 30 × k -point grids have opposite effects on the band gapand exciton energies leading to fortuitous error cancella-tion. Our calculations applying the truncated Coulombinteraction and 45 × k points give G W @LDA gapsof 2 .
77 eV (direct) and 2 .
58 eV (indirect) and binding en-ergy of the lowest exciton of 0 .
55 eV. This places the low-est exciton at ∼ . . W @LDA approximation or the use of static screeningin the BSE. VI. ACKNOWLEDGEMENT
KST acknowledges support from the Danish Coun-cil for Independent Research’s Sapere Aude Programthrough grant no. 11-1051390. The Center for Nanos-tructured Graphene (CNG) is sponsored by the DanishNational Research Foundation, Project DNRF58.0 ∗ Electronic address: [email protected] † Electronic address: [email protected] K. S. Novoselov, D. Jiang, F. Schedin, T. J. Booth, V. V.Khotkevich, S. V. Morozov, and A. K. Geim, Proc. Natl.Acad. Sci. U.S.A. , 10451 (2005). B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti,and A. Kis, Nature Nanotechnology , 147 (2011). K. Kaasbjerg, K. S. Thygesen, and K. W. Jacobsen,PHYSICAL REVIEW B , 115317 (2012). H. Shi, H. Pan, Y.-W. Zhang, and B. I. Yakobsen, Phys.Rev. B , 155304 (2013). J. Yan, K. W. Jacobsen, and K. S. Thygesen, PHYSICALREVIEW B , 045208 (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, et al., SCIENCE , 1311(2013). 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). M. Bernardi, M. Palummo, and J. C. Grossman, NanoLett. (to be published). Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman,and M. S. Strano, Nature Nanotechnology , 699 (2012). M. V. Bollinger, J. Lauritsen, K. W. Jacobsen, J. K.Nørskov, S. Helveg, and F. Besenbacher, Phys. Rev. Lett. , 196803 (2001). T. F. Jaramillo, K. P. Jørgensen, J. Bonde, J. H. Nielsen,S. Horch, and I. Chorkendorff, Science , 100 (2007). X. Zong, H. Yan, G. Wu, G. Ma, F. Wen, L. Wang, andC. Li, J. Am. Chem. Soc. , 7176 (2008). R. F. Frindt and A. D. Yoffe, Proc. R. Soc. Lond. A ,69 (1963). B. L. Evans and P. A. Young, Proc. R. Soc. Lond. A ,402 (1965). J. A. Wilson and A. D. Yoffe, Adv. Phys. , 193 (1969). R. A. Neville and B. L. Evans, Phys. Stat. Sol. B , 597(1976). G. L. Frey, S. Elani, M. Homyonfer, Y. Feldman, andR. Tenne, Phys. Rev. B , 6666 (1998). J. P. W. Newcomer and G. A. Samara, J. Appl. Phys. ,7934 (1997). R. Coehoorn, C. Haas, and R. A. de Groot, Phys. Rev. B , 6203 (1987). T. Cheiwchanchamnangij and W. R. L. Lambrecht, Phys.Rev. B , 205302 (2012). H.-P. Komsa and A. V. Krasheninnikov, Phys. Rev. B ,241201 (2012). C. Ataca and S. Ciraci, The Journal ofPhysical Chemistry C , 13303 (2011), http://pubs.acs.org/doi/pdf/10.1021/jp2000442, URL http://pubs.acs.org/doi/abs/10.1021/jp2000442 . A. Ramasubramaniam, Phys. Rev. B , 115409 (2012). A. Molina-S´anchez, D. Sangalli, K. Hummer, A. Marini,and L. Wirtz, Phys. Rev. B , 045412 (2013), URL http://link.aps.org/doi/10.1103/PhysRevB.88.045412 . H. J. Conley, B. Wang, J. I. Ziegler, R. F. Haglund, S. T.Pantelides, and K. I. Bolotin, Nano Lett. , 3626 (2013). L. V. Keldysh, JETP Lett. , 658 (1979). A. Molina-Sanchez and L. Wirtz, Phys. Rev. B , 155413(2011). J. Enkovaara, C. Rostgaard, J. J. Mortensen, J. Chen,M. Duak, L. Ferrighi, J. Gavnholt, C. Glinsvad,V. Haikola, H. A. Hansen, et al., Journal of Physics: Con-densed Matter , 253202 (2010), URL http://stacks.iop.org/0953-8984/22/i=25/a=253202 . C. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, andA. Rubio, Phys. Rev. B , 205119 (2006). S. Ismail-Beigi, Phys. Rev. B , 233103 (2006). F. H¨user, T. Olsen, and K. S. Thygesen, Phys. Rev. B ,235132 (2013), URL http://link.aps.org/doi/10.1103/PhysRevB.87.235132 . D. Y., Y. Wang, J. Ni, L. Shi, S. Shi, and W. Tang, PhysicaB , 2254 (2011). J. Yan, J. J. Mortensen, K. W. Jacobsen, and K. S.Thygesen, Phys. Rev. B , 245122 (2011), URL http://link.aps.org/doi/10.1103/PhysRevB.83.245122 . S. ADLER, PHYSICAL REVIEW , 413+ (1962), ISSN0031-899X. N. Wiser, Phys. Rev. , 62 (1963), URL http://link.aps.org/doi/10.1103/PhysRev.129.62 . P. Cudazzo, I. V. Tokatly, and A. Rubio, Phys. Rev. B ,085406 (2011), URL http://link.aps.org/doi/10.1103/PhysRevB.84.085406 . T. C. Berkelbach, M. S. Hybertsen, and D. R. Reichman,Phys. Rev. B , 045318 (2013), URL http://link.aps.org/doi/10.1103/PhysRevB.88.045318 . E. E. Salpeter and H. A. Bethe, Phys. Rev. , 1232(1951). G. Strinati, Phys. Rev. B , 5718 (1984). L. Wirtz, A. Marini, and A. Rubio, Phys. Rev. Lett. ,126104 (2006). A. Marini, C. Hogan, M. Gr¨uning, and D. Varsano, Comp.Phys. Comm. , 1392 (2009). The gpaw code is available as a part of the CAMPOS soft-ware: This is most easily seen by performing a 3D Fourier trans-formation of δn , then multiplying by 1 /q2