Local field effects in ultrafast light-matter interaction measured by pump-probe spectroscopy of monolayer MoSe$_{\boldsymbol 2}$

Using a novel approach to ultrafast resonant pump-probe spectroscopy we investigate the spectral shape and dynamics of absorption features related to the A exciton in an hBN/MoSe$_2$/hBN van der Waals heterostructure. While in a pure two-level system a pump-probe experiment measures the occupation or the polarization dynamics, depending on the time ordering of the pulse pair, in the transition metal dichalcogenide (TMD) system both quantities get thoroughly mixed by strong exciton-exciton interaction. We find that for short positive delays the spectral lines experience pronounced changes in their shape and energy and they relax to the original situation on a picosecond time scale. For negative delays distinctive spectral oscillations appear indicating the first-time observation of perturbed free induction decay for a TMD system. The comparison between co-circular and cross-circular excitation schemes further allows us to investigate the rapid inter-valley scattering. By considering a three-level system as a minimal model including the local field effect, excitation induced dephasing and scattering between the excited states we explain all phenomena observed in the experiment with excellent consistency. Our handy model can be even further reduced to two levels in the case of a co-circular excitation, for which we derive analytic expressions to describe the detected signals. This allows us to trace back the spectral shapes and shifts to the impact of local field effect and excitation induced dephasing thus fully reproducing the complex behavior of the observed effects.


Introduction
Absorption and emission of light from semiconducting transition metal dichalcogenide (TMD) monolayers (MLs) is governed by tightly bound excitons [1,2]. The surface of a TMD single layer can be protected by embedding it between sheets of high quality, atomically flat hexagonal boron nitride (hBN). Such shielding prevents from aging effects due to erosive chemistry at ambient conditions [3,4]. Moreover, the hetero-structuring can flatten the TMD monolayers, eradicating disorder due to wrinkling and strain variations across macroscopic distances of many microns [5,6]. As a result, exciton transitions in modern hBN/TMD/hBN heterostructures reach line widths close to the homogeneous limit in the range of several milli-electronvolts (meV) at cryogenic temperatures. Such energies correspond to exciton dynamics on the sub-picosecond time scale as confirmed by nonlinear spectroscopy [7]. Femto-second (fs) multi-pulse spectroscopy is therefore required in order to investigate exciton dynamics in such systems.
Here, we perform resonant fs pump-probe measurements on the hBN/MoSe 2 (ML)/hBN heterostructure sample shown in Fig. 1(a). Considering co-and cross-circularly polarized excitations we study population decay and inter-valley scattering rates. We analyze the shape of the absorption features close to temporal overlap of the laser pulses in the regime of strong optical excitation. Characteristic energy shifts and line shapes give insight into the impact of local field effects and excitation induced dephasing (EID). Importantly, when probing the coherence by a probe-pump sequence, we detect characteristic spectral oscillations on the exciton line reaching over few tens of meV. This is usually attributed to the extinction of the exciton polarization transient by the second pulse [8][9][10] or by EID [11,12]. Here, it is additionally related to a rapid change of the character of the pump-probe signal accompanied by a frequency shift introduced by the local field effect [13]. This study is the first reported observation of this effect for excitons in TMDs. Our data are simulated by a few-level model including local field effects and EID, reaching excellent agreement with the experiment.

Sample and experiment
Our hBN/MoSe 2 (ML)/hBN heterostructure (hBN thickness bottom = 85 nm and top = 20 nm) was fabricated using a triaxial micromanipulator by a PDMS dry transfer method [14]. The used layers were exfoliated from bulk crystals with a standard micro-mechanical exfoliation technique. The sample was characterized by micro-photoluminescence (micro-PL) mapping at cryogenic temperature of = 5 K with continuous wave excitation at = 630 nm. The measured spectra in Fig. 1(b) display all characteristic features for this type of heterostructure found in the literature [4]. The PL in the bottom panel consist of a pronounced neutral exciton line (X) and a charged exciton transition (trion, T) shifted to lower energies by Δ XT ≈ 30 meV. Statistical analysis of the PL-map was performed by fitting the neutral exciton line with a standard Gaussian spectral shape. The histogram of exciton emission energies in Fig. 1(c) forms a smooth, peaked distribution between 1625 meV and 1650 meV. This effect is caused by the remaining local strain [6] of the crystal lattice induced in the stamping process as well as changes in the dielectric environment influenced by the presence of air bubbles trapped between the different layers of the heterostructure stack. These inhomogeneities are also reflected by the distribution of the exciton's spectral widths in Fig. 1(d), where we find that only a small fraction of measured spectra show a minimal width around 5 meV. Compared to previous measurements reporting widths of 2 meV [4,6,7] we can conclude that inhomogeneous broadening could still contribute to the optical response for the narrowest spectra. The further studies are done on selected spots exhibiting low linewidth broadening (5 meV) and spatial homogeneity as extracted from the PL mapping. In particular, we aim to reduce the influence of inhomogeneous broadening to avoid additional shortening of the exciton dephasing time and to stay within the temporal resolution of our experiment [15].
The applied pump-probe technique is a convenient tool to measure the dynamics of charge carriers and excitons in various semiconductor nanostructures [16,17], that employs a pair of fs laser pulses. Usually, it consists of exciting the sample with a strong optical pulse, called the pump, followed by a weaker pulse, called the probe. The intensity and spectral lineshape of the differential reflection signal launched by the probe are sensitive to the modifications induced in the medium by the pump pulse. Dynamics of such variations are monitored by varying the time interval between pump and probe, called the delay . The temporal resolution in our setup is only limited by the duration of the pulses Δ = 21 fs (FWHM 50 fs), measured via intensity autocorrelation (see Supporting Information). The pulses themselves are generated in a 76 MHz oscillator pumped with a Ti:Sapphire laser and are tuned to = 760 nm, which corresponds to the neutral exciton transition at a temperature of 5 K. We apply a novel approach to the pump-probe experiment in a micrometer resolution setup. It is performed in back-reflection geometry with efficient spatial separation between the pump and the probe beam. The laser pulses are focused on the sample with a 4 mm focal length, high-NA aspheric lens with 5 mm total diameter. The diffraction limit of the used laser beams corresponds to a spot diameter of = 3.8 µm. After the reflection from the sample the probe beam is directed into an imaging spectrometer and recorded on a CCD camera. Spatial separation between the parallel pump and probe laser beams allows for a high degree of extinction of the pump signal (in the order of 10 4 ) on the CCD, while maintaining their overlap on the sample. This makes it possible to perform measurements with spectrally degenerate laser pulses that have the same energy and polarization. Furthermore, this method is a complementary alternative to technically more involved coherent detection via optical heterodyning employed in recent experiments [5,7]. The reached micrometer resolution is also important for the case of heterostructures with a significant inhomogeneity across the sample where it allows for isolating spots of good, homogeneous optical quality. While recent developments in epitaxial techniques partly solve this issue [18], still a notable amount of scientific research relies on the preparation of samples by the lift-off method. Finally, small residual interference originating from the cross-talk between the pump and the probe is removed by periodically changing the optical path of the pump beam and averaging the signal over time. The relative change in the optical path is generated by a mirror located on a piezo element that oscillates with a frequency of about 30 Hz. The total spatial shift corresponds approximately to the laser wavelength, thus allowing for filtering out interference without deteriorating the temporal resolution. The experiments are performed at a temperature of = 5 K. The reflectance measured at the area investigated in the pump-probe study is shown in Fig. 1(b) (top panel). From this signal we calculate the imaginary part of the optical susceptibility Im( ) via a Kramers-Kronig transform according to Ref. [19] in the middle panel.

Theory
In the experiment we are resonantly exciting the lowest exciton states, i.e., the A excitons in the K and K' valley, with in general differently circularly polarized laser pulses. To treat this system theoretically, a minimum of three states is required as schematically shown in Fig. 2(a). The ground state |0⟩ has no exciton and the two excited states |±⟩ with the same energy ℏ 0 have an exciton in the K or the K' valley, respectively. Each circular polarization orientation addresses one of the two excitons, e.g., − -polarized light excites the |−⟩ exciton and + -polarized light the |+⟩ exciton. Further we include scattering between the two excitons, which leads to a transition of the occupation of the excited states. In practice we consider the following equations of motion to describe the microscopic polarizations + = ⟨|0⟩⟨+|⟩ and − = ⟨|0⟩⟨−|⟩ of the two excitons and the respective occupations + = ⟨|+⟩⟨+|⟩ and The excitons' decay rate Γ and dephasing rate are the same for both excitons, and the transition rate between the exciton occupations is . Especially the phenomenological dephasing rate may include different microscopic sources like fluctuations of the environment or phonon scattering which are for simplicity not further specified in this context. The process of occupation transfer leading to valley depolarization is a known phenomenon in TMD materials [20][21][22]. While the reason for the occupation transfer is still under debate, we choose here the simplest possible model which is a rate equation balancing the occupations in both valleys. Among the proposed valley depolarization mechanisms are scattering processes with phonons [23] and exciton-exciton interaction [24]. Consistently, we choose for the simulations a scattering rate in the sub-ps range. We always consider > Γ/2 + /2 to ensure that the density matrix is positive-definite. Note that since the states |+⟩ and |−⟩ are only coupled via scattering processes no coherence ⟨|+⟩ ⟨−|⟩ is created. The external optical driving is given by the electric field of the pulses ± ( ) translating here to instantaneous Rabi frequencies Ω ± ( ) = ± · ± ( )/ℏ, where ± are the dipole moments, allowing to address the two excited states individually. In addition we consider a local field effect parametrized by and an excitation induced dephasing (EID) by . We want to remark that the local field term results from a combination of different microscopic phenomena [25]: band gap renormalization, reduction of the exciton binding energy due to Pauli blocking and exciton-exciton interaction. Following the motivation of this effect in Ref. [26] the local field coupling should depend on the emission polarization of the considered exciton. Therefore, we only consider the local field coupling within a given valley. The EID stems from exciton-exciton scattering events and depends therefore on the total exciton density + + − . Note that here and are real-valued, while in Ref. [13] the two quantities were summed in a complex-valued .
To simulate the experimentally detected pump-probe signal we numerically calculate the system's dynamics for a two-pulse excitation with where the maximum of the probe pulse defines the time = 0. The delay between the pulses is and the pulses have phases 1 and 2 , pulse areas 1 and 2 and the same duration Δ . Both pulses are resonant to the exciton transition which is additionally renormalized by the local field strength in this model [13]. We identify the pulse with index 1 as the pump pulse and pulse 2 as the probe, meaning 1 ≫ 2 in agreement with the experiment. A positive delay > 0 indicates that the probe pulse arrives after the pump pulse, while the inverse ordering is described by negative delays. To isolate the pump-probe signal we filter the exciton polarization after the probe pulse 2 ( , ) with respect to the phase of the probe pulse, which characterizes the propagation direction 2 of this pulse in the experiment, via From this we simply retrieve the pump-probe spectrum via [12,27,28] where the tilde marks the Fourier transformed signal from the time to the spectral domain, defined as︀ In the following part of this theory section we consider a co-circular excitation which optically addresses only the states |0⟩ and |+⟩. In Ref. [13] we have shown that for a two-level system (2LS) the dynamics can be calculated analytically and brought into a compact form in the limit of ultrafast laser pulses. To carry out an equivalent analysis here we will reduce the three-level system (3LS) to a 2LS. To do so in principle we just have to disregard the exciton scattering rate choosing = 0. By doing so we have to take care of the other rates in the model. The scattering from |+⟩ into |−⟩ leads to a dynamical loss of population of the optically driven exciton. Therefore, we introduce the effective decay rate Γ ′ to compensate for the in general time-dependent influence of ̸ = 0. When comparing the ultrafast pulse limit in the 2LS with finite pulse durations in the 3LS this loss of occupation results in slightly different occupations for the same pulse area directly after the pulse. A detailed discussion of this effect for the present parameters is given in the Supporting Information. The phenomenological dephasing rate is considered to stay the same because it already sums up all dephasing processes in the 3LS. In the following we set = + , = + , and Ω( ) = Ω + ( ) and the reduced equations of motion read In the limit of ultrafast laser pulses, treated as delta functions in time, the full system's dynamics can be solved analytically as shown in Ref. [13]. Based on the analytical expression derived in the following, we can explain the lineshapes and extract approximations for energy shifts and dynamics that we will later compare to our numerical simulations for non-vanishing pulse durations and the experimental findings.
In the pump-probe experiment we have to distinguish between positive and negative delays. The usual signal generated for positive delays, where the stronger pump pulse arrives before the weaker probe pulse, is schematically shown in Fig. 3(a). Note that the schematic depicts occupation (yellow) and polarization dynamics (green) stemming from an actual numerical simulation for a given realization of phase combination and pulse areas. It also shows a respective pump-probe signal dynamics pp in red. In the limit of ultrafast pulses we get in the first order of the probe pulse 2 , which combines local field and EID in one complex quantity. The pump-probe polarization pp consists of two parts: The first term in the curly brackets stems from the linear polarization generated by the second pulse which automatically carries the phase 2 and also appears in the pure 2LS without local field and EID. The second term appears because of a phase mixing of the polarization generated by the pump pulse ∼ 1 with the phase-difference of pump and probe ∼ ( 2− 1) due to local field and EID as explained in Ref. [13] corresponding here to the diffraction by the transient grating generated due to the different propagation directions of the two pulses. Because both admixtures in this latter term are damped by dephasing, this contribution to pp decays with twice the dephasing 2 and with the EID in regarding the delay dynamics. We note that the Fourier transform of the probe pulse in Eq. (4) is a constant that is proportional to the pulse area 2 in the delta pulse limit. With this the pump-probe signal is given by This can be evaluated analytically in the case that the decay is slow in comparison to the investigated timescale, i.e., Then the pump-probe spectrum reads with the effective dephasing˜and transition frequency˜resulting from the EID and local field respectively, defined as˜( The signal consists of three contributions. In the pure 2LS ( = = 0,˜= ,˜0 = 0 ) only the first term would remain, while the other two stem from EID ∼ and local field ∼ , respectively. The first term has a Lorentzian lineshape of width˜and is centered around˜0. The EID contribution (third term) exhibits a minimum directly at and two symmetric maxima around =˜0. The local field contribution (second term) is asymmetric and changes its sign from negative at <˜0 to positive at >˜0. Such dispersive line shapes are known to appear in few-level systems that exhibit additional internal dynamics [28]. We find that all contributions from EID and local field vanish quickly for > 0. The overall amplitude is damped with 2 and via ( ), while˜and˜0 are additionally damped with the decay rate Γ ′ . Therefore, the spectral width and the spectral maximum relax exponentially towards and 0 − , respectively.
At exact pulse overlap for = 0 + the spectrum simplifies to This result already shows that for increasing pump pulse areas 1 the dispersive, asymmetric line shape becomes more pronounced while the symmetric Lorentzian gets weaker. At the same time the spectrum gets broader and its center experiences a blue shift. Such energy shifts are well known from quantum dots [29][30][31]. All these effects are at least of second order in the pump pulse area 1 . As will be discussed in detail below, and cannot be determined independent from 1 as long as the lowest order of the optical fields dominates. Our experiments are performed in this regime and therefore we will use 2 1 and 2 1 as fitting parameters to reproduce the experimental results. A more details discussion of possible values for and is given in Ref. [13], where a connection to the parameters retrieved from microscopically derived models [25,32] is given.
Next we discuss the case of negative delays. Here, the probe pulse arrives first. The polarization is then simply given by the free decay of the polarization If the delay is large enough that the entire polarization is decayed, the spectrum is simply given by a Lorentzian of width centered around 0 − . However, as schematically shown in Fig. 3(b) for smaller negative delays the dynamics are interrupted by the arrival of the pump pulse. Before the pump the signal dynamics follow the polarization created by the probe pulse. The appearing contributions from local field and EID change the pump-probe polarization dynamics drastically and it reads All in all, the polarization (i) is instantaneously reduced because the pump pulse redistributes the excitonic wave function, (ii) changes its oscillation frequency due to the local field contribution in (clearly seen in Fig. 3(b)), and (iii) dephases faster due to the EID entering in . This rapid change of the temporal evolution results in spectral oscillations similar to those well known from quantum wells [33][34][35] and quantum dots [8,29,30], that were also reported in multi-wave coherent control experiments on individual excitons [9]. Note, that for a pure 2LS only effect (i) is present. However, by including the local field effect not only the amplitude of the signal is rapidly interrupted by the pump pulse but also the oscillation frequency and damping of the signal changes which leads to a more involved origin of the spectral oscillations.

Results and Discussion
We first perform the co-circularly polarized (⟲⟲) pump-probe experiment and vary the delay between the two pulses from = −2 ps to 2 ps. The retrieved measured spectra are plotted in Fig. 4(a) as a function of the delay . The corresponding simulation in the 2LS model is plotted in the same way in Fig. 4(b), where we considered a pulse duration of Δ = 21 fs in agreement with autocorrelation measurements shown in the Supporting Information. In addition we choose ℏ 2 1 = ℏ 2 1 = 11.7 meV to reach the overall excellent agreement with the measurement. The other parameters are fitted as independently as possible to the experiments as explained in more detail in the following. Note, that the simulations have to be performed numerically when considering non-vanishing pulse durations. Nevertheless, we will in the following use the derived equations in the ultrafast pulse limit to qualitatively explain the found spectral dynamics.
Before discussing all details of the signal's dynamics we directly see that the simulation almost perfectly reproduces the measured data. Therefore, we can analyze both plots simultaneously. As already explained in the Theory section the signal behaves entirely different for positive delays, where the probe pulse comes after the pump, and negative delays, where the signal is created already after the first arriving pulse. Starting from large negative delays ≈ −2 ps the spectrum is given by a single peak. As can be seen in Eq. (11) the width of this peak is given by the dephasing rate and is determined to = 4 ps −1 . It is therefore fitted to the experiment independently from the other parameters. Although the fitted Lorentzian agrees well with the measured spectrum, small contributions to the line width stemming from sample inhomogeneity are compensated by the choice of . We define the center of the peak as natural exciton transition energy X , which corresponds to X = ℏ( 0 − ) in the model. For < 0 and approaching zero the spectrum begins to develop characteristic spectral oscillations. Such features are well known for pump-probe spectra on quantum wells [33][34][35] or quantum dots [8,29,30], where they originate from a sudden decrease of the polarization, an increased EID by the pump pulse (arriving second) or a Coulomb-induced spectral shift, like in our case. This introduces an asymmetry in the spectrum that is observed in the shape of the oscillations. Having a close look at the spectra for ≲ 0 we indeed find that it splits into two peaks, which is not expected for a normal 2LS without local field effect. Focusing now on positive delays, starting at = 0 the spectrum again consists of a single pronounced maximum which is shifted to larger energies compared to X . For increasing delays this maximum moves back to its original energy within approximately 2 ps. At the same time we see that the intensity starts with a maximum and is slightly quenched for a short time interval during this relaxation process, both in experiment and theory. The entire spectral dynamics are quite involved, especially for negative delays where spectral oscillations build up. Although the two color plots look alike we try to find a more quantitative comparison between experiment and theory for the entire delay scan. For that purpose in Fig. 4(c) we plot the positions of local maxima in the spectra. The found positions are plotted as blue crosses for the measurement and as red circles for the simulation. We find that the traces of local spectral maxima match perfectly, which demonstrates the high accuracy of our model. Considering the spectral relaxation for > 0 we again have a look at the analytical Eqs. (9) for ultrafast pulses and find that the signal shift decreases exponentially with the effective decay rate Γ ′ . Therefore, we use this decay to determine the rate to Γ ′ = 1.6 ps −1 again independently from the other system parameters. More details on this population relaxation are given in the Supporting Information. The found strong similarity between the co-circular measurement and the calculation in the 2LS show that it is a justified approximation for this excitation scheme. In the Supporting Information we also directly compare the numerical and analytical simulations in the 2LS.
Next, we move to the cross-circularly polarized (⟲⟳) pump-probe experiment, which is depicted in Fig. 4(d) in the same way as the co-polarized one. The fact that we detect spectral oscillations and a spectral shift of the signal maximum around = 0 shows that we have to consider the entire 3LS in Fig. 2(a). A cross-polarized pulse sequence cannot be described in a 2LS model. The corresponding simulation is depicted in Fig. 4(e) and the direct comparison of the local spectral maxima is given in Fig. 4(f). To achieve, once more, a remarkable agreement between simulation and measurement we determine the exciton decay to Γ = 0.6 ps −1 , the dephasing to = 4 ps −1 (same as for ⟲⟲) and adjusted the pulse area slightly such that the parameters are now given by ℏ 2 1 = ℏ 2 1 = 13 meV. The reason for the changed pulse area is at least partly due to a renormalization of the pulse area when changing the decay channels between 2LS and 3LS as discussed in more detail in the Supporting Information. To find the best agreement, we set the inter-valley scattering rate to = 4.4 ps −1 . As expected, we find that the effective decay rate in the 2LS is larger than the one in the 3LS. While in principle the same features as in the co-polarized case are found, they are significantly less pronounced. When the occupation is transferred between the valleys, the maximal occupation in the not-pumped valley does not exceed half of the maximal occupation in the pumped valley. Therefore the energy shifts are expected to be much smaller in cross-polarized excitations. One remarkable qualitative difference is found for small positive delays around the maximal energy shift of the spectral maximum, where the cross-polarized spectrum is strongly suppressed. The reason is that the phase-mixing processes between pump and probe, that enter in Eq. (9) via the combination of 2 − 1 and 1 , are not possible because the two pulses are orthogonally polarized. Therefore, the probe polarization is simply experiencing an additional oscillation due to the local field from the occupation already transferred into the probed valley. The signal is significantly damped due to the dephasing and the EID which depends on the total occupation. This total occupation is preserved under inter-valley scattering and decays only with the rate Γ and the damping influence of the EID is therefore maximal for = 0. This finally results in the significant suppression of the signal at small .
Although a cross-polarized excitation in principle allows for a biexciton creation, we do not find pronounced spectral features in the measured pump-probe spectra characteristic for a biexcitonic transition. On the one hand, the reason might be that the biexciton binding energy of several tens of meV [36,37] is large enough, such that the respective optical transition is not efficiently driven by the laser pulse centered at the exciton line. On the other hand, calculations for WSe 2 within a microscopic model show that the biexciton contribution is strongly broadened and mainly consists of a weak shoulder on the low-energy side of the exciton line [25].
To confirm the consistency of our model we performed the simulations for the co-circularly polarized excitation also in the 3LS model considering the same system parameters as for the cross-circular excitation. By slightly adjusting the pump pulse areas to ℏ 2 1 = ℏ 2 1 = 13 meV we achieved an equally good agreement with the measurement as in the 2LS model depicted in Figs. 4(a) -(c). The corresponding simulations in the 3LS for co-polarized excitation are shown in the Supporting Information.
As shown in Eq. (10) we expect that the influence of local field and EID manifest in the particular spectral shape at = 0 in the 2LS. In addition we find that their influence can be controlled by changing the pump pulse area 1 . In Fig. 5(a), we present typical pump-probe spectra with co-circular polarization (⟲⟲) obtained at pulse overlap, i.e., a delay of = 0, for increasing intensities of the pump pulse from bottom to top. The energy axis is shifted by X , i.e., to the peak energy at the largest negative delay of = −2 ps and all depicted curves are normalized in amplitude to the first spectrum (blue). The dark lines show the measured data and the light ones the simulations based on the 3LS (a corresponding simulation in the 2LS is shown in the Supporting Information). Starting with the smallest considered excitation power we find a single nearly symmetric peak that is slightly shifted to energies larger than X . With increasing pulse powers we find four striking changes of the spectra: (i) The peak maximum moves to larger energies, reaching a shift of approximately 6.5 meV for a pump power of 600 µW, (ii) the peak intensity shrinks by nearly 40% for the largest power, (iii) the peak gets significantly broader, and (iv) the spectrum develops an increasingly pronounced dispersive feature that even reaches negative values for the two largest pulse powers.
Comparing the measured spectra with the simulations we obviously find a strong similarity for all considered pump powers. The fitted system parameters are Γ = 0.6 ps −1 , = 4.4 ps −1 (both the same as in Fig. 4(d)-(f)), and = 3 ps −1 . In the Supporting Information we directly compare Fig. 5(a) with Eq. (10) to additionally confirm that our analytic calculations in the delta-pulse limit accurately describe the numerical simulations with a pulse duration of Δ = 21 fs. There, we find only slight deviations in the local field induced energy shift between the analytic and the numeric results. The considered pulse areas in the simulation are listed in Fig. 5(a) next to each spectrum. Starting from the lowest pulse area (0) 1 corresponding to 100 µW the higher areas grow according to the increase in pulse power in the experiment. As the pulse area is proportional to the electric field , the pulse power scales with its intensity 2 .
To extract a quantitative measure from the set of spectra we consider finding (i) and determine the energy max of the pronounced maximum in the spectrum by fitting a Lorentzian locally around the maximum and plot this quantity as a function of excitation density in Fig. 5(c). The blue crosses give the experimental and the red circles the theoretical data. As expected from the prediction in Eq. (10) we find a linear shift of the spectral resonance for small pulse intensities ∼ 2 1 . At this point it is important to note again that the influence of local field and EID enter the model in the lowest order of the pump field via 2 1 and 2 1 , respectively. Therefore, when operating in this regime of the optical fields, shown by the linear fit in Fig. 5(c), it is not possible to determine the actual pulse area and , independently. This means that in the simulations the choice of a larger 1 can be compensated by smaller and . However, we can determine the product of pulse power and local field factor, respectively EID strength, to Also the other findings (ii) -(iv) can be traced back to results from our theoretical model in Eq. (10). (ii) The fading of the signal strength stems from the term ∼ − sin( 1 /2) reducing the amplitude of the Lorentzian contribution of the spectrum due to the EID. (iii) The same effect leads to the increasing width of the spectrum in¯[Eq. (10)] such that the EID approximately preserves the integrated intensity. Considering Eq. (10) we find that the total intensity decreases ∼ 2 1 in the lowest order, as the integral over the last two terms vanishes. (iv) The appearance of a minimum in the spectrum can directly be traced back to the local field effect contributing with a dispersive feature in Eq. (10). These findings of the influence of local field and EID depending on the excitation density are in agreement with Refs. [25,32].
In Fig. 5(b) we perform the same pump power analysis for the cross-circular excitation (⟲⟳) where we consider ℏ  Fig. 4 we find that the energy shifts are less pronounced than in the co-polarized case. Also no clear minimum develops in the spectrum at increased excitation powers. However, the spectral line still significantly broadens, which shows that the influence of the EID remains nearly unperturbed. This demonstrates that the effect is only depending on the total exciton density as we have considered in our model in Eq. (1). At the same time the local field depends on the polarization orientation of the excitons. Therefore, the induced energy shift is much less pronounced in a cross-polarized excitation scenario because the probed valley has to get occupied to develop a local field coupling. This finding is again summed up in Fig. 5(c) where the energy shifts for ⟲⟳ exhibit nearly half the slope of ⟲⟲.

Conclusions and Outlook
In summary, we have studied the shape of the absorption features in the hBN/MoSe2/hBN heterostructure in the regime of ultrafast resonant excitation. All characteristic optical signatures related to exciton dynamics that appear in the pump-probe experiments were fully reproduced by the applied local-field model. By combining experiment and theory we have demonstrated that for pulse delays shorter than the exciton lifetime the A exciton response shows a blue shift of a few meV, which we could trace back to the impact of the local field effect. In addition, by investigating spectral line widths and amplitudes we could study the influence of excitation induced dephasing effects. It was found that both effects approximately have the same strength and also influence the appearance of spectral oscillations for an inverted pulse ordering. When moving from a co-circularly polarized excitation scheme to cross-circular excitation we had access to the inter-valley scattering rate which we found to be in the range of a few ps −1 .
Our results are in line with previous nonlinear spectroscopy studies of TMD systems and further strengthen the potential of this technique to explore ultrafast dynamics in layered materials. A handy fewlevel model, which even allows to develop analytical expressions in some special cases, explains the observed features at least qualitatively, thus offering insight into the physics behind the spectral dynamics. Therefore, we have a useful tool for example to analyze higher excited states like biexcitons or the fundamental impact of external magnetic fields in the context of local field coupling and excitation induced dephasing. Such external perturbations lead to distinct shifts of the exciton energies in opposite valleys, which have an impact on inter-and intra-valley scattering mechanisms. In the forthcoming work, we will explore the impact of such high exciton density effects on coherent optical response in more advanced experimental and theoretical configurations, for example offered in four-wave mixing spectroscopy.  These are compared to cooresponding simulations in the two-level system (2LS) in Fig. 1(b), which are the same as in Fig. 4(b) in the main text. Regarding the spectral oscillations for negative delays and the energy shift and relaxation for positive delays, which are highlighted in Fig. 1(c), we find a good agreement. This shows that the experiment is well reproduced in both models. However, we find one slight difference in the 3LS compared to the 2LS: The amplitude of the signal recovers significantly slower for positive delays after it is quenched for ≳ 0. This leads to an overall better agreement with the experiment in Fig. 4(a) in the main text. The reason for this is that the decay rate in the 3LS is smaller than the effective one in the 2LS Γ < Γ ′ . Given that the scattering does not reduce the total occupation and that the EID depends on this total exciton occupation, the influence of the EID decays much slower in the 3LS. Therefore the signal recovery takes much longer in the 3LS. To demonstrate that also the pump-probe spectra at = 0 in co-circular polarization can be well reproduced in the 2LS Fig. 2(a) shows the respective results. Panel (b) is the same as Fig. 5(a) in the main text. To reach the excellent agreement between simulation (bright) and measurement (dark) we determined the system parameters to Γ ′ = 1.6 ps −1 , = 3 ps −1 , and ℏ (0) 1 2 = ℏ (0) 1 2 = 2.2 meV. We again find that the effective decay rate has to be increased with respect to Γ in the 3LS to compensate for the missing inter-valley scattering channel. The energy shift as a function of the pump intensity is depicted in Fig. 2(c) and shows an excellent agreement between the 2LS, the 3LS simulations and the experiment.

Simulations in the ultrafast pulse limit
In Fig. 3 we directly compare the numerically simulated spectra at pulse overlap, i.e., = 0, in the 2LS with the ones from the analytical derivations. The solid lines stem from the numerical simulation and are the same as in Fig. 5(a) in the main text. The dashed lines are the respective results in the ultrafast pulse limit in Eq. (11) in the main text. In Fig. 3(a) we choose exactly the same parameters for the numerical and the analytical calculations.
Overall we find that the spectral shapes agree very well for each considered pump pulse area. But we find that each spectrum in the delta-pulse limit shows a larger shift to higher energies. There are two reasons for this: (i) In the analytic result we entirely disregard the exciton decay (Γ ′ = 0), which leads to a larger occupation and therefore a larger local field induced shift. (ii) During the excitation in the numerical calculation the system already dephases, which results in reduced polarization and occupation compared to the analytic result. This additionally results in an increased energy shift as explained before. To confirm this in the numerical simulations in Fig. 3(b) be set Γ ′ = 0 and choose a renormalized pulse area to compensate for the dephasing during the pulse (see following section for more details). Now we indeed find an almost perfect agreement between the two sets of simulations. The remaining deviations stem from the still non-vanishing pulse durations of Δ = 21 fs in the numerical treatment. If this value is reduced to approximately 1 fs both simulations agree perfectly (not shown here).

Impact of decay channels
As briefly discussed in the Theory section of the main text, the occupation after a single pulse depends on the considered system. As a reference we choose the pure 2LS without any decay excited with an ultrafast pulse where the occupation is given by = sin 2 ( /2), which is shown as black solid line in Fig. 4(a). When considering a non-vanishing pulse duration of Δ = 21 fs and adding the dephasing rate and the effective decay rate Γ ′ from Fig. 5(a) in the main text we have additional contributions reducing the occupation after the pulse. This leads to the solid blue curve which clearly exhibits smaller occupations for a given pulse area. The occupation is further reduced when local field and EID from Fig. 5(a) are included in the red curve. Finally going to the 3LS by including the inter-valley scattering , and choose Γ from Fig. 5(b) (main text) instead of Γ ′ we get the green dashed curve with the smallest occupations. This discussion explains why we needed to choose smaller pulse areas in the 2LS than in the 3LS because be needed to reach approximately the same occupation in the different systems. It also shows why the simulations in the ultrafast pulse limit in Fig. 3 exhibit a larger energy shift because we chose the same pulse area as for the non-vanishing pulse duration. To again embark on the interplay between the intervalley scattering and the decay rates Γ (′) in Fig. 4(b) we plot the occupation dynamics after a pulse excitation at = 0. The dashed blue curve depicts the 2LS with effective Γ ′ , while the solid red and green curve show + and − in the 3LS, respectively. The laser pulse only addresses the |+⟩ exciton in the 3LS, which after the excitation rapidly scatters into the |−⟩ exciton. After the two occupations are balanced the decay slows down significantly. As explained in the Theory section of the main text the loss rate changes over time due to the scattering process. To compensate for this time dependence in the 2LS we choose the effective decay rate Γ ′ such that the mono-exponential decay is slower at the beginning and faster at the end of the relaxation as seen in Fig. 4(b).

Estimation of carrier concentrations
One of the consequences of the spatial separation of pump probe laser beams in the setup is the misalignment of the pump beam with the optical axis of the lens. This in turn results in a larger excitation spot on the sample, which on the one hand fully encompasses the probing area, but on the other hand makes it difficult to estimate the actual density of photogenerated excitons. In order to calibrate the concentration of carriers during the pumpprobe measurement we measured the reflection from the sample with a single probing beam as a function of its power. The path of the probe laser beam lies on the optical axis of the lens allowing for good control over the size and shape of the laser spot on the surface of the sample. The final results presented in Fig. 5(c) show a pronounced blueshift of the exciton resonance that increases with the power of the laser. From this we can extract how it depends on the concentration. To do so we calculate the photogenerated carriers via = photon /( 2 ), where is the total absorption of the ML, photon = / avg is the number of photons per pulse with the laser power , the laser repetition period , the average energy of photons avg , and the laser spot radius .
To calculate the absorption one needs to look at the overlap of the ML absorption coefficient and the fs pulse spectrum. In order to take into account contributions originating from light interfering between different layers of our sample we simulate the reflection by a transfer matrix method (TMM) similarly to previous works [1,2]. In this approach the total reflection from the heterostructure is given by the ratio of transfer matrix elements = | 21 / 11 | 2 . Here, is the product of all successive interface and layer matrices of the full heterostructure. The propagation of light on a single interface is given by where 1 and 2 are the refraction coefficients of neighboring materials 1 and 2. The propagation within a layer is described by where is the wavevector in vacuum and the refraction index and the thickness of layer . In the simulation we use the following values for the refractive indexes of the heterostructure materials: n hBN = 2.1, n SiO2 = 1.54, n Si = 3.9 [3][4][5]. In Fig. 5(a) we present the measured reflectance contrast from the heterostructure (blue) and the TMM simulation (red). The fitting was done via the resonance parameters included as an imaginary addition to the permittivity function of an isolated optical transition where is the amplitude, the resonance width, and 0 the resonance energy. The obtained resonance shape shown in Fig. 5(b) bares close resemblance to the imaginary susceptibility derived from a Kramers-Kronig transformation used in the main text. The absorption coefficient function ( ) = 4 Im ( ) 2 is then used in order to calculate the total absorption of our laser while considering its overlap with the resonance. Finally, = (6.3 ± 0.5)% is used to calculate the density of photogenerated carriers in the reflection measurement. In Fig. 5(d) we show the exciton resonance as a function of carrier density. By fitting a linear function Δ = to the data in the low excitation regime we retrieve = (0.9 ± 0.7) × 10 −12 meV cm 2 which is a value equivalent to what can be found in other works for similar MoSe 2 heterostructures and falls within an order of magnitude to a theoretical estimation of = 5 × 10 −12 meV cm 2 [1,6]. Finally, by considering the blueshift in the pump-probe measurements we can translate the measured pumping power into the density of photogenerated excitons (see Fig. 5 in the main text). For the data presented in Fig. 5 an average pumping power of 100 µW corresponded to = 10 12 /cm 2 .

Determination of the pulse duration
The temporal resolution of the pump-probe experiment is given by the laser duration, which we determine by an autocorrelation measurement. In order to take into account possible dispersion effects the probing point for the autocorrelation was chosen right before entering the cryostat after the beams have passed through all optical elements in the setup. The signal intensity is presented in Fig. 6 and fitted with a gaussian function of standard deviation Δ int = 30 fs. To retrieve the correct pulse duration for the electric field in Eq. (2) in the main text we have to scale fitted value for the intensity via Δ = Δ int √ 2 ≈ 21 fs and a corresponding full width at half maximum (FWHM) of 50 fs.