Estimation of the single-particle band gap and exciton binding energy in two dimensional insulators: a modified G0W0-BSE method approach

In this paper we present an alternative G0W0-BSE procedure, suitable for calculation of the quasi-particle and optical properties in 2D semiconductors. The method completely excludes the spurious Coulomb interaction with 2D crystal replicas. The calculated band gap energies of hexagonal boron nitride (hBN), MoS2 and MoTe2 monolayers are in good agreement with other theoretical results. The 2D Bethe–Salpeter equation is derived and reduced to a 2D-hydrogen Schrödinger equation in which enter the G0W0 band gap, DFT effective masses, and RPA screened Coulomb interaction. This formulation is applied to the problems of determining exciton binding energies and estimating the quasiparticle band gap in hBN, as well as in some transition-metal dichalcogenides. A semiclassical procedure is used in the limit of high polarizability λ in order to obtain the analytical expression for exciton binding energies.


Introduction
For the last decade, the condensed matter physics has been dominated by the experimental and theoretical investigations of two-dimensional materials. One class of such materials are the direct gap two-dimensional insulators. The best known examples are monolayers (ML) of the members of transition-metal dichalcogenides (TMDs) [1,2] and the hexagonal boron nitride (hBN) monolayer. Many interesting phenomena are attributed to TMDs and hBN. Heterostructures with hBN improve graphene transport properties, enhancing its application in electronics [3,4]. By applying uniaxial or biaxial strain, the optical properties of single-layer TMDs can be tuned [5,6]. Single-layer TMD (MoS 2 or WS 2 ) transistors are synthesized and tested for application in digital electronics [7][8][9]. Finally, TMDs can be potentially applied in optoelectronics as photodetectors [10] or as photodiodes and photovoltaic devices [11][12][13].
Although hBN and TMDs appear as simple honeycomb crystals with just two or three atoms per unit cell, theoretical calculations of their electronic structures and optical properties tend to be complex. Specifically, the very inefficient screening causes the long-range Coulomb correlations to prevail and the DFT methodology systematically fails in estimating the single particle band gap. In addition, the strong (weakly screened) Coulomb interaction gives rise to excitons with large binding energies.
The well-established ab initio methodology which includes quasiparticle corrections (GW method) and solving the Bethe-Salpeter equation (BSE, the so-called GW-BSE method) is capable of giving a very accurate single particle band gap, excitonic binding energies, and oscillator strengths [14][15][16][17][18][19][20]. The semi-analytic study of optical properties in hBN/substrate systems in the framework of Bethe-Salpeter and Wannier-Keldysh approach were done in reference [21]. The authors provided very extensive investigation of excitons, exciton-polaritons, absorption and spontaneous radiative decay of excitons in hBN on SiO 2 or graphite surfaces. The authors also presented the analytic expressions for hBN/surface absorption spectra and quasi-particle corrections of hBN/surface band gaps. They suggested how their semi-analytic approach

Ab initio studies of hBN-ML quasi-particle and optical properties
In this section the modified q2D-G 0 W 0 method adapted for the calculation of quasiparticle band gaps in 2D crystals is presented. The equation of motion is formulated for the 2D electron-hole operator which is analogous to 2D-BSE. This equation contains q2D-G 0 W 0 corrected band gap and RPA 2D screened Coulomb interaction w. Finally, it is shown that this methodology produces a fairly accurate band-gap and exciton binding energies.

The q2D-G 0 W 0 method
It is well known that the LDA underestimates the semiconducting band gap. Therefore, in order to obtain the accurate exciton energy, providing quasiparticle corrections is of crucial importance. The calculation of the exchange-correlation self-energy within G 0 W 0 aproximation, introduced a long time ago by Hedin, Strinati and Louie [31][32][33][34] and more recently by Louie, Thygesen and others [20,24,35], is applied to calculate various crystal quasiparticle band structures. Here we briefly present some modifications which should be done to the standard G 0 W 0 approximation when used in calculating the band gaps of 2D crystals. The G 0 W 0 exchange-correlation (XC) self-energy corresponding to the Bloch state (n, k) is where the exchange self-energy is and the correlation self-energy is Here f n k is the Fermi-Dirac distribution, G 0 nk (ω) = 1/ ω − E n k + iη sgn(E n k − E F ) is the one-particle Green's function at T = 0 and the charge vertices are The τ = 0 + and η = 0 + are infinitesimally small positive numbers. The induced part of dynamically screened Coulomb interaction (which enters into (2.3)) can be written as The time-ordered response function is the solution of the Dyson equation where the RPA time-ordered irreducible polarizability is explicitly given bỹ . (2.7) The Coulomb interaction kernel is is a 3D momentum transfer, r = (ρ, z) is a 3D position vector, G = (G , G z ) are the reciprocal lattice wave vectors and V is the normalization volume. The equations (2.1)-(2.7) represent the standard G 0 W 0 procedure for the electronic self-energy calculation in 3D crystals. The question now is how to transform these formulas in a way that they can be applied to the calculation of the electron self-energy in an isolated 2D crystal. Suppose that we want to calculate the electron self-energy in the i-th hBN monolayer which is part of periodically stacked non-overlapping hBN layers, as sketched in figure 1(a). The charge density fluctuations ρ nk,mk+q created by electron-electron scattering or by electron-hole transitions |nk → |mk + q in ith hBN layer are represented by the blue dots while the red squares represent the same in the surrounding hBN layers ( j = i). A blue dot can interact with a blue dot in the same layer via Coulomb interaction (2.8), as well as with red squares in the surrounding hBN layers. The corresponding correlation self-energy (2.3) is represented by the Feynman diagrams in figure 1(b). It consists of the process (I) where the electron propagator in ith layer interacts with the polarization in the same layer, and the processes (II-IV) where it interacts with the polarization in the surrounding layers. Since we are interested in how the electronic state in the ith layer is influenced only by the polarization of the same layer, the processes (II-IV) should be neglected and only the process (I) should be retained. However, the process (I) still contains the influence of the surrounding ( j = i) hBN layers. The time-ordered RPA response function χ T , represented by the Feynman diagrams in figure 2(a), consists of processes where the 'bubble' diagrams in the ith layer interact with the 'bubble' diagrams in the same layer (processes (i-iii)), but also of processes where the 'bubble' diagrams in the ith layer interact with the 'bubble' diagrams in the surrounding layers (e.g. process (iv)). Therefore the processes of type (iv) should be neglected and the Dyson equation for χ T takes the form shown by Feynman diagrams figure 2(b). Finally, this results in the complete annulment of the interaction with the surrounding hBN layers. Technically, the processes (II-IV) in figure 1(b) and processes of type (iv) in figure 2(a) can be neglected in a way that the formulation (2.1)-(2.7) is retained, however, the Coulomb interaction kernel (2.8) is modified as follows [44,45] v 3D where Also, the following transformation must be done q, k → Q, K (2.10) where Q = (q x , q y ) and K = (K x , K y ) are the wave-vectors parallel to the x-y plane. The 3D Brillouin zone q summation should be replaced by 2D summation as 11) and the same is true for the k summation in (2.7). Here S is the normalization surface and L is the unit cell constant in the perpendicular (z) direction, as sketched in figure 1(a). An example of how to transform the Dyson equation for 3D or bulk response function (2.6) into the Dyson equation for q2D response function which represents the charge density fluctuations in an isolated 2D crystal is presented in appendix A. There, it is clearly demonstrated that the only modification which should be done in (2.6) is the transformation (2.9). The quasiparticle energy of the Bloch state (n, K) is calculated by solving the equation for ω. Here V XC nK represents the LDA exchange-correlation energy of Bloch state |nK . For non-overlapping layers, in the exchange term (2.2) the charge density fluctuation ρ nk,mk+q always interacts within the same layer through the bare Coulomb interaction. This means that the expression (2.2) can still be used without further modifications. The only difference between a Coulomb kernel (2.9) and a Coulomb truncation (equation 16) in reference [46] or (equation 43) in reference [24] is that in equation (2.9) the Coulomb interaction is 'closed' within z ∈ [−L/2, L/2], i.e. the charge at z = L/2 cannot interact with any other charge, even at z = L/2 + 0 + , while the Coulomb interaction (equation (16) or (43)) is such that the charge at z = L/2 can interact with all of the charges in the region z ∈ [L/2 − R, L/2 + R]. This requires more layers of vacuum if one wants to avoid interaction with neighbouring cells. The detailed description how (2.9) is derived is given in reference [45]. The truncated interaction (2.9) enables that in a periodic superlattice arrangement it is sufficient to choose inter-layer spacing such that electronic densities in adjacent layers do not overlap. Any increase in this distance will not change the result. So the interaction (2.9) allows for supercells with a very few layers of vacuum, which is another important advantage over the truncated interaction proposed in reference [46].
The introduction of the Coulomb kernel (2.9) easily overcomes the most critical issue of the 3D periodicity: how to calculate the long-range image potential at any point z produced by an individual slab, located at −L/2 < z < L/2, which is a part of the supper-lattice. This issue is also solved in reference [47]. Here, in appendix A is demonstrated that if the Dyson equation for the response function χ T (z, z ) is solved in the restricted real space −L/2 < z, z < L/2, the density-density Coulomb interaction remains limited within the boundaries [−L/2, L/2]. This results in χ T (z, z ) representing the dynamical response of an individual slab without any admixture of other slabs. Consequently, for the induced Coulomb interaction the other layers becomes 'invisible', i.e. it properly describes the induced interaction in point z produced by charge in z where z and z are no longer limited to −L/2 < z, z < L/2. Here the bare Coulomb interaction v 2D (Q, z, z ) is given by (A.4). Moreover, the integration of static potential W ind (Q, ω = 0, z, z) over parallel wave vector Q provides the correct image potential at any point z, as presented in detail in reference [48]. A very similar argumentation applies to calculation of the correlation self-energy (2.3) with the induced interaction (2.5), which after transformation v → v 2D , represents exactly the Fourier transform of (2.13).

Implementation of the crystal local field effects in perpendicular z direction
The crystal local field effects CLFE can be explained as dispersivity (or microscopic distribution) of induced potential W ind or induced density ρ ind on the unit cell scale. For example, if a unit charge is placed somewhere in the crystal it will induce a charge density at macroscopic scale but also variatioins of microscopic charge density within the unit cell. If the CLFE are excluded then the variations of microscopic charge density are neglected and the macroscopic ones remain. From the point of view of the reciprocal space the exclusion of CLFE results that in matricesχ 0 , χ T and W ind we retain just the G = G = 0 component. In bulk crystals the exclusion of the CLF is sometimes justified, especially in long wavelength limit Q → 0. However, in the low dimensional crystals (0D, 1D or 2D crystals) the the strong dispersivity (fast variation) of the induced charge density in the direction of the broken periodicity (e.g. in 2D crystals it would be perpendicular to crystal plane) requires inclusion of the CLFE, at least in the direction of the broken periodicity. For these reasons in this calculation we retain the CLFE in the perpendicular z direction and neglect them in the parallel direction. From the point of view of the reciprocal space this means that in all equations (2.1)-(2.9) we put G = G = 0, but we keep the z component different from zero, i.e. G z = 0 and G z = 0. By doing so we capture the important microscopic variations of the induced density in z direction and ignore the not so important microscopic dispersivity in parallel (x-y) directions. This reduces the dimension of matricesχ 0 , χ T and W ind , significantly reducing the memory allocation requirements and substantially accelerating the computation. At the same time the accuracy of the calculation remains satisfactory. The inclusion or exclusion of the CLFE does not affect the definition of the truncated interaction (2.9) and it remains unchanged. For example, if the CLFE are included just in z direction, then G = G = 0 and the irreducible polarisability and the response function matrices becomeχ 0 respectively. This implies that in the matrix Dyson equation for the same form as (2.9) with the exception that G = G = 0.
In the optical limit (Q 1/a, 1/L) a q2D crystal can be considered a 2D sheet crystal. In this case we can average the dynamical response in the z direction In this way the strong dispersivity in z direction counterintuitively resulted in exclusion of the crystal local field effects in the z direction (G z = G z = 0). This works well in the Q 1/a, 1/L limit, but for example in the correlation self energy (2.3) short wavelength contributions (Q ≈ 1/a, 1/L) or the dispersivity of the induced interaction W ind within the unit cell become important and the 2D model is no longer adequate. So when solving the BSE we use the 2D limit (G = G = 0) but when we provide the quasiparticle GW corrections we include the CLFE in the direction in which they are strongest, the z direction, so we put G = G = 0 and leave G z = 0, G z = 0.

Calculation of quasi-particle band gap
At the DFT stage of the calculation the Kohn-Sham (KS) wave functions φ * nk and energy levels E nK , i.e. the band structure of a hBN-ML, are determined using the plane-wave DFT code QUANTUM ESPRESSO (QE) [36]. The core-electron interaction is approximated by the norm-conserving pseudopotentials [37], and the exchange correlation (XC) potential by the LDA Perdew-Zunger (PZ) functional [38]. For the hBN-ML primitive cell constant, a = 4.746a 0 (a 0 is the Bohr radius) is used and the superlattice constant in the z direction is L = 23.73a 0 . The ground state electronic densities of the hBN-ML are calculated using the 12 × 12 × 1 Monkhorst-Pack k-point mesh [39] of the first Brillouin zone (BZ). For the plane-wave cut-off energy we used 60 Ry (816 eV). The RPA time-ordered irreducible polarizability (2.7) is calculated using 201 × 201 × 1k-point mesh sampling. This k-point mesh sampling enables the minimum transfer wave vector q min = 0.0076 a −1 0 . The damping parameter used is η = 50 meV and the temperature is k B T = 10 meV. The band summation is performed over 40 bands, which proved sufficient for the proper description of the electronic excitations up to 30 eV. The crystal local field effects (CLFE) are included only in the perpendicular (z) direction and CLFE energy cut-off is chosen to be 10 Ry (136 eV). The Q-point mesh used in self-energy calculation (2.1)-(2.7) is 21 × 21 × 1. In the exchange energy calculation (2.2) the 3D CLFE are included and CLFE energy cut off is is chosen to be 20 Ry (272 eV).
The hBN-ML is a direct gap insulator with the conduction band minimum and the valence band maximum located at the K point of the Brillouin zone. The DFT band gap obtained in this calculation is Δ DFT = 4.5 eV and after the q2D-G 0 W 0 quasiparticle correction (2.1)-(2.12) it is increased to Δ GW = 7.7 eV. The value of the q2D-G 0 W 0 band gap is slightly overestimated compared to the value of about 7.4 eV proposed in previous calculations [23,24,40,41]. The slight difference may arise from the fact that in the correlation self-energy (2.3) the CLFE are included only in the z direction. On the other hand, it seems that the number of 40 bands is well converged, considering that after using 80 bands the gap was reduced by only 20 meV. However, it should be noted that the number of bands and CLFE energy cut-off are not independent [20,24]. Thus, even if the number of bands is well converged, the low CLFE cut-off may be responsible for the disagreement. Figure 3 shows the hBN-ML band structure obtained using LDA. The red dots represent the q2D-G 0 W 0 correction of highest occupied (or valence) B(2p z ) and the lowest unoccupied (or conduction) N(2p z ) bands in Γ, K and M points of Brillouin zone. It is worth noticing that after G 0 W 0 correction the hBN possess the indirect band gap of 6.2 eV which slightly underestimate the values 6.58 eV reported in reference [24]. However, the q2D-G 0 W 0 correction causes a mostly equal shift upwards of the conduction band which is also in accordance with previous, more detailed calculations [23,24,40,41]. This allows us to use the scissor operator approach, i.e. at K point, the parabolic conduction band will be shifted only relative to parabolic valence band for the value of q2D-G 0 W 0 band-gap. The LDA effective masses of the conduction and the valence bands at K point are m * c = 0.8m e and m * v = 0.75m e , respectively.

Optical properties of q2D
One of the ways to describe the optical properties of a semiconducting system is to solve the 3D Bethe-Salpeter equation (BSE) for the electron-hole propagator. In the leading order of Coulomb interaction the BSE consists of an RPA term (creation and annihilation of electron-hole pairs), a ladder term (interaction between excited electron and hole) and of quasiparticle G 0 W 0 corrections (renormalization of the electron and hole self-energies).
In the very first moment of the optical absorption process, the electron and the hole are created and they propagate for a short time without any scattering. Considering that the photon wavelength is usually much larger than q2D crystal thickness (λ = 1/|Q| L), such pure creation or annihilation processes can be described using an intrinsically 2D model. However, in a very short time after their creation, the electron and the hole being to interact with each other via the screened Coulomb interaction w and scatter with other excitations and impurities in the crystal. The short range and the local character of these scattering processes (they occur for larger wave vectors |Q | ≈ 1/a) requires the microscopic treatment of the screened Coulomb interaction, i.e. its spatial dispersivity in all three dimensions becomes important. Because of that, the DFT-RPA screened Coulomb interaction w, implemented in the G 0 W 0 self-energy and in the ladder term, should be calculated at a high level of accuracy, including a large number of unoccupied bands, as well as the 3D crystal-local field effects. For this reason, the 2D model is adequate at the RPA stage of calculation but inadequate when G 0 W 0 and ladder corrections are included. However, below we still use the 2D approach because the short-range self-energy corrections are already implemented in band structure through the scissor operator approach. Moreover, considering the weak dependence of the ladder term on the number of bands, it is sufficient to include in the calculation only the scattering within one valence (v) and one conduction (c) band. Additionally, in the ladder term we restrict the calculation to only the strongest contribution Q ≈ 0 which finally reduces to a 2D model.

Equation of motion for the electron-hole propagator
The dynamics of charge density fluctuation at wave vector Q can be described by a Heisenberg equation for the charge density operatorˆ (Q) = where we omit writing the spin index. The Hamiltonian of the system is given by and electron-electron interaction is given by Using the Wick theorem [42], the solution of the equation (2.15), to the first order in the Coulomb interaction v q , can be written as which are already improved by the q2D-G 0 W 0 self-energy. The first and the second term on the right-hand side represent RPA and ladder contributions, respectively. In the Q ≈ 0 limit the interband charge vertex ρ vK,cK+Q ∝ Q while the intraband is ρ cK,cK+Q ≈ 1. Hence the RPA term goes linearly with Q and it can be put to zero. Secondly, this is the reason why we can neglect the 'off-resonance' scattering processesˆ cv K ,K +Q in the higher order contribution to the ladder part in equation (2.19). Such 'off-resonance' processes that requires interband v ↔ c transitions are neglected here in accordance with the so called the Tamm-Dancoff approximation [43]. We can also restrict our consideration to transitions which are only in narrow area around the K point of the Brillouin zone. Approximating the valence and the conduction bands by parabolas with effective masses m * v and m * c respectively, the left-hand side in (2.19) becomes is the reduced mass and the band gap is Δ. By applying these approximations and defining the exciton energies Ω relative to the bottom of the conduction band, i.e. Ω = ω − Δ, In the above expression we have put Q = K − K. Main contributions to the sum in (2.21) come from Q ≈ 0 due to the singular behavior of the v Q ∼ 1/Q . In this regime the charge vertices can be approximated as ρ vK+Q ,vK ≈ 1 and equation (2.21) finally gives The particular set of higher order contributions to the Coulomb interaction in (2.22) can be provided by replacing the bare Coulomb interaction by a dynamically screened Coulomb interaction v Q → w Q (ω).
Subsequently, the equation (2.22) is solved in the cases of bare and screened Coulomb potential.

Inserting (2.27) into (2.25) the 2D dynamically screened Coulomb potential becomes
which after a Fourier transformation to the direct space becomes the Keldysh potential [49][50][51] W(r) = e 2 π 2λ Y 0 (r/λ) − N 0 (r/λ) . (2.29) Here Y 0 (x) and N 0 (x) represent Struve and Neumann functions, respectively. Considering the wide hBN-ML band gap, the static approximation is also valid in the dynamic limit. The λ(ω = 0) is derived such that we first calculate the static (ω ≈ 0) 2D dielectric function (2.26) for different wave vectors Q. In the interval Q < 0.2 a.u. it behaves exactly as (2.27) such that the screening length λ(ω = 0) is derived by fitting (2.26) with (2.27) and extracting the slope parameter. Due to large hBN band-gap, the dynamical 2D dielectric function can be approximated as (2.27) for λ(ω = 0) up to ω ≈ 4 eV. This is the assumption we use in a further estimation of the optical properties of hBN. From (2.27), it is obvious that λ(ω = 0) = λ can also be interpreted as a screening length. The DFT screening length of hBN-ML is λ = 12.6a 0 [44]. However, given that λ strongly depends on the band-gap energy (λ decreases as the band-gap increases), in further calculations we use the value obtained using q2D-G 0 W 0 corrected band-gap, which is λ = 10a 0 .

Wannier model
Here we briefly mention the Wannier procedure for solving the equation ( The eigenvalues are well-known 2D hydrogen-like energies [52] Ω n = − μ m e 1Ry (n − 1/2) 2 , (2.32) where each state denoted by the quantum number n is (n − 1) times degenerated (not including the valley degeneracy which brings an extra factor of 2), since = 0, 1, . . . , n − 1. By inserting the reduced mass μ = 0.35m e in (2.32) the hBN-ML ground state exciton energy can be estimated as Ω 1 = −1.4 Ry (−19 eV). Thus obtained binding energy is highly overestimated because in (2.31) a bare Coulomb interaction is used. The binding energies will be more realistic when the interaction is replaced by a screened Coulomb interaction (2.29).

Exciton energies and spatial extent
In this section we present the exciton binding energies obtained by solving the Schrödinger equation (2.31) for the screened potential (2.29) where the screening length λ = 10a 0 is obtained from DFT-RPA calculations. The exciton binding energies given in table 1 are presented in terms of two quantum numbers (n r , ). The radial quantum number n r gives the number of nodes in the radial part of the wave function and is the orbital quantum number. The states denoted by the = 0, 1, 2, . . . are labeled s, p, d, . . . This is an adequate choice since the absolute square of the angular part of the wave function behaves like |Φ(ϕ)| 2 ∼ cos 2 ( ϕ), thus resembling to the 2D projections of the 3D atomic hydrogen orbitals. Moreover, the principal quantum number n = 1 + n r + can be introduced. Then the energy states labeled Ω(n r , ) can be equally labeled Ω n . For example, a state Ω(n r = 1, = 1) is equivalent to the state Ω 3p , etc. All energy states having the same n are given in the same colour in table 1.
The replacement of the bare Coulomb interaction by screened one (even wide band-gap hBN screening is not very efficient) causes strong reduction of the ground state exciton binding energy, from 19 eV (obtained from equation (2.32)) to Ω 1s = −2.53 eV, as can be seen in table 1.
This result are similar to results recently obtained by combined ab initio and tight-binding Wannier calculation [23] although the discrepancies can be easily traced to the different values of the effective masses and especially to larger screening length (λ ≈ 20a 0 ) used in reference [23]. Figure 5 shows the exciton (excitation) energies in hBN-ML how they would appear in, for example, an optical absorption experiment. The excitation energies are calculated using formula ω = Δ + Ω(n r , ) where first ten binding energies (or energy levels) are taken from the table 1. The energy of 1s exciton ω ≈ 5.14 eV agrees very well with energy of q ≈ 0 singlet exciton ω ≈ 5.3 eV obtained by solving full ab initio Bethe-Salpeter equation [17]. Also we obtain excellent agreement with optical absorption [53,54] or photoluminescence (PL) [55] experimental result, which reported ω 1s ≈ 6 eV.
The states with higher have lower energy for the same n. This can be seen from figure 5 in the case of n = 3 series Ω 3d < Ω 3p < Ω 3s . This energy ordering is experimentally observed in the two-photon absorption experiments on tungsten disulphide [56]. The mean exciton radius, defined as the average electron-hole separation in the state ψ n , is calculated as r n = ψ n |r|ψ n . In the ground state r 1s ≈ 6 a 0 which is approximately equal to the unit cell dimension. The mean exciton radius decreases with Ω n . For example, r 2s ≈ 22 a 0 and r 2p ≈ 15 a 0 , while for the highest calculated energy level r 4s ≈ 75 a 0 . In the figure 4 the screened potential W(r) is plotted as a function of the screening length λ. As λ increases, the logarithmic nature of the potential becomes more apparent for small electron-hole separations. This can be seen by taking the two opposite limits of the expression (2.29) where γ is Euler-Mascherion constant. The shape of the above potential suggests the reason why the states within the same shell n differ in energy. The radial part of the wave function of the states with lower is located closer to the origin where it is governed by a weaker logarithmic potential (3.2). Therefore, the binding energy is smaller and the mean radius is larger compared to the states with larger values of which are under the influence of stronger bare Coulomb potential (3.1) and thus have larger binding energy and smaller spatial extension. Table 1. The first ten exciton energy levels in eV. The states within the same shell are given in the same colour: n = 1, 2, 3, 4 in maroon, red, orange and olive green, respectively.  Here, it is also suggestive to see how the ground state exciton binding energy Ω 1s depends on the screening length λ. This dependence is shown in the insert of figure 4. In the limit λ → 0 the Ω 1s (λ) = −19 eV as obtained from equation (2.32). However, with the increase of the screening length, the exciton binding energy rapidly decreases such that for λ GW = 10a 0 it is Ω 1s (λ) = −2.53 eV, as already discussed in table 1. In the opposite limit (λ → ∞) a saturation of the Ω 1s (λ) can be seen. The analytical approximation of Ω 1s (λ → ∞) will be considered in the last section.
Observing figure 5, some conclusions can be made regarding the interaction of excitons with phonons, impurities, and electrons, changing their appearance in the absorption spectrum. Due to these interactions, the exciton signals will be broader and slightly shifted in energy. This causes an overlap between the individual exciton signals that are close in energy, to the point that they can even be joined with the single particle threshold. This would imply the indistinguishability of partial contributions originating from the single particle excitations and the excitons in the optical absorption spectra.

The single particle gap estimation
The results obtained in section 3.1 can be used to estimate the hBN-ML band gap which is experimentally still not observed and ab initio results are not yet in consensus, considering that the values between 7 and 9 eV are reported [14,15]. The band gap Δ can be estimated by finding a value Δ for which the theoretically obtained exciton energy ω 1s is equal to the experimental value ω This equation can be solved numerically by calculating the binding energy Ω 1s as a function of Δ which is then inserted in (3.3). However, here we take into account that the binding energy Ω 1s changes slowly for screening lengths in the range of λ = 10a 0 ± 2a 0 which essentially covers the values of the static screening length calculated for the bandgap between 4.5 eV and 10 eV, as shown in figure 4 insert. Hence, for a reliable estimate we can choose the screening length value λ = 10a 0 obtained from the G 0 W 0 band structure. Moreover, the same value for exciton binding energy Ω 1s (λ) = −2.53 eV can be used. In optical absorption or PL experiments on quasi hBN-ML [53][54][55] the exciton ground state signal appears at the energy ω exp 1s ≈ 6 eV. Therefore using this data and (3.3) the quasiparticle band gap is estimated as Δ ≈ 8.5 eV, which overestimate the q2D-G 0 W 0 band gap Δ GW = 7.7 eV, as presented in section 2.1.

The limit of high polarizability-EBK procedure
Here we present the analytical form for the ground state exciton energy Ω 1s in the limit of high screening length (λ → ∞). As the ab initio results predict, hBN-ML can not be considered a material where the above condition applies, since λ = 10a 0 . However, other two-dimensional materials, such as transition-metal dichalcogenides, have a much larger λ. We have also performed the DFT calculations for molybdenum disulfide (MoS 2 -ML) and molybdenum ditelluride (MoTe 2 -ML) and obtained λ = 76a 0 and 120a 0 and μ = 0.25m e and 0.3m e , respectively. The relatively large λ for these two materials originates from their spatial structure where a transition metal plane lies between the planes of two chalcogenide atoms. This geometrical coordination reduces the electronic hopping elements between the atomic orbitals of the neighbouring atoms [57], therefore causing the flattening of the electron bands. Smaller dispersivity of the electron bands combined with the smaller band gap, such as Δ DFT ≈ 1.8 eV in the case of MoS 2 -ML and Δ DFT ≈ 1.2 eV in the case of MoTe 2 -ML, give a static screening length almost an order of magnitude larger compared to λ = 10a 0 in hBN-ML. This will certainly be responsible for the logarithmic behaviour of the screened potential over a sizeable electron-hole spatial extension, as can be seen in figure 4, making the logarithmic potential (3.2) adequate for determining the ground state exciton binding energy.
In this case the semiclasical EBK approach can be applied to calculate low-lying ( = 0) exiton energy levels. The general EBK approximation [58] asserts that in the spherical symmetric problem the phase integral of the radial impulse is quantized as Here r 0 is the classical turning point, n r = 0 is the radial quantum number for the ground state energy case, and z = 2 is the Maslov index [58], which gives the number of classical turning points. This leads to the implicit expression for the exciton ground state energy Ωλ/e 2 e Ωλ/e 2 − Ωλ/e 2 = π e γ z/4 32μe 2 λ . (3.5) In the case of a large screening length (λ/e 2 → ∞), when the exponential function dominates and the error function is Erf(x → ∞) = 1, equation (3.5) can be simplified, which finally leads to the analytical expression for the ground state exciton binding energy The similar expression has been found as the limiting solution of the Schrödinger equation for the logarithmic potential [59], while the results of reference [60] predict Ω 1s (λ) ≈ −(3/4π)e 2 /λ behaviour in the same limit suggesting that they neglected the logarithmic term. The expression (3.6) gives a fairly accurate description of Ω 1s (λ) in λ → ∞ limit, but it can be improved by setting z = 3.5. For that purpose we introduce the dimensionless values λ and μ trough λ = λ a 0 and μ = μ m e and by putting z = 3.5 in the expression (3.6) we have (3.7) Figure 6 shows the comparison between the ground state energy Ω 1s (λ) obtained using the static screened potential (2.29) (solid lines) and using the analytical form (3.7) (dashed lines), for two different reduced masses μ. It is evident that the agreement between the (solid and dashed) curves starts only for large values of λ. Using formula (3.3) outlined in section 3.2, together with formula (3.6) and the experimentally determined exciton energies ω exp 1s (MoS 2 ) ≈ 1.9 eV and ω exp 1s (MoTe 2 ) ≈ 1.2 eV [61][62][63], it is possible to estimate the single particle band gap of MoS 2 and MoTe 2 monolayers. Inserting the calculated polarizabilities and the reduced masses into (3.6) and using (3.3) gives Δ MoS 2 ≈ 2.5 eV and  a The values in parentheses are taken from reference [24]. b The values in parentheses are taken from reference [20] c The values in parentheses are taken from reference [25].
Δ MoTe 2 ≈ 1.6 eV. The DFT polarizabilities, the experimental exciton energies and the estimated quasiparticle band-gap for the three studied 2D crystals are summarized in table 2. It is worth noticing that the estimated band-gaps Δ agree satisfactorily well with other BSE-GW calculations [18][19][20][25][26][27] considering that the relative error is not larger than 15%, while the method used to derive Δ is very simple. The q2D-G 0 W 0 band-gaps Δ GW in hBN, MoS 2 and MoTe 2 are also compared with other G 0 W 0 calculations. The hBN direct gap (K → K) 7.7 eV slightly overestimates and the indirect gap (K → Γ) 6.2 eV (shown in figure 3) slightly underestimates the values 7.37 eV and 6.58 eV, respectively, reported in reference [24]. The MoS 2 band gap of 2.6 eV slightly underestimates the value 2.78 eV reported in reference [20]. As in reference [20], there is a problem here concerning the band gap in the M point which is 4.8 eV instead of 3.9 eV obtained by using 10 4 bands and a 500 eV CLFE cut-off. However, the high-resolution STM spectroscopy measurements [64,65] show that MoS 2 monolayer band-gap is in the range 2.2-2.4 eV. This substantial reduction is probably a result of MoS 2 monolayer being deposited on graphite which due to metallic screening reduces the band gap. Finally the MoTe 2 band-gap 1.9 eV overestimates the other results which are in range 1.72-1.77 eV [25][26][27]. The disagreement is probably because here the CLFE are included only in the perpendicular (z) direction which simultaneously results in the band gap converging with fewer bands.

Conclusion
The alternative q2D-G 0 W 0 -BSE method was applied to calculate the quasiparticle band-gap and the exciton binding energies in hBN, MoS 2 and MoTe 2 . The q2D-G 0 W 0 -BSE method is complementary to existing high performance GW-BSE methodology because it is computationally less demanding (at the G 0 W 0 stage the CLFE are included only in the direction of broken periodic symmetry), it excludes the spurious Coulomb interaction and takes advantage of the reduction of BSE to its 2D counterpart. The obtained results for hBN, MoS 2 and MoTe 2 single-particle band-gap and exciton binding energies agree well with other theoretical results. Using the calculated ground state exciton binding energies Ω 1s and the experimental exciton (excitation) energies ω 1s the realistic value of the single particle band gap Δ in hBN, MoS 2 and MoTe 2 monolayers are estimated and compared with other results. An analytical expression for the ground state exciton binding energy Ω 1s is obtained in the high polarizability limit using the EBK procedure. This was shown to be valid for the family of TMDs and therefore can be a very useful tool for experimentalists. where