Critical assessment of $G_0W_0$ calculations for 2D materials: the example of monolayer MoS$_2$

Two-dimensional (2D) materials combine many fascinating properties that make them more interesting than their three-dimensional counterparts for a variety of applications. For example, 2D materials exhibit stronger electron-phonon and electron-hole interactions, and their energy gaps and effective carrier masses can be easily tuned. Surprisingly, published band gaps of several 2D materials obtained with the $GW$ approach, the state-of-the-art in electronic-structure calculations, are quite scattered. The details of these calculations, such as the underlying geometry, the starting point, the inclusion of spin-orbit coupling, and the treatment of the Coulomb potential can critically determine how accurate the results are. Taking monolayer MoS$_2$ as a representative material, we employ the linearized augmented planewave + local orbital method to systematically investigate how all these aspects affect the quality of $G_0W_0$ calculations, and also provide a summary of literature data. We conclude that the best overall agreement with experiments and coupled-cluster calculations is found for $G_0W_0$ results with HSE06 as a starting point including spin-orbit coupling, a truncated Coulomb potential, and an analytical treatment of the singularity at $q=0$.


I. INTRODUCTION
The isolation of graphene in 2004 can be regarded as a milestone in materials science that triggered a novel field of research, namely atomically thin 2D materials [1].Compared to their 3D counterparts, 2D materials have a higher surface-to-volume ratio, making them ideal candidates for catalysts and sensors [2,3].Due to the confinement of electrons, holes, phonons, and photons in the 2D plane, the electronic, thermal, and optical properties of 2D materials present unusual features not found in their 3D counterparts [4][5][6].For instance, their electronic structure -especially band gaps -can be easily adjusted by acting on the vertical quantum confinement through e.g., the number of atomic layers, or external perturbations, such as an external electric field, and strain [7,8].The sensitivity to strain, i.e., to structural details, implies that 2D materials also exhibit strong electron-phonon coupling [8].In addition, exciton binding energies are significantly larger than in 3D materials, and they can be tuned by the dielectric environment, e.g., by encapsulation or deposition on substrates [9][10][11].All these characteristics make them outstanding components in novel applications for electronics and optoelectronics [12][13][14][15][16][17].
For a deep understanding of 2D materials, an accurate description of their band structure is a must.Many-body perturbation theory within the GW approach has become the stateof-the-art for ab initio electronic-structure calculations of materials.In this sense, many studies have employed GW to investigate the electronic properties of 2D materials .
Surprisingly, as illustrated in Fig. 1 for monolayer MoS 2 , they show a wide dispersion in the fundamental band gap.The same has been found for a number of 2D materials that have been extensively studied in the last years.Results for MoS 2 , MoSe 2 , MoTe 2 , WS 2 , WSe 2 , BN, and phosphorene are summarized in the Appendix.In the extreme cases of MoS 2 , WS 2 , WSe 2 , and BN, the calculated band gaps are scattered between 2.31-2.97,2.43-3.19,1.70-2.89,and 6.00-7.74eV, respectively; in the worst case, the deviation (ratio between largest and smallest values) is as much as 61% .Moreover, for some materials, such as for MoS 2 , MoTe 2 , WS 2 , WSe 2 , and BN, not even the gap character is uniquely obtained -being direct or indirect, depending on the details of the calculation.@LDA @LDA, SOC @LDA, trunc.@LDA, SOC, trunc.@PBE @PBE, SOC @PBE, trunc @PBE, SOC, trunc @HSE @HSE, SOC FIG. 1. Literature results for the G 0 W 0 energy gap of MoS 2 at the K point of the Brillouin zone as a function of the in-plane lattice parameter a.They are obtained with and without SOC (filled and open symbols, respectively), as well as with and without truncation of the Coulomb potential (triangles and circles, respectively), using LDA, PBE, and HSE as starting points (green, orange and blue, respectively).Table VIII summarizes these data.
Many factors contribute to this confusing and unsatisfactory situation: • In various works, different geometries have been adopted.In this context, it must be said that the properties of 2D materials are highly sensitive to structural parameters [18,19,66].Small changes in the lattice constant a already have a large impact on the energy gap, as seen in Fig. 1 and Tables VIII-XIV.Moreover, often the lattice parameter alone is not sufficient to unambiguously determine the structure of a 2D material.For instance, phosphorene is characterized by four structural parameters; transition metal dichalcogenides require, besides a, the distance between two chalcogens (for MoS 2 , d SS as depicted in Fig. 2).These "other" structural parameters have a notable effect on the electronic properties as well [66,67].Unfortunately, in several studies, only the lattice parameter is reported, which prevents not only a fair comparison between published results but also reproducibility.
• A second reason can be attributed to the various ways of performing G 0 W 0 calculations.First, there is the well-known starting problem [68][69][70][71][72][73][74][75][76].Then, especially for 2D materials, G 0 W 0 energy gaps converge very slowly with respect to the vacuum thickness and the number of k-points.Even slabs with a vacuum layer of 60 Å together with a 33 × 33 × 1 k-grid have been shown to be insufficient to obtain fully converged results [28].However, by truncating the Coulomb potential, convergence can be achieved with a reasonable amount of vacuum [49].The number of k-points can be drastically reduced by an analytic treatment of the q = 0 singularity of the dielectric screening or by using nonuniform k-grids [24,38].
• Last but not least, also spin-orbit coupling (SOC) plays an important role in many cases.Besides decreasing the size of the fundamental gap mainly through a splitting of the valence band [77], in some 2D materials and for certain geometries and methods, disregarding or including this effect may change its character from indirect to direct or vice versa [66].
In this manuscript, we address all these issues and provide a benchmark data set of DFT and G 0 W 0 calculations, taking monolayer MoS 2 as a representative 2D material.Due to its unique properties, it can be considered the most important 2D material after graphene.MoS 2 exhibits high electron mobility [14,78]; moderate SOC that can be exploited in spinand valleytronics [79][80][81][82][83][84][85]; a direct fundamental band gap with intermediately strong exciton binding, which is suitable for (opto)electronic devices operating at room temperature [14,29,31,38,78,86].For these reasons, there are many experimental and theoretical works in the literature that investigate MoS 2 , allowing for a better comparison with our results.
We employ the linearized augmented planewave + local orbital (in short LAPW+LO) method as implemented in the exciting code.LAPW+LO is known to achieve ultimate precision for solving the Kohn-Sham (KS) equations of density functional theory (DFT) [87] and high-level GW results [88].Besides the local and semilocal DFT functionals LDA [89] and PBE [90,91] respectively, we include HSE06 [92][93][94] both for geometry optimization and as a starting point for G 0 W 0 .So far, HSE06 has not often been used for such calculations of MoS 2 [20,22,29,33], and to the best of our knowledge, neither a Coulomb truncation nor an adequate treatment of the singularity at q = 0 was applied.For brevity, hereafter, we will refer to HSE06 as HSE.In our G 0 W 0 calculations, we truncate the Coulomb potential [28,49,95], and apply a special analytical treatment for the q = 0 singularity [24].Moreover, we investigate the role of SOC at all levels.We carefully evaluate the impact of all these elements and conclude what leads to the most reliable electronic structure of this important material.Besides a detailed analysis of energy gaps, we address effective masses and spinorbit splittings.

A. Linearized augmented planewave + local orbital methods
The full-potential all-electron code exciting [96] implements the family of LAPW+LO methods.In this framework, the unit cell is divided into non-overlapping muffin-tin (MT) spheres, centered at the atomic positions, and the interstitial region in between the spheres.exciting treats all electrons in a calculation by distinguishing between core and valence/semi-core states.For core electrons, assumed to be confined inside the respective MT sphere, the KS potential is employed to solve the Dirac equation which captures relativistic effects including SOC.The KS wavefunctions |Ψ nk ⟩ for valence and semicore states, characterized by band index n and wavevector k, are expanded in terms of augmented planewaves, |ϕ G+k ⟩, and local orbitals, |ϕ γ ⟩, as |ϕ G+k ⟩ are constructed by augmenting each planewave with reciprocal lattice vector G, living in the interstitial region, by a linear combination of atomic-like functions inside the MT spheres.In contrast, the LOs |ϕ γ ⟩ are non-zero only inside a specific MT sphere.They are used for reducing the linearization error [96,97], for the description of semicore states, as well as for improving the basis set for unoccupied states [87,88].The quality of the basis set can be systematically improved by increasing the number of augmented planewaves (controlled in exciting by the dimensionless parameter rgkmax) and by introducing more LOs [87,96].With all these features, exciting can be regarded as a reference code not only for solving the KS equation [98], where it is capable of reaching micro-Hartree precision [87], but also for G 0 W 0 calculations [88].
In the G 0 W 0 approximation, one takes a set of KS eigenvalues {ε nk } and eigenfunctions {Ψ nk } as a reference, and evaluates first-order quasi-particle (QP) corrections to the KS eigenvalues in first-order perturbation theory as where Z nk is the renormalization factor, V xc is the exchange-correlation (xc) potential, and Σ is the self-energy.The latter is given as the convolution with G being the single-particle Green function and W the screened Coulomb potential.
In this non-selfconsistent method, the quality of the QP eigenvalues may depend critically on the starting point.In many cases, LDA and PBE have been proven to be good starting points for G 0 W 0 , leading to QP energies that agree well with experiments [73,99,100].

C. Coulomb truncation
In calculations with periodic boundary conditions, for treating 2D systems, a sufficient amount of vacuum is required to avoid spurious interaction between the replica along the outof-plane direction.Local and semilocal density functionals have a (nonphysical) asymptotic decay much faster than 1/r, facilitating convergence of unoccupied states -and thus KS gaps-with respect to the vacuum size.In G 0 W 0 , the 1/r decay of the self-energy complicates this task.Specifically for MoS 2 , even a vacuum layer with 60 Å thickness turned out not being sufficient to obtain a fully converged band gap [28,49].Truncating the Coulomb potential along the out-of-plane direction z, however, leads to well-converged G 0 W 0 band gaps with a considerably smaller vacuum size [28,49].
Here, we truncate the Coulomb potential with a step function along z.Setting the cutoff length to L/2, where L is the size of the supercell along z (Fig. 2), the truncated Coulomb potential can be written in a planewave basis as [95]: where Q = q + G, and q is a vector in the first Brillouin zone (BZ).
D. Treatment of the q = 0 singularity On the down-side, truncating the Coulomb interaction slows down the convergence in terms of k-points because of the non-smooth behavior of the dielectric function around the singularity at q = 0 [23, 24,41].To bypass this problem, we follow an analytical treatment of W c , the correlation part of W , close to the singularity as proposed in Ref. [24].Without this treatment, the correlation part of the self-energy Σ c (ω) can be written as [107] ⟨Ψ where εnk = ϵ nk + i η sign(E F − ϵ nk ) and E F the Fermi energy.M i nm (k, q) are the expansion coefficients of the mixed-product basis, an auxiliary basis to represent products of KS wavefunctions.Like LAPWs and LOs, they have distinct characteristics in the MT spheres and interstitial region [107,108].To treat the q = 0 case separately, the corresponding term in Eq. ( 5) is replaced by where Ω 0 is a small region around q = 0. Analytical expressions for W c ij (q, ω ′ ) in the limit q → 0 [24] are then employed to calculate the integral in Eq. ( 6).

E. Spin-orbit coupling
In this study, SOC is treated via the second-variational (SV) scheme [109].In LDA and PBE calculations, the conventional SV approach is utilized [96,[110][111][112]: First, the scalarrelativistic problem, i.e., omitting SOC, is solved.A subset of the resulting eigenvectors is then used as basis set for addressing the full problem.The number of eigenvectors is a convergence parameter.In this work, the SOC term is evaluated with the zero-order regular approximation (ZORA) [113,114].
A ground-state calculation with HSE is performed via a nested loop [115]: In the outer loop, the nonlocal exchange is computed for a subset of KS wavefunctions, using a mixedproduct basis.The inner loop solves the generalized KS equations self-consistently, where only the local part of the effective potential is updated in each step.Within this inner loop, the SV scheme is applied self-consistently to incorporate SOC.The corresponding term, evaluated within the ZORA, is based on PBE, which is justified by the minimal contribution due to the gradient of the nonlocal potential [116][117][118].
G 0 W 0 calculations with SOC are performed in two steps on top of ground-state calculations that include SOC.First, the QP energies are computed as explained in Sec.II B, using the scalar-relativistic KS eigenvalues and eigenvectors.In the second step, the obtained QP energies are used, together with the SV KS eigenvectors, to evaluate SOC through the diagonalization of the SV Hamiltonian.This approach is sufficient in the case of MoS 2 since SOC does not cause band inversion [108,118].

III. COMPUTATIONAL DETAILS
In our calculations, we employ the all-electron full-potential code exciting [96].The only exception is for obtaining the HSE equilibrium geometry, where FHI-aims [119,120] is used, since so far exciting lacks geometry relaxation with hybrid functionals.Even though exciting and FHI-aims implement very different basis sets to expand the KS wavefunctions, the two codes have been shown to be among those with the best mutual agreement [98].

Moreover, a comparison of energy gaps and geometry relaxations for MoS 2 confirms this finding (see Appendix A).
In all calculations, the unit-cell size L along the out-of-plane direction (Fig. 2) is set to 30 bohr.Different flavors of xc functionals are applied, namely LDA, PBE, and HSE.
In the latter, we use the typical parameters [94], i.e., a mixing factor of α = 0.25 and a screening range of ω = 0.11 bohr −1 .To determine the respective equilibrium geometries, we relax the atomic positions until the total force on each ion is smaller than 10 µHa/bohr.
For these geometries, the electronic structure is calculated with all three functionals, with and without SOC, giving rise to a set of 18 calculations.These calculations are followed by G 0 W 0 calculations, taking the respective DFT solutions as starting points.
The dimensionless parameter rgkmax that controls the exciting basis-set size is set to 8.
In the LDA and PBE calculations, we use a 30×30×1 k-grid.In HSE and G 0 W 0 calculations, we employ 400 empty states and an 18 × 18 × 1 k-grid.In G 0 W 0 , the correlation part of the self-energy is evaluated with 32 frequency points along the imaginary axis, and then analytic continuation to the real axis is carried out by means of Pade's approximant.For the bare Coulomb potential, we use a 2D cutoff [95] combined with a special treatment of the q = 0 singularity as described in Section II D. Furthermore, we carefully determine the minimal set of LOs, sufficient to converge at least the lowest 400 KS states.This is achieved with 2 and 6 LOs for sulfur s and p states, respectively, as well as 3, 6, and 10 LOs for molybdenum s, p, and d states, respectively.Overall, we estimate a numerical precision of 50-100 meV in the energy gaps obtained with our G 0 W 0 calculations.To determine effective masses, we use parabolic fits within a range of 0.05 Å−1 around the valence-band maximum (VBM) and the conduction-band minimum (CBm) at the K point of the BZ.As will be shown further below, depending on the respective case, K can host either global or local extrema.

A. Ground-state geometries
The geometry of MoS 2 , depicted in Fig. 2, is determined by the in-plane lattice parameter a and the distance between sulfur atoms, d SS .In Table I

B. Electronic structure
Table II summarizes the energy gaps obtained with different functionals for the different geometries.We consider here the direct gap at the K point (E g (KK)) as well as the indirect gaps between Γ and K (E g (ΓK)) and between K and T (E g (KT)).For the definition of the T point, see Fig. 2. For each geometry and methodology, the bold font highlights the fundamental gap.For the calculations that include SOC, Fig. 3 displays the energy gaps given in Table II with respect to the lattice parameter.In the DFT calculations, regardless of the calculation method, E g (KT) (squares) shows a weak dependence on the geometry.
The fundamental gap obtained with LDA, PBE, and HSE, is always direct at K. In G 0 W 0 , the fundamental gap is E g (KT) for the structures with smaller lattice parameter andE g (KK) for larger lattice constants.
Experimental results for bulk MoS 2 are also provided (c [122]).For the definition of d SS , d MoS and θ, see Fig. 2. two functionals also agree on the location of the VBM and the CBm.For both, the band gap is direct if SOC is included, and indirect otherwise for the PBE and HSE geometries.

This work
In contrast, HSEgives a direct gap, regardless of whether SOC is considered.
Also G 0 W 0 @LDA and G 0 W 0 @PBE are very close to each other, the largest difference being 0.04 eV.This can be attributed to the similarity between LDA and PBE when the same geometry is adopted.However, the similarity between G 0 W 0 @LDA and G 0 W 0 @PBE results is material dependent, as observed in other works [126][127][128][129].
For a given geometry, the locations of the VBM and the CBm are independent of the starting point, with the only exception being the HSE geometry when SOC is included.
When comparing the three geometries, we encounter three different scenarios.First, for the LDA geometry, the fundamental gap changes from a direct KS gap at K to an indirect QP gap (between K and T), independent of the starting point.Second, for the PBE and HSE geometries, when SOC is disregarded, the indirect gap ΓK obtained with LDA and PBE becomes direct and located at K upon applying G 0 W 0 .Third, for the HSE geometry, and SOC being included we observe an indirect band gap for G 0 W 0 @LDA while it is direct for PBE and HSE as starting points.This can be understood in terms of the small differences or not (N).The experimental gap is 2.6 eV (direct at K) [86].Note that this value is corrected for the zero-point renormalization energy of 75 meV [125].between the KK and KT gaps, ∆ KT , which are 0.01 eV, 0.05 eV, and 0.15 eV for G 0 W 0 @LDA, G 0 W 0 @PBE and G 0 W 0 @HSE, respectively.Including SOC, splits the conduction band state at T (K) by ∼ 0.07 eV (∼ 3 meV), decreasing ∆ KT by ∼ 0.03 eV.This is enough to make ∆ KT negative for G 0 W 0 @LDA, but not for G 0 W 0 @PBE, and G 0 W 0 @HSE.A more detailed discussion about ∆ KT can be found in Section IV D 1.
Figure 4 shows the band structures obtained for the HSE geometry including SOC.As expected, the differences between G 0 W 0 @HSE and HSE bands are less pronounced than [eV] DFT those between G 0 W 0 @LDA (G 0 W 0 @PBE) and LDA (PBE) bands.For all three starting points, the G 0 W 0 corrections are not uniform over all k-points, i.e., a simple scissors approximation is, strictly speaking, not applicable.We will explore this in a more quantitative fashion further below.The SOC splitting in the valence band is zero at the Γ point and increases toward the K point.The SOC effect on the conduction bands is much less pronounced.These observations are in agreement with other theoretical [25,29,[130][131][132][133] and experimental works [86,134].
The impact of the self-energy correction on the energy gaps, shown in Fig. 5 for the case when SOC is included.Clearly, ∆E g is more significant for (semi)local DFT starting points than for HSE.Interestingly, for a given starting point, ∆E g hardly depends on the geometry.For LDA and PBE as the starting points, the ranges of ∆E g (KK), ∆E g (ΓK), and ∆E g (KT) are 0.84-0.90,1.0-1.1, and 0.61-0.65 eV, respectively.
The dependence is even weaker for G 0 W 0 @HSE with values of 0.64-0.66,0.74-0.78,and 0.39-0.40eV, respectively.Very similar results are observed when SOC is disregarded.

C. Spin-orbit splittings and effective masses
In Table III, we report for the HSE structure the spin-orbit splitting ∆ val (∆ cond ) at the K point for the highest occupied (lowest unoccupied) band.Effective electron (hole) masses m e (m h ) calculated at the K point along different directions are shown as well.We observe that neither the spin-orbit splittings nor the effective masses are very sensitive to the geometry (see Table XV for more details).

LDA PBE HSE
The effective hole mass, m h , exhibits minor variations, not larger than 0.04m 0 , when going from KΓ to KM.This is in line with other calculations [27,142].The measured value for freestanding MoS 2 is m h = 0.43±0.02[141] and, apart from LDA and PBE, all theoretical results show excellent agreement.The electron mass, m e , is more isotropic than m h .For LDA, it differs by at most 0.02m 0 between KΓ and KM.Our values are in line with other calculated results [18,21,27,34,131,142,143]. The measured counterpart of 0.67 ± 0.08 [140], is significantly larger than the calculated value reported here and in other theoretical works [18,21,27,34,131,142,143]. As discussed in Ref. [34], the difference could originate from the heavy doping of the measured sample, which may introduce metallic screening.

Comparison with experiment
The experimental band gap for free-standing MoS 2 , determined by photocurrent spectroscopy, is 2.5 eV [86].In order to compare with experiment, it is important to account for the zero-point renormalization energy of 75 meV [125].This means that the theoretical value computed without this correction should be 2.575 ∼ = 2.6 eV to match its measured counterpart.For our discussion here, we consider the HSE and PBE geometries, which are closer to experiment than that obtained by LDA.For the following assessment, we refer to the values in Table II.For the structure optimized with HSE, G 0 W 0 performed on LDA, PBE, and HSE as starting points gives E g (KK) of 2.52, 2.52, and 2.76 eV, respectively, i.e., G 0 W 0 @LDA and G 0 W 0 @PBE underestimate the measured value by about 0.08 eV, whereas G 0 W 0 @HSE overestimates it by 0.16 eV.However, even though G 0 W 0 @LDA agrees better with experiment than G 0 W 0 @HSE, it erroneously predicts an indirect band gap E g (KT) which is 0.02 eV smaller than E g (KK).G 0 W 0 @PBE shows the best agreement with experiment and also predicts the gap to be direct.Interestingly, considering the PBE geometry, as done in several works [19, 21-24, 28, 31-33, 35, 40-43, 143], G 0 W 0 @HSE, giving a direct band gap of E g (KK) = 2.68 eV, agrees best with experiment, the deviation being 0.08 eV only.With LDA and PBE as starting points, the calculated G 0 W 0 band gap is direct, but 0.16 and 0.15 eV, respectively, below the experimental value.(a [131,140]; b [131,144]).All values are given in eV.Other relevant aspects of the band structure concern relative energy differences, in particular ∆ KT (introduced above) as well as the maximum energy at the Γ point wrt the VBM at the K point, ∆ ΓK = E g (KK) − E g (ΓK) (see Fig. 4).Experimentally, ∆ KT is expected to be ≳ 60 meV [131,140] and ∆ ΓK ≈ 140 meV [131,144].Taking our calculations with SOC at the HSE geometry, (see Table IV), the G 0 W 0 @HSE value of 0.12 eV reproduces ∆ KT best.
On the other hand, HSE satisfies ∆ ΓK best with a value of 0.12 eV (0.02 eV smaller than in experiment), while the values obtained with the other methods differ from experiment by 0.09 eV (G 0 W 0 @HSE), -0.09 eV (LDA and PBE), 0.10 eV (G 0 W 0 @PBE), and 0.11 eV (G 0 W 0 @LDA), respectively.At the PBE geometry, G 0 W 0 @HSE and G 0 W 0 @PBE give the same value for ∆ ΓK .With an overestimation of 0.07 eV, it is closer to experiment than the value at the HSE geometry.Again, HSE is the only starting point for which G 0 W 0 predicts ∆ KT in agreement with experiment.
In summary, considering the band gap as well as the energy differences ∆ KT and ∆ ΓK , we conclude that at the PBE geometry, G 0 W 0 @HSE including SOC shows the best overall agreement with experimental data.Also for the HSE structure, HSE is the best starting point, with results that are overall only slightly worse.Overall, G 0 W 0 @HSE at the HSE geometry can be considered more appropriate, since only one xc functional is needed for providing decent results for both, the geometry and the electronic properties and thus the most consistent picture.Also for other materials, HSE has been found to be a superior G 0 W 0 starting point [76,[145][146][147] compared to LDA and PBE.For such materials with intermediate band gaps, [71,72,[101][102][103][104][105][106], this functional better justifies the perturbative self-energy correction [68,73,145].Figure 5 confirms this for MoS 2 .

Comparison with theoretical works
By employing coupled-cluster calculations including singles and doubles (CCSD) [148,149], Pulkin et al. obtained energy gaps of E g (ΓK) = 2.93 eV and E g (KK) = 3.00 eV with an error bar of ±0.05 eV [143] for the PBE geometry of Ref. [133]; SOC was not included.
For our PBE geometry and also omitting SOC, the G 0 W 0 @HSE results are the ones closest to these values, with E g (ΓK) differing by 0.04 eV and E g (KK) by 0.24 eV.
For a fair comparison with other published G 0 W 0 values with LDA and PBE as starting points, we restrict ourselves here to results obtained by using a Coulomb truncation in combination with either a special treatment of the q = 0 singularity or a nonuniform k-grid sampling, since these methods ensure well converged gaps.In Ref. [24], disregarding SOC and adopting the PBE geometry (a = 3.184 Å, d SS = 3.127 Å), a direct band gap of 2.54 eV was reported for G 0 W 0 @PBE which is very close to ours (E g (KK) = 2.52 eV), i.e., differing by only 0.02 eV.Including SOC and the thus slightly changed PBE geometry (a = 3.18 Å, d SS = 3.13 Å), a G 0 W 0 @LDA value of 2.48 eV was obtained in Ref. [23]; at basically the same geometry (differences in the order of 10 −3 Å), our results of E g (KK) = 2.44 eV is only 0.04 eV smaller.
For a lattice parameter of 3.15 Å, Rasmussen et al. calculated a G 0 W 0 @PBE band gap of 2.64 eV without SOC [24].For the same lattice constant, but including SOC, Qiu et al. reported a G 0 W 0 @LDA band gap of [38] of E g (KK) = 2.59 eV with the plasmon-pole model and E g (KK) = 2.45 eV with the contour deformation method.In our case, the structure optimized with HSE (a = 3.160 Å) is closest to a = 3.15 Å.For this structure, without including SOC, we compute a G 0 W 0 @PBE band gap of E g (KK) = 2.60 eV, which agrees quite well with the one by [24], differing by less than 0.04 eV.When we include SOC, we obtain E g (KK) = 2.52 eV with G 0 W 0 @LDA, although at this geometry, we obtain an indirect gap that is 18 meV smaller than E g (KK).This is in good agreement with Ref. [38], with a difference of 0.07 meV only.
For the experimental geometry and neglecting SOC, Hüser et al. [28] reported values of E g (KT) = 2.58 eV and E g (KK) = 2.77 eV for G 0 W 0 @LDA.In our case, at the HSE geometry, we obtain E g (KT) = 2.61 eV and E g (KK) = 2.60 eV.As the HSE geometry is close to experiment, we may attribute the discrepancies mainly to the different underlying KS states.Indeed, at the LDA level, the energy gaps in Ref. [28] are E g (KK) = 1.77 eV and E g (ΓK) = 1.83 eV [28], while ours are E g (KK) = 1.73 eV and E g (ΓK) = 1.71 eV.The values for ∆E g , however, compare fairly well (∆E g (KK) = 1.00 eV, ∆E g (KT) =0.6-0.7 eV in Ref. [28], compared to ∆E g (KK) = 0.87 eV, ∆E g (KT) = 0.63 eV in the present work).
When it comes to G 0 W 0 @HSE, there are only a few results for MoS 2 in the literature, neither obtained with Coulomb truncation nor by any special treatment of the q = 0 singularity.For MoS 2 , these two aspects lead to opposite effects, competing with each other when converging band gaps with respect to the vacuum size and the number of k-points [24]: Neglecting them, band gaps increase when the vacuum layer is enlarged, whereas denser k-grids make them decrease.Hence, due to fortunate error cancellation, an insufficient vacuum length combined with a coarse k-grid may lead to G 0 W 0 band gaps that agree well with those obtained in a highly converged situation [24,29].In Ref. [33], using the PBE geometry and taking SOC into account, a KK gap of 2.66 eV was reported.With 15 Å of vacuum, a 6 × 6 × 1 k-grid, and adopting the PBE geometry, in Ref. [22], band gaps of 2.05 and 2.82 eV at the HSE and G 0 W 0 levels, respectively, have been obtained.The HSE band gap agrees quite well with ours (2.04 eV), whereas our G 0 W 0 @HSE06 gap is 0.14 eV smaller.The band gap of 2.72 eV reported in Ref. [29] is based on the experimental lattice parameter of 3.16 Å, a 12 × 12 × 1 k-points grid, and a vacuum layer of 17 Å, and includes SOC effects.The authors state to have chosen these settings to take advantage of error cancellation in the band gap [29], and, surprisingly, our band gap of 2.76 eV obtained with G 0 W 0 @HSE at the HSE geometry (a = 3.160 Å) agrees quite well.

V. CONCLUSIONS
We have employed the LAPW+LO method to provide a set of benchmark G 0 W 0 calculations of the electronic structure of two-dimensional MoS 2 .We have addressed the impact of geometry, SOC, and DFT starting point on the energy gaps, spin-orbit splittings, and effective masses.We find that the self-energy corrections to the band gaps hardly depend on the adopted geometry.As could be expected, employing LDA and PBE as starting points does not make a significant difference when the same structure is used.The best agreement with experimental results is achieved by G 0 W 0 @HSE at either the HSE or PBE geometry, considering SOC.The spin-splittings obtained with all methods agree well with experimental results.This also holds true for the effective hole mass, using either HSE or G 0 W 0 on top of any of the considered starting points (LDA, PBE, and HSE).In line with other theoretical works, we highlight the importance of a Coulomb truncation and an adequate treatment of the Coulomb singularity around q = 0 as being fundamental for high-quality calculations.
Our findings are expected to be valid for other two-dimensional materials as well.    is provided in the respective reference, only the functional used for relaxation is shown.
SP stands for the G 0 W 0 starting point, KS refers to the specific pseudopotential method (NC means norm-conserving, PAW stands for the projector augmented wave method).
FI indicates the method for frequency integration, where AC stands for imaginary frequencies and analytic continuation; CD for the contour-deformation approach; PP for plasmon-pole model; RA for full-frequency integration along the real axis.The symbols * , ⋄ , and , we list these structural parameters as obtained with LDA, PBE, and HSE, and include the Mo-S bond length d MoS and the angle θ between Mo and S atoms as well.As expected, LDA underestimates the lattice spacing, PBE slightly overestimates it, and HSE shows the best performance with respect to experiment.All three xc functionals underestimate the S-S bond length, PBE being

FIG. 2 .
FIG. 2. Top view (left) and side view (middle) of the MoS 2 slab geometry, determined by the inplane lattice constant a and the distance between sulfur atoms, d SS .The Mo-S bond length d MoS and the S-Mo-S angle θ are shown as well.L is the unit-cell size along the out-of-plane direction z, including a vacuum layer.High-symmetry points and paths used in the band-structure plots (Fig. 4) are highlighted in the BZ (right panel).

4 FIG. 3 .
FIG. 3. Calculated energy gaps of MoS 2 as a function of the lattice parameter a.The values for KK, ΓK, and KT are represented as circles, triangles, and squares, respectively.Dotted (solid) lines stand for DFT (G 0 W 0 ) results; all include SOC.Note that these values only indirectly reflect the S-S distance d SS .

FIG. 4 .FIG. 5 .
FIG. 4. Band structures including SOC obtained for the HSE geometry.The parameters ∆ KT and ∆ ΓK are discussed in Section IV D 1.

ACKNOWLEDGMENTSFIG. 6 .
FIG. 6. Dependence of d SS , E g (KK), and E g (ΓK) on the lattice constant a of MoS 2 obtained with exciting (straight lines) and FHI-aims (dots) with different xc functionals.

TABLE II .
Energy gaps (in eV) obtained for the 18 different cases considered.For each case, the fundamental gap is highlighted in bold.The second column indicates whether SOC is included (Y)

TABLE IV .
∆ KT and ∆ ΓK from G 0 W 0 calculations including SOC, compared to measured values

TABLE V .
Maximum absolute difference in d SS , E g (KK), and E g (ΓK) between results obtained with exciting and FHI-aims.Figure 6 depicts for LDA, PBE, and HSE the dependence of d SS , E g (KK), and E g (ΓK) on the lattice parameter a.For LDA and PBE, comparison of exciting (straight lines)and FHI-aims (dots) reveals excellent agreement.TableVpresents the maximum absolute differences of these quantities which are smaller than 1.0 × 10 −3 Å for d SS and smaller than 5.2 meV for the energy gaps.From this comparison we expect a similar level of agreement for HSE.

TABLE VI .
Convergence behavior of the energy gaps E g (KK) and E g (ΓK) with the k-grid (n k × n k × 1) with and without using the special treatment of the q = 0 singularity (indicated by the Y and N in the first column).

TABLE VII .
Extrapolated G 0 W 0 energy gaps depending on the range [n k,i , n k,f ] of the data used in the fit (B1) with and without using the special treatment of the q = 0 singularity (indicated by the Y and N in the first column).The fitting coefficient E g (∞), is shown in the table as a function of the used data range [n k,i , n k,f ].The results obtained with the analytic treatment are overall less sensitive to this range, rapidly converging to 2.50 eV and 2.65 eV for KK and ΓK, respectively.In contrast, not using the analytic scheme, shows slow convergence to 2.54 eV and 2.69 eV, respectively.For the analytic treatment, the energy gaps from the 18 × 18 × 1 mesh deviate from E g (∞) by only 20-30 meV.Without the special treatment, even a k-grid of 30 × 30 × 1 leads to errors of 130-140 meV.To obtain a precision of 20-30 meV in this case, we estimate that k-grids in the range 66 × 66 × 1 to 84 × 84 × 1 may be necessary.

TABLE VIII :
G 0 W 0 energy gaps from the literature, indicating the used supercell geometry, the method / approximation, and computational parameters.If no value for d SS • mean, respectively, analytical treatment of the q = 0 singularity, nonuniform neck subsampling, and Sternheimer G 0 W 0 .The symbol † indicates the use of Coulomb truncation; ⊕ refers to extrapolation of the gap assuming a 1/L behavior with the vacuum size.

TABLE IX :
Same as Table VIII, but for MoSe 2 .d SeSe is the distance between two Se atoms, analogous to d SS .

TABLE X :
Same as Table VIII but for MoTe 2 .d TeTe is the distance between two Se atoms, analogous to d SS .

TABLE XI :
Same as Table VIII but for WS 2 .

TABLE XII :
Same as Table VIII but for WSe 2 .

TABLE XIII :
Same as Table VIII but for BN.G stands for atom-centered Gaussian orbitals and US for ultra-soft pseudopotentials.

TABLE XIV :
Same as Table VIII but for phosphorene.aand b are the two in-plane lattice parameters; d refers to the layer thickness; and θ is the in-plane P-P-P angle.Appendix D: Spin-orbit splittings and effective masses for different geometries g [eV] Ref. a × b [ Å2 ] d [ Å] θ [°] L [ Å] SP KS SOC FI k-grid Misc.ΓΓ Table XV provides spin-orbit splittings in the valence and conduction band at the K point for the LDA, PBE, and HSE geometries together with effective electron and hole masses at the K point along the directions KΓ and KM.∆ val and ∆ cond are not very sensitive to the geometry in the sense that, for a given methodology, ∆ val (∆ cond ) varies less than 5 meV (1 meV) between the different geometries.A similar trend is also observed for the effective masses.

TABLE XV .
Spin-orbit splittings in the valence (conduction) band ∆ val (∆ cond ), in meV, evaluated at the K for the various geometries considered in this work as well as effective electron (hole) masses m h (m e ), in units of m 0 , at the K point along different directions (in parentheses).GeometryLDA PBE HSE G 0 W 0 @LDA G 0 W 0 @PBE G 0 W 0 @HSE