# Band structure and excitonic properties of WSe2 in the isolated monolayer
limit in an all-electron approach

Niloufar Dadkhah and Walter R. L. Lambrecht  Department of Physics, Case
Western Reserve University, 10900 Euclid Avenue, Cleveland, OH-44106-7079

###### Abstract

A study is presented of the electronic band structure and optical absorption
spectrum of monolayer WSe2 using an all-electron quasiparticle self-consistent
$GW$ approach, QS$G\hat{W}$, in which the screened Coulomb interaction
$\hat{W}$ is calculated including ladder diagrams representing electron-hole
interaction. The Bethe-Salpeter Equation is used to calculate both the
screened Coulomb interaction $\hat{W}$ in the quasiparticle band structure and
the imaginary part of the macroscopic dielectric function. The convergence of
the quasiparticle band gap and lowest exciton peak position is studied as
function of the separation of the monolayers when using periodic boundary
conditions. The quasiparticle gap is found to vary as $1/d$ with $d$ the size
of the vacuum separation, while the excitonic lowest peak reaches convergence
much faster. The nature of the exciton spectrum is analyzed and shows several
excitonic peaks below the quasiparticle gap when a sufficient number of k
points is used. They are found to be in good agreement with prior work and
experiment after adding spin-orbit coupling corrections and can be explained
in the context of the Wannier-Mott theory adapted to 2D.

##  I Introduction

It is well known that 2D materials show strong excitonic effects thanks to the
reduced screening in two dimensions [1, 2] and the intrinsic differences
between the 2D and 3D Coulomb problem [3]. The lower dimensionality also has a
strong effect on the quasiparticle self-energy [4, 5]. Transition metal
dichalcogenides are a well studied family of 2D materials which can be
isolated in monolayer form. In reality these monolayers are always placed on a
substrate and often also capped with additional layers such as hBN. This
provides additional screening of the surroundings of the monolayer.
Determining the true isolated monolayer properties of a freestanding layer
experimentally is thus a challenging task. On the other hand, computationally,
determining the isolated limit is also challenging because almost all
calculation rely on periodic boundary conditions and thus the convergence as
function of interlayer spacing needs to be studied carefully. In this paper,
we pick WSe2 as an example to study this question.

Many-body-perturbation theory (MBPT) provides currently the most accurate
approach to the electronic band structure in the form of Hedin’s $GW$
approach[6, 7] where $G$ is the one-electron Green’s function and $W$ the
screened Coulomb interaction. Likewise the Bethe Salpeter Equation (BSE)
approach [8, 9, 10, 11, 12, 13, 14] provides an accurate approach for optical
absorption including local field and electron-hole interaction effects. Most
of the implementations of these methods start from density functional theory
(DFT) in the local density approximation (LDA) or generalized gradient
approximation (GGA) and make use of pseudopotentials and plane wave basis
sets. Here we use an all-electron implementation based on the full-potential
linearized muffin-tin orbital (FP-LMTO) method [15, 16], which uses an
auxiliary basis set consisting of interstitial plane waves and products of
partial waves [17] inside the spheres to represent two-point quantities, such
as the bare and screened Coulomb potentials, $v,W$, the polarization
propagator $P$ and the inverse dielectric response function
$\varepsilon^{-1}$.

Different levels of the $GW$ approximation are in use. The most commonly used
method is the $G_{0}W_{0}$ approximation in which perturbation theory is used
to find the correction due to the difference between the $GW$ self-energy
operator and the DFT exchange correlation potential used as zeroth order.
Self-consistent $GW$ in the sense conceived by Hedin [6] does not necessarily
improve the results because it appears to require one to include vertex
corrections and go beyond the $GW$ approximation. A successful way to make the
results of $GW$ independent of starting point is the quasiparticle self-
consistent $GW$ approach (QS$GW$)[18, 15]. Nonetheless, this method
overestimates the quasiparticle gaps because it underestimates screening.
Typically, it underestimates dielectric constants of semiconductors by about
20 % as illustrated for example in Fig. 1 in Ref. 19. To overcome this
problem, several approaches are possible. One can include electron-hole
interactions in the screening of $W$ by means of a time-dependent density
functional approach by including a suitable exchange-correlation kernel [20,
21]. The approach we use here was recently introduced by Cunningham et al.
[22, 23, 24] and includes ladder diagrams in the screening of $W$ via a BSE
approach for $W({\bf q},\omega=0)$. It can also be viewed as introducing a
vertex correction in the Hedin equation for $W$. No vertex corrections are
included in the self-energy $\Sigma$ but this is justified within the QS$GW$
approach by the approximate cancellation in the ${\bf q}\rightarrow 0$,
$\omega\rightarrow 0$ limits of the $Z$ factors in the vertex $\Gamma\propto
1/Z$ and the coherent part of the Green’s function $G=ZG^{0}+\tilde{G}$ [24].
Here $Z=(1-\partial\Sigma/\partial\omega)^{-1}$ is the renormalization factor
from the quasiparticlized $G^{0}$ to the dynamic $G$ Green’s functions and
$\tilde{G}$ is the incoherent part of the Green’s function.

In this paper, we apply this all-electron approach to the 2D material WSe2 and
study the convergence of various aspects of the method, such as number of
bands included in the BSE approach, QS$GW$ without and with ladder diagrams,
which we call QS$G\hat{W}$, and most importantly the dependence on the size of
the vacuum region. We compare with previous results and with experimental work
and discuss the nature of the excitons. Our results confirm the main
conclusion of Komsa et al. [5] in a study of MoS2 that the long-range effects
of the self-energy in $GW$ lead to a slow $1/d$ convergence of the
quasiparticle gap with $d$ the size of the vacuum region but with a similar
and compensating slow convergence of the exciton binding energy, resulting in
a much faster convergence of the optical gap corresponding to the lowest
allowed exciton peak. Besides the behavior as function of distance, we show
that the convergence of the excitons as function of the k mesh included in the
BSE calculations is crucial to obtain good agreement with the exciton gap. We
study the spectrum of exciton levels below the quasiparticle gap, including
dark excitons and explain their localization in k space and real space in
relation to the 2D Wannier-Mott theory. Spin-orbit coupling effects are
included in the bands and added as a-posteriori corrections to the excitons.

##  II Computational method

The band structure calculations carried out here use the Questaal package,
which implements DFT and MBPT using a linearized muffin-tin-orbital (LMTO)
basis set [16]. The details of the QS$GW$ approach and its justification
compared to fully self-consistent $GW$ in the Hedin set of equations, can be
found in Ref. 15. The recent implementation of the BSE in this code is
documented in Refs 22, 23, 24. Here we summarize our convergence parameters
for its application to WSe2. We use a double ($\kappa$,$R_{sm}$) basis set
$spdfgspdf$ on W, with $5p$ semicore-states and high-lying 6d states treated
as local orbitals. Here $\kappa^{2}$ is the Hankel function kinetic energy and
$R_{sm}$ is the smoothing radius of the smoothed Hankel function envelope
function. For Se, we use a $spdfspd$ basis set. The quasiparticle self-energy
matrix is calculated up to a cut-off of 2.5 Ry and is approximated by an
average diagonal value above 2.0 Ry. Convergence is studied as function of the
${\bf k}$ mesh on which the self energy is calculated and subsequently
interpolated to a finer mesh or along the symmetry lines by Fourier
transforming back and forth to the real space LMTO basis set. We find that a
$9\times 9\times 3$ set is converged to within 0.1 eV. Somewhat surprisingly,
even for a vacuum region as large as 35 Å we need more than one ${\bf k}$
division in the short direction of the reciprocal unit cell perpendicular to
the layers in order to obtain a clear $1/d$ behavior of the quasiparticle gap.
This reflects the long-range nature of the self-energy, as will be discussed
in more detail later. For the BSE calculations, and the nature of the
excitons, on the other hand, a finer in-plane mesh is important. The in-plane
lattice constant is chosen as $a=3.32$ Å in the $P\bar{6}m2$ space group,
which corresponds to an AA stacking of the trigonal prism unit cell in each
layer. The stacking is irrelevant since we focus on the non-interacting
monolayer limit. The $c$ lattice constant is varied to study the convergence.
The in-plane lattice constant is taken from Materials Project [25, 26] and was
optimized in the generalized gradient approximation (GGA). The convergence
with $d$, the distance between the layers, is a major object of our study and
discussed in Sec. III

##  III Results

###  III.1 Quasiparticle bands and band gaps

Fig. 1: Band structure of WSe2 in the LDA (blue dashed), QS$GW$ (red dotted),
and QS$G\hat{W}$ (green solid). The bands are referred to the valence band
maximum within each method.

First, we explore the electronic structure at different levels of theory (DFT,
QS$GW$ and QS$G\hat{W}$) for several interlayer distances $d$. The bandgap is
found to be direct with the VBM (valence band maximum) and CBM (conduction
band minimum) at points $K$ and $-K$. Fig. 1 shows the DFT, QS$GW$ and
QS$G\hat{W}$ band structure for $d=25$ Å specified by blue, red, green,
respectively. The band gap changes from 1.53 eV at LDA to 2.95 eV and 2.71 eV
in QS$GW$ and QS$G\hat{W}$, respectively using the converged $9\times 9\times
3$ mesh. In QS$G\hat{W}$, including the electron-hole coupling in terms of the
ladder diagrams in the calculation of the polarization propagator $P$ reduces
$W=(1-Pv)^{-1}v$ to $\hat{W}$ and hence the self-energy $\Sigma=iG\hat{W}$ and
the gap are also reduced. We note that
$(E_{g}^{\text{QS}G\hat{W}}-E_{g}^{\text{LDA}})/(E_{g}^{\text{QS}GW}-E_{g}^{\text{LDA}})\simeq
0.83$, which shows that the reduction of the gap is in good agreement with the
empirical approach of hybrid QS$GW$ (h-QS$GW$), which proposes to reduce the
QS$GW$ self-energy by a universal factor of $\sim$0.8. This approach was found
to lead to good agreement with experimental band gaps for several materials
[27, 19].

Fig. 2: Band gaps in QS$GW$, QS$G\hat{W}$ and QS$G\hat{W}$+BSE (optical gap)
as a function of the inverse interlayer distance using a $9\times 9\times 3$ k
mesh at $d=25$ Å. Dashed lines show the linear fitting to estimate the band
gap at $1/d=0$

To attain the strict monolayer limit, we change the interlayer distance from
17.5 Å up to 27 Å and examine the band gaps using different $\mathbf{k}$-mesh
samplings. Our results show that it is only for $9\times 9\times 3$ and
$8\times 8\times 3$ k mesh that we can see a clear linear trend for the gap as
a function of $1/d$. One expects a $1/d$ slow convergence because of the long-
range nature of the screened Coulomb interaction [5]. Fig. 2 shows the band
gaps versus $1/d$ for the $9\times 9\times 3$ k mesh. We note that although
the cell is quite large in the direction perpendicular to the plane, the long-
range nature of the self energy requires us to take more than one k point in
that direction to achieve convergence. As shown in Fig. 2, QS$GW$ and
QS$G\hat{W}$ band gaps converge quite slowly with $d$, such that in the limit
for infinite $d$, the extrapolated band gaps are around 3.5 eV and 3.2 eV,
respectively which are about 0.7 eV larger than the corresponding gaps at
$d=17.5$ Å.

###  III.2 Optical response

Fig. 3: (a)-(b) Imaginary part of the macroscopic dielectric function tensor
$\epsilon_{2}(\omega)$ for $x$ polarization within the IPA and BSE,
respectively; shown for four interlayer distances; (c) Oscillator Strengths
along $x$ associated with the excitonic states from the BSE on a logarithmic
scale for interlayer distance 27 Å; (d) Number of excitonic states (on
intervals of 0.1 eV) normalized to 1.

In this section, we turn to the optical absorption, which is closely related
to the imaginary part of the dielectric function $\varepsilon_{2}(\omega)$.
Figure 3(a,b) show $\varepsilon_{2}(\omega)$ for several interlayer distances
as calculated in the Independent Particle Approximation (IPA) and in the
Bethe-Salpeter Equation (BSE). The former omits both local field and electron-
hole interaction effects and corresponds to vertical transitions between the
bands. These calculations have been done using the QS$G\hat{W}$ bands as input
and the photon polarization was chosen along the $x$-direction. By symmetry
the macroscopic dielectric function is isotropic in the $xy$-plane and we are
primarily interested in the in-plane polarization. For the BSE calculations
shown in 3(b), we include 6 valence and 6 conduction bands in the two-particle
Hamiltonian on a $9\times 9\times 3$ k mesh. We use the Tamm-Dancoff
approximation which has been shown to be valid in most materials [28]. The
spectrum is broadened by including an imaginary part which linearly increases
with energy. Although there is an overall red-shift of the spectrum, one can
see that the difference between the IPA and BSE does not amount to a simple
shift toward lower energies but BSE significantly modifies the shape of the
spectrum. Moreover, in BSE a sharp peak occurs below the quasiparticle gap,
which corresponds to the first bright exciton. The position of the various
peaks in $\varepsilon_{2}(\omega)$ does not exhibit significant variation with
respect to the interlayer distance. In particular, unlike the QS$GW$ gaps
which change strongly with the size of the vacuum region, the position of the
first peak or the optical gap stays practically constant for this k mesh. This
indicates an increasing binding energy with size of the vacuum region, as
expected, given the long range nature and decreasing screening of the screened
Coulomb interaction $W$ with size of vacuum region. This is shown also in Fig.
2, which shows that the optical gap is nearly constant as function of
interlayer distance within the range of distances shown. In fact, it appears
to slightly decrease with layer distance, but this variation is within the
error bar of the calculations. It may be due to a slightly decreasing
convergence of the basis set for larger interlayer distances. Furthermore, as
we discuss in the next section, the exciton peak is sensitive to the k-mesh
sampling. In the result shown in Fig. 2 the mesh is kept fixed at the same
mesh as used in the QS$G\hat{W}$ band calculation, and provides only a lower
limit of the optical or exciton gap.

Fig. 3(c) shows the oscillator strengths for all eigenvalues of the two-
particle Hamiltonian on a log-scale and Fig. 3(d) shows the density of
excitonic states (normalized to 1), both at $d=27$ Å. The difference between
the density of excitonic states and the actual $\varepsilon_{2}(\omega)$
indicates the modulation of the spectrum by the dipole matrix elements. Fig.
3(c) shows that the matrix elements can vary by several orders of magnitude.
Only the levels with highest oscillator strength in (c) show up as peaks in
(b). The oscillator strengths here are not normalized and thus only the
relative intensity between different levels has physical meaning.

###  III.3 Convergence of exciton gap

To converge the energy of the first exciton peak, it is necessary to use a
finer k mesh, since it is localized in k space. On the other hand, as we will
show in III.5, for the lowest excitons, we can restrict the active space of
bands to fewer valence and conduction bands included in the BSE. Fig. 4 shows
the exciton gaps as function of $1/N_{k}^{2}$, where $N_{k}$ is the total
number of k in the 2D mesh and for a distance $d=20$ Å, which is close to the
converged monolayer limit and including 4 valence and 4 conduction bands. For
the direction normal to the plane, three k points are used. In fact, as
mentioned earlier, from Fig. 2 the BSE calculated lowest exciton peak does not
vary significantly with $d$. The optical gap for the isolated monolayer limit
and for $N_{k}\rightarrow\infty$ is thus found to be $\sim$2.1 eV, whereas the
quasiparticle gap in that limit was found to be 3.2 eV, implying an exciton
binding energy of 1.1 eV.

Fig. 4: Energy of the lowest exciton as a function of the inverse of the
number of k points squared for $d=20$ Å.

###  III.4 Discussion of exciton spectrum in relation to prior theory work and
experiment

Fig. 5 shows $\varepsilon_{2}(\omega)$ in the BSE as in Fig. 3(b), for $d=25$
Å, but using a finer k mesh of $30\times 30\times 3$ with $N_{v}=N_{c}=1$.
This finer mesh provides us with a better picture of the profile of the
spectrum below and just above the quasiparticle gap. We can readily notice
that compared to Fig. 3(b), more weight is shifted to the first peak, making
it the strongest optical excitation.

Fig. 5: Imaginary part of the macroscopic dielectric function tensor
$\epsilon_{2}(\omega)$ for $x$ polarization within the BSE using one valence
and one conduction band on $30\times 30\times 3$ k mesh. Label X shows the
ground state exciton without SOC, know as $1s$. Adding SOC, splits X exciton
into two excitons commonly known as excitons A and B, which are schematically
shown as magenta A and B labels. X′ shows the second bright exciton, known as
$2s$. A third prominent peak also appears, which is labeled as X′′. Here, the
spectrum is shifted so that the energy of the first peak coincides with the
converged value of 2.1 eV.

Table 1 shows some of the experimental and theoretical results for the energy
of A exciton (first bright), binding energy of the A exciton and the energy
splitting between A and A′(second bright). To compare our results here with
experiment and other literature work, it is important to include spin-orbit
coupling (SOC) effects. Questaal’s BSE code does not yet allow for the
inclusion of SOC; however, we can include spin-orbit coupling perturbatively
in our QS$G\hat{W}$ calculations. In other words, we add the SOC Hamiltonian
to the QS$G\hat{W}$ Hamiltonian matrix in the LMTO basis set and rediagonalize
it. We don’t just evaluate its effect on each separate eigenstate but we do
not include SOC effects in the calculation of the self-energy and do not use a
full spin-dependent Green’s function in the $GW$ steps. Adding SOC splits the
lowest exciton peak into two commonly known as A and B excitons which results
from the SOC splitting of the highest valence and lowest conduction bands at
and around point $K$ [29]. Our results show that at point $K$, the valence
band splitting is about 445 meV with the VBM-B going 225 meV down and VBM-A
219 meV up in energy. However, the CBM, which has a much smaller splitting
also shifts down for both spin directions by about 0.035 eV. Including these
a-posteriori SOC effects, the A exciton in our results is placed at $\sim$1.8
eV which is in excellent with results from [30]. These authors include SOC
fully as well as perturbatively and it is with the latter that we best agree.
Most of the experimental results also place the A exciton at about 1.7-1.8 eV,
whereas some older theoretical calculations [31] find it near 1.5 eV. In Fig.
5 we label our exciton peaks as X, X′ and X′′ to indicate that they do not
include SOC. Each would split into an A and B peak when adding SOC.

Our calculated binding energy of $\sim$1 eV is slightly higher than previously
reported. Ramasubramaniam [31] obtains a gap of 2.42 eV at K, but his
calculation includes spin-orbit coupling (SOC). If we compare with the average
of the SOC split states, his quasiparticle gap would be 2.67 eV but this is at
a interlayer separation $d=15$ Å. Furthermore, this author obtained a lower
indirect $K$-$\Lambda$ gap from the valence band at $K$ to a secondary valley
in the conduction band along the $\Gamma$-$K$ axis $\Lambda$. In Fig. 2 the
smallest distance we considered is 17.5 Å and it already gives a slightly
higher gap. Our extrapolated quasiparticle gap for infinite distance is
significantly higher. Ramasubramaniam’s lowest exciton peak occurs at 1.52 eV
for the A exciton and 2.00 eV for the B exciton, showing that this splitting
simply corresponds to the corresponding valence band SOC splitting at K. His
exciton binding energy is thus reported as 0.9 eV. He used a k mesh of
$6\times 6\times 1$. We can see that our exciton binding energy is only
slightly higher than his while our optical gap and quasiparticle gap are both
significantly higher but in better agreement with [30] and experiment.

Moreover, Fig. 5 shows that between the lowest exciton and the quasiparticle
gap, other peaks occur and by examining the oscillator strengths, we see that
there are also several dark excitons in this region as shown in Fig. 3(c). We
identify the next bright exciton which is slightly lower in intensity but
clearly visible in the $\varepsilon_{2}(\omega)$ in Fig. 5 as the X′ exciton,
which in literature is usually associated with either a $2s$ or $2p$ type
exciton using the hydrogenic Rydberg series nomenclature. Whether $2s$ or $2p$
are showing up in experiment depends on the experimental probe used: two-
photon absorption reveals the $2p$ excitons, while linear absorption or
luminescence reveals the one-photon absorption. This arises from the different
matrix elements and symmetry selection rules. We should note that this
classification does not imply a strict hydrogenic Rydberg series, because the
system is 2D rather than 3D and the screened Coulomb potential does not behave
simply as $1/\varepsilon r$ because the screening in 2D is strongly distance-
dependent [32]. Furthermore, a $2p$ type exciton having an envelope function
which has nodes would be dark in the $\varepsilon_{2}(\omega)$ spectrum but
accessible in two-photon absorption experiments. Adding SOC will split this
peak into two peaks known as A′ and B′. Following our previous discussion of
including SOC effects, We expect the splitting between A and A′ to stay almost
the same as X and X′ since we associate them to $s$-like states. We call the
next peak X′′, which is close to X′ in energy but relatively higher in
intensity. We also caution that the details of the exciton spectrum shown in
Fig. 5, such as peak splittings and positions depend strongly on the
convergence with k points. However, qualitatively Fig. 5 is in good agreement
with the results of [30] without SOC and shows the same relative intensities
for X, X′ and X′′. In particular, our X-X′ splitting is about 255 meV which is
very close to the reported value of 215 meV by [30] for a $39\times 39\times
1$ k-mesh sampling. This value is higher than values in Table 1, which the
authors contribute to the substrate effects used in experiments. These authors
also report a value of 140 meV using a coarser k mesh, which is closer to
experimental results.

We conclude that both convergence as function of k mesh and vacuum thickness
plays an important role. This was also concluded by Qiu et al. [33] in the
case of MoS2. These authors in fact dealt with the monolayer separation by
using a truncated Coulomb potential and used an even finer k mesh. Their use
of a truncated Coulomb interaction allows for an easier convergence to
isolated monolayers and they compare the screening dielectric constant with
the strict 2D Keldysh limit [1].

Table 1: A exciton energy, A exciton binding energy and A-A′ energy splitting of monolayer WSe2 as reported in several experimental and theoretical studies. | Method | Optical gap (eV) /Eb (meV) | A-A′ splitting (meV)
---|---|---|---
hBN/ML-WSe2/hBN [34] | Magneto-PL | 1.727/170 | 131
hBN/ML-WSe2/hBN [35] | Magneto-PL | 1.712/172 | 131
ML-WSe2/Si/SiO2 [36] | Linear absorption and two-photon PL | 1.65/370 | 160
hBN/ML-WSe2/hBN [37] | Magnetoabsorption | 1.723/161 | 130
hBN/ML-WSe2/hBN [38] | Gate-dependent reﬂection and PL |  $\sim$1.72/- | 120
hBN/ML-WSe2/hBN [39] | Magneto-photocurrent | 1.725/168.6 | 128.6
hBN/ML-WSe2/hBN [40] | EELS | 1.697/- | -
ML-WSe2/Si/SiO2 [41] | PL and reﬂectance | 1.744/- | -
ML-WSe2/PC [42] | Absorption |  $\sim$1.66/710 | 200
ML-WSe2 [43] | Momentum-resolved EELS | 1.69/650 | -
ML-WSe2/Si/SiO2 [44] | 1- and 2-photon PL excitation | 1.75/600$\pm$200 | 150
ML-WSe2/SiO2 [45] | Reflectance | 1.67/- | -
hBN/ML-WSe2/hBN [46] | PL |  $\sim$1.75/- | -
ML-WSe2/hBN/n-doped Si [47] | Photoemission | 1.73/390$\pm$150 | -
hBN/ML-WSe2/hBN [48] | Magnetoabsorption | 1.722/- | -
ML-WSe2[31] |  $GW$-BSE | 1.52/900 | -
ML-WSe2[30] |  $GW$-BSE |  $\sim$2.1/596 | -
ML-WSe2[30] |  $GW$-BSE-Pert. SOC |  $\sim$1.8/533 | 140
ML-WSe2 [30] |  $GW$-BSE-SOC |  $\sim$1.7/550 | 140
ML-WSe2 [49] | LDA+$GdW$+BSE+SOC | 1.76/- | -
ML-WSe2 [31] |  $G_{0}W_{0}$+BSE+SOC | 1.52 |
ML-WSe2 [44] |  $GW_{0}$+BSE |  $\sim$2 |
ML-WSe2 [50] |  $GW$+BSE | 1.78/620 | -

###  III.5 Exciton wave function analysis

####  III.5.1 Overview

To better understand the nature of excitons, we can analyze the k-space and
real-space distribution of the first few excitons. First, we calculate
$A_{\lambda}^{v,c}(\mathbf{k})=\innerproduct{v,c,\mathbf{k}}{\Psi_{\lambda}^{Ex}}$
for a particular band pair $\\{v,c\\}$ along several high-symmetry paths in
the BZ where $\lambda$ denotes the excitonic eigenstate. Here we focus on the
first few excitons by separately considering the ranges of energies which
include only $\lambda=1,2$ (degenerate-bright), $\lambda=3$ (dark),
$\lambda=4$ (dark), $\lambda=5,6$ (degenerate-semi-bright), $\lambda=7,8$
(degenerate-semi-bright) and $\lambda=21,22$ (degenerate-bright). We should
note that bright excitons $\lambda=1,2$, $\lambda=5,6$ and $\lambda=7,8$
correspond to X, X′ and X′′ peaks in Fig. 5, respectively.

In this section we use $N_{v}=1$, $N_{c}=1$, which then allows us to use a
finer k mesh. This will be justified in Sec. III.5.2. However, we have checked
that they give essentially the same results for the exciton region below the
quasiparticle gap as a $N_{v}=4$, $N_{c}=4$ for the smaller k meshes used. For
a study of the overall $\varepsilon_{2}(\omega)$ up to higher energies in the
continuum, it is important to include more bands but this limits the number of
k points we can afford. When focusing on the excitons, the opposite is true
and the important convergence parameter is the number of k points.

Fig. 6 in section III.5.2 shows the band weight of various excitons. We use
the weight $W_{\lambda}^{c}({\bf k})=\sum_{v}|A^{vc}_{\lambda}({\bf k})|^{2}$
at each k by the size of the colored circle in the conduction bands $c{\bf k}$
and likewise $W_{\lambda}^{v}=\sum_{c}|A^{vc}_{\lambda}({\bf k})|^{2}$ for the
weight on the valence bands. The colors only distinguish between different
bands but have no physical meaning.

Next, we can investigate the $A_{\lambda}^{v,c}({\mathbf{k}})$ for a
particular exciton and band pair as function k on a mesh of k points in Sec,
III.5.3. These exciton expansion coefficients, in fact are the Fourier
transforms of a slowly varying envelope function within the Wannier-Mott ${\bf
k}\cdot{\bf p}$ theory of excitons [33]. This slowly varying envelope function
would be a hydrogenic function in case of a 3D isotropic system and the full
real-space exciton function is then modulated by the product of the Bloch
functions at the valence and conduction band edges at the ${\bf k}$ point
where the exciton is centered. Here they are 2D analogs and somewhat distorted
by the anisotropies. In 2D the angular behavior of the envelope function is
described by the 2D Laplacian which has solutions of the form
$\cos{m\phi},\sin{m\phi}$ and $m=0,1,2,\dots$ are referred to as $s,p,d,\dots$
in analogy with the 3D case. The general solution $e^{im(\phi-\phi_{0})}$
contains an arbitrary phase, resulting from the overall arbitrary phase of the
eigenstates of the two-particle Hamiltonian diagonalization, but this only
affects the orientation of these exciton wavefunctions in k space. This is
valid as long as the slowly varying potential for which the slowly varying
envelope function is a solution of the effective Schrödinger equation, is
indeed axially symmetric. This is true for k points close to the point $K$ in
the region contributing to the exciton. Beyond a certain distance in k space,
one might expect a trigonal warping. The real-space excitons at short distance
thus will not show a circular but rather hexagonal shape. Conversely, the less
strongly bound excited excitons are more spread in real space and hence more
localized in k space and hence in the region where circular symmetry is valid.
The main advantage of this analysis is that it shows us the sign changes of
the envelope functions as function of angle and hence the angular patterns,
which can explain why certain excitons are dark.

Finally, we will examine the excitons fully in real space in Sec. III.5.4. In
that case, we look at $|\Psi({\bf r}_{h},{\bf r}_{e})|^{2}$ and we can either
fix the hole at ${\bf r}_{h}$ or the electron at ${\bf r}_{e}$ and then plot
the probability distribution of the other particle. Since the valence and
conduction bands, from which the excitons are derived, are all mainly W-$d$
like, either choice is equally good and we choose to fix the hole. These will
be shown in Sec. III.5.4. These will show both the large scale envelope
function behavior and the local atomic orbital character.

Before we present our exciton analysis results let us mention that our
discussion of bright and dark excitons here is limited to singlet excitons.
Excitons can also be dark due to the spin structure as was discussed for WSe2
and MoS2 by Deilman et al. [49]. The darkness of excitons considered here is
rather related to the excited nature of the excitons in the Rydberg series,
and their slowly varying envelope function in a Wannier-Mott exciton model.

We also note that all bright excitons are doubly degenerate and all dark
excitons are singly degenerate. The double degeneracy of the first bright
exciton results from the equal contributions at $K$ and $-K$. These are
related by time reversal symmetry and hence, there is a Kramers degeneracy.
Their double degeneracy can also be understood from the symmetry analysis of
Song et al. [51]. In fact, the exciton is expected to be localized near the W
atom since the top valence and bottom conduction bands are both dominated by W
orbitals. The excitons can thus be classified according to the symmetry of the
W-site which is the full symmetry $D_{3h}$ of the crystal. Because we focus on
excitons with polarization in the plane, this immediately tells us that all
bright excitons will be doubly degenerate and correspond to the $E^{\prime}$
representation. Dark excitons which are symmetric with respect to the
horizontal mirror plane (as we find to be the case for all the ones we
examined) must thus correspond to the $A_{1}^{\prime}$ or $A_{2}^{\prime}$
irreducible representations which are both non-degenerate. The distinction
between the two is whether they are even or odd with respect to the vertical
mirror planes and or 2-fold rotations in the plane.

####  III.5.2 Bandweights

Fig. 6 shows the band weights for various excitons. We can clearly see that
almost all the contribution comes from the pair of the highest valence and
lowest conduction bands focused on and around the point $K$ in the Brillouin
zone and decaying as we move away from $K$ in an anisotropic way. For example
the conduction band has a larger mass along the $KM$ than the $K\Gamma$
direction and the exciton spreads out further along the $KM$ direction. We can
explain this in terms of the Wannier-Mott hydrogenic model of the exciton.
First, indeed, the localization in k space implies a large spread in real
space consistent with the Wannier-Mott type exciton. Second, within this model
the effective Bohr radius of the exciton is inversely proportional to the band
effective mass and hence gives a smaller Bohr radius in real space for the
$MK$ direction which corresponds to a larger spread in k space. Figures 6(b)
and 6(c) demonstrate the weights for the second and third excitons in energy,
which are $|A_{\lambda=3}^{vc\mathbf{k}}|$ and
$|A_{\lambda=4}^{vc\mathbf{k}}|$, respectively. For these two dark and non-
degenerate excitons, we can clearly see that almost all the contribution still
comes from only the pair of the highest valence and lowest conduction bands
focused around the point $K$ in the Brillouin zone but now with a zero value
at point $K$ itself in clear contrast to the first bright exciton where the
contribution at point $K$ is the largest, as shown in Fig. 6(a). Fig. 6(d),
shows the weight
$|A_{\lambda=5}^{vc\mathbf{k}}|,|A_{\lambda=6}^{vc\mathbf{k}}|$ which
corresponds to fourth exciton in energy, doubly degenerate and semi-bright,
i.e. less bright than the first bright exciton. This exciton has
characteristics in between those of the previous bright and dark excitons,
with an overall distribution like the dark excitons but nonzero contribution
at point $K$. Fig. 6(e) corresponds to fifth exciton in energy which resembles
the first bright exciton in oscillator strength but less extended around $K$.
Fig. 6 (f) shows the excitonic weights for an exciton with higher energy
which, in its k-space distribution, resembles the dark excitons even though it
has a strong oscillator strength. In k space we see it has larger dip at and
around point $K$.

Fig. 6: Weights of exciton wave function contributed by different bands
(a)–(f) for six excitonic eigenvalues. The size of the colored circles
indicates the exciton weight $|A^{\lambda}_{vc{\bf k}}|^{2}$ for a given
exciton $\lambda$. The colors have no meaning and only serve to distinguish
different bands.

####  III.5.3 k-space analysis

To explain the behaviors noted in Sec. III.5.2, we need to fully examine the
k-space relations of the $A_{\lambda}^{vc}({\bf k})$ instead of only looking
at their absolute values along high-symmetry lines. Because the bright
excitons are doubly degenerate, any linear combination of the two eigenstates
is also an eigenstate and this makes it difficult to visualize them. Thus, in
this section we focus on the dark excitons.

In Fig. 7 we show the real, imaginary and absolute values of the envelope
function $A_{\lambda}^{vc}({\bf k})$ for the first two dark excitons over the
2D k space in one unit cell as a color plot. This requires an even finer k
mesh but we can now restrict the $v$ and $c$ to just the valence and
conduction band edges. We use here a $30\times 30$ mesh. In spite of the
rather coarse pixelation, we can see that around both the points $K$ and
$K^{\prime}=-K$ there is a three fold symmetry with alternating positive and
negative signs as we go around the central point $K$ or $K^{\prime}$. In
between we may see a white line of pixels indicating a node. This indicates a
$cos(3\phi)$ type of behavior for the envelope function. There appears to be a
slight offset of the nodal lines from the symmetry directions but this may
arise from an arbitrary overall phase of the two particle eigenstates. It
indicates that if the matrix elements of the velocity operator in the plane
are fully symmetric around $K$ as we expect, then the exciton will be dark by
this modulation of the three fold symmetry envelope function. Furthermore we
can see that for the $\lambda=3$ exciton the signs of the k-space envelope
function are opposite near $K$ and $K^{\prime}$ or more precisely the function
is odd with respect to a vertical mirror plane going through the $\Gamma-M$
line while for $\lambda=4$ the signs are the same. So, we conclude that
exciton 3 can be described as an $A_{1}^{\prime}$ exciton while $\lambda=3$ is
an $A_{2}^{\prime}$ exciton. But both become dark by the three-fold symmetry
envelope function.

As already discussed in Sec. III.5.1 a $\cos{(m\phi)}$ behavior is expected
but the question arises why we only see a $m=3$ or $f$-like envelope function
for these lowest energy dark excitons. States with lower $m$ are indeed
expected but would be more spread out in real space (the higher the anguar
momentum the more localized by the centrifugal term (of the form
$m^{2}/r^{2}$) and hence more localized in momentum space; therefore, in order
to access $p$-like or $d$-like states originating from and around point $K$,
as reported in [33] for MoS2, we would need a much finer k mesh to capture the
fine structure around this point. With the k mesh used in this study, we only
access points in momentum space which are not close ”enough” to $K$
responsible for $p$\- or $d$-like states.

Our analysis in terms of the envelope function resembles the approach used by
Qiu et al. [33] for MoS2. However, they used a special dual mesh approach,
which allowed them to use k meshes as large as $300\times 300$. With our 10
time smaller divisions we still obtain a rather pixelated view of the k space
behavior. Therefore we can unfortunately not obtain the nodal structure near
each point $K$ as accurate. The effective k-spacing is $(2\pi/a)/30$ is about
0.06 Å-1 but the size of the exciton spread itself is of order 1 Å-1 and thus
we don’t have sufficient resolution to fully see the fine structure which
would reveal the Rydberg like excited states of the exciton. The number of
excitons we can resolve also depends strongly on the fineness of the k mesh
and as mentioned before does not allow us to capture the $p$ or $d$-like
excitons.

In fact, we did notice that the number or excitons depends on the size of the
k-mesh and went up to $35\times 35\times 3$ with $N_{v}=1$, $N_{c}=1$ but were
still not able to obtain the $p$\- or $d$-like excitons.

Fig. 7: Real, imaginary and absolute values of the envelope function
$A^{\lambda}_{vc{\bf k}}$ for the first two dark excitons over the
$\mathbf{k}$ space. First and second rows belong to $\lambda=3$ and
$\lambda=4$, respectively. High symmetry points are labeled as magenta in the
absolute value panel as reference. Fig. 8: Real-space probability density of
the first three bright excitons in energy and one higher energy exciton, shown
for a section that goes through the tungsten atoms perpendicular (Top) and
parallel (bottom) to the out-of-plane axis, respectively. In the bottom row of
each panel, the section is situated in the middle of the suppercell. For the
color bar limits, the same minimum (cyan) and maximum (yellow) values are used
in all three maps. Values smaller (larger) than this range are shown by the
saturation colors cyan (yellow). Fig. 9: Real-space probability density of the
first two dark excitons in energy, shown for a section that goes through the
tungsten atoms perpendicular (Top) and parallel (bottom) to the out-of-plane
axis, respectively. In the bottom row of each panel, the section is situated
in the middle of the supercell. For the color bar limits, the same minimum
(cyan) and maximum (yellow) values are used in both maps. Values of
$|\Psi|^{2}$ smaller (larger) than this range are shown by the saturation
color cyan (yellow).

####  III.5.4 Real-space analysis

Next we turn our attention to the real space spread of the exciton wave
functions. In Figures 8 and 9 we visualize the exciton wave functions in real
space as a colormap. The interlayer distance of 25 Å is used for this
calculation.

Fig. 8(a) shows the spatial distribution of the lowest energy exciton on a
section going through the W atom, perpendicular to the c-axis (top) and
perpendicular to the a-axis (bottom). As we can see, the overall shape
resembles the hexagonal symmetry of the structure rotated by 60∘ such that the
highest weight (yellow) lies on a hexagonal region extending from two to four
unit cells away from the hole. Looking closer, we can see that there is a
3-fold symmetry. Also, around each W atom, a dumbbell-like distribution of the
electron density can be seen, which originates from the $d$-orbitals of the
individual W atoms. This can best be seen in the edge-on projection, which
clearly shows $d$-orbital shapes. The orientation of this dumbbell-like
distribution changes locally going from one edge to the other. For the two
dark excitons we see a different character shown in Fig. 9. In 9(a), the
distribution is the highest around W while in 9(b), it is almost zero around
W. Furthermore, in Fig. 9(a) we see clear nodal lines at 120∘ from each other,
reminiscent of the $\cos{(3\phi)}$ pattern in the k-space figure. In Fig. 9(a)
we still see regions of low values in rings around each atom but still showing
some three-fold symmetry. However, their overall sizes in the real space are
pretty close to each other and somewhat smaller than the previous bright
exciton. On the other hand, for the next bright exciton shown in Fig. 8(b) is
spread out substantially farther. As the exciton energy approaches the gap and
the exciton binding energy is reduced, the exciton spreads further in real
space.

Interestingly, even for the lowest energy bright exciton the envelope function
is not at its maximum in the center but is rather ring shaped with a dip in
the very center. In the second bright exciton, shown in Fig. 8(b), the central
region is distinctly triangular shaped corresponds to a larger central region
with very low value. The spatial extent of this exciton is significantly
larger than the first bright exciton, consistent with the $2s$ character.
These features can be explained in terms of the 2D envelope functions which
correspond to the electron-hole attraction potential. Because of the 2D
screening being very distant dependent, it behaves like a weak logarithmic
potential $\ln{(r)}$ near the origin and as a unscreened $1/r$ potential far
away [2]. The solutions for the 2D Schrödinger equation with a logarithmic
potential were studied by Atabek, Deutsch and Lavaud [52]. Because of the weak
logarithmic potential the $m=0$ (or $s$-like states) do not start as $r^{0}$
or a constant but more like $r$ as can be seen in Fig. 4 of [52]. They go
through a single maximum at finite $r$ consistent with the ring shaped
envelope we see. For the $2s$ state, there is a node but the onset starts not
at $r=0$ but at finite $r$. The stronger unscreened Coulomb interaction at
large distance means that excited states like $2s$ states are pushed out to
larger $r$. This is consistent with the large triangular hole in this exciton
wave function. The triangular shape at close distance in real space results
from the triangular warping of the bands involved in k space at large k. Even
for the lowest exciton, the center region is triangular rather than circular.
The next bright exciton, shown in Fig. 8(c) shows very similar behavior as the
first bright exciton. We recall that the doubly degenerate exciton labeled 5,6
is the X′ exciton while 7,8 corresponds to the X′′ exciton. (We here call them
X instead of A,B because our calculation does not include spin-orbit
coupling). It is somewhat surprising that X′′ and X are so similar in spatial
extent. However, we can see from Fig. 8(a,c) that they have different
probability density. We picked the minimum and maximum values mapped to the
color scale the same for both to facilitate a an absolute comparison. They
also have somewhat different k-space localization: X′′ falls of faster with
$k$ away from the band edge $K$ than X. Finally, the bright exciton with
largest exciton energy (Fig. 8(d) is seen to have a wider spread and a clear
radial node structure.

Overall, the real space and reciprocal space exciton wave functions show an
interesting variety of behaviors which can largely be understood in terms of
the 2D problem with a Keldysh type distant dependent screened Coulomb
potential.

##  IV Conclusions

In this paper we presented a study of the band structure and excitons in
monolayer WSe2 using an all-electron muffin-tin-orbital implementation of the
QS$GW$ and BSE approaches. The calculations went beyond the usual IPA and RPA
based $GW$ approach by including the electron-hole interaction effects in the
screening of $W$, which is then called $\hat{W}$. The quasiparticle gaps were
found to converge slowly as function of interlayer distance, as $1/d$ and the
replacement of $W$ by $\hat{W}$ essentially provides a constant downward shift
of the gap. The exciton spectrum was found in excellent agreement with prior
theory work [30] and various experiments when SOC corrections of the bands are
added a-posteriori. As expected, the excitons show a smaller variation as
function of interlayer distance, as a result of the compensation of the
exciton binding energy and quasiparticle gap self-energy shift which are both
proportional to $\hat{W}$. Converged exciton eigenvalues were shown to require
a very fine in-plane ${\bf k}$ mesh in order to obtain agreement with
experimental exciton gaps. The nature of dark and bright excitons was analyzed
in terms of their reciprocal space and real space envelope functions and
showed distinct behaviors. First, all the excitons examined were derived
almost exclusively from the top valence and lowest conduction bands near $K$.
But the dark ones had zero value at $K$ itself. The darkness of two of the
excitons was found to result from the three-fold symmetry nodal pattern of the
envelope function but these two dark excitons still show distinct behavior in
real space. The exciton symmetries were analyzed in the context of the
Wannier-Mott theory adjusted for 2D. The three-fold symmetry of the first dark
excitons found here indicates $f$-like excitons. We expect that there are also
$p$\- and $d$-like envelope functions in a Rydberg like series, but the reason
why the lower angular momenta $p$\- or $d$-like excitons could not be revealed
is related to the limitations in the number of k points we could include in
the BSE. A finer mesh near the band extrema at $K$ and $-K$ is required to
capture these excitons which are expected to be more localized in k space.
Among the bright excitons, we were able to identify the A′ exciton and
identify it with a $2s$ like exciton state. Our present study provides a
detailed look at how the Mott-Wannier model explains qualitatively the exciton
series. However, it shows that the strongly distance dependent screening of
the Coulomb interaction between electron and hole modifies the spatial extent
and shape of the exciton envelope functions. These features emerge naturally
from the first-principles BSE approach even though we do not explicitly solve
for a Wannier-Mott slowly varying envelope function.

###### Acknowledgements.

This work was supported by the U.S. Department of Energy Basic Energy Sciences
(DOE-BES) under Grant No. DE-SC0008933. Calculations made use of the High
Performance Computing Resource in the Core Facility for Advanced Research
Computing at Case Western Reserve University.

## References

  * Keldysh [1979] L. V. Keldysh, Coulomb interaction in thin semiconductor and semimetal films, JETP Lett. 29, 658 (1979).
  * Cudazzo  _et al._ [2011] P. Cudazzo, I. V. Tokatly, and A. Rubio, Dielectric screening in two-dimensional insulators: Implications for excitonic and impurity states in graphane, Phys. Rev. B 84, 085406 (2011).
  * Shinada and Sugano [1966] M. Shinada and S. Sugano, Interband Optical Transitions in Extremely Anisotropic Semiconductors. I. Bound and Unbound Exciton Absorption, Journal of the Physical Society of Japan 21, 1936 (1966).
  * Delerue  _et al._ [2003] C. Delerue, G. Allan, and M. Lannoo, Dimensionality-Dependent Self-Energy Corrections and Exchange-Correlation Potential in Semiconductor Nanostructures, Phys. Rev. Lett. 90, 076803 (2003).
  * Komsa and Krasheninnikov [2012] H.-P. Komsa and A. V. Krasheninnikov, Effects of confinement and environment on the electronic structure and exciton binding energy of MoS2 from first principles, Phys. Rev. B 86, 241201 (2012).
  * Hedin [1965] L. Hedin, New method for calculating the one-particle green’s function with application to the electron-gas problem, Phys. Rev. 139, A796 (1965).
  * Hedin and Lundqvist [1969] L. Hedin and S. Lundqvist, Effects of electron-electron and electron-phonon interactions on the one-electron states of solids, in  _Solid State Physics, Advanced in Research and Applications_ , Vol. 23, edited by F. Seitz, D. Turnbull, and H. Ehrenreich (Academic Press, New York, 1969) pp. 1–181.
  * Strinati [1988] G. Strinati, Application of the green’s functions method to the study of the optical properties of semiconductors, La Rivista del Nuovo Cimento (1978-1999) 11, 1 (1988).
  * Hanke [1978] W. Hanke, Dielectric theory of elementary excitations in crystals, Advances in Physics 27, 287 (1978).
  * Hanke and Sham [1980] W. Hanke and L. J. Sham, Many-particle effects in the optical spectrum of a semiconductor, Phys. Rev. B 21, 4656 (1980).
  * Albrecht  _et al._ [1998] S. Albrecht, L. Reining, R. Del Sole, and G. Onida, Ab initio calculation of excitonic effects in the optical spectra of semiconductors, Phys. Rev. Lett. 80, 4510 (1998).
  * Rohlfing and Louie [1998] M. Rohlfing and S. G. Louie, Electron-Hole Excitations in Semiconductors and Insulators, Phys. Rev. Lett. 81, 2312 (1998).
  * Benedict  _et al._ [1998] L. X. Benedict, E. L. Shirley, and R. B. Bohn, Optical absorption of insulators and the electron-hole interaction: An ab initio calculation, Phys. Rev. Lett. 80, 4514 (1998).
  * Onida  _et al._ [2002] G. Onida, L. Reining, and A. Rubio, Electronic excitations: density-functional versus many-body Green’s-function approaches, Rev. Mod. Phys. 74, 601 (2002).
  * Kotani  _et al._ [2007] T. Kotani, M. van Schilfgaarde, and S. V. Faleev, Quasiparticle self-consistent GW method: A basis for the independent-particle approximation, Phys.Rev. B 76, 165106 (2007).
  * Pashov  _et al._ [2019] D. Pashov, S. Acharya, W. R. Lambrecht, J. Jackson, K. D. Belashchenko, A. Chantis, F. Jamet, and M. van Schilfgaarde, Questaal: A package of electronic structure methods based on the linear muffin-tin orbital technique, Computer Physics Communications , 107065 (2019).
  * Aryasetiawan and Gunnarsson [1994] F. Aryasetiawan and O. Gunnarsson, Product-basis method for calculating dielectric matrices, Phys. Rev. B 49, 16214 (1994).
  * van Schilfgaarde  _et al._ [2006] M. van Schilfgaarde, T. Kotani, and S. Faleev, Quasiparticle Self-Consistent $GW$ Theory, Phys. Rev. Lett. 96, 226402 (2006).
  * Bhandari  _et al._ [2018] C. Bhandari, M. van Schilfgaarde, T. Kotani, and W. R. L. Lambrecht, All-electron quasiparticle self-consistent $\mathit{GW}$ band structures for ${\mathrm{SrTiO}}_{3}$ including lattice polarization corrections in different phases, Phys. Rev. Materials 2, 013807 (2018).
  * Shishkin  _et al._ [2007] M. Shishkin, M. Marsman, and G. Kresse, Accurate Quasiparticle Spectra from Self-Consistent GW Calculations with Vertex Corrections, Phys. Rev. Lett. 99, 246403 (2007).
  * Chen and Pasquarello [2015] W. Chen and A. Pasquarello, Accurate band gaps of extended systems via efficient vertex corrections in $GW$, Phys. Rev. B 92, 041115 (2015).
  * Cunningham  _et al._ [2018] B. Cunningham, M. Grüning, P. Azarhoosh, D. Pashov, and M. van Schilfgaarde, Effect of ladder diagrams on optical absorption spectra in a quasiparticle self-consistent $\mathit{GW}$ framework, Phys. Rev. Materials 2, 034603 (2018).
  * Cunningham  _et al._ [2023] B. Cunningham, M. Grüning, D. Pashov, and M. van Schilfgaarde, $\mathrm{QS}G\hat{W}$: Quasiparticle self-consistent $GW$ with ladder diagrams in $W$, Phys. Rev. B 108, 165104 (2023).
  * Radha  _et al._ [2021] S. K. Radha, W. R. L. Lambrecht, B. Cunningham, M. Grüning, D. Pashov, and M. van Schilfgaarde, Optical response and band structure of ${\mathrm{LiCoO}}_{2}$ including electron-hole interaction effects, Phys. Rev. B 104, 115120 (2021).
  * Jain  _et al._ [2013] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. A. Persson, Commentary: The Materials Project: A materials genome approach to accelerating materials innovation, APL Materials 1, 011002 (2013).
  * [26] Materials Project, mp-1023936.
  * Deguchi  _et al._ [2016] D. Deguchi, K. Sato, H. Kino, and T. Kotani, Accurate energy bands calculated by the hybrid quasiparticle self-consistent gw method implemented in the ecalj package, Japanese Journal of Applied Physics 55, 051201 (2016).
  * Sander  _et al._ [2015] T. Sander, E. Maggio, and G. Kresse, Beyond the Tamm-Dancoff approximation for extended systems using exact diagonalization, Phys. Rev. B 92, 045209 (2015).
  * Cheiwchanchamnangij and Lambrecht [2012] T. Cheiwchanchamnangij and W. R. L. Lambrecht, Quasiparticle band structure calculation of monolayer, bilayer, and bulk MoS2, Phys. Rev. B 85, 205302 (2012).
  * Marsili  _et al._ [2021] M. Marsili, A. Molina-Sánchez, M. Palummo, D. Sangalli, and A. Marini, Spinorial formulation of the $gw$-bse equations and spin properties of excitons in two-dimensional transition metal dichalcogenides, Phys. Rev. B 103, 155152 (2021).
  * Ramasubramaniam [2012] A. Ramasubramaniam, Large excitonic effects in monolayers of molybdenum and tungsten dichalcogenides, Phys. Rev. B 86, 115409 (2012).
  * Qiu  _et al._ [2013] D. Y. Qiu, F. H. da Jornada, and S. G. Louie, Optical Spectrum of ${\mathrm{MoS}}_{2}$: Many-Body Effects and Diversity of Exciton States, Phys. Rev. Lett. 111, 216805 (2013).
  * Qiu  _et al._ [2016] D. Y. Qiu, F. H. da Jornada, and S. G. Louie, Screening and many-body effects in two-dimensional crystals: Monolayer ${\mathrm{MoS}}_{2}$, Phys. Rev. B 93, 235435 (2016).
  * Chen  _et al._ [2019] S.-Y. Chen, Z. Lu, T. Goldstein, J. Tong, A. Chaves, J. Kunstmann, L. S. R. Cavalcante, T. Woźniak, G. Seifert, D. R. Reichman, T. Taniguchi, K. Watanabe, D. Smirnov, and J. Yan, Luminescent Emission of Excited Rydberg Excitons from Monolayer WSe2, Nano Letters 19, 2464 (2019).
  * Liu  _et al._ [2019] E. Liu, J. van Baren, T. Taniguchi, K. Watanabe, Y.-C. Chang, and C. H. Lui, Magnetophotoluminescence of exciton rydberg states in monolayer $\mathrm{WS}{\mathrm{e}}_{2}$, Phys. Rev. B 99, 205420 (2019).
  * He  _et al._ [2014] K. He, N. Kumar, L. Zhao, Z. Wang, K. F. Mak, H. Zhao, and J. Shan, Tightly Bound Excitons in Monolayer ${\mathrm{WSe}}_{2}$, Phys. Rev. Lett. 113, 026803 (2014).
  * Stier  _et al._ [2018] A. V. Stier, N. P. Wilson, K. A. Velizhanin, J. Kono, X. Xu, and S. A. Crooker, Magnetooptics of exciton rydberg states in a monolayer semiconductor, Phys. Rev. Lett. 120, 057405 (2018).
  * Liu  _et al._ [2021] E. Liu, J. van Baren, Z. Lu, T. Taniguchi, K. Watanabe, D. Smirnov, Y.-C. Chang, and C. H. Lui, Exciton-polaron Rydberg states in monolayer MoSe2 and WSe2, Nature Communications 12, 6131 (2021).
  * Wang  _et al._ [2020] T. Wang, Z. Li, Y. Li, Z. Lu, S. Miao, Z. Lian, Y. Meng, M. Blei, T. Taniguchi, K. Watanabe, S. Tongay, D. Smirnov, C. Zhang, and S.-F. Shi, Giant Valley-Polarized Rydberg Excitons in Monolayer WSe2 Revealed by Magneto-photocurrent Spectroscopy, Nano Letters 20, 7635 (2020).
  * Woo  _et al._ [2023] S. Y. Woo, A. Zobelli, R. Schneider, A. Arora, J. A. Preuß, B. J. Carey, S. Michaelis de Vasconcellos, M. Palummo, R. Bratschitsch, and L. H. G. Tizei, Excitonic absorption signatures of twisted bilayer WSe2 by electron energy-loss spectroscopy, Phys. Rev. B 107, 155429 (2023).
  * Arora  _et al._ [2015] A. Arora, M. Koperski, K. Nogajewski, J. Marcus, C. Faugeras, and M. Potemski, Excitonic resonances in thin films of WSe2: from monolayer to bulk material, Nanoscale 7, 10421 (2015).
  * Schmidt  _et al._ [2016] R. Schmidt, I. Niehues, R. Schneider, M. Drüppel, T. Deilmann, M. Rohlfing, S. M. de Vasconcellos, A. Castellanos-Gomez, and R. Bratschitsch, Reversible uniaxial strain tuning in atomically thin WSe2, 2D Materials 3, 021011 (2016).
  * Hong  _et al._ [2020] J. Hong, R. Senga, T. Pichler, and K. Suenaga, Probing Exciton Dispersions of Freestanding Monolayer ${\mathrm{WSe}}_{2}$ by Momentum-Resolved Electron Energy-Loss Spectroscopy, Phys. Rev. Lett. 124, 087401 (2020).
  * Wang  _et al._ [2015] G. Wang, X. Marie, I. Gerber, T. Amand, D. Lagarde, L. Bouet, M. Vidal, A. Balocchi, and B. Urbaszek, Giant Enhancement of the Optical Second-Harmonic Emission of ${\mathrm{WSe}}_{2}$ Monolayers by Laser Excitation at Exciton Resonances, Phys. Rev. Lett. 114, 097403 (2015).
  * Li  _et al._ [2014] Y. Li, A. Chernikov, X. Zhang, A. Rigosi, H. M. Hill, A. M. van der Zande, D. A. Chenet, E.-M. Shih, J. Hone, and T. F. Heinz, Measurement of the optical dielectric function of monolayer transition-metal dichalcogenides: ${\mathrm{MoS}}_{2}$, $\mathrm{Mo}\mathrm{S}{\mathrm{e}}_{2}$, ${\mathrm{WS}}_{2}$, and $\mathrm{WS}{\mathrm{e}}_{2}$, Phys. Rev. B 90, 205422 (2014).
  * Wierzbowski  _et al._ [2017] J. Wierzbowski, J. Klein, F. Sigger, C. Straubinger, M. Kremser, T. Taniguchi, K. Watanabe, U. Wurstbauer, A. W. Holleitner, M. Kaniber, K. Müller, and J. J. Finley, Direct exciton emission from atomically thin transition metal dichalcogenide heterostructures near the lifetime limit, Scientific Reports 7, 12383 (2017).
  * Madéo  _et al._ [2020] J. Madéo, M. K. L. Man, C. Sahoo, M. Campbell, V. Pareek, E. L. Wong, A. Al-Mahboob, N. S. Chan, A. Karmakar, B. M. K. Mariserla, X. Li, T. F. Heinz, T. Cao, and K. M. Dani, Directly visualizing the momentum-forbidden dark excitons and their dynamics in atomically thin semiconductors, Science 370, 1199 (2020).
  * Li  _et al._ [2020] J. Li, M. Goryca, N. P. Wilson, A. V. Stier, X. Xu, and S. A. Crooker, Spontaneous valley polarization of interacting carriers in a monolayer semiconductor, Phys. Rev. Lett. 125, 147602 (2020).
  * Deilmann and Thygesen [2017] T. Deilmann and K. S. Thygesen, Dark excitations in monolayer transition metal dichalcogenides, Phys. Rev. B 96, 201113 (2017).
  * Deilmann and Thygesen [2019] T. Deilmann and K. S. Thygesen, Finite-momentum exciton landscape in mono- and bilayer transition metal dichalcogenides, 2D Materials 6, 035003 (2019).
  * Song and Dery [2013] Y. Song and H. Dery, Transport Theory of Monolayer Transition-Metal Dichalcogenides through Symmetry, Phys. Rev. Lett. 111, 026601 (2013).
  * Atabek  _et al._ [1974] O. Atabek, C. Deutsch, and M. Lavaud, Schrödinger equation for the two-dimensional Coulomb potential, Phys. Rev. A 9, 2617 (1974).