Unraveling energy and charge transfer in type-II van der Waals heterostructures

Recent experiments observed significant energy transfer in type-II van der Waals (vdW) heterostructures, such as WS2/MoSe2, which is surprising due to their staggered band alignment and weak spectral overlap. In this work, we carry out first-principles calculations to shed light on energy and charge transfer in WS2/MoSe2 heterostructure. Incorporating excitonic effect in nonadiabatic electronic dynamics, our first-principles calculations uncover a two-step process in competing energy and charge transfer, unravel their relative efficiencies and explore the means to control their competition. While both Dexter and Förster mechanisms can be responsible for energy transfer, they are shown to operate at different conditions. The excitonic effect is revealed to drive ultrafast energy and charge transfer in type-II WS2/MoSe2 heterostructure. Our work provides a comprehensive picture of exciton dynamics in vdW heterostructures and paves the way for rational design of novel vdW heterostructures for optoelectronic and photovoltaic applications.

In vdW heterostructures, both charge and energy transfer can take place and compete with each other 3,23 . While charge transfer is crucial in applications such as photovoltaics 24 , photocatalysis 25 , and photodetectors 26 , energy transfer is essential to light-emitting diodes 27 , sensors 28 , and lasers 29 , etc. As charge and energy transfer compete with each other, one of them has to be promoted while the other suppressed for a given application. While charge transfer is believed to be a dominant process in type-II heterojunctions [30][31][32] , energy transfer is generally expected to take place in type-I heterojunctions [33][34][35][36] . However, recent experiments have observed significant energy transfer in type-II TMD heterojunctions [37][38][39] . Kozawa et al. 37 were among the first to report energy transfer in type-II MoSe 2 /WS 2 heterostructure on the basis of enhanced photoluminescence (PL) in MoSe 2 layer under an excitation resonant with the optical band gap of WS 2 layer. The energy transfer from WS 2 to MoSe 2 is efficient and ultrafast (<1 ps) at low temperatures (5-200 K) 37 . Moreover, the energy transfer can be enhanced by separating MoSe 2 and WS 2 with insulating h-BN layers 37 . Energy transfer was also observed in type-II heterostructures of MoS 2 /WS 2 separated by insulating layers of h-BN 38,39 . Although the experimental evidence of energy transfer in type-II heterostructures is clear, the underlying mechanism remains elusive. Based on temperature dependence of energy transfer rates, Kozawa et al. 37 proposed that the energy transfer follows Forster mechanism 40 , which involves excitons in WS 2 exciting higher-order excitons in MoSe 2 via a direct dipole-dipole interaction. However, low-energy excitons in TMD heterostructures tend to be non-emissive dark excitons whose dipole-dipole interactions are normally too weak to yield efficient Forster energy transfer. Based on this argument, Dexter mechanism 41 , consisted of simultaneous transfer of electron and hole, was proposed to be responsible for ultrafast energy transfer in type-I WSe 2 /MoTe 2 heterostructure 36 . Compounding to the confusion, the origin of efficient energy transfer in MoSe 2 /WS 2 heterostructure is poorly understood; it is not clear what is the driving force of simultaneous electron and hole transfer in the staggered type-II band alignment. This is particularly puzzling given weak temperature dependences of energy and charge transfer rates 37 , which rules out thermal fluctuations as the driving force. Therefore, it is of fundamental importance to elucidate the mechanism and unravel the competition of energy and charge transfer in type-II vdW heterostructures. It is also of timely, practical interest to explore exciton dynamics and seek control of energy and charge transfer in vdW heterostructures. However, despite scientific importance and practical interest, little work has been devoted to the competition of energy and charge transfer in vdW heterostructures, especially from a theoretical perspective. To the best of our knowledge, no energy transfer in TMD heterostructures has been studied by first principles. Thus, in this work we aim to fill the knowledge gap by providing a first-principles picture of energy and charge transfer in MoSe 2 /WS 2 heterostructure.
The main challenge in first-principles studies is that excitonic effect must be included because the strong Coulomb interactions between electrons and holes have profound effects on the dynamics of energy and charge transfer in TMD heterostructures 3,42 . The conventional first-principles approach to capture the excitonic effect is GW-Bethe-Salpeter equation (GW-BSE) method [43][44][45] based on many-body perturbation theory. However, it is prohibitively expensive to perform dynamics simulations with the GW-BSE method. To overcome this challenge, we employ a recently developed first-principles method based on the linearresponse time-dependent density functional theory (LR-TDDFT) 46,47 , which is computationally much more expeditious than the GW-BSE method. To unravel the competing dynamics of energy and charge transfer, we combine the nonadiabatic molecular dynamics (NAMD) with LR-TDDFT in this work. Specifically, we first perform ab initio Born-Oppenheimer molecular dynamics (BOMD) simulations for the MoSe 2 /WS 2 heterostructure. At each BOMD timestep, the excitation energies and many-body wavefunctions are determined with LR-TDDFT calculations using optimally tuned, screened, and range-separated hybrid (OT-SRSH) exchange-correlation (XC) functionals [48][49][50][51] implemented in conjunction with plane waves and pseudopotentials [52][53][54][55][56][57][58][59][60] . The time-dependent excitonic wavefunctions are then expanded in terms of these many-body wavefunctions 61 , with expansion coefficients computed from the fewest-switches-surface-hopping (FSSH) NAMD simulations [62][63][64] . Each coefficient represents the probability amplitude that a given exciton makes a transition from an initial state to a specified excitonic state. The phonon-assisted excitonic transitions are captured by the nonadiabatic coupling matrix elements which depend on the time-dependent wavefunctions and excitation energies of the system 61 . In sharp contrast to local and semi-local XC functionals, the OT-SRSH functionals can reproduce the correct long-range electron-electron and electron-hole interactions in solids with appropriate screening 51 , which enables the LR-TDDFT to capture the excitonic effect in extended systems. The details of the LR-TDDFT and NAMD methods can be found in Supplementary Methods. In particular, comparisons with available experimental and GW-BSE values of the band gaps and exciton binding energies are included in Supplementary Table 1. The combined LR-TDDFT and NAMD framework has been used to study exciton dynamics in WS 2 /MoS 2 heterostructures 42 , ZnO/polymer heterostructures 65 , conjugated polymers 61 , and small molecules 66,67 for photovoltaic applications.

RESULTS
Static properties of MoSe 2 /WS 2 heterobilayer The single-particle band structure calculated using the conventional hybrid XC functional (HSE06) 68 combined with SOC for the hexagonal unit cell is displayed in Fig. 1b. For the equilibrium structure (d = 3.3 Å), our calculations confirm the type-II band alignment in MoS 2 /WSe 2 heterostructure 37 , where the conduction band minimum (CBM) and the valence band maximum (VBM) reside on WS 2 (in deep red) and MoSe 2 (in deep blue) layer, respectively. Based on a band unfolding analysis, we can identify the original K-valley states of the constituent monolayers in the band structure of the heterostructure, and label them by black circles for K@WS 2 and red triangles for K@MoSe 2 in Fig. 1b. The conduction band energy of K@MoSe 2 is only 42 meV higher than that of K@WS 2 , which suggests a small band offset for the electrons at K-valley and agrees with the previous results 69, 70 . There is a hybridized valence band (~100 meV below the VBM), highlighted in the dashed box with a light blue color, which arises from the interlayer coupling. As the interlayer distance is increased to d = 8.3 Å, this hybridized band drops in energy and disappears from the energy window ( Supplementary Fig. 3), which has consequences as discussed below. Note that the single-particle band structure shown here only serves the purpose to define the relevant excitonic transitions. The energies of these excitons are actually determined based on LR-TDDFT with the OT-SRSH functional.
We next perform the LR-TDDFT calculations with the OT-SRSH functional to examine the excitonic properties at zero temperature. In the TDDFT and subsequent BOMD (and NAMD) calculations, the larger orthogonal supercell (the red dashed box in Fig.  1a) is used instead. In this work, we focus on singlet excitons as they have lower energies and are the subject of the experimental studies which we will compare to. Specifically, five low-energy singlet excitons are considered, and their charge densities (see Supplementary information for the definition) are shown in Fig.  1b: (i)S 1 is an interlayer exciton comprised primarily of an electron from CBM@WS 2 and a hole from VBM@MoSe 2 ;S 1 is the lowestenergy exciton in the heterostructure (1.64 eV); (ii) S Mo is an intralayer exciton with both its electron and hole from K@MoSe 2 , corresponding to the A exciton of MoSe 2 monolayer. S Mo is the second lowest-energy exciton (1.65 eV) and responsible for energy transfer in the heterostructure; (iii-iv)S 2 andS 3 are two interlayer excitons whose holes are spread over both layers while electrons are localized at WS 2 and MoSe 2 monolayer, respectively. The energy ofS 2 andS 3 is 1.75 and 1.80 eV, respectively; and (v) S 0 is an intralayer exciton with 100% of its electron and 80% of its hole localized at WS 2 . S 0 (2.01 eV) corresponds to the A exciton in WS 2 monolayer and also the initial excitation (or pump) in the pumpprobe experiments 37 . To distinguish the intra-and inter-layer excitons, we label the interlayer excitons with a tilde (~) on the top. Owing to the atomic-scale interlayer distance (~3 Å), even the intralayer excitons have small amounts of electrons or holes spilled across the layers, dictated by the uncertainty principle. The calculated oscillator strength for the interlayer excitons,S 1 ,S 2 ,S 3 is 0.005, 0.028, 0.03, respectively, and these excitons are considered as dark excitons. In contrast, the oscillator strength for the intralayer excitons S Mo and S 0 is 0.74 and 0.15, respectively, and thus they are bright excitons with stronger absorption and radiative emission.

Exciton dynamics from many-body NAMD simulations
In the following, we perform NAMD simulations at 300 K to examine the decay dynamics of the initial exciton (S 0 ) and as well as the formation dynamics of the interlayer excitons (S 1 ,S 2 ,S 3 ) and the intralayer exciton S Mo . The transitions from S 0 toS 1 ,S 2 ,S 3 contribute to charge transfer while the transition from S 0 to S Mo is considered as energy transfer. As shown in Fig. 2a, there is an ultrafast decay of the initial exciton S 0 , whose population decreases drastically in~11 fs. During this time, the average energy of the excitons is only slightly reduced as shown Fig. 2c, suggesting ultrafast transfer of S 0 to intermediate "hot" excitons. These hot excitons subsequently cool down and relax to lowenergy excitons at a longer timescale. In 500 fs, about 15% of the initial S 0 excitons decay into the intralayer exciton (S Mo ) via energy transfer, and~30% of the initial excitons decay into the interlayer excitons (S 1 ,S 2 ,S 3 ) via charge transfer. The detailed dynamics for the low-energy excitons are shown in the Fig. 2b. The remaining initial excitons transfer to intermediate excitons with approximate 32% of interlayer character and 18% of intralayer character. The dominant exciton species contributing to charge transfer is the interlayer excitonS 1 while the dominant species for energy transfer is the intralayer exciton S Mo , thanks to their low energies. The bright intralayer exciton S Mo would eventually radiatively recombine to yield the enhanced PL intensity observed in the experiment 37 . The two-step dynamics, i.e., the ultrafast decay to intermediate hot excitons followed by a slower relaxation to lowenergy excitons, is also observed at 50 K as shown in Fig. 2d. This two-step process appears general in TMD vdW heterostructures 36,42,71 . From Fig. 2d, two characteristic timescales (~48 and 200 fs) can be identified, with the former (~48 fs) corresponding to the ultrafast decay of S 0 and the latter (~200 fs) being where S 0 decay curve intercepts the energy/charge transfer curves. Both match well to the experimentally observed upper (15-45 fs) and lower (~200 fs) bounds for the timescale of energy/charge transfer at low temperatures (5-200 K) 37 . Also consistent with the experiment 37 , the charge and energy transfer rates show a weak dependence on temperature (50 K vs. 300 K). Interestingly, the energy decay is much greater at 300 K than at 50 K as shown in Fig. 2d, suggesting much stronger electron-phonon coupling at 300 K, leading to a faster decay time (11 fs) at 300 K. We have also performed NAMD calculations at 50 K with the MoSe 2 A exciton (S Mo ) as the initial state. The results are shown in Supplementary  Fig. 4; similar sub-50 fs ultrafast charge transfer (43 fs) from S Mo tõ S 1 is observed.
To shed light on competing energy and charge transfer, we examine the energy variation of the excitons (S 1 , S Mo, S M , and S 0 ) at 50 K as shown in Fig. 3a. S M represent "hot" intermediate states whose energy range is also highlighted in Fig. 3b. Due to thermal fluctuations, the degeneracy of the excitonic states at 0 K is broken, and as a result, four groups of energetic curves are formed, each represented by a different color. The energy of the initial exciton S 0 overlaps strongly with that of S M , yielding ultrafast decay (~50 fs) from S 0 to S M via adiabatic transitions. The intermediate hot excitons (S M ) subsequently relax to the lowenergy excitons (e.g.,S 1 and S Mo ) via nonadiabatic processes. The much longer relaxation process is due to the substantial energy difference between the intermediate excitons (S M ) and the low-energy excitons (S 1 and S Mo ). On the other hand, the energy overlap between the intralayer exciton S Mo and the interlayer excitonS 1 signifies the competing energy and charge transfer in MoSe 2 /WS 2 . This is in contrast to other type-II TMD heterostructures (e.g., MoS 2 /WS 2 42 ) with larger band offsets. Much less energy transfer was observed in MoS 2 /WS 2 42 due to a larger energy gap between the corresponding intralayer and interlayer excitons.
Since the average energy ofS 1 is lower than S Mo , charge transfer is expected to be more efficient than energy transfer as observed in Fig. 2. On the other hand, being the lowest-energy exciton,S 1 does not have broad energy overlap with other excitons, thus phonon-assisted transitions are limited forS 1 , whose energy exhibits a stronger harmonic (or quasi-periodic) variation as compared to other excitons. In Fig. 4, we present the Fourier transform of time-dependent energy variation of S 0 ,S 1 , and S Mo . It is found that S Mo is driven primarily by phonons around 200 cm −1 , corresponding to out-of-plane A 1g and in-plane E 1g optical modes in MoSe 2 layer. In contrast,S 1 is strongly coupled to phonons around 400 cm −1 , corresponding to A 1g optical mode in WS 2 layer. In addition to A 1g mode in WS 2 , S 0 also couples to lowfrequency acoustic modes. Selective phonon coupling ofS 1 and S Mo suggests that the competition of energy and charge transfer may be tuned by exciting different phonons. Finally, the smaller energy variation of S Mo (compared toS 1 and S 0 ) shown in Fig. 3a can be attributed to the lower energy phonons that it couples to.
To elucidate the energy and charge transfer mechanism, we examine the charge density evolution of the initial exciton S 0 during the first 50 fs of the NAMD simulation at 50 K. Specifically, one of the most probable NAMD trajectories is identified and along which the charge density of S 0 is calculated, shown in Fig.  3c. A direct and ultrafast charge transfer-with both electron and hole-between the two layers is observed. Specifically, a simultaneously transfer of the electron and the hole takes place at 30 fs, indicating that the Dexter mechanism is at play.  Table 1), thus their binding energies are greater than 0.17 eV. We have also calculated the oscillator strengths of the intermediate excitons, shown in Fig.  3b. They are found to have much smaller oscillator strengths compared to the bright exciton S Mo (dashed oval). This suggests that the direct dipole-diploe interactions between them may be too weak to yield ultrafast and efficient energy transfer. Thus, Forster energy transfer is not believed to play an important role in MoSe 2 /WS 2 bilayer. The fact that Dexter mechanism is primarily responsible for energy transfer in MoSe 2 /WS 2 should not be surprising because ultrafast charge transfer (<150 fs) of either electron or hole has been routinely observed in TMD heterostructures with type-II band alignment 30,32,[71][72][73] , thanks to their small interlayer distances (<1 nm).

Single-particle NAMD simulations
To uncover the excitonic effect on energy transfer in MoSe 2 /WS 2 bilayer, we perform a contrasting NAMD simulation in the absence of many-body excitonic interaction. Specifically, the ground state DFT-as opposed to the excited state LR-TDDFT-calculations are carried out at each BOMD timestep. The time-dependent singleparticle Kohn-Sham (KS) orbitals and energies, as opposed to the many-body wavefunctions and exciton energies, are obtained from these ground state DFT calculations. The time-dependent single-particle wavefunctions of "excited" electrons or holes are then expanded in terms of these KS orbitals and the expansion coefficients represent the probability amplitude of the "excited" electron or hole occupying one of the KS orbitals, which can be determined from FSSH-NAMD simulations. The excited electron initially occupies the lowest K-valley-conduction band of WS 2 , while the hole occupies the highest K-valley valence band of WS 2 , corresponding to the resonant excitation (A exciton) of WS 2 in the experiment 37 . The population dynamics of the "excited" hole and electron are displayed in Fig. 5a, b, respectively. It is found that the "excited" hole decays to intermediate valence states in~50 fs and about 20% of the holes are transferred to the K-valley of MoSe 2 in 1 ps. However, the "excited" electron in WS 2 cannot transfer to MoSe 2 , but remains at the K-valley of WS 2 during the entire NAMD simulation time of 1 ps. The forbidden electron transfer is due to the energy gap between the K-valley states of the electron at WS 2 and MoSe 2 , as shown in Fig. 5c. The thermal fluctuation at 50 K is not sufficient to overcome the energy gap. Therefore, in the absence of the excitonic effect, one would conclude that while charge transfer (via holes) is possible, energy transfer is forbidden in WS 2 /MoSe 2 at low temperatures; the latter of course is contradictory to the experiment 37 . These results have several important implications: (1) They demonstrate that the excitonic effect is the origin of energy transfer in type-II TMD heterostructures. The strong electron-hole interaction renders the overlap in energy between the interlayer excitonS 1 and intralayer exciton S Mo , thus promotes efficient energy transfer. In other words, the electron-hole binding provides the driving force for the electron to overcome the energy barrier. The aid from thermal fluctuations is not critical for energy transfer, which explains the weak temperature dependence of energy transfer in WS 2 /MoSe 2 .
(2) NAMD simulations based on the single-particle description are not justified in TMD heterostructures even if they may yield seemingly correct numbers. (3) The commonly used energy alignment diagrams for TMD heterostructures can be misleading because they are based on single-particle energy levels of monolayers. Without considering the excitonic effect, one may draw a wrong conclusion that electron transfer cannot take place in type-II WS 2 /MoSe 2 , thus rule out the Dexter mechanism erroeneously 37 .
Tuning the competition of energy and charge transfer It is of significant interest to be able to tune the competition of energy and charge transfer in vdW heterostructures. In the following, we explore two means of tuning energy and charge transfer in vdW heterostructures-the application of an electric field and altering interlayer distance. First, we apply a vertical electric field of pointing from WS 2 to MoSe 2 and carry out the same LR-TDDFT calculations. The exciton energies at 0 K under the electric field are shown in the last rows of Table 1. Since the intralayer excitons (S 0 and S Mo ) have in-plane dipoles, the Stark interaction is negligible, thus their energies are barely changed. The interlayer excitonsS 2 andS 3 with delocalized holes in both layers also exhibit minor changes (~0.05 eV) in energy under the electric field. In contrast, the interlayer excitonS 1 has a larger vertical dipole moment, hence its energy increases by 0.11 eV under the electric field. We next perform NAMD simulations to examine the energy and charge transfer dynamics under the electric field. As shown in Fig. 2e, the similar two-step dynamics is observed. However, the relaxation process is slowed down considerably in the presence of the electric field because the energy gaps between S 0 and the intermediate excitons are increased by the electric field. More importantly,~20% of the initial excitons decay into the intralayer exciton S Mo via energy transfer, doubling the energy transfer efficiency in the absence of the electric field. As shown in Fig. 2f, the population of the intralayer exciton S Mo increases substantially at the expense of the interlayer excitonS 1: Thus, the electric gating appears to be an effective means to modulating the competition of energy and charge transfer. Our results are consistent with the experimental observation 74 that the electric field pointing from WS 2 to MoSe 2 can enhance energy transfer in WS 2 /MoSe 2 heterostructure.
Lastly, we examine how interlayer distance may affect charge and energy transfer in vdW heterostructures 37,38 . Specifically, we determine the energies of low-lying excitons in WS 2 /MoSe 2 heterostructure with enlarged interlayer distances d. Note that d = 25.3 Å is comparable to the experimental interlayer distance (6 nm) in WS 2 /hBN/MoSe 2 heterostructure 37 . As shown in Table 1, the energies of the interlayer excitons (S 2 ,S 3 ) increase drastically, exceeding the energy of S 0 , due to the change in the hole energy. As alluded to earlier, the hybridized valence band shown in Supplementary Fig. 3 drops its energy as the interlayer distance d is increased. In contrast, the energies of the intralayer excitons (S Mo and S 0 ) are barely affected since their electrons and holes are at the same layer. Hence, as the interlayer distance d increases, charge transfer (via the interlayer excitons) is suppressed and as a result, energy transfer is expected to increase. For the same reason, Dexter mechanism involving electron-hole transfer is switched off and Forster mechanism becomes active for energy transfer at large distances (>1 nm) 40,75 . In other words, both the energy transfer efficiency and mechanism are dependent on interlayer distance. Once again, our results agree with the experimental observation that the PL intensity of MoSe 2 A exciton is increased in WS 2 /hBN/MoSe 2 compared to WS 2 /MoSe 2 at the same excitation energy of 2.0 eV 37 .

DISCUSSION
Before closing, we discuss two points relevant to the experiments. In the work of Kozawa et al. 37 , MoSe 2 and WS 2 monolayers were randomly stacked and as a result no moiré excitons were observed. Accordingly, we do not consider moiré excitons in the present work. In addition, the screening effect due to the quartz substrate is not considered in present study. Previous firstprinciples calculations 76 suggested that the effect of substrate screening is less important in TMDs. In particular, the reduction in exciton binding energy due to the substrate screening is roughly offset by the renormalization of the quasiparticle band gap, thus the optical gap (or the exciton energy) remains unchanged. In other words, the energy differences between the excitons are not expected to change substantially and the same for their dynamics. We also believe that the energy transfer mechanisms remain effective, although their relative contributions may be altered somewhat by the substrate screening.
In summary, we provide a first-principles and comprehensive description of energy and charge transfer in type-II WS 2 /MoSe 2 heterostructure. Focused on the dynamics of low-energy interlayer and intralayer excitons, we reveal that the charge and energy transfer proceeds by an ultrafast decay (10-50 fs) of the initial excitation to intermediate hot excitons, followed by a slower relaxation (100-200 fs) to the lowest-energy excitons. Under   Fig. 5 Results from single-particle NAMD simulations. Time-dependent population of the "excited" hole (a) and electron (b) at 50 K from the single-particle NAMD simulations without considering the excitonic effect. c Time evolution of the KS energy levels for the "excited" electron at K-valley of WS 2 (Red) and MoSe 2 (blue) layer from the BOMD simulation at 50 K.  The energies (eV) of the excitons with different interlayer distances d and with the applied electric field (last row). equilibrium conditions, charge transfer is~2 times more efficient than energy transfer. However, the energy transfer efficiency can be doubled by applying an electric field of 0.1 V/Å pointing from WS 2 to MoSe 2 . One can further boost the energy transfer efficiency by increasing the interlayer distance. The prevalence of low-energy dark excitons precludes Forster mechanism from playing an important role in energy transfer and Dexter mechanism is found to be responsible for energy transfer in the equilibrium heterostructure. However, when the interlayer distance is increased, the direct exchange of the electron and hole across the interface is suppressed, hence the Dexter mechanism is no longer effective, and Forster energy transfer mechanism takes over. While the dynamics of energy and charge transfer may be enhanced by coupling to different phonon modes, they exhibit a weak temperature dependence. Instead, the ultrafast and efficient energy and charge transfer in WS 2 /MoSe 2 heterostructure is driven primarily by the excitonic effect. Complementary to experiments, our first-principles framework represents a powerful tool in rational design of novel vdW heterostructures for optoelectronic and photovoltaic applications.

Computational models and details
As shown in Fig. 1a, the MoSe 2 /WS 2 heterobilayer is modeled by a hexagonal unit cell (blue dashed box) containing 75 atoms with 16.1°m isorientation and a negligible (0.06%) lattice mismatch between the two constituent monolayers. In order to capture the thermodynamics accurately in the ab initio BOMD simulations, we also adopted an orthogonal supercell (red dashed box) containing 450 atoms in the ab initio MD and the corresponding LR-TDDFT calculations. A vacuum slab of 15 Å is included in the supercell to remove spurious interactions between the periodic images. After an initial structural relaxation, ab initio BOMD simulations are performed to heat up the heterobilayer to a given temperature following the velocity rescaling protocol. The temperature is then kept for 2 ps to reach the thermal equilibrium, and finally a microcanonical BOMD production run is carried out for 2 ps with a timestep of 1 fs. The SOC is included in the calculation of the exciton energies and wavefunctions via the first-order perturbation theory 77 . The time-dependent population of each exciton is determined by averaging 2 × 10 4 NAMD trajectories. The oscillator strengths of the excitons are calculated by a recently developed LR-TDDFT method 78 based on a stochastic formulation of range-separated hybrid functionals. The rangeseparation parameters used in the stochastic LR-TDDFT calculations are the same as those in the (TD)DFT-OT-SRSH method. The structural relaxation and BOMD simulations are carried out with the Vienna ab initio Simulation Package (VASP) 79 , using the projector-augmented-wave potentials 80 with an energy cutoff of 400 eV and the generalized gradient approximation in the Perdew-Burke-Ernzerhof (PBE) 81 form for the XC functional. The semiempirical DFT-D2 method 82 is adopted to account for the vdW interaction between the monolayers.

DATA AVAILABILITY
The data that support this work are available in the article and Supplementary information file. Further raw data are available from the corresponding author (G.L.) upon request.