Exciton Relaxation Cascade in two-dimensional Transition Metal Dichalcogenides

Monolayers of transition metal dichalcogenides (TMDs) are characterized by an extraordinarily strong Coulomb interaction giving rise to tightly bound excitons with binding energies of hundreds of meV. Excitons dominate the optical response as well as the ultrafast dynamics in TMDs. As a result, a microscopic understanding of exciton dynamics is the key for a technological application of these materials. In spite of this immense importance, elementary processes guiding the formation and relaxation of excitons after optical excitation of an electron-hole plasma has remained unexplored to a large extent. Here, we provide a fully quantum mechanical description of momentum- and energy-resolved exciton dynamics in monolayer molybdenum diselenide (MoSe2) including optical excitation, formation of excitons, radiative recombination as well as phonon-induced cascade-like relaxation down to the excitonic ground state. Based on the gained insights, we reveal experimentally measurable features in pump-probe spectra providing evidence for the exciton relaxation cascade.


Exciton Relaxation Cascade in two-dimensional Transition Metal Dichalcogenides
Samuel Brem 1 , Malte Selig 2 , Gunnar Berghaeuser 1 & Ermin Malic 1 Monolayers of transition metal dichalcogenides (TMDs) are characterized by an extraordinarily strong Coulomb interaction giving rise to tightly bound excitons with binding energies of hundreds of meV. Excitons dominate the optical response as well as the ultrafast dynamics in TMDs. As a result, a microscopic understanding of exciton dynamics is the key for a technological application of these materials. In spite of this immense importance, elementary processes guiding the formation and relaxation of excitons after optical excitation of an electron-hole plasma has remained unexplored to a large extent. Here, we provide a fully quantum mechanical description of momentum-and energyresolved exciton dynamics in monolayer molybdenum diselenide (MoSe 2 ) including optical excitation, formation of excitons, radiative recombination as well as phonon-induced cascade-like relaxation down to the excitonic ground state. Based on the gained insights, we reveal experimentally measurable features in pump-probe spectra providing evidence for the exciton relaxation cascade.
Atomically thin semiconductors, such as transition metal dichalcogenides (TMDs), have revolutionized research in optics and electronics [1][2][3][4][5] . Alongside with the tremendous potential for innovative technologies, their optical and electronic properties further allow to study fundamental many-particle processes in an unprecedented scope. In particular, the reduced dimensionality of these materials leads to a weak dielectric screening, resulting in an optical response dominated by the emergence of excitons, i.e. quasi-particles comprised of Coulomb-bound electron-hole pairs [6][7][8][9][10] . Equivalent to hydrogen orbitals, excitons exhibit discrete internal degrees of freedom 11,12 , which makes them appear as a semiconductor analogon of atoms in a gas, but with properties which can be easily externally tuned [13][14][15][16] by e.g. changing the dielectric environment. The microscopic understanding and controlled manipulation of the internal structure of excitons could enable an entire new class of light emitters and absorbers.
Although, the analogy to atoms might help to illustrate some exciton properties, the formation and relaxation of these collective excitations are exclusively guided by many-particle effects [17][18][19] . Despite the continuous progress made in 2D material research during the last years, one of the key processes, the formation of bound excitons out of a quasi-free electron-hole plasma 20,21 , has not yet been theoretically investigated in TMDs. In standard luminescence experiments the semiconductor is excited with frequencies above the single particle band gap. Subsequently, the excited carriers relax via the emission of phonons and emit light corresponding to the energy of the excitonic ground state, which is far below the band gap. Recent experiments 22,23 revealed that the corresponding energy dissipation mechanisms is highly efficient in TMD monolayers, yielding a formation of 1 s excitons after off-resonant excitation on a picosecond timescale. The underlying microscopic mechanism enabling this ultrafast relaxation of excitons has not been investigated in literature yet.
In this work, we provide a fully quantum mechanical description giving access to the momentum-, energy-and time-resolved kinetics of electron-hole correlations. In particular, our model allows to track the phonon-induced relaxation of excited electron-hole pairs into their excitonic ground state after an optical excitation with frequencies far above the 1 s resonance. We apply the developed approach to map the relaxation dynamics of excitons in monolayer MoSe 2 on a SiO 2 substrate after an optical excitation 20 meV below the single particle band gap. We predict an ultrafast relaxation into the 1 s ground state with a relaxation time of about 1.1 ps. Further, we use the gained insights into the momentum-and energy-resolved exciton dynamics to calculate experimentally accessible ultrafast pump-probe spectra. Here, we reveal that the transient occupation of energetically higher excitonic

Microscopic Model
To illustrate the formation process of bound excitons after off-resonant excitation, we transfer the single-particle band structure into a two-particle dispersion, which is strongly modified by the Coulomb-attraction between electrons and holes. Figure 1 schematically illustrates the excitonic band structure and phonon-assisted relaxation processes in molybdenum-based TMD monolayers. We treat valence and conduction band in effective mass approximation, which allows to mathematically separate relative and center-of-mass motion of electron-hole pairs. The relative motion is determined by the mutual Coulomb attraction, yielding a Rydberg-type series of exciton eigenenergies 11,12 denoted with 1 s, 2 p, 2 s, etc. in Fig. 1. Additionally, excitons as a whole can move like free particles through the host crystal with a center-of-mass momentum Q. Therefore each excitonic state exhibits a parabolic dispersion, whose curvature is given by the sum of the effective electron and hole mass.
Apart from the energetically bound exciton states, there is also a continuum of scattering states with energies above the single-particle band gap, representing quasi-free electron-hole plasma states. Moreover, excitons can be composed of electrons and holes located either at the same valley (KK excitons) or at different valleys of the Brillouin zone [24][25][26][27][28][29][30] . Due to valley dependent effective masses, the intervalley excitons exhibit different dispersions and binding energies and further have an energy minimum at Q ≠ 0, which is illustrated in Fig. 1 for KΛ excitons. Here the hole is located at the K-valley whereas the electron resides at the Λ-valley. Note that excitons with Q ≠ 0 are dark, i.e. they can neither be excited optically nor decay radiatively, since photons provide only a negligible momentum transfer. In addition to this so called momentum-forbidden dark states, there are also angular-momentum-forbidden (relative motion with angular momentum l ≠ 0) 31 and spin-forbidden dark states (electrons and holes with different spin) [32][33][34] . In order to obtain access to the complex scattering mechanisms guiding the phonon-induced exciton relaxation, we neglect the exchange Coulomb interaction, which has been reported to induce a center-of-mass dependent energy renormalization and intervalley coupling 24 . This mechanism will lead to quantitative corrections of our predictions, but does not affect the principal mechanisms and typical timescales guiding the relaxation of hot excitons into their groundstate.
In this work, we model the situation, where a pump pulse with an energy close to the single-particle band gap excites electron-hole pairs at the scattering continuum at Q = 0, cf. yellow arrow in Fig. 1. These hot electron-hole pairs subsequently dissipate their excess energy via emission of phonons, performing a cascade-like relaxation through energetically lower lying exciton orbitals, cf. orange arrows in Fig. 1. Due to this process, the excited quasi-free electron-hole plasma cools down and forms a population of bound 1 s excitons. In the following we present a theoretical model, which predicts the above described mechanism based on a well established many-body quantum theory and allows to track the full relaxation kinetics of excitons after off-resonant excitation. We use an equation-of-motion approach to derive the dynamics of a coupled electron, phonon and photon system 18,19,35,36 . We solve the Heisenberg equation for interband coherences where we have introduced operators which create (annihilate) valence band electrons † v k (v k ) and conduction band electrons † c k (c k ) with momentum k. To obtain access to the relaxation dynamics of the system after an optical excitation, we take into account electron-phonon, electron-photon and electron-hole interactions.
In contrast to conventional semiconductors, the effective Coulomb interaction between charge carriers in TMDs leads to a strong coupling of interband coherences at different momenta, yielding an excitonic eigenspectrum for interband transitions. To calculate the energies μ E 0 and wave functions Φ μ k of these excitonic eigenmodes μ, we solve the Wannier equation 18,36-39 Figure 1. Intra-excitonic relaxation cascade in MoSe 2 . A pump pulse with an energy close to the single-particle band gap excites quasi-free electron-hole pairs at the scattering continuum (yellow arrow). The hot electronhole pairs dissipate energy via a sequence of phonon emissions, performing a cascade-like relaxation through energetically lower lying exciton states (orange arrows) including indirect KΛ excitons.
1 is given by the effective masses of valence-(m v ) and conduction band (m c ), while E g denotes the single particle band gap. For these material parameters we use values from ab initio calculations (PBE) provided in ref. 40 . Moreover, we take into account the dielectric environment of the two-dimensional monolayer by applying a Rytova-Keldish potential V q for the Coulomb interaction 39,41,42 . Based on the solutions of the Wannier equation we transform all electronic observables to the excitonic basis, by expanding their momentum dependence in terms of exciton wave functions. Therewith we get access to the excitonic polarisation = ∑ Φ μ μ ⁎ P P k k k (coherent excitons) and the exciton occupations (incoherent excitons) of the state μ and center-of-mass momentum Q 19,36,39 . Here, the mass weights v account for the different curvatures in valence and conduction band. In excitonic basis the dynamics of the system is given by The first equation reflects the evolution of coherent excitons, which are initially excited by the optical pump pulse . Here, M k denotes the optical matrix element for interband transitions 39 , whereas e 0 and m 0 are the elementary charge and the electron rest mass, respectively. The optically generated polarizations decay radiatively with the rate γ µ rad 26,43 as well as via non-radiative dephasing due to the interaction with phonons and form incoherent excitons μ N Q that are described in the second equation. The probability to scatter from the state (μ, Q) to the (v, Q′) is given by The appearing sum over + and − accounts for phonon emission and absorption, respectively, while ω λ q and λ n q stand for the phonon frequency and the occupation in the mode λ and momentum q. The phonon dispersion is treated in Einstein and Debye approximation for optical and acoustic modes respectively, with energies and sound velocities taken from ref. 44 . The delta function accounts for the energy conservation during the scattering processes. We have applied the Markov approximation for carrier-phonon triplets and therefore only allow strictly energy conserving scattering processes 19,36 . Finally, the exciton-phonon scattering cross section is given by the matrix element, where λ g c v q k / are the electron-phonon matrix elements in conduction/valence band at momentum k, regarding the phonon mode λ and momentum transfer q. Here we use the deformation potential approach for acoustic phonons (LA,TA) and the approximation for the Fröhlich interaction for optical phonons (LO,TO) deduced from DFPT calculations in ref. 44 . Therewith we take into account intravalley scattering within the K and Λ valley as well as intervalley scattering allowing for a transfer of excitons between both valleys, cf. Fig. 1.
While the interaction with phonons leads to a decay of the excitonic polarisation P μ in Eq. (2), the exact same terms drive the formation of incoherent exciton populations μ N Q , cf. first term in Eq. (3) 19,28 . This reflects the so called polarization-to-population transfer 17,18 , in which coherent excitons generated at Q = 0 lose their phase information and obtain a center-of-mass momentum due to the scattering with a phonon. Due to the Boltzmann-like scattering terms in Eq. (3), the incoherent excitons in turn redistribute along all exciton states via emission and absorption of phonons until a thermal equilibrium is reached. Additionally, exciton populations with vanishing center-of-mass momenta can recombine radiatively with the rate γ μ rad , which is calculated in analogy to ref. 26 .
The approach described above is restricted to the case of low excitation powers. When exciting high carrier densities the exciton-phonon scattering rates as well as the Wannier equation would need to be modified by phase space filling factors, which incorporate the fermionic substructure of excitons via Pauli-blocking. Moreover, when increasing the excitation density exciton-exciton scattering will significantly influence the relaxation cascade by offering additional scattering and recombination channels such as exciton-exciton annihilation. The quantum mechanical description of collisions of two-particle clusters is theoretically highly challenging but will be addressed in future studies. In this work we focus on the thermalization and equilibration guided by phonon-scattering, which is the dominant mechanism at low and moderate excitation densities, where excitons are well described as non-interacting bosons.

Exciton Relaxation Dynamics
The numerical solution of Eqs (2) and (3)  following, we describe the exciton dynamics in MoSe 2 monolayers on SiO 2 at room temperature after an exemplary optical excitation 20 meV below the single particle band gap, which corresponds to the eigenenergy of the 5 s-state, cf. yellow arrow in Fig. 1. Furthermore, we focus on the occupation dynamics of the so called A-series, i.e. excitons composed of electrons and holes located at the lowest conduction and highest valence band at the K-point, respectively. Moreover, we also include momentum-forbidden dark KΛ excitonic states, cf. Fig. 1. We take into account the 30 lowest lying direct as well as indirect exciton states up to the 12 s-, 10 p-and 9 d-state. This allows us to track the exciton dynamics in high resolution. Figure 2 shows the dynamics of the overall number of excitons (corresponding to the momentum-integrated exciton occupation μ N Q ) together with a color-coded decomposition into the appearing exciton states. Note that the different exciton states were added up successively, so that the colors need to be read similar to a bar chart, representing the fractional contribution of the particular state to the total number of excitons. Further, to obtain a good overview, we have summarized the very dense lying states with energies E > E 2s . The investigated TMD material is excited at t = 300 fs by a 100 fs long Gauss pulse with a central frequency matching the 5 s-eigenenergy (yellow shaded region in Fig. 2), which leads to the rapid generation of electron-hole pairs (black line). The optical field induces a polarisation, which corresponds to coherent excitons P μ , cf. Eq. (2). The formation of incoherent excitons simultaneously occurs via the so called polarisation-to-population transfer. Here, the interaction with phonons leads to a strong dephasing of the excitonic polarisation P μ generating incoherent exciton populations. Hence, the total number of excitons, increases rapidly during the pump pulse and shortly afterwards, cf. the black line in Fig. 2. After about 1 ps the remaining polarization has completely decayed into incoherent populations and the total number of excitons stays almost constant due to the large radiative recombination times in the 100-picosecond range. The fast decay of exciton populations on a picosecond timescale measured in experiments 22,23 , predominantly stems from defect assisted non-radiative recombination. Since these processes strongly depend on the sample quality and the used substrate, we have not included any disorder mediated mechanisms in our simulation and focus on the intrinsic phonon-guided and radiative relaxation. However for a direct comparison with experimentally measured dynamics, an overall decay of all exciton populations due to non-radiative recombination needs to be included. The color-coded decomposition into different exciton states shows that during the first few 100 fs only states > 2 s are formed (mainly the resonantly driven 5 s-and surrounding states), cf. red fraction in Fig. 2. Those higher excitonic states now decay into lower lying states on a picosecond timescale via emission of optical and acoustic phonons. After about 1.5 ps a 1/e-fraction of the total number of excited pairs has decayed into the 1 s ground, which is in good agreement with the experimentally measured ultrafast appearance of spectral resonances stemming from 1 s excitons on a picosecond timescale 22,23 .
Since the scattering probability crucially depends on the overlap of excitonic wave functions of the involved states, scattering between excitons with the same angular momentum is most favorable. The optically excited (bright) s-states preferably decay into the lowest (since in momentum space broadest) s-exciton, given that its dispersion can be reached under emission of a phonon. An other important aspect, is that the 2 s-state very effectively couples to the momentum indirect 1 s-KΛ, which in turn effectively couples back to the direct 1 s ground state. As a result, the indirect 1 s-KΛ-state effectively acts as a catalyst, rapidly promoting excitons from the 2 s to the 1 s-state under emission of two Λ-phonons, as illustrated in the right relaxation path in Fig. 1.
In total, the cascade ns → 2 s → 1 s-KΛ → 1 s leads to an ultrafast relaxation of band gap near electron-hole pairs into the excitonic ground state. Thereby, the relaxation cascade mainly occurs via s-states, since the overlaps of wave functions with different angular momentum are small due to phase cancellations. However, due to the large number of densely spaced p-and d-type states, there is also a notable amount of excitons scattered into those states. Due to their small scattering cross sections these anisotropic states slow down the overall relaxation process by trapping excitons in excited states.
The main features of the described relaxation cascade are illustrated in more detail in Fig. 3. Here, we present momentum-and energy-resolved snapshots of exciton occupations for the most relevant states 1 s, 1 s-KΛ, 2 p, 2 s, 3 s, 4 s and 5 s. We have chosen three times illustrating (a) formation, (b) relaxation and (c) thermalization  N μ (t). While the black line shows the overall number of incoherent excitons, the surface below the black line is divided into colored areas, whose surface ratio represents the relative fraction of the respective exciton state. After the optical excitation close to the single-particle band gap at t = 300 fs (yellow area) most electron-hole pairs occupy exciton states above the 2 s exciton. Those excited excitons decay into lower lying states via the emission of phonons and at approximately 1.5 ps a 1/e-fraction of the excited pairs has decayed into the 1s-ground state.
ScienTific REPORTS | (2018) 8:8238 | DOI:10.1038/s41598-018-25906-7 of excitons after a close-to-bandgap excitation. Figure 3(a) shows the exciton density at t = 250 fs. Here the pump pulse has just started to induce polarizations, so we can directly observe the optical formation process. In all bright exciton bands, we observe sharp peaks appearing at small momenta, which is characteristic for the polarisation-to-population transfer. Coherent excitons are optically generated within the light cone (Q ≈ 0) and are transferred to incoherent populations mainly via absorption of acoustic phonons. The sharp peaks appear at the intersection of phonon dispersion and exciton parabola. The inverse process (emission of acoustic phonons) leads to the accumulation of incoherent excitons at Q = 0, which explains the threefold features at small momenta. On the 2 s parabola, we also see several stripes at larger momenta, which result from the absorption or emission of an optical phonon. Figure 3(b) shows the exciton distribution at t = 400 fs shortly after the pump pulse has past the sample. Here, we still see the pronounced threefold structure at low momenta, but it is now much broader due to interactions with acoustic phonons. Additionally, we see a manifold of peaks at large momenta, resulting from the replication of the sharp peaks at low momenta due to absorption and emission of optical phonons. Especially in the 3 s and 2 s band we observe a very broad distribution of excitons, which follows from the decay of energetically higher occupations. Notably, the 2 p-state is almost empty, since it only couples weakly to the s-states due to the symmetry mismatch.
Finally, Fig. 3(c) shows the almost relaxed exciton occupation at t = 2 ps. The sharp density peaks stemming from the optical formation have been softened due to acoustic phonons and we observe smooth, almost thermalized exciton distributions. While excitons in the 1 s ground state are already in thermal equilibrium with the phonon bath giving rise to a Boltzman distribution with a finite occupation of the KΛ state, there is still a large amount of excitons trapped in long-living excited states. Especially exciton states with non-vanishing angular momentum are now occupied, which is exemplary shown for the 2 p state.
Apart from the exciton dynamics shown in Figs 2 and 3, we have also performed calculations for MoSe 2 on a diamond substrate. As recently shown in a joint experiment-theory study 16 , excitonic wavefunctions and binding energies are strongly influenced by the dielectric environment of the layer. We find that the binding energy of the groundstate decreases from about 440 meV on fused silica (ε sub = 3.9) to about 330 meV for a diamond substrate (ε sub = 5.7). Surprisingly, regardless of the decreased excess energy, the overall relaxation timescale (1/e-rise time of 1 s occupation) is about 200 fs slower on the diamond substrate. This observation is a consequence of a change in the exciton-phonon coupling, which depends on the overlap of exciton wavefunctions in momentum space, cf. Eq. (5). An increased screening of the Coulomb interaction manifests as larger exciton radius in real space, which in turn leads to different overlaps of the shifted wavefunctions in Eq. (5). While intravalley scattering with acoustic phonons becomes slightly enhanced, the scattering with optical modes (as well as intervalley acoustic modes) becomes less efficient on a diamond substrate. Since the exciton cascade is mainly driven by the emission of optical phonons, an increased screening slows down the relaxation of hot excitons.

Ultrafast Pump-Probe-Spectroscopy
After presenting a microscopic view on the relaxation of excitons in monolayer MoSe 2 , we now discuss the possibility to experimentally measure the predicted cascade-like intra-exciton dynamics. Since the investigated processes occur on a picosecond timescale, the best way to capture the dynamics is to use ultrafast sequences of pump and probe pulses. As the investigated quasi-particles are bosons, they do not induce a measurable bleaching of interband absorption via Pauli blocking. However, excitons can absorb/emit light by performing transitions between internal energy levels, as long as they are dipole-allowed, such as transitions from s-to p-states. Those transitions can also occur at center-of-mass momenta larger than zero, since the exciton just travels vertically between the parabolas of two states. Similar to atomic spectroscopy, the absorbed light of a certain frequency is directly proportional to the amount of excitons in initial and final state. The induced optical susceptibility Δχ after the pump pulse is governed by the difference in exciton occupations ν μ N N , Since the response at a certain frequency directly reflects the involved exciton occupations, the probing of intra-excitonic transition after an interband pump pulse allows to indirectly measure the relaxation dynamics on a femtosecond timescale 16,22,23,45,46 . To identify experimentally measurable traces of the predicted exciton cascade, we calculate the time-dependent intra-excitonic response function based on the exciton dynamics discussed in the previous section. Figure 4 shows the induced absorption as function of the probe energy for four different times. The colored spectra result from the dynamics presented in the last section (room temperature), while the black curves show the induced absorption at 77 K. For short delay times after the pump pulse at t = 0.5 ps and t = 1 ps (red and orange curve) we observe several absorption peaks for energies below 150 meV. They predominantly stem from high energy excitons, which occupy s-states above the 2 s and perform transitions to energetically close p-states. These spectral features provide direct evidence for the transient but ultrafast occupation of intermediate states during the relaxation from the optically excited close-to-bandgap states. Moreover, we predict that transient population inversions between s-and p-states arise during the relaxation process, since the phonon scattering preferably occurs between states with the same orbital symmetry. Since the broadening γ vμ of the absorption lines is determined by the exciton-phonon scattering rates, cf. Eq. (6), it strongly depends on phonon occupations and therefore on the temperature. Thus, experiments at low temperatures provide much narrower lines allowing to distinguish the contributions from close lying transitions. In particular, for the spectra measured at 77 K we observe pronounced gain signals (i.e. negative absorption) e.g. stemming from the population inversion between 2 s and 2 p at about 40 meV.
For larger delay times (blue curves in Fig. 4) we observe a strong increase of the 1s-np resonances stemming from excitons in the ground state. The rise time of this signal provides direct access to the overall relaxation timescale. At low temperatures, we can directly observe the transient occupation of the dark intervalley KΛ exciton during the relaxation process. For short delay times, we predict a clear shoulder on the high energy side of the 1s-2p peak stemming from the corresponding transition in the Λ valley, which disappears after equilibration of the system. At t = 4 ps the spectrum is dominated by transitions starting from the 1s-state, reflecting the fully relaxed exciton distribution. In our study we have focused on the dynamics and optical response of the A exciton series. In experiments the optical pumping close to the band gap will simultaneously also inject B excitons, which will undergo an equivalent relaxation process. Since phonon scattering does not allow for spin-flips the relaxation cascades of A and B excitons will occur independently. However, the measured pump probe spectra will be a superposition of responses stemming from both spin configurations. The two spin bands corresponding to A and B resonances exhibit very similar effective masses, so that the transition energies of both series are only slightly shifted (5 meV for 1s-2p transition). This results in an apparent broadening of the shown resonance peaks. In total, Fig. 4 shows that the temporal evolution of the pump-probe spectrum represents a powerful tool to map the ultrafast and cascade-like intra-excitonic relaxation dynamics. Here, upcoming and disappearing peaks can be directly related to a redistribution of excitons along the Rydberg-like series of excitonic states. So far, optical pump infrared/THz probe experiments have focused on the dynamics of the 1s-2p transition allowing to capture the formation of the groundstate population. Our simulation of the exciton dynamics and the corresponding optical response shows that when extended to a broad frequency range, this experimental technique further allows to study details of the relaxation cascade. In particular, experiments carried out at low temperatures will allow to clearly separate resonances stemming from close lying exciton transitions.

Conclusion
We have presented a microscopic view on phonon-driven ultrafast relaxation of excited electron-hole pairs into bound excitons. In particular, we have investigated the dynamics of A-excitons in MoSe 2 monolayers on a SiO 2 substrate after an optical excitation 20 meV below the single particle band gap. Based on the calculated momentum-and energy-resolved exciton dynamics, we have revealed symmetry dependencies within the excitonic relaxation cascade and predict a formation of 1 s excitons with a time constant of about 1.1 ps. Finally, we have presented signatures of intra-exciton dynamics in absorption spectra, which are accessible in ultrafast pump-probe experiments. For pump-probe delay times smaller than 2 ps, we predict pronounced multi-peak features for energies below 150 meV. For low temperatures, we find a clear gain signal resulting from a transient population inversion between 2s and 2p states and a pronounced high energy shoulder of the 1s-2p peak stemming from the intermediate occupation of dark intervalley KΛ states. The presented work reveals new microscopic insights into the ultrafast intra-excitonic dynamics in atomically thin 2D materials and provides a concrete recipe for future experimental studies.