Exciton landscape in van der Waals heterostructures

van der Waals heterostructures consisting of vertically stacked transition-metal dichalcogenides (TMDs) exhibit a rich landscape of bright and dark intra- and interlayer excitons. In spite of a growing literature in this field of research, the type of excitons dominating optical spectra in different van der Waals heterostructures has not yet been well established. The spectral position of exciton states depends strongly on the strength of hybridization and energy renormalization due to the periodic moir\'e potential. Combining exciton density-matrix formalism and density-functional theory, we shed light on the exciton landscape in TMD homo- and heterobilayers at different stackings. This allows us to identify on a microscopic footing the energetically lowest-lying exciton state for each material and stacking. Furthermore, we disentangle the contribution of hybridization and layer polarization-induced alignment shifts of dark and bright excitons in photoluminescence spectra. By revealing the exciton landscape in van der Waals heterostructures, our work provides the basis for further studies of the optical, dynamical, and transport properties of this technologically promising class of nanomaterials.

In this work, we combine the density matrix formalism with density functional theory (DFT) calculations to reveal the exciton landscape for untwisted vertically stacked TMD homo-and heterobilayers. In particular, we identify the lowest lying exciton state for different TMD homo-and heterobilayers at different stackings. Furthermore, we disentangle the spectral shifts of different exciton types (Fig. 1c) stemming from hybridization and energy renormalization due to a charge polarization of the two layers (resulting in a periodic moiré * joakim.hagel@chalmers.se potential). We find a strong impact of hybridization on momentum-dark excitons, which can turn them into the energetically lowest states and the material into an indirect semiconductor. We also analyze the optical footprint of the exciton landscape by calculating photoluminescence (PL) spectra and taking into account direct and phonon-assisted exciton recombination. The latter allows an indirect visualization of dark excitons via the emergence of phonon sidebands [28]. This allows us to identify the microscopic origin of resonances appearing in PL spectra of different TMD homo-and heterobilayers at different stackings.

II. THEORY
We exploit the exciton density matrix formalism [29][30][31] with input from DFT to construct a material specific fully microscopic framework for modeling the exciton energy landscape in van der Waals heterostructures. The goal is to disentangle the microscopic contributions from hybridization and layer polarization-induced alignment shifts, to the final energies of hybrid excitons, as illustrated in Fig. 1.
We start by formulating a many-particle Hamilton operator in second quantization. The interaction-free electronic Hamiltonian for vertically stacked TMDs consists of three parts H = H 0 + H ren + H T , where H 0 is the stacking-independent kinetic energy. The renormalization Hamiltonian H ren describes the stacking-dependent layer polarization-induced alignment shift [12,32]. In the following, we refer to this as the alignment shift. Finally, H T is the hybridization Hamiltonian taking into account the overlapping electronic wavefunctions giving rise to hybrid exciton states. In the basis of monolayer

eigenstates, the bilayer Hamiltonian thus reads
where l/l are layer indices while α = (λ, ξ) is a compound index with ξ denoting the valley and λ = (c, v) the conduction and the valence band, respectively. The goal of this framework is to explicitly take into account dark intervalley exciton states [33]. Here, we have introduced the electronic creation (annihilation) operators a † (a). Furthermore,Ẽ α lk (S) = E α lk + ∆ε λ l (S) includes the stacking-independent kinetic energy of the monolayers (E α lk ) and the alignment shift (∆ε λ l (S)). Finally, T α ll (S) is the tunneling matrix element, where the stacking dependence (S) results from the varying interlayer distance and orbital symmetry of the Bloch functions involved. For the K point, the conduction and valence band Bloch waves are composed of d-orbitals [32,34,35]. Using the angular symmetry of d-orbitals, the tunneling around the three equivalent K points in an untwisted structure can be written as [32] T α ll (S) = 1 3 where C (n) 3 denotes the three-fold rotation operator and D(S) the stacking dependent lateral displacement between the layers. Furthermore, t α (S) describes the tunneling strength and τ is a prefactor which equals 1 for the K and −1 for the K point. The factor 1/3 reflects that the three equivalent K points have the same tunneling strength. For H-type structures, this expression obtains the additional phase e i(−1) l 2π/3 for the valance band [32].
Due to symmetry, tunneling can only occur for certain stacking configurations. To demonstrate this effect, we first consider tunneling around the K point for R-type structures (schematically shown in Fig. 3c) and apply the corresponding lateral displacement vectors D(S) for each high-symmetry stacking. We find that the tunneling term vanishes for R M h and R X h around the K point [36]. For an H-type structure one of the two layers is inverted (rotated by 180 • ) relative to the R-structure. Here ) stacking. In these cases, Eq. 2 simply becomes T α ll (S) = t α (S). For tunneling around the Γ and Λ points, there are no equivalent points within the first Brillouin zone (BZ) and thus the tunneling matrix element is only given by the tunneling strength T α ll (S) = t α (S) [32]. The latter can be extracted from DFT calculations of the bilayer band structure by considering the Hamiltonian in Eq. 1 as a 2 × 2 matrix with respect to the layer index l/l . Here, the diagonal components are given byẼ α lk (S) and the off-diagonal terms correspond to the tunneling T α ll (S). The eigenvalues of this matrix are given by the avoidedcrossing formula where ∆ α k (S) = E α 1k (S)− E α 2k (S) is the spectral difference between the monolayer energies as extracted from DFT calculations shifted by the layer polarization-induced alignment potential, sometimes referred to as the ferroelectric potential [37] (Fig. 1a). Furthermore, |t α (S)| is the tunneling strength (see Eq. 2) and E α ±,k (S) are the hybrid energies corresponding to the bilayer eigenenergies extracted from DFT calculations. As illustrated in Fig. 1b, the hybridization is particularly prominent at the Γ point due to the strong tunneling probability of holes. In contrast to the bands at the K point, which are mostly composed of d-orbitals localized around the metal atom, the bands at the Γ point also have a significant contribution from chalcogen atoms, thus considerably increasing the interlayer overlap [22,35,38]. Exploiting Eq. 3, the tunneling strength can be calculated for each band, valley, and stacking, yielding where ∆E α k (S) = E α +,k (S) − E α −,k (S) is the spectral difference between the hybridized electronic states.
To obtain a microscopic picture of the energy landscape in van der Waals heterostructure, we performed DFT calculations. First, we relax the investigated TMD homo-and heterobilayers using VASP [39][40][41]. The vdW-DF-cx functional [42] is used to account for the van der Waals interactions of the bilayers. The plane wave energy cutoff is 500 eV and the BZ is sampled using a 18 × 18 × 1 k point mesh. The electron density of the relaxed structures is subsequently evaluated using the local density approximation in GPAW [43,44] with a grid spacing of 0.15Å. The alignment shift ∆ε λ l (S) stems from a charge transfer induced dipole moment [45]. Its origin, the electrostatic potential is calculated by solving the Poisson equation. The alignment shifts can then be found within the project augmented wave formalism (Appendix A). Exploiting the computed monolayer energies (E α lk ), alignment shifts (∆ε λ l (S)) and the hybrid bilayer energies (E α ±,k (S)), the tunneling strength can be calculated from Eq. 4. The values are summarized for several homo-and heterobilayers in the appendices.
Having obtained both the alignment shifts and the tunneling-induced hybridization, we investigate now the impact of these microscopic processes on the final exciton energy. To achieve this goal, we transform the Hamiltonian in Eq. 1 into an exciton basis resulting in [12,30,32], Here, L = (l e , l h ) and ξ = (ξ e , ξ h ) are compound layer and valley indices, the first labeling the two intra-and interlayer excitons respectively and the second labeling the valley configuration. The index S denotes the stacking and Q the center-of-mass momentum. Furthermore, X † LξQ (X LξQ ) is the exciton creation (annihilation) operator and T ξ LL (S) is the exciton tunneling matrix element, which takes into account the tunneling strength and exciton form factors. The energy E ξS LQ is the exciton energy associated with the binding energy, the stackingindependent monolayer band edges and the alignment shift. The completely stacking-independent part of this energy, i.e. without considering the alignment shifts, yields the unperturbed exciton energies (the dashed lines in Fig. 1c). The interlayer and intralayer exciton binding energies are microscopically obtained by solving the Wannier equation , where the strong Coulomb interaction between electrons and holes is taken into account. Furthermore, the screening of the Coulomb potential is modeled with the generalized Keldysh screening [46]. The change in interlayer distance does in principle affect the screening, but since the variation in distance is small, this effect has an negligible impact on the binding energy.
The four different intra-and interlayer exciton states couple to each other through hybridization and split into four hybrid exciton states. The final hybrid exciton operator can be written in terms of a linear combination of the composing intra-and interlayer excitons i.e.
where η is the hybrid exciton quantum number. The mixing coefficient C ξη L (Q) fulfills the eigenvalue problem Here, E ξ ηQ (S) denotes the final hybrid exciton energies. Equation (6) provides microscopic access to the exciton landscape in different TMD homo-and heterobilayers, allowing us in particular to identify the energetically lowest exciton state for different stackings.

A. Exciton energy landscape
The final hybrid exciton energies are obtained by solving the eigenvalue problem in Eq. 6. We use the calculated tunneling strengths and alignment shifts, where the effective masses and valley separations are described in a monolayer basis [47]. This also yields the respective degree of hybridization for each exciton valley (K-K, K-K , K-Λ, Γ-K, Γ-K , Γ-Λ). As an exemplary material, we first investigate the untwisted MoS 2 -WS 2 heterostructure on a SiO 2 substrate. The latter is taken into account via a substrate screening of the Coulomb potential.
In figure 2, the lowest lying exciton state for all layer configurations at the R M h -stacking is displayed. It is evident that the Γ-K exciton is clearly the energetically lowest state. It is located about 300 meV below the lowest bright A exciton (intralayer K-K in MoS 2 ). Its purple color indicates that the Γ-K exciton is strongly hybridized reflecting a strong hole tunneling. The lowest lying K-K state is an interlayer exciton (red color). Figure 3a illustrates the lowest lying exciton state for each exciton valley as a function of stacking, while Fig. 3b shows the corresponding degree of hybridization. Energies and coefficients were computed for the three high- Fig. 3a,b), whereas intermediate displacements have been interpolated with the fit function given in Appendix A. We predict that the Γ-K( ) excitons are the lowest lying states regardless of stacking. Note that the Γ-K exciton is about 1 to 5 meV lower than the Γ-K state, this is mainly due to the slightly increased effective mass at the K point compared to that of the K point. Furthermore, we find that the lowest K-K transition is given by interlayer excitons, except for the R X h -stacking, where the intralayer exciton gives the largest contribution to the lowest lying K-K exciton. The K-Λ exciton, although very strongly hybridized, is located about 100 meV above the K-K exciton for all stackings. This is mostly due to the large energy difference between the Λ and the K point band edges in MoS 2 . The Γ-Λ exciton interestingly exchanges place with the K-K exciton, at R h h compared to R M/X h . This can be ascribed to the change in interlayer distance, which strongly affects the hybridization and has been implicitly taken into account in the energies extracted from our DFT calculations as in the latter we allowed for relaxation of the interlayer spacing.
By turning on and off certain parts of the Hamiltonian in Eq. 5) and consequently modifying the eigenvalue problem in Eq. 6, we can now disentangle the relative microscopic effects of hybridization and alignment shifts. perturbed interlayer (intralayer) exciton and the filled red bar, the lowest lying interlayer (intralayer) exciton after applying the alignment shifts. The two purple bars display the energetic position for the two final hybrid exciton states including the effect of hybridization. We first focus on the microscopic contributions of hybridization and alignment shifts for the K-K exciton (Fig. 4a). We find that the interlayer exciton is the energetically lowest state at R h h and R M h stacking (red bar). This is due to the large band offset between the layers (Fig. 1a) and the alignment red-shift. Rather unexpectedly at R X h , the intralayer exciton becomes the energetically lowest state (blue bar). Here, the metal atoms of the upper layer (WS 2 ) are on top of the chalcogen atoms of the bottom layer (MoS 2 ) (Fig. 3c). This variation in atomic configuration leads to a change in direction of the layer polarization field, thus changing the nature of the energy shift relative to R h h . As a result, the lowest lying interlayer exciton becomes blue-shifted -opposite to the red-shift for the R M h stacking (solid vs hatched red bar), where the chalcogen atoms of the upper layer (WS 2 ) are on top of metal atoms of the bottom layer (MoS 2 ) (Fig. 3c). This blue-shift is large enough to compensate for the otherwise large band offset between the layers, thus making the intralayer exciton the lowest lying state. However, this effect is difficult to observe in PL experiments due to the small oscillator strength of the interlayer exciton and the overall dominant nature of Γ-K excitons in MoS 2 -WS 2 (Fig. 3a). In addition, we observe that the K-K exciton is only very weakly affected by hybridization at R h h and not at all at R M/X h stacking (compare the purple bars with red and blue bar, respectively). The vanishing tunneling at R M/X h is enforced by the orbital-symmetry of the states at the K point (Eq. 2). As a result, the variation in the energy for the K-K exciton is only due to the alignment shift. The energy corresponds to the spectral distance (in meV) to the lowest lying A exciton (intralayer), in the corresponding stacking and material. For excitons denoted with K( ), the K -related exciton lies about 1 to 5 meV lower. The notation (I) indicates that it is an interlayer exciton.
Considering the situation for the Γ-K exciton in Fig. 4b, we immediately observe the significant role of hybridization that is reflected in the large red-shift of the dark purple bars. This is expected due to the strong tunneling of holes at the Γ point. The strong variation in energy between different stackings for the hybrid Γ-K exciton mainly results from the varying interlayer distance, which has a strong impact on the tunneling strength (see Eq. 4). However, also the stacking-dependent layer polarization field changes the energetic distance of the initial intra-and interlayer states and therefore indirectly influences the hybridization process.
We performed calculations for multiple TMD homoand heterobilayers mapping out their entire exciton energy landscape. The results are summarized in Table I including the lowest lying exciton for each material and stacking and its energetic distance to the bright K-K exciton. In the MoS 2 -WS 2 heterostructure, we predict the momentum-dark Γ-K( ) exciton to be the lowest regardless of stacking. The situation is completely different in MoSe 2 -WSe 2 , where the bright K-K exciton is the lowest in R-type stacking, the K-Λ exciton in H h h and the K-K exciton in H M/X h -stacking. Interestingly, the K-Λ exciton becomes the lowest due to the increased tunneling strength with the varying interlayer distance, while the K-K exciton is the lowest due to the nature of the H-type stacking, where the spin ordering changes in one of the layers. This means that the K-K exciton in Htype structures corresponds to K-K excitons in R-type stackings [36]. For TMD homobilayers, R M h and R X h are equivalent stackings just mirrored in the out-of-plane direction. Similarly to MoS 2 -WS 2 , the Γ-K( ) exciton is also the lowest state in MoS 2 -MoS 2 bilayers due to the strong tunneling of holes around the Γ point. In contrast, in WSe 2 -WSe 2 the K-Λ exciton is the lowest reflecting the efficient tunneling of electrons at the Λ point, whereas the Γ-valley cannot compete due to the large energy separation to the K point in the monolayer case. Interestingly, we find a different lowest exciton in WS 2 -WS 2 depending on the stacking: K-Λ exciton at the R h h -stacking and Γ-Λ exciton at R M h and R X h , where the interlayer distance is reduced resulting in a larger red-shift of the Γ-Λ exciton due to the hybridization.

B. Photoluminescence
So far we have been focusing on exciton energies in TMD homo-and heterobilayers revealing a rich landscape of various direct and indirect excitons. Now we investigate the optical response of these excitons by calculating PL spectra. In particular, we take into account direct and indirect recombination of excitons, the latter driven by emission and absorption of optical and acoustical phonons. To this end, we exploit the generalized Elliot formula for phonon-assisted PL derived in Refs. [28,32] reading where η(η ) is the hybrid exciton quantum number, σ the polarization of the photon, ξ(ξ ) the valley index, q the involved phonon momentum, j the phonon mode and ± denotes phonon absorption (+) and emission (−). The first part of the equation describes the direct recombination of bright excitons within the light cone at the K point. Here, L is a Cauchy-Lorentz distribution and M ξη σ is the exciton-photon matrix element determining the oscillator strength of the hybrid excitons. Furthermore, E ξ ηQ (S) is the hybrid exciton energy as calculated by the eigenvalue problem given in Eq. 6. The radiative and non-radiative broadening are described by γ ξ η and Γ ξ η , respectively. Since this work mainly focuses on the spectral position of the peaks, we account for these phenomenologically (Appendix D).
After optical excitation, excitons distribute throughout different valleys according to the Boltzmann distribution N ξ ηq . By emitting or absorbing a phonon, momentumdark excitons can become optically visible via (indirect) phonon-assisted recombination or in pump-probe spectra [48]. This higher-order process is determined by the second part of Eq. 7, where D ξ η q ξηj0 is the exciton-phonon matrix element that takes into account the electronphonon coupling deformation potentials (Appendix D) and (d) R X h stacking. The green and the red shaded areas indicate the phonon sidebands of the Γ-K exciton and the K-K intralayer exciton, respectively. Phonon sidebands stemming from emission (−) and absorption (+) of optical (Op) and acoustical (Ac) phonons are labeled accordingly. [49,50]. Furthermore, Ω qj denotes the phonon energy andñ ± qj = 1/2 ∓ 1/2 + n B (Ω qj ) the phonon occupation, which is given by the Bose-Einstein distribution n B (Ω qj ). In this equation we have neglected the unlikely higher-order process involving two phonons, which would be needed to scatter the Γ-Λ exciton to the bright K-K state.
By evaluating Eq. 7 numerically for different TMD bilayers at various temperatures, we reveal the optical footprint of the exciton landscape for different materials at different stackings. Figure 5a shows the temperaturedependent PL in the exemplary case of MoS 2 -WS 2 in R h h stacking. Here the position of the A exciton is fixed to experimental values [23]. PL spectra from other TMD homo-and heterobilayers for different stackings are shown in the appendices. We find that at higher temperatures, the bright K-K exciton dominates the PL, while at lower temperatures, phonon sidebands from the Γ-K exciton become pronounced. They stem from emission of acoustical and optical phonons and are clearly separated at very low temperatures, while they overlap at intermediate temperatures reflecting the increased non-radiative broadening. In the temperature range around 150 K, we find the emergence of additional peaks stemming from phonon absorption. We can also see a clear asymmetrical broadening of the peaks at higher temperatures. This is due to the Boltzmann nature of the exciton distribution and is characteristic for the appearance of phonon sidebands.
To further understand the nature of the multi-peaked PL, we show the PL at 70 K at different stacking in Fig. 5b-d. We identify the peak with the highest intensity as a phonon sideband stemming from the emission of an optical phonon. The next highest peak can be traced back to the emission of an acoustical phonon. These peaks have about equal contribution from their respective transverse and longitudinal components. There is a slight variation in the relative intensity of these peaks with stacking. This is due to the weak change in the degree of hybridization (Fig. 4b), which means that the phonon contributions from each layer will be weighted differently. A small peak stemming from the acoustical absorption phonon mode is also visible at higher energies. The PL at the three different stackings is qualitatively the same. We observe a clear stacking-dependent red-shift of the phonon sidebands reflecting different hybridization of the Γ-K exciton (Table I).
Our results are in good agreement with experiments, where two predicted phonon sidebands have been observed about 300 meV below the bright K-K exciton in MoS 2 -WS 2 [23]. This corresponds well to the calculated PL at R M h -and R X h -stacking, which is the energetically most favorable R-type stackings. The relative intensity of the two predicted phonon sidebands favoring the lower one agrees well with the experimental observation. An indirect dark exciton has also been observed in MoS 2 -MoS 2 [24,51] and WSe 2 -WSe 2 [27,51] located about 300 meV and 200 meV below the bright exciton, respectively -in a very good agreement with our results summarized in Table I. Similar observations have been made also for WS 2 -WS 2 including a peak roughly 200 meV below the K-K exciton [26]. This is in good agreement with the predicted phonon sideband from the K-Λ exciton in our calculations. Overall, our work provides microscopic insight into the exciton landscape and the optical signatures for a variety of TMD bilayers. In particular, the energy landscapes obtained for different lateral displacements play an important role for the understanding of the spatially dependent stacking domains and resulting moiré potentials [52,53] in twisted van der Waals bilayers.

IV. CONCLUSION
We have presented a material specific and fully microscopic model revealing the exciton landscape in TMD homo-and heterobilayers at different high-symmetry stackings. We identified the energetically lowest lying exciton state for each respective material and stacking. By combining the exciton density matrix formalism and density functional theory, we have gained microscopic insights into different contributions determining the hybrid exciton states including tunneling-induced hybridization and layer polarization-induced alignment shifts. Finally, we predict the optical footprint of the calculated exciton landscape in photoluminescence spectra including direct and phonon-assisted indirect exciton recombination processes. This allows us to uncover the exciton type behind the most prominent resonances for different TMD bilayers at different stackings. Furthermore, revealing the energy landscape at different stackings gives us insight into the expected behavior of the moiré potential at highsymmetry points in a twisted structure. Presenting the exciton landscape and its impact on optical spectra, our work can guide future theoretical and experimental studies in the growing field of van der Waals heterostructures.  Table III for R-type structures and in Table IV for H-type structures. The indices k and ξ have been dropped from the alignment shift ∆ε λ l (S), since the latter is a valley-independent energy renormalization of the entire band structure. The intermediate displacements have been interpolated with a continuous fit function given by [12], where E λ l is the stacking-independent monolayer energy and σ l = (−1) l for H-type stacking, otherwise σ l = 1. Furthermore, α λ l and β λ l determine the stackingdependent shift. These parameters are fitted to the band structure for the K point, obtained from first-principle calculationsẼ λK l0 (S) = E α l0 + ∆ε λ l (S), where the momentum k = 0, i.e the minimum of the parabolic dispersion. For the correct energetic position of other valleys, one has to take into account the spectral valley separation and spin orbit splitting Ref. [47] listed in Table II. The fitted parameters for all materials can be found in Tables V and VI. The tunneling strength as expressed in Eq. 4 in the main part of the manuscript is given in Table VII for different homo-and heterobilayers.

APPENDIX B: EXCITON HAMILTONIAN
We transform the Hamilton operator given in Eq. 1 into an exciton basis [12,32]. For this purpose, we expand the electronic creation and annihilation operators in pair operators P † ik,jk = a c † ik a v jk in the low density limit (with i = (ξ i , l i ) as a compound index) yielding This allows us to expand the Hamiltonian in center-ofmass coordinates, using the exciton operators [30,32] yielding Here, µ is the exciton quantum number, which we set to 1s, since we are focusing on the lowest lying exciton states. Furthermore, Ψ µ ij (k) is the exciton wavefunction as obtained from the solution of the Wannier equation [46]. Here, we have introduced . This consequently results in an exciton Hamiltonian where Q = k − k is the center-of-mass momentum, ξ = (ξ e , ξ h ) the exciton valley index and L(L) a compound index L = (l e , l h ) replacing the generic indices i, j. Here, we have assumed an untwisted structure and thus no momentum transfer is involved. The tunneling matrix element is given by The matrix elements for the conduction and valance band T c(v)ξ ll (S) are described in the main part of the manuscript. Finally, the exciton form factors F ξ LL are given by By expanding the exciton Hamiltonian into hybrid exciton operators, we can derive an eigenvalue problem, which diagonalizes the Hamiltonian and yields the final hybrid exciton energies as shown in the main text.

APPENDIX C: EXCITON ENERGY LANDSCAPE
Here, we discuss the exciton energy landscape for the four homobilayers MoS 2 , WS 2 ,WSe 2 , MoSe 2 and the heterostructure MoSe 2 -WSe 2 that have not been covered in the main text. Figure 6a shows the exciton landscape for MoS 2 − MoS 2 for the three R-type high-symmetry stackings. Similar to the heterostructure MoS 2 − WS 2 discussed in the main text, in this homobilayer Γ-K exciton is the energetically lowest state for all stackings. The Γ-Λ exciton is also significantly red-shifted, but due to the large valley separation of the Λ point, it still remains above the Γ-K exciton.
The exciton landscape for the MoSe 2 −MoSe 2 homobilayer is illustrated in Figure 6b. Here, we can see that the
K-K exciton is the lowest at R h h , while K-Λ becomes the lowest at R M h and R X h . This reflects the increased tunneling probability of electrons when the interlayer distance decreases at the latter two stacking configurations.
In the WSe 2 − WSe 2 homobilayer shown in Fig. 6c, we can instead see that the K-Λ exciton is the lowest state. This is mostly due to the very strong tunneling of electrons around the Λ point. In WSe 2 , the Γ point is very far away from the global valence band maximum and even though it has a larger tunneling strength than at the Λ point (see Table IV), the associated Γ-K/Λ excitons remain above K-Λ states. Due to the significant spinorbit splitting in the conduction band at the K point, we have also included the K -related excitons. Figure 6d illustrates the stacking-dependent exciton energies for the WS 2 − WS 2 homobilayer. Here, we find that Γ-Λ excitons are the energetically lowest states. In this material, the Γ (Λ) point of the monolayer is energetically much closer to the valence (conduction) band edge at the K point. Consequently, considering the strong hybridization in these valles, the Γ-Λ exciton becomes by far the lowest lying exciton in R M/X h -stacking. However, we observe that at R h h stacking, the K-Λ is the lowest state instead. This is due to the increased interlayer distance at these stackings.
For the four homobilayers discussed above, we can see no change between the R M h and R X h -stacking, since these stackings are mirror-symmetric for a homobilayer. We Photoluminescence spectra for the MoSe2-WSe2 heterostructure at 70 K. also observe a very strong hybrid nature for the low lying dark excitons (Γ-K, Γ-Λ, K-Λ) in comparison to the heterostructure discussed in the main part. This is mostly due to the lack of a band offset in homobilayers, which consequently reduces the detuning of the intra-and interlayer exciton and thus enhances the layer mixing.
In figure 7 the exciton landscape for the heterostructure MoSe 2 −WSe 2 is shown. In contrast to MoS 2 −WS 2 , which is discussed in the main part of the manuscript, the K-K exciton is the lowest lying exciton, regardless of stacking. This is mainly due to the large valley separation for Γ (see Table II) and the large band offset between the different materials. We can also see that the K-Λ exciton lies very close to the bright K-K exciton, especially at R M h and R X h stackings.

APPENDIX D: PHOTOLUMINESCENCE SPECTRA
As described in the main part of the manuscript, the photoluminescence (PL) spectra are calculated by considering both the resonant emission of photons from excitons and phonon-assisted recombination of excitons as a higher-order process. The exciton-photon matrix element which governs the direct photoemission reads [28,32], where M Lξ σ = κ rad L k Ψ ξ L (k) is the oscillator strength and κ rad L the square root of the radiative decay rate [5]. Furthermore, δ ξe,ξ h ensures that only bright excitons are taken into account in the completely resonant case.
For the indirect, phonon-assisted exciton recombination, the exciton-phonon matrix element D ξ η q ξηj0 plays the crucial role and is given by [28,32] Here, C ξη L (q) is the mixing coefficient and it holds C ξη L (0+ q) ≈ C ξη L (0) for untwisted structures, since the mixing coefficient only weakly vary due to similar masses of intraand interlayer excitons. The transferred valley momentum ∆ξ = (ξ e − ξ h − ξ e + ξ h ) ensures the correct transformation between the globally defined phonon momentum and the exciton/electron momenta defined in valley local coordinates. The selection rule δ l,l enforces that electron-phonon scattering takes place within one monolayer and that no tunneling occurs during this process. This assumption holds, as the interlayer atomic interaction is given by the van der Waals interaction, which is very weak in comparison to the in-plane atomic bonds. As a result, it is highly unlikely for phonons to tunnel during this process. The exciton-phonon matrix elements G c(v)j ξL,ξ L (q) for the conduction and valance band are given by [28,32], Here, g c(v)j ξξ (q) are the electron-phonon coupling deformation potentials, calculated via ab-inito calculations in [49,50]. The intravalley scattering with acoustic phonons is approximated as linear in momentum, otherwise the coupling is given by a constant. Similarly the phonon dispersion in Eq. 7 in the main text is approximated as a constant (according to the Einstein model), while the long range acoustic phonons (intravalley scattering) are approximated as linear in momentum (Debye model). The values for the phonon dispersion is also calculated via ab-inito calculations in [49,50].
The radiative γ ξ η ≈ 1 meV and non-radiative broad-ening Γ ξ η of the exciton resonances are treated phenomenologically. The non-radiative broadening is approximated with a linear temperature dependence in accordance with Ref. [32,54] with Γ Bright η (T ) = 5.0 meV + 0.05 · 10 −3 (meV/K) T and Γ Dark η = Γ Bright η /2.0. The non-radiative broadening of dark excitons is significantly smaller, since they are lower in energy and therefore have less relaxation channels for scattering with phonons.
In Fig. 8, the Pl spectra for the heterostructure MoSe 2 -WSe 2 is illustrated. Here we can see that the spectra is dominated by the bright K-K exciton at the R h h and R X h stacking. At R h h we also observe a small peak from the A exciton. A special case arises when considering R M hstacking: The lowest lying exciton is the interlayer K-K exciton, but due to its weak oscillator strength, the PL is dominated by the K-Λ exciton instead. Figure 9a-h show the PL spectra for the four homobilayers (MoS 2 -MoS 2 , WSe 2 -WSe 2 , WS 2 -WS 2 , MoSe 2 -MoSe 2 ) for different stackings. The spectrally lowest Γ-K and K-Λ resonances are marked green and yellow, respectively. We observe that these momentum-dark excitons clearly dominate the PL at 70 K for MoS 2 -MoS 2 , WSe 2 -WSe 2 and WS 2 -WS 2 . In the case of WS 2 -WS 2 , the actually lowest Γ-Λ exciton is not visible, since here an unlikely two-phonon process is required for the indirect exciton recombination. In the case of MoSe 2 -MoSe 2 , the bright K-K exciton dominates with a low-energy shoulder stemming from the momentum dark Γ-K and K-Λ excitons at R M h and R X h stackings.