Theory of Excitons in Atomically Thin Semiconductors: Tight-Binding Approach.

Atomically thin semiconductors from the transition metal dichalcogenide family are materials in which the optical response is dominated by strongly bound excitonic complexes. Here, we present a theory of excitons in two-dimensional semiconductors using a tight-binding model of the electronic structure. In the first part, we review extensive literature on 2D van der Waals materials, with particular focus on their optical response from both experimental and theoretical points of view. In the second part, we discuss our ab initio calculations of the electronic structure of MoS2, representative of a wide class of materials, and review our minimal tight-binding model, which reproduces low-energy physics around the Fermi level and, at the same time, allows for the understanding of their electronic structure. Next, we describe how electron-hole pair excitations from the mean-field-level ground state are constructed. The electron-electron interactions mix the electron-hole pair excitations, resulting in excitonic wave functions and energies obtained by solving the Bethe-Salpeter equation. This is enabled by the efficient computation of the Coulomb matrix elements optimized for two-dimensional crystals. Next, we discuss non-local screening in various geometries usually used in experiments. We conclude with a discussion of the fine structure and excited excitonic spectra. In particular, we discuss the effect of band nesting on the exciton fine structure; Coulomb interactions; and the topology of the wave functions, screening and dielectric environment. Finally, we follow by adding another layer and discuss excitons in heterostructures built from two-dimensional semiconductors.

1. General Overview

1.1. 2D van der Waals Materials

The fact that many crystals have a layered form with strong intralayer and weak interlayer bonds is well known [1,2,3]. Renewed widespread interest in such crystals reignited around 2005 due to the realization that high-quality few-layer crystals can be obtained via a mechanical exfoliation process [4,5,6]. Immediately, it was realized that other materials can be obtained using this technique [7], i.e., insulators, semiconductors, metals, and superconductors. More recently, large-scale data mining approaches established that from about experimentally known crystals, approximately 1800 should be exfoliable and thermodynamically stable [8]. Around 40 of them belong to the group called transition metal dichalcogenides (TMDs). In this family, contrary to graphene, a single layer is built from three planes of atoms: two of them are built out of chalcogen X, surrounding a symmetrically metal M plane (X-M-X structure). Usually, they have a trigonal (1T), hexagonal (2H), or rhombohedral (3R) polymorphic structure [9]. Around half of them in bulk form are metallic, and the other half are semiconducting [10], sharing the common theme that when thinned down, the electronic band gap increases. In several cases, the transition from metal to semiconductor is predicted. The semiconducting TMDs that sparked the greatest attention due to their promising optoelectronic properties [11] were compounds with formulas MoS, MoSe, MoTe, WS, and WSe.

What is worth noting, is that the properties of the monolayer crystals depend not only on their chemical composition and structural phase, but also on the choice of the substrate. When the first samples of TMDs were obtained, the standard procedure after their exfoliation was to place them on a properly prepared SiO/Si substrate, which for a certain thickness of the SiO layer enhanced the contrast in optical microscopy [12]. It was also understood that optical quality depends heavily on the surrounding dielectric environment [13,14], which posed the question as to whether SiO is truly the best substrate for optical studies. Currently, it is clear that the substrate on which the subtle optical features of TMDs can be most easily observed is hBN [15,16], for which, e.g., exciton lines approach ∼2 meV widths.

1.2. Optical Properties of TMDs and Their Heterostructures

One of the results that sparked widespread attention was the observation that MoS, when thinned down to a single layer, exhibits significantly stronger photoluminescence (PL) than -layer crystals [17,18,19]. A similar trend was then confirmed in MoSe [20], WS [21], WSe [20,21], and MoTe [22,23] semiconductors with varying magnitudes of differences between samples and with different numbers of layers. This phenomenon was related to the indirect–direct band gap transition, caused by different renormalizations of indirect () band gaps compared to the much smaller renormalizations of the direct gap when samples are thinned down. This trend was confirmed theoretically by ab initio calculations [24,25] and experimentally by photoemission spectroscopy for valence bands [26].

The first inkling that strong PL is dominated by an excitonic (X) transition comes from the comparison of the PL peak energy position to the known spectral positions of the excitonic transitions in bulk MoS at K point [17,27]. The generation of excitons in a given valley is possible via specific selection rules, which couple circularly polarized light with excitation in a given valley [28,29,30,31,32,33,34,35], (see Figure 1 for schematic picture). The valley coherence of such excitons has been probed in many experiments [36,37,38,39,40]. The spontaneous choice of valley after excitation with unpolarized light has also been observed [41], pointing towards a valley-polarized ground state [42,43,44]. Excited exciton states were measured in WSe [45,46] and WS [47] prior to molybdenum-based TMDs due to the larger separation between A and B excitons (∼0.4 eV), so that the B exciton was not “masking” the excited A exciton states. Complementary PL excitation spectroscopy allowed for the measurement of excited states in MoS, MoSe, WS, and WSe [48,49,50]. The excited B exciton series is more challenging to measure [46,49] due to overlap with the wide PL signal from so-called C excitons. All of those experiments confirmed that the spectrum of s-like states deviates strongly from the 2D hydrogen-like series. Further experiments in magnetic fields reached up to excited states [51,52,53,54,55], addressing also the complicated problem of the diamagnetic shifts of exciton lines resulting from three contributions (see Figure 1), as discussed theoretically in one of our other works [56]. The exciton series consists not only of bright s-series but also of two types of dark states—spin-forbidden and excited excitons states with symmetry different from ‘s’. The former one can be probed by PL dynamics [57,58], different emission configurations [59,60,61,62,63,64], the different temperature dependence of emissions [65], or by using a tilted magnetic field [60,66,67,68,69]. The latter ones are usually probed using two-photon spectroscopy [70,71]. Dark p states can also be studied in pump–probe experiments [72,73,74,75,76,77,78], where excitons are generated and transition from to states, which can be measured by a probe terahertz beam [79]. Additionally, excitons with finite center-of-mass momentum, which require phonon-assisted excitation due to momentum conservation can also be probed, revealing the complicated landscape of so-called momentum-indirect excitons [80,81,82,83].

Figure 1

Schematic structure of spin-split bands in monolayer MX/ valleys. Respective A/B excitonic transitions are shown as dashed arrows. Additional solid arrows denote different contributions to the and valleys’ responses to magnetic fields. The solid blue and red arrows show the bare electron Zeeman contribution in magnetic field ; the green—the atomic orbital Landé contribution and the black ones—the valley Zeeman contribution.

In addition to strongly bound excitons, TMDs exhibit another type of bound complexes in their photoresponse, namely trions [84,85,86,87]. In the simplified picture, this complex consists of three particles, i.e., two electrons + one hole (X) or two holes + one electron (X) for negatively and positively charged types, respectively. In GaAs quantum wells, the rule of thumb is that trion-binding energy is an order of magnitude smaller than that for exciton, and the same rule seems to work in TMDs. Because the binding energy of the exciton is on the order of hundreds of meV in TMDs, the trion-binding energy is given in tens of meV [88,89,90,91,92,93,94,95,96,97]. The existence and properties of such complexes depend heavily on the presence of excess carriers [13,36,91,98,99,100,101,102,103,104,105,106,107,108], and their stability has been shown up to room temperatures [88,97]. For higher doping regimes, a picture of an exciton interacting with the Fermi sea has been proposed [43,109]. Interestingly, just as in the exciton case, trions exhibit a fine structure [107,108,110,111,112,113,114,115,116], involving both dark/bright states and also an excited trion states [117]. Moreover, a strong electron–phonon interaction allows to couple trions and excitons by an up-conversion process [118,119], in which low-energy photon excitation (exciting trion) induces an optical emission from a higher energy state (exciton). It is worth noting that larger optical complexes have also been experimentally studied [104,110,120,121,122,123,124,125,126,127,128,129,130,131,132]. Their interpretation is, however, non-trivial due to the various possible couplings with other quasiparticles, e.g., related to excess charge carriers [86], collective excitations of electron gas [133,134], and polarons [43,135].

Atomically thin crystals offer a novel route to engineering material properties due to the possibility of assembling them into configurations not present in nature [136]. Starting from proof-of-concept double layers of graphene with hBN layers between them [137], it has been shown that graphene can act as an atomically sharp contact to TMDs [138]. Recently, one of the greatest improvements of the optical quality of semiconducting crystals was achieved by placing them in-between thin hBN crystals [16,139]. Different band gaps and band edge positions allowed for creating analogues of type I and type II heterostructures [140,141]. In such systems, we can distinguish two types of excitons: intralayer excitons, localised in one of the layers building the heterostructure, and interlayer excitons (sometimes called indirect), built by electrons and holes from distinct layers. Interlayer excitons have been observed [141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158], showing the electrical tunability of both PL and excitons energy [156,159,160,161]. Order-of-magnitude-longer lifetimes than intralayer excitons were observed due to the spatial separation of electron and hole [142,158,162]. Unlike interlayer excitons in GaAs heterostructures, for TMD heterostructures, we deal with interlayer excitons stable at room temperatures [160,163,164,165]. Larger complexes with carriers residing in different layers have already been proposed and some evidence of such behavior has been found [156,161,166,167,168,169]. Long interlayer exciton lifetimes naturally led to the proposal of exciton Bose condensates [164,165], with the first evidence in the MoSe-WSe system [170]. It is worth noting that semiconductor heterostructures offer new means of control by the twist angle, an area that has been actively pursued in recent years [148,159,161,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192].

1.3. Review of Theory of Correlated Optical Excitations

The successful description of bound excitonic states requires an analysis of two-particle correlations, described by the Bethe–Salpeter equation (BSE) [193,194,195,196,197,198,199,200]. Currently, several ab initio codes implement the DFT + GW + Bethe–Salpeter framework [201,202,203,204,205,206,207,208]. First-principles Green’s function theory of optical excitations can also be applied to materials with reduced dimensionality. Many studies of excitons in 2D semiconductors within DFT + GW + BSE have emerged [25,70,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225]. Although common difficulties with the numerical aspects of excitons in TMDs [211,221] are known, those studies correctly identified that excitonic effects in TMDs should be large [25,209,210] due to the strong electron–hole interaction (weak screening) compared to bulk crystals. The inclusion of the spin-orbit coupling (SOC) confirmed the general double-peak structure of the optical response (A and B excitons), associated with the valence band split-off [226], as well as the fine structure related to the spin arrangement in the conduction band [227] and momentum-dark excitons [228]. The first ab initio calculations of the excited states of exciton series predicted large deviations from the non-hydrogenic series, the existence of – splitting in the first shell [70], and the topological splitting of the states [78]. Non-analytic (Dirac-like) dispersion of exciton-excited branches due to inter-valley coupling has been predicted [218,219]. The mixing of A and B excitons via the inter-valley exchange interaction has shown that understanding them in terms of “Ising” excitons may not always be satisfactory [225]. Interestingly, the DFT + GW + Bethe–Salpeter methodology has been recently extended to study charged exciton complexes [117,229,230,231].

Large differences between the binding energies of excitons on different substrates have been early established [215,217,220]. Interestingly, MoS put on the two different, most popular substrates, SiO and hBN, which have almost exactly the same energy in the optical gap. This is surprising because it is known that the binding energies of excitons on those substrates should be significantly different. Another confusing result is that the absolute peak position (G0W0 + BSE), almost identical to the experimental, is calculated for MoS in vacuum. The similarity of these values is coincidental and stems from the cancellation of self-energy and excitonic renormalizations, which are true even when the dielectric environment (vacuum, SiO, hBN) is vastly different. Recent GW calculations of MoS on hBN [232] help to rationalize this, showing the actual position of the GW-normalized “free particle” band gap.

Similarly to the case of monolayers, many theoretical studies for TMD heterostructures emerged as well [223,226,233,234,235,236,237,238,239,240,241,242]. Bright interlayer excitons were predicted [236] with A and B exciton structures [223,226,233,236,237,238] analogous to the monolayer. It is worth noting that exciton binding energy is strongly correlated with the stacking and emergence of moiré patterns [233,237,238,239,240,241]. Furthermore, both the temperature dependence of the exciton fine structure [223,236,238] and the electrical control of exciton energy [234] were analyzed. Additionally, the spectrum of excited states was predicted [223] and the exciton lifetime was determined [241,243].

Due to many computational challenges inside DFT + GW + Bethe–Salpeter theory, simplified approaches to the excitonic problem have been studied, including various levels of approximations, i.e., tight-binding [79,80,82,83,244,245,246,247,248,249,250,251], k·p (massive Dirac fermion and beyond) [28,92,218,252,253,254,255,256,257,258,259], and effective mass approximations [47,53,139,260,261,262,263,264,265,266,267,268,269,270,271]. Those approximations routinely use DFT as a staring point for the calculation of the effective mass of carriers and then simplify somehow their interaction, e.g., assuming only Rytova–Keldysh screening [272,273] without form factors coming from the Bloch part of the electron and hole Bloch wave functions. These methods, being conceptually and numerically traceable, allowed for many fascinating advances in understanding exciton properties in TMDs. For the case of the heterostructure, k·p models were also developed, allowing studies of intra- and interlayer exciton properties [77,274,275,276,277,278,279,279]. Further simplification of the k·p methodology to the massive Dirac fermion model was used to analyze intra- and interlayer exciton properties [77] and the Bose–Einstein condensate of excitons [277,278,280].

Due to complicated numerical and conceptual challenges inside DFT + GW + BSE theory and the somehow approximate nature of the simplified methods mentioned above, one may turn to a tight-binding methodology, offering an optimal trade-off between numerical tractability and the ability to capture essential physical ingredients in the excitonic problem. Tight-binding studies [79,80,82,244,248,249,251,252,281] allowed, e.g., to track exciton fine-structure evolution dependence on electron doping, the exchange interaction, and magnetic and electric fields. Theoretical results on excitonic landscapes have been shown [83,247], together with results confirming the topological splitting of states in the exciton spectrum due to the effect of the Berry’s curvature. In our own work [281], as further described in more detail, we studied the role of band nesting, screening, and wave function topology on exciton fine structure.

2. Electronic Structure of MX Semiconductors

2.1. Ab Initio Insight into Electronic Structure

In the following subsection, we review our understanding of TMD monolayers based on ab initio techniques that we used in several recent works [24,56,281,282]. Single-layer MX crystals in a 2H phase are arranged in a trigonal prismatic structure. The unit cell consists of three atoms, one metal and two chalcogens, as shown in Figure 2. From the top, the arrangement of the atoms reminds one of graphene. From the side, one can notice that the atoms are actually organized in three planes, one (central) metal plane and two chalcogen planes shifted by . The distance between the metal atom and the central position between the two chalcogens is denoted by . We define the primitive vectors of the real space lattice as and . The lattice constant can be written as and . Following the commonly practiced ab initio procedure, described in detail elsewhere [282], the lattice constants calculated using the PBE exchange-correlation functional for MoS are Å. Reciprocal lattice vectors satisfy the relation . This gives a 2D set of four equations for . The solution yields and . We note that such choice of real and reciprocal space-primitive vectors gives the following position of the hexagonal Brillouin zone (BZ) corner at , which is called the K point, just as in graphene. All other K points are obtained by the successive application of the rotation. We note that in half of the distance between the K points and the point there are so-called Q points (alternatively [251] or [268] points). Halfway between the two nearest K points lies the M points.

Figure 2

(a) Top view of the structure of MX in the 2H phase: metal atoms are denoted by blue dots, and chalcogen by red ones. (b) Side view of MX, showing that the atoms are organized into three layers, central metal and two chalcogen, with structural constants parametrized by and .

Now, we describe the general features of the DFT band structure, focusing on a representative example of MoS. In Figure 3, one can observe 11 bands around the Fermi level (set to 0). The direct band gap at the K point without SOC is equal to eV. We note that both the valence band (VB) and the conduction band (CB) have secondary extrema, at the point in the VB and at the Q points in the CB. Next, we study the general properties of the DFT Kohn–Sham wave functions. First, we choose three spheres around one Mo and two S atoms, and calculate how much of the wave function is localized inside each sphere. We note that 11 bands around the Fermi level have wave functions that are in general localized on both Mo and S atoms; however, there is a clear trend that VB and bands above are mainly localized on Mo atoms, and VB-1 and lower bands on S atoms. DFT wave functions inside spheres can be next projected onto Slater orbitals. We first confirm that the largest overlap is achieved, in case of MoS, with localized Slater, single-zeta-basis-consistent [283,284] with orbitals with a principal quantum number on Mo and on S. For other compounds, e.g., WSe, there is a clearer and larger numerical overlap with orbitals with quantum numbers for W and for Se, confirming that such projection is a reliable tool in studying symmetry and the orbital composition properties of DFT wave functions. In the next step, summarized in Figure 3, for each wave vector and energy E, the wave functions are projected onto symmetric and anti-symmetric orbitals with respect to the metal plane. One can observe that VB and CB are symmetric across most of the Brillouin zone, with the exception of parts close to the point in the CB. This ordering of bands and their symmetry is general for all MX-family compounds. In the next step, we study the orbital-resolved decomposition of wave functions. The results confirm that orbital composition close to the K point in the VB bands is a combination of and orbitals. Contrary to VB, in the CB, the state coupled to the state dominates. The higher symmetric conduction band (CB + 1) () is coupled to the state. Coupling between anti-symmetric d and p states is also obtained, in which the state is coupled to anti-symmetric and to both anti-symmetric and states.

Figure 3

Color-mapped localization of a given k-resolved eigenenergy on Mo and S spheres and symmetry of eigenvalues across the Brillouin zone. Circles (crossed rectangles) denote symmetric (anti-symmetric) orbitals with respect to the metal plane.

Due to the heavy nature of the atoms in the MX family, one can expect that a relativistic spin-orbit interaction may influence the band structure. Here, we discuss briefly the results of the PBE + SOC calculation for MoS. The main effect of SOC is visible in the VB at K point, where bands are spin-split by meV. It is worth mentioning that this is the smallest splitting in the VB at the K point in MX. Significantly smaller splitting is observed in the CB at K point in MoS ( meV). This less-pronounced order-of-magnitude effect comes from the fact that the majority of () d-like orbitals in the CB have the quantum number , therefore splitting must come from the admixture of p-orbitals () from chalcogen atoms. What is also important is that there is a significant spin splitting in the second minimum of the CB at Q point, which is larger than at K point, pushing the relative distance between the CB minimum at K at Q point () close to small values, making the CB almost degenerate, e.g., in WSe. Let us also note that by symmetry there is no spin splitting along the line in the BZ.

2.2. Minimal Tight-Binding Hamiltonian

In the following subsection, we review our construction of a tight-binding model for TMDs monolayers [56,282]. Building on our DFT analysis presented in the previous subsection, we conclude that there are several common features in the band structures of analyzed MX semiconductors. They have a direct band-gap at K points, opposite to their -layer form, which are indirect-gap semiconductors [24]. These materials also have secondary minima in the conduction band, localized close to the Q points, and secondary maximum in the valence band at the point. Analyzing orbital compositions, we conclude that CB and VB can be described by orbitals, even with respect to the metal plane. This result is consistent with several other works [285,286,287,288,289,290,291,292,293,294,295]; therefore, we can start building a tight-binding model of those materials from the orbitals contributing most to the band structure around the Fermi level.

Our conclusion is that the minimal tight-binding model has to include orbitals that are even with respect to the inversion symmetry about the z plane of the metal atoms (). Motivated by ab initio results, at least d—orbitals from metals—and symmetric p-orbitals from the top and bottom sulfur atoms must be taken into account. We note only that for the d orbitals of metals the situation is clear, namely orbitals with a given and are centered around M atoms. However, the orbital construction for two chalcogens must be performed with care. Because we define the so-called dimer orbitals, which are centered around the same plane as metal atoms, we begin with upper (U, with respect to the dimer center) and lower (L, )) p-orbitals with quantum numbers , and define dimer orbitals as a proper combination of those two, as described in detailed elsewhere [56,282].

In the next step, we construct the standard linear combination of the atomic orbitals for all of the orbitals discussed above. First, we note that due to our dimer orbital construction, it makes sense now to talk about the sublattices A (metal positions, inside the unit cell) and B (chalcogen dimer centers, ) on the real-space hexagonal lattice. Because we aim to understand the fundamental properties of those materials in terms of graphene physics, we focus on the sublattice A–sublattice B interaction. First, let us try to understand the physics of why the MX materials considered are semiconductors instead of semimetals, such as graphene. To achieve that, let us analyze the tunneling matrix element between the central A atom with the potential and the three nearest neighbor B atoms at positions . One can notice, that exactly at the K point, this formula gives where are standard Slater–Koster integrals. For such combinations of and quantum numbers, that tunneling matrix element is non-zero . This is different from the same tunneling matrix element in graphene, in which only , orbitals play a role and give . Non-vanishing tunneling will, therefore, open a gap at the K point, resulting in a gapped, massive Dirac fermion dispersion instead of the Dirac point at K. Non-zero interaction between different and orbitals at K point leads to the following pairs of orbitals that are coupled: , , and . Similar analysis at point leads to the selection rule . The resulting couplings are then , , and . Finally, at point, we obtain a different scheme of couplings: , , and . We note that all those couplings explain the couplings detected in DFT.

In next step, we move to a discussion of the tight-binding Hamiltonian. The standard procedure outlined by Slater and Koster (SK) [296] is used. All details of this straightforward derivation are discussed elsewhere [56,282]. We note only that because in our basis we are using complex orbitals, SK formulas cannot be used directly and we had to implement them into a more complicated fashion, which is hidden in our notation behind functions , depending on the SK parameters. The next-nearest neighbor Hamiltonian can be written in block form as: where the matrix describing the metal–metal next-nearest neighbor interactions is given by and the corresponding matrix describing the X–X dimer interactions is given by Finally, metal dimer tunneling is described by The functions f and g, depending on the wave vector , combine the proper sum of the plane wave functions to the nearest and next-nearest neighbors. They are listed in Ref. [56]. As discussed in one of our works [56], we find that the next-nearest neighbor Hamiltonian is a minimal model that correctly describes the band gap across the whole BZ and is able to quantitatively reproduce the orbital compositions of the VB and CB. We note that the last column (except the diagonal element) in the Hamiltonian has a sign opposite to the one used in Ref. [56] due to two different conventions of complex spherical harmonics that could be used (Condon–Shortley phase). Here, complex and real orbitals are related by .

After the derivation of our tight-biding Hamiltonian, we turn to the problem of fitting the SK parameters. Our goal is to obtain dispersion for even bands from TB as close as possible to even bands obtained and detected using the ab initio method, as described above. We note that in principle it is possible to calculate SK parameters directly from a first-principles calculation [297]; however, the resulting electronic dispersions are not satisfactory, neither for graphene [298], nor for MoS [289], especially for unoccupied electronic states (CB and higher bands). Those parameters could be in principle treated as a starting point for our analysis; however, we take another route. The general problem with SK parameters is that even if we fix the structural parameters and , we still end up with a 10-dimensional, highly non-linear optimization problem, depending on the energies , , , , , , , , , and . We have tested various schemes to choose those parameters to fit the dispersion of such a model to the DFT [282]. At the end of the day, we used the one closest to the Monte Carlo philosophy of random probing of multidimensional parameter space.

The results of two fitting procedures are presented in Figure 4. We note that the fit presented in the left of Figure 4 was found by assuming equal weights for all k points and all even bands; therefore, it is called “best all bands”. Additionally to reproducing the overall energies of even bands well, this parametrization reproduces VB very well. To obtain the best transition energy (right panel of Figure 4), a much more complicated procedure was used, first, increasing the weights for VB and CB on the entire line, then, in subsequent sweeps, the weights around the K point in the VB and CB and the Q point in the CB were further increased. Thanks to those weighting procedures, both the VB and CB along the line are reproduced very well.

Figure 4

Left: TB dispersion obtained after optimizing the SK parameters to reproduce all even DFT bands. Right: TB dispersion optimized to reproduce the transition energy between the VB and CB. We note that the former reproduces the VB very well, while the latter one reproduces the CB very well, especially on the line.

To include spin-orbit coupling in our tight-binding model, we analyze the matrix elements of the operator. It turns out that only non-zero spin-orbit operator matrix elements are diagonal in our basis and read (for spin-up) The full Hamiltonian with SOC can be written therefore as: The parameters and have to be chosen in order to reproduce the spin splitting of the bands. This choice in general changes the dispersion of the bands due to the modification of the diagonal parts of the Hamiltonian and, in principle, requires additional fitting. We checked that setting eV and reproduces the correct splittings in the VB (0.148 eV) and CB (0.003 eV). We note that the order-of-magnitude-higher value of stems from a different contribution of orbitals to CB (order of ).

To better understand how SOC affects the band structure across the BZ, let us analyze spin-split bands along the direction. Our choice of parameters reproduces spin splitting in the VB and CB at K points well and our model catches spin inversion between the and points. Interestingly, the general feature of the whole MX family is the spin inversion of bands in the CB close to the K point, taking place between K and Q points. This feature is better visible when the lowest spin-split band is shown throughout the BZ, as shown in Figure 5a,b. For example, in the point, the bottom of the CB has a spin orientation the same as the top of the VB. However, the region of the BZ where this property holds is very small (red region around in Figure 5b), and quickly, the other spin becomes lower in energy. This situation changes again approximately half-way between the K and Q points, around which the spin is again oriented in the same way as at the point in the VB. Interestingly, in both Mo and W-based TMDs, the spin at Q point is always oriented the same way as in the VB at K point, irrespective of the spin ordering change in W-based materials. The same spin orientation between the VB at K and Q at the CB means that all momentum-indirect excitons with momentum will have a spin-“bright” configuration and, when activated, e.g., by phonons, they should be optically detectable, in contrast to the spin-forbidden lowest excitons that have a small momentum around K points in tungsten-based TMDs.

Figure 5

(a) Spin texture of the CB, showing the spin orientation of the lowest CB spin-split band. (b) Corresponding VB spin texture.

In the next step, we discuss the low-energy Hamiltonian around K point. We begin by noticing that at K point, the top of the VB is built almost solely from the orbital, while the bottom of the CB is built from a combination of and orbitals. However, assuming the low-energy basis as and , and expanding and functions around K point (), one can immediately end up with a massive Dirac fermion (mDF) Hamiltonian [299]: The best parametrization for dispersion including up to 1/4 distance between the K and Q point we found is eV, Å, eV. The simple mDF model defined by Equation (8) gives very similar results of dispersion compared to models with further corrections, such as trigonal warping, and when applied to excitonic calculations there is no qualitative and very little quantitative difference. We note that mDF can be even further reduced to the parabolic (effective mass) model. This can be performed by noting that eigenergies of mDF are: Keeping only the first order of , we end up with The effective electron mass is given as . We note that the choice of parameters depends heavily on how large a portion of the BZ we want to fit as closely as possible. One can also note that the top of the VB and the bottom of the CB are described better if different effective masses of electrons () and holes () are taken. For MoS, when using the parabolic model, we take and .

As a summary, in Figure 6, we present the ab initio result and different low-energy models of the VB and CB: the tight-binding, massive Dirac fermion, and parabolic models. One can observe that all of them describe well the neighborhood of K point; however, by only using the TB model it is possible to obtain a second minimum at the Q point in the CB and a correct second maximum of the VB at the point. Approximately 10% of the line from the K point is properly described by both mDF and parabolic models, with the mDF model being better for the CB description. We note also that both mDF and parabolic models can be extended to include spin; however, in both cases it is necessary to find a different parametrization for two spin species (e.g., two effective electron and hole masses and two gap parameters for the parabolic model). Such spin-dependent low-energy model parameters for MX semiconductors are also available in the literature [300].

Figure 6

Comparison between dispersion models along the line. DFT dispersion is denoted by black circles, TB by blue rectangles, massive Dirac fermion by green diamonds, and parabolic (effective mass) model by red triangles. The corresponding connecting lines are shown as a guide to the eye.

3. Tight-Binding Theory of Optical Excitations

3.1. Bethe–Salpeter Equation

The simplest picture of light absorption by semiconducting material can be understood in terms of the transitions of carriers from the valence to the conduction band due to photon excitation with energy , changing the angular momentum by when circularly polarized light is used. In MX semiconductors, the situation becomes more complicated, because dipole transitions between d orbitals ( in the VB and in the CB) and negligible p–d transitions cannot explain that circularly polarized light excites carriers within one valley. The solution to this problem comes from a realization that the velocity matrix elements inside the transition matrix elements in general have two contributions [301]: dipole transitions between localized orbitals and terms related to electron hopping between lattice sites. This hopping contribution sets the phase of the velocity matrix element between the VB and CB Bloch wave functions and defines the optical selection rules in the gapped “chiral” fermion systems [302,303].

Because it is possible to selectively excite carriers from the band edge of one valley (e.g., ), in theoretical investigations it is useful to think about two non-equivalent parts of the hexagonal BZ. A unique association of k points that belong to one or the other valley has to be performed, as shown in Figure 7. Starting from cutting “wedges” around three equivalent (related to each other by reciprocal lattice vector G translations) K points, as shown in Figure 7a, one can move respective wedges to one neighborhood of the K point, creating a triangle around it, as in Figure 7b. We note that a similarly constructed triangle for the valley has to be rotated by symmetry, and both triangles for and valleys placed next to each other create a rhomboidal BZ equivalent to a full hexagonal BZ.

Figure 7

(a) Choice of valley on the whole BZ. (b) Construction of the valley around a single K point.

Exciton properties can be calculated from the configuration-interaction approach to interacting electrons [304,305], particularly useful in studies of nanostructures, e.g., quantum dots [306,307,308,309,310,311,312,313,314,315,316,317,318,319,320,321,322] or electron gas in quantum wells in strong magnetic fields [323,324]. This approach has only recently been utilized to study excitons in reciprocal space [230,281].

Now, we give a brief summary of steps performed for the derivation of the Bethe–Salpeter equation for the bound electron–hole problem, shown schematically in Figure 8a. We note that the complete derivation is presented in detail in Refs. [281,282]. We assume that we have only one valence and one conduction band. Next, to construct the ground state, we fill all states in the valence band, as shown in Figure 8b. Then, the exciton state is formed as a linear combination of excitations with coefficients being complex electron–hole amplitudes. Exciton states can be calculated using the standard eigenvalue problem , where is an interacting excitonic Hamiltonian. After calculating the matrix elements of this Hamiltonian using the Wick theorem, neglecting ground state energy correction due to interactions and incorporating electron and hole self-energies into dispersion energies and , we obtain the Bethe–Salpeter equation for the exciton (for center-of-mass momentum ) In this equation, the summation over the states is understood as over all vectors in the first BZ, the same as the number of unit cells in the crystal. Electron–hole interaction matrix elements are generally defined as:

Figure 8

(a) The single electron–single hole picture (exciton in effective mass approximation) in which interaction creates a spectrum of bound states. (b) Exact picture where the “hole” is created by exciting the electron from the filled ground state in the VB. The exciton is then constructed as a coherent superposition of all possible excitations for a given center-of-mass momentum interacting via the Coulomb interaction.

We note that in our Equation (11), we expand the exciton wave function only in electron–hole pair excitations from the VB to CB and not reverse (just as Equation (16) in Ref. [200]). Such an expansion is called the Tamm–Dancoff approximation. If the interaction matrix elements are not screened, Equation (11) is called time-dependent Hartree–Fock. In our work, we always used screened versions of matrix elements (see discussion below); therefore, in our case, formally, we should be using the name simplified Bethe–Salpeter (simplified in the sense of using an approximate, not self-consistent, screening). In our approach, based on the configuration interaction picture, it is natural to use the same screening for both direct and exchange electron–hole terms, producing results that are consistent with other similar works [230].

We note that in general the issue of the screening of the exchange electron–hole interaction is a long-standing problem, with several works discussing its various aspects (see, e.g., references in Ref. [230]). It is understood from a formal point of view that the iterative solution for the two-particle Green’s function ( is the non-interacting Green’s function describing the independent propagation of electron and hole, K is interaction kernel) does not allow the screening of the exchange electron–hole part of K, as discussed in Ref. [325]. However, there are reasons to screen it just as one direct stemming, e.g., from the finite Hilbert space of the problem when the Bethe–Salpeter equation is solved in practice. The screening of exchange in DFT + GW + Bethe–Salpeter has recently been discussed in the context of excitons in MoS monolayers in Refs. [326,327]. In an alternative approach to this problem, in paper [328], some of us expanded the systematically excited states of a graphene quantum dot in multi-electron–hole pair excitations. Figure 7 shows the evolution of ground and excited states as a function of the number of excited pairs. Figure 9 shows the singlet–triplet splitting. The main result is obtained by including an extra two pairs which do the screening, i.e., produce an effective interaction for singlet and triplet excitations.

Figure 9

Graphical representation of two types of interaction between an electron and a hole: (a) direct process (b) exchange process.

3.2. Coulomb Matrix Elements

Let us identify the type of matrix elements within the summation over in Equation (11). Note that they are written in electron-only language. As shown in Figure 9, the first element (with the minus sign) describes the process in which electrons in conduction bands in state and electrons in valence bands in state scatter via the Coulomb interaction to electrons in state in conduction bands and electrons in state in valence bands. This description is equivalent to the electron–hole pair scattering from state to . The second process in Equation (11) describes electrons starting as previously in the state in the conduction band and the state in the valence band and scattering to the same and , but changing band indices to valence and conduction, respectively. We identify the first process as an attractive direct electron–hole interaction, and the second process as a repulsive-exchange electron–hole interaction.

Next, we analyze direct electron–hole Coulomb matrix elements. Substituting the Bloch form of the wave functions and utilizing the fact that functions of coordinates in plane can be analyzed in reciprocal space (interaction and co-densities), we use their Fourier components as a general strategy to re-group six-dimensional integrals. The final expression for the direct matrix element (with coefficient resulting from sum to integral transition) is where and interaction form factor is given by: Pair densities can be evaluated numerically for every coordinate z by using the explicit form of the Bloch wave functions, constructed using localized Slater orbitals , as described before, as We note that further discussions of the details of the properties of the direct Coulomb interaction, e.g., the behavior of pair densities, convergence issues, matrix elements symmetries, complex phase properties, and descriptions of their inter-valley behavior, are described elsewhere [281,282].

Analogously, the exchange Coulomb matrix elements are expressed by with form factors There are several differences between direct and exchange electron–hole interactions. As can be seen from Equation (11), direct interaction comes with a negative sing (electron–hole attraction) and an exchange one comes with a positive sign (electron–hole repulsion). Contrary to direct matrix elements, there is no overall dependence of magnitude of those elements; therefore, potentially exchange interaction can be large whenever direct matrix elements are small due to the large distance. We note also that interaction form factors consist of a pair of densities diagonal in band indices () and off-diagonal with wave vector indices (), while for exchange-interaction form factors, the situation is opposite, i.e., there is off-diagonal dependence on band indices () and diagonal wave vector dependence (). Due to those properties, we found that direct matrix elements are generally complex numbers, while exchange matrix elements are real. Additionally, due to dependence only on the diagonal wave vector, the exchange matrix element can be computed much faster than the direct matrix elements. At this point, we note also that at first glance there is singularity in . This singularity is generally problematic for 3D crystals and different methods of dealing with it have been discussed in the literature [329]. However, for a 2D crystal, the singular term is equal to 0 and can be excluded from the summation over G vectors in Equation (16).

We note that the Bethe–Salpeter equation defined in Equation (11) is gauge-invariant, i.e., it does not depend on the gauge choice of wave functions, which can be created arbitrarily for every band and k point. On the other hand, the numerical phase affects the phase of the exciton coefficients. This arbitrariness of phase choice may influence the apparent symmetry of the excitonic ground state [253]. One of the solutions to this problem is to study only observable quantities, such as the imaginary part of the dielectric function , in which complex exciton amplitudes are multiplied by velocity operator matrix elements, which are also constructed from wave functions. Keeping track of the phase in both should be enough to obtain a gauge-independent answer. However, sometimes one is interested in studying excitonic wave functions themselves. To obtain ground excitonic states with symmetry, various numerical approaches are used; however, their details are rarely discussed in the literature [200,253]. In our procedure, we follow the idea introduced by Rohlfing and Louie [200], in which the global phase of tight-binding wave functions is chosen in such a way that the sum of imaginary parts is 0, i.e., . In the next step, we rotate the global phase to obtain phase 0 in the orbital, which means that the phase of the second TB coefficient is set to zero (). The second step breaks the first property, although we found that it is necessary to numerically obtain excitons in the and valleys, the wave functions of which have properties, as expected from time-reversal symmetry arguments.

Just as for the exchange Coulomb interaction discussed above, it is easy to understand that the diagonal term of the direct electron–hole interaction in BSE given by Equation (11) is singular at and . Renormalization due to this singularity has to be included in simulations on the finite lattice, otherwise numerical results give vastly different results when compared with theoretical predictions. Dealing with this singularity is connected with the discretization of the BZ associated with a single valley. Below, we discuss the rectangular discretization of the lattice. Confining our discussion to the case of for a moment, one of the methods that allow to integrate out the singularity is to note that form factor and make an approximation that the exciton wave function inside the box centered around the point (,) takes a constant value. This approximation is naturally more and more exact with an increasing number of points (and a decreasing area associated with each k point), into which the BZ is discretized. This allows us to approximate the diagonal of the BSE interaction kernel as We note that this result assumes a constant, static screening of the electron–hole interaction and in principle should be recalculated when the k-dependent screening model is used. Additionally, the term is the leading one; however, summation over G vectors introduces further corrections into the singular term. We checked numerically that both effects introduce contributions that are at least two orders of magnitude smaller than the expression given in Equation (18) and are neglected in further studies. We note that the constant is calculated for the rectangular lattice with an area of BZ around the single point given by , and it should be changed to for a rhomboidal lattice with an analogous rhombus area given by under the same approximations as for the rectangular discretization scheme.

In the next step, we build a systematic theory of approximation for the direct electron–hole matrix elements. The calculation of the interaction form factors is a major bottleneck of both ab initio and TB calculations due to the necessity of calculating them for all combinations of k and , and the additional summation over reciprocal lattice vectors G. The easiest solution is to assume all form factors to be one (their highest possible value, exact for and ) and note that the highest value-entering sum over G vectors comes from Gs minimizing distance. We checked that this approximation is useful around the K point, for which the neighborhood-form factor absolute value deviates from one rather slowly. The approximation described above allows us to reproduce numerically analytical solutions for a simple model of the exciton with an electron/hole dispersion in a parabolic approximation. Its main deficiency is that it does not include any effects related to the orbital composition of the bands (no Bloch function effect) and is purely real, which always gives degenerate states in the exciton spectrum with the same exciton angular momentum quantum numbers. To motivate a way around this problem, let us first discuss which parts of the full, complex form factor influence mainly its value. For concreteness, let us discuss, for example, the difference between a form factor describing the scattering between one chosen k point to some point ( is assumed to be small). One can quickly deduce that the effect of on the matrix element is complicated: it affects (1) the denominator ; (2) tight-binding coefficients inside ; (3) exponent values depending on z, z′; and (4) details of the in-plane integration of Slater orbitals and the summation over unit cells. We checked numerically, implementing in code the selective turning on/off of all the above contributions so that, actually, the first two corrections (denominator and TB coefficient) yield a very good approximation to the matrix elements and the last two (exponent and details of integration) do not contribute too much. This motivated us to extract TB coefficients from form factors, giving an expression formally equivalent to Equation (14) with an implicit summation over sublattices () and orbitals (), which can be written as and the quantity we call “orbital form factor” is given by Those orbital form factors depend on analogues of product densities, now in the form that does not depend on TB coefficients. This method, being equivalent to Equation (14), is not faster; however, it helps to realize that the coefficient can be calculated very fast from the TB model and we can take only the orbital form factor at the limit. This approximation is conceptually equivalent to treating the long-range Coulomb interaction as not dependent on the details of the orbital structure and taking the maximal absolute value of the orbital form factor (=1). Interestingly, taking this limit also produces a phase , related to non-zero G and the position of atoms inside unit cell: The final equation for the direct form factor is, therefore, simplified to This expression is an extension of results presented in Refs. [245,247,250,330] to include summation over non-zero G vectors. We note that the phase rotation presented in Equation (22) is either 1 or rotation equal to .

As is known from the literature and from our own experience, the computation of the excitonic spectrum is a numerically challenging task. To lower computational complexity even further, let us discuss now how excitonic fine-structure calculations are performed in practice. First, we discretize Equation (11), neglecting the electron–hole exchange interaction, since it is much weaker than the direct electron–hole interaction. Additionally, we choose in summation over wave vectors only those associated with one valley, as shown in Figure 7. Remembering about singular terms as discussed previously, the Bethe–Salpeter-like equation for one valley takes the form This equation represents a dense, Hermitian matrix that is diagonalized numerically. The primary convergence parameter is the number of k vectors into which a single valley was discretized. Details of computations and convergence studies are presented in our works [281,282].

A further step in fine-structure calculations is to add spin splitting to both valence and conduction bands. Spin-splitting modifies the electron–hole energy difference . To calculate matrix elements, we choose to use spinless wave functions due to their negligible dependence on spin. By this method, we are able to obtain a fine structure in one valley, i.e., both bright and dark A and B excitonic series. The calculation of fine structure in one valley automatically gives fine structure in the other one due to the symmetry of energies for spin-flipped configurations of excitons. We also found that with a proper gauge of matrix elements, the property of the matrix elements between valleys is satisfied. Implementing this symmetry in BSE, one can formally prove that exciton wave functions have to be related to each other as . We checked numerically, performing separate, full calculations in and valleys that our implementation satisfies these properties, giving within numerical precision the same excitonic spectrum and phases of excitonic states related by the mirror symmetry of k vector and complex conjugation.

3.3. Screening of Coulomb Interactions

An additional complication in the realistic description of the Coulomb electron–hole interaction stems from the screening of the interaction by other carriers. As noted early in studies of graphene [331], electron–electron interaction screening in 2D crystals behaves differently than in 3D crystals. Therefore, it is necessary to use a more involved screening model. The most well known model that captures the major linear dependence of screening in k-space was derived in the physics of thin dielectric slabs and is called the Rytova–Keldysh (R.-K.) model [272,273,331]. This approximation has been used in numerous works [53,260,265,269,270,271,331,332]. In this theory, the bare Coulomb matrix elements are divided by a dielectric function Counterintuitively, the static screening part depends not on the material itself, but on the surrounding material’s dielectric properties . This expression comes from the theory of a dielectric slab with a finite width d surrounded by two semi-infinite dielectrics with electric permeability . In the numerical results for MX discussed in the next sections, we consider two specific cases of the dielectric environment. In the first step, we analyze MoS (with “effective” thickness Å) on top of the bulk SiO () crystal [97,333], see Figure 10a for reference. We assume that from the top, the monolayer is surrounded by a vacuum; therefore, we take . Next, we consider the case of both MoS and MoSe encapsulated by hBN (as shown in Figure 10b), where we take [54].

Figure 10

(a) Schematic picture of the dielectric environment of the MoS monolayer on the SiO substrate. (b) Slab model of the MoSe/WSe heterostructure. Each of the MX layers with the width , respectively, is described by the dielectric constant .

As is known from the literature [334], more advanced models of screening do not significantly affect exciton spectra. This is related to the property that the R.-K. model describes well—screening from ab initio and only for large k-space distance (when matrix elements are small due to dependence) there is substantial difference between “correct” and approximated screening models.

The realistic description of the Coulomb electron–hole interaction becomes even more complicated in the case of heterostructures [278,280,335,336,337]. Furthermore, the complex form of the system leads to different formulas of dielectric functions for intra- and interlayer electron–electron interactions [335]. The Coulomb interaction potential for MX bi-layers can be written in the form analogous to Equation (24) as Following Ref. [335], analysing the result of Poisson’s equation in the limit of , we can rewrite as and In the above, R.-K. screening formulas and can be treated as an effective polarizability estimated from experimentally known -state energies of both intra- and interlayer excitons. One can expect that the interlayer screening model leads to the renormalization of the exciton spectrum, similarly to the monolayer case [281]. Interestingly, the interlayer R.-K. model gives a more precise dielectric screening description compared to the monolayer case [335]. We note that when we discuss the numerical results for the MoSe/WSe heterostructure encapsulated by hBN, we use [54].

4. Mechanisms of the Renormalization of the X Spectrum

In the following section, several mechanisms affecting excitonic fine structures are briefly discussed, with in-depth discussion already presented elsewhere [281]. Importantly, we note that our main motivation to construct the presented theory was understanding the effects of the multivalley (K and Q points) band structure of MX and the topology of the wave functions on the excitonic fine structure. Instead of focusing on agreement with ab initio + GW + Bethe–Salpeter, we note that our approach offers greater tunability, allowing to understand various physical mechanisms (e.g., role of dispersion, effective mass, strength of interactions, topology of wave functions, and screening mechanisms, etc.) contributing to the complicated problem of the excitonic spectrum. The understanding of these aspects is necessary to build a theory of optical response, e.g., in quantum dots, twisted bi-layer systems, or in the presence of magnetic fields, problems that are usually not easy to solve in the ab initio + GW + Bethe–Salpeter framework. We start with the effect of the band structure on excitonic levels. How different approximations to conduction and valence band energies affect the electron–hole pair energy is summarized in Figure 6. Direct transition energies enter as diagonal to the Bethe–Salpeter equation. In the parabolic-effective mass approximation with simplified bare Coulomb interactions and static screening, the excitonic spectrum gives the binding energy of the first state equal to −4 Ry and the degenerate second shell of three states: , , , whose energies are Ry. When we change the dispersion model to a massive Dirac fermion, keeping interaction and screening as previously, we notice that the binding energy of the state lowers to ≈−5.5 Ry. This behavior can be understood as an increased “average” effective mass, i.e., carriers in the massive Dirac fermion model on average taken over the BZ are described better by a higher effective mass than in the parabolic model. This naturally leads to stronger binding, corresponding to a lower energy of the state. We also observe a larger separation between the first and second shells and a small breaking of degeneracy between and states within the second shell. Finally, when a tight-binding dispersion model is used, the effects described for the massive DF model become even more pronounced. The energy of the state lowers as much as to −10 Ry, there is large renormalization of the – states’ energy differences, and states with exciton angular momentum () and are no longer degenerate. A similar effect is observed for higher shells. In addition to the “average lowering” of the effective mass process, there is a pronounced contribution coming from the existence of Q points that can be observed in excitonic wave functions. The breaking of the degeneracy of the and state for is also observed, together with the breaking of the degeneracy of states with different L. We conclude that due to existence of three Q points around K point in the single valley, the full rotational symmetry of s-like states is broken and therefore those states react more strongly than others to changes in electron–hole energy dispersion from the parabolic to the TB.

As discussed in the previous section, constant screening does not faithfully describe the actual screening of interactions in 2D semiconductors. The simplest correction, capturing the non-locality of screening via its dependence can be modeled by Rytova–Keldysh theory and parametrized by polarizability . The excitonic spectrum is heavily renormalized for TB model dispersion and the interaction when static screening is switched from the homogeneous one to the R.-K. model. Non-local screening has an opposite effect than changing dispersion from parabolic to TB, and state energy rises from Ry back to approximately Ry. Additionally, split and excited exciton states drastically change their energies, reversing the order of arrangement of excitonic states.

In addition to the effects of dispersion and screening that renormalizes the 2D exciton spectrum, the direct electron–hole interaction form factors have to be taken into account. Focusing on complex values of the form factors, the result of a calculation on a 7000 k-point grid (the largest we were able to study with the full effect of wave functions) that allows for a reliable discussion of the first and second excitonic shell is presented in Figure 11. Generally speaking, because its value is ≤1, its averaged effect translates into lowering the value of binding energy and pushing excitonic shells toward each other, as shown in Figure 11a. Due to the strong effect of tight-binding wave functions, it is understood that with realistic form factors, different parameters of screening have to be chosen, for example, as shown for static screening in Figure 11b. We conclude that the collective effect of the renormalization of the 2D Rydberg series by dispersion, non-local interaction screening, and carrier wave functions is to make s-like series look like a “more than 3D” exciton. It means that even though the exciton is physically confined to the 2D plane of MX, its excited state series resembles more the 3D excitonic Rydberg ladder of states. We note that non-hydrogenic Rydberg series is usually explained as an effect of non-local screening [47]. We find that it is not only related to screening, but also depends heavily on dispersion, especially secondary minima in the CB at Q points, and wave function contributions affecting the electron–hole interaction.

Figure 11

The first two shells of the excitonic spectrum with full tight-binding direct interaction form factors. (a) Effect of the form factors compared to approximation. (b) Effect of different static screening on the spectrum. On both (a,b), tight-binding energies of electron and hole are used.

As mentioned before, taking into account electron and hole wave functions affects the exciton fine structure by contributing to the renormalization of s-like states. In addition to that, another interesting effect occurs. Because the direct electron–hole interaction depends on wave functions, its value is in general complex and phases of matrix elements cannot be set to zero by any gauge transformation. The role of a complex interaction phase manifests itself in two effects: first, it slightly renormalizes positions of s-like states, however, this effect is small. On the other hand, a clearly visible effect is that in the second shell, normally degenerate p-like states become split in energy and mix, forming states. More generally, due to complex interactions, all states with a non-zero exciton momentum L also mix and split (e.g., , , etc.).

The splitting of states with non-zero L reminds of the situation when a magnetic field is applied. In our case there is no magnetic field but instead there is a “geometric” field resulting from the topology of wave functions. This orbital moment, described by the gauge invariant Berry connection, has to come from the properties of wave functions. How to trace this Berry’s connection effect on the electron–hole interaction has been shown elsewhere [282,338,339].

Now, let us focus for a moment on the effect of spin splitting on the exciton fine structure. The largest effect, related to the splitting of bands in the VB, is A-B exciton splitting. For meV we obtain A-B 1s excitonic states split by ≈125 meV for and the full TB model. A much more subtle effect, introducing the splitting of the A exciton state to spin bright (same spin arrangement in the VB and CB) and spin dark (opposite spins) is related to interplay between spin splitting in the CB and the different dispersion of bands, which can be understood approximately as different effective masses for carriers with different spins. Even though the single-particle arrangements of bands point to a “bright smaller than dark” arrangement, different effective masses cause an inversion of excitonic states and the spin-dark state becomes a ground excitonic state. We conclude that this situation happens in MoS, which has very small meV. This conclusion is consistent with recent GW-BSE calculations [219,230] and some experiments [67].

Up to this point, we neglected generally small electron–hole exchange interactions. This interaction, however, controls the splitting of dark and bright states, because it affects only excitons built from an electron and hole with the same spin. It increases the energy of the exciton as a result of quantum mechanical electron–hole repulsion, in opposite to direct interaction. Starting with spin-degenerate states to make our analysis more transparent, when no exchange interaction is present, dark and bright states are degenerate. When an unscreened electron–hole interaction is turned on, the bright-state energy increases (the binding energy is lowered) up to ≈20 meV, depending on the choice and details of the screening model. When homogenous screening or R.-K. models are applied also to exchange interactions, this value is significantly lowered. We conclude that in all cases, the trend is such that the dark excitonic state in MoS has the lowest energy, adding up to a similar conclusion from the spin splitting discussion above. We note that in the DFT + GW + Bethe–Salpeter approach [219], the electron–hole exchange interaction is always taken as the bare, unscreened one. In our approach, derived from the CI approach, both direct and exchange interactions are treated on an equal footing and should be screened in the same way. This issue has been discussed in the literature [230,325], but in our view further studies are necessary to understand the source of this discrepancy. On the other hand, irrespective to details of screening, we conclude that in MoS, the lowest in the energy excitonic state should be dark due to spin and exchange effects.

In next step, we analyze the excitonic spectrum of the MoS monolayer using the combination of the minimal tight-binding model, Bethe–Salpeter equation, and simplified Coulomb interactions theory, with form factors given by Equation (22) and the Rytova–Keldysh screening model. The exciton fine structure has been determined in two cases, namely the MoS monolayer on the SiO substrate and the hBN encapsulated MoS, respectively. We restrict ourselves to the first three shells (). For MoS on the SiO and MoS encapsulated in hBN we set the polarizability to correctly determine the experimentally known exciton ground state [54,247]. The 1s-state has been taken as eV and eV, so that we obtained for the SiO substrate and for the hBN surrounding.

Figure 12 shows the excitonic spectrum determined for the MoS monolayer. We restricted ourselves to (first three shells) only. For both types of surroundings, the state is characterized by a smaller binding energy than the states. Moreover, we deal with generic topological effects manifested in the exciton spectrum resulting from the properties of electron and hole Bloch wave functions, especially the Berry curvature. Inclusion of the effect of wave functions on the exciton spectrum results in the topological splitting of , , and states.

Figure 12

Exciton fine structure for the MoS monolayer (a) on the SiO substrate and (b) encapsulated with hBN, in a full TB model with complex electron–hole interaction included. Results restricted to the first three shells. The topological splitting of , , and states in the excitonic spectrum of the MoS layer is presented.

The excitonic spectrum for the MoS monolayer on the SiO substrate, shown in Figure 12a, is characterized by the exciton binding energy meV. The energy splitting between the 1s and 2s states has been determined to be meV. In the hBN-encapsulated MoS case presented in Figure 12b, both the excitonic binding energy meV and the determined – energy splitting meV are smaller than for the MoS on SiO.

We now turn to the effect of the second atomic layer on the exciton spectrum. Using the combination of effective mass approximation and the Bethe–Salpeter equation, including the interlayer Rytova–Keldysh screening [335], we now analyze the interlayer exciton spectrum for MoSe/WSe encapsulated in hBN [340]. For various spin-split combinations of bands, we study the fine structure for excitons with zero total momentum (at K point). Excitonic states built from electrons and holes from distinct layers present a rich spectrum of different types of optical transitions. In Figure 13 we consider interlayer A/B/ exciton series of MoSe/WSe, where A denotes the transition between WSe and MoSe spin-up, B transition MoSe–WSe spin-down, transition MoSe–WSe spin-up, and transition WSe–MoSe spin-down, respectively. We restrict ourselves to , , and states. The parameter was set in order to correctly determine the experimentally known 1s-state interlayer exciton energies [142]. We take the ground state binding energy of the interlayer exciton as meV [142] and obtain .

Figure 13

Interlayer exciton fine structure for the MoSe/WSe heterostructure encapsulated with hBN in the EMA, restricted to , , and states. Interlayer A/B/ exciton types have been distinguished, where A denotes the transition WSe–MoSe spin-up, B transition MoSe–WSe spin-down, transition MoSe–WSe spin-up, and transition WSe–MoSe spin-down, respectively.

In Figure 13, we show the excitonic spectrum of the MoSe/WSe heterostructure, restricted to 1s, 2p, and 2s states only. Considering all interlayer exciton types, namely A/B/ excitons, we can analyze both the A-B / and – splittings. Interestingly, the values for all interlayer 1s-2s splittings have been determined as meV, while meV and meV. Including in our studies various spin-split combinations of bands, considering not only the energetically lowest optical transitions (A excitons), we predict rich interlayer exciton series resulting with additional peaks due to the presence of all A/B/ exciton types. More work is needed to understand intra- and interlayer excitons in 2D heterostructures.

5. Summary

In summary, in this review we describe the excitonic problem in transition metal dichalcogenide semiconductors. We discussed building an ab inito-based tight-binding model that captures all the important features of the electronic structures of these materials. Then, we describe our theory of exciton, focusing on issues related to the evaluation of the interaction matrix elements. Then, we presented how such a theory can be used to understand the physics of the excitonic spectrum in both mono- and bi-layer heterostructures. We discussed the effect of band nesting on the exciton fine structure; Coulomb interactions; and the topology of wave functions, screening, and the dielectric environment. Finally, we followed by adding another layer and discussed excitons in heterostructures built from two-dimensional semiconductors. We hope that this review will be helpful to people who are interested in entering the fascinating field of the optical properties of low-dimensional semiconductors.