Ultrafast Coherent Exciton Couplings and Many-Body Interactions in Monolayer WS2

Transition metal dichalcogenides (TMDs) are quantum confined systems with interesting optoelectronic properties, governed by Coulomb interactions in the monolayer (1L) limit, where strongly bound excitons provide a sensitive probe for many-body interactions. Here, we use two-dimensional electronic spectroscopy (2DES) to investigate many-body interactions and their dynamics in 1L-WS2 at room temperature and with sub-10 fs time resolution. Our data reveal coherent interactions between the strongly detuned A and B exciton states in 1L-WS2. Pronounced ultrafast oscillations of the transient optical response of the B exciton are the signature of a coherent 50 meV coupling and coherent population oscillations between the two exciton states. Supported by microscopic semiconductor Bloch equation simulations, these coherent dynamics are rationalized in terms of Dexter-like interactions. Our work sheds light on the role of coherent exciton couplings and many-body interactions in the ultrafast temporal evolution of spin and valley states in TMDs.


Experimental setup
Pump-probe and two-dimensional electronic spectroscopy (2DES) experiments are performed using a home-built, non-collinear optical parametric amplifier (NOPA).This is pumped by a fiber amplifier (Tangerine V2, Amplitude Systèmes) operating at 175 kHz and generating 260-fs pulses (Full-Width at Half-Maximum, FWHM) at 1030 nm.Fig. S1 displays the NOPA layout, which is adapted from Ref. 1.A small fraction of the 1030-nm pump beam is directly used to generate a supercontinuum in a 4-mm-thick YAG crystal (yellow beam in Fig. S1).The remaining portion of the pump beam is converted to its third harmonic (TH) via a set of two BBO crystals.The first crystal frequency-doubles the pump beam (cut angle  = 23.4°,thickness 3 mm) to generate the second harmonic (SH, green beam in Fig. S1).Then the SH is combined in a second BBO crystal with the remaining fundamental light ( = 32.5°,thickness 2 mm) to generate the TH (blue beam in Fig. S1).The YAG seed and TH pump are focused into another BBO crystal ( = 37.0°, thickness 1 mm) under an angle of 6.5° to drive a non-collinear optical parametric amplification process.In order to allow for broadband amplification, the chirped white light supercontinuum pulse is compressed to match the duration of the pump pulse using a pair of chirped mirrors (DCM9, Laser Quantum).The amplified seed beam (red beam in Fig. S1) is pre-compressed to compensate dispersion using another pair of chirped mirrors (DCM9, Laser Quantum).
The 2DES setup used for the experiments was discussed in Ref. 2. The NOPA output is split into a pump and a probe beam via a beam splitter.Afterwards, the pump beam goes through a mechanical chopper system (MC2000B, Thorlabs) with a custom-made wheel (500 slots) synchronized with the laser source.
Chopping is performed at a modulation frequency of 43.75 kHz, one quarter of the repetition rate of the laser.1c).The NOPA pulses are used in a home-built 2DES setup (right). 2 Key components of the 2DES setup are: (i) a passively phase-stabilized TWINS interferometer, (ii) a fast line camera detector operating at 87. 5

kHz and (iii) a home-built high-repetition-rate mechanical chopper.
Using an in-line interferometer based on birefringent wedges (TWINS, Translating-Wedge-Based Identical Pulses eNcoding System) 3 we then generate a pair of phase-stable pump pulses with tunable pulse separation (coherence time ) to record 2DES data.A small fraction of the pump light is sent to a photo diode behind the interferometer that records a field autocorrelation for each coherence time scan.The relative delay between the pump-pulse pair and the probe can be set via a retro reflector on a motorized translation stage (M126.DG1, Physik Instrumente).This delay between the second pump pulse and the probe is denoted as the waiting time .The pump and probe beams are focused onto the sample under a small few-degree angle using an off-axis parabolic mirror (OAP, reflected focal length of 50 mm).Both are vertically polarized.Their spot sizes are characterized with a CMOS camera to be 12x21 µm² for the probe and 19x22 µm² for the pump.The reflected beams are collimated using the same OAP.The probe beam is spectrally dispersed in a grating monochromator (Acton SP2150i, Princeton Instruments) and detected with an attached fast (up to 126 kHz) and highly sensitive line camera (Aviiva EM4, e2v).Using the chopper wheel to periodically modulate the pump pulses in pairs of two (denoted as on or off) and the fast line camera (operating at 87.5 kHz), we record differential reflectivity spectra: with a rate of 43.75 kHz as a function of the coherence time , waiting time  and detection energy   , measured using the grating monochromator.In a pump-probe measurement,  = 0 fs and differential spectra are recorded as a function of .For a 2DES measurement, a  scan is performed for each T, and differential spectra are recorded on the fly.By performing a Fourier transform of the recorded interferogram along  and taking the real part of this signal we obtain absorptive 2DES maps: These are energy-energy maps of the differential reflectivity as a function of the excitation energy   , and detection energy   , recorded for each .
A cross correlation second harmonic frequency-resolved optical gating (SH-FROG) measurement between pump ( = 0 fs) and probe at the sample position, taken with a 10-µm BBO crystal, is in Fig. S2, giving a ~9fs experimental pulse duration at the sample position.All experiments are performed at room temperature (RT).

Figure S2: SH-FROG measurement of the cross correlation between pump and probe pulses at the sample position.
The retrieved pulse duration is 9 fs.

Sample preparation
5][6] 1L-WS2 flakes are prepared by micro-mechanical exfoliation from bulk 2H-WS2 (HQ Graphene source) on Nitto Denko tape, 7 and then exfoliated again on a polydimethylsiloxane (PDMS) stamp placed on a glass slide for inspection under an optical microscope.Optical contrast is optimized to identify 1L-WS2 prior to dry transfer. 8Selected flakes are aligned and stamped on the desired location with a micro-manipulator at 60°C, before increasing the temperature to 80°C, so the flakes detach from the PDMS and adhere preferentially to the substrate.Raman and photoluminescence (PL) spectroscopy (Fig S3) confirm that the investigated flake with a size of ~50x70 µm² is a 1L-WS2.The Raman spectrum shows the main Raman modes at ~418.8±0.2 cm -1 (A ' 1), ~356.8±0.2 cm -1 (E') and ~352.1±0.2 cm -1 (2LAM(A)) and 8 additional Raman modes that are expected for 1L-WS2 under 514 nm laser excitation as discussed in Ref.
9. The difference in Raman shift between the A ' 1 and the E' peak of ~62 cm -1 is also consistent with 1L-WS2. 9Fig.S3b shows a PL emission at ∼616.1 nm, also consistent with 1L-WS2. 9

Fluence study
To ensure that all experiments are performed well below the Mott transition 10 and within the  (3) regime, a fluence study is conducted.Fig. S4 presents pump-probe measurements ranging from 5 to 243 µJ/cm².Within this almost two orders of magnitude variation in fluence, no significant changes to the pump-probe maps can be observed (Figs.S4a-c).A closer inspection of the lineshape and dynamics (Figs.S4d-e) reveals minor fluence-dependent changes.By increasing the incident pump power, the A exciton resonance broadens.This is consistent with our interpretation that many-body effects due to the increase in excitation density, such as EID, 11 dominate the optical nonlinearity of 1L-WS2.In addition, the dynamics of the A exciton shows a flattening of the initial decay with increased fluence, but no changes of the early, ~100-fs dynamics.No indications of exciton-exciton annihilation processes or changes to the delayed rise of the pump-probe signal are seen.The most striking effect of the increase in fluence is shown in Fig. S4f.The integrated pump-probe signal at the A exciton resonance for  > 150 fs has a strong deviation from a linear behavior for fluences >30 µJ/cm², and even a saturation of the signal strength at ~200 µJ/cm².This might reflect the influence of many-body effects, such as EID, 12 rather than going beyond the  (3) regime.
The experiments in the main manuscript are performed at a pump fluence of ~30 µJ/cm² and a probe fluence of ~20 µJ/cm², with the latter being sufficiently weak to not affect the observed pump-probe signal.
The ultrafast coherent oscillations observed at the B exciton resonance in Fig. 2c are also seen up to the maximum fluence of 243 µJ/cm² (Fig. S4c).To estimate the excitation density, we take into account the laser spectrum and the linear reflectivity of the sample (Fig. 1c).We can thus estimate the number of photons remaining in the sample.By taking the photon energy and pump fluence we estimate an excitation density of ~1.7 ⋅ 10 12 cm -2 .

Homogeneous and inhomogeneous broadening
As seen from the experimental 2DES maps (Fig. 3), the sample shows signatures of an inhomogeneous broadening contribution to the exciton linewidths, especially noticeable for the diagonal A exciton peak.
Using the 2DES data we can quantify the homogeneous and inhomogeneous contribution to the lineshape using the diagonal and crossdiagonal peak profile. 13Fig.S5 presents these crosscuts for the diagonal A exciton feature for  = 0 fs.A Gaussian fit to the experimental profiles (black and blue lines) yields FWHM of ~16 and 33 meV for the crossdiagonal and diagonal curves, respectively.The crossdiagonal lineshape is strongly affected by EID, thus the value resulting from the fit represents a lower limit for the homogeneous broadening.While the crossdiagonal profile usually is a marker for the homogeneous linewidth and therefore the dephasing rate   of the exciton, the diagonal profile marks ensemble effects due to the finite focus size such as local strain 13 or disorder 11 in 1L-TMDs.Therefore, the sample does exhibit a finite amount of inhomogeneous broadening.We account for the effect of experimentally observed inhomogeneity in exciton energy in the pump-probe and 2DES simulations (see section 5), using a pure dephasing time for the  = ,  excitons:  2, * = 1/  (S3) with dephasing rate   .We use  2, * = 50 fs and  2, * = 20 fs to reach good agreement with the experimental data.This corresponds to ℏ  = 13.2 meV and ℏ  = 32.9meV, in agreement with theoretical predictions. 14

Simulation of pump-probe and 2DES data
To support the interpretation of our experimental pump-probe and 2DES data, we compute the timedependent density matrix  ̂ of our system upon interaction with the electric fields of up to three ultrashort laser pulses in order to calculate the sample polarization in the time-domain. 15For this, we numerically solve the master equation in Lindblad form: 16,17 to obtain the evolution of the density matrix ρ ̂ under the total Hamiltonian  ̂=  ̂ +  ̂int (), given as the sum of the system Hamiltonian  ̂ and the light-matter interaction Hamiltonian  ̂int ().System-bath interactions, such as population relaxation and dephasing phenomena, are accounted for via Lindblad operators  ̂. 16,17  this non-perturbative approach, [18][19][20] the system interaction with the coherent laser pulses can be described in the dipolar approximation, 21,22 via the interaction Hamiltonian  ̂ () = −Ê() that couples the transition dipole moment operator ̂ with the electric field: 23 This total electric field consists of up to three Gaussian laser pulses (pump 1, pump 2 and probe) with amplitude  0, , pulse duration  = 5 fs and frequency   = 2.15 eV/ℏ.Additionally, the phases   can be set for phase cycling, and we can shift the pulses in time via   ′ .We label the delay between the two pump pulses as =  2 ′ −  1 ′ , and the delay between the second pump pulse and the probe pulse as  =   ′ −  2 ′ .The time  is the detection time, defined to be zero at the arrival time of the probe pulse, and is used for the numerical integration of Eq. (S4).From this numerical integration we obtain the timedependent density matrix  ̂(), which allows us to calculate the time-domain polarization: 21,22 () = (ρ ̂()).
(S6) as the trace of the expectation value of the transition dipole moment operator.To isolate the signals corresponding to the linear and nonlinear data measured in the experiment, we calculate both linear  (1) () (pump off) and total polarization   (, , ) (pump on).From these, the linear and total susceptibilities can be derived as: 24 (1) (  ) = 1  0 ℱ( (1) ())/ℱ(  ()) (S7) via Fourier transforms along the detection time using the vacuum dielectric constant  0 .We thus obtain   .The total susceptibility includes all possible interactions with the three laser pulses, thus also signals not measured in the experiment, since they are not radiated into the probe direction.9][20] Here,   (, , ) is repeatedly calculated for 4 phase settings , , 3 2 ], keeping   = 0.By averaging these, we obtain    (, ,   ) and can isolate the experimentally accessible nonlinear susceptibility 25   (, ,   ) =    (, ,   ) −  (1) (  ) (S9) by removing the linear contribution from the total phase-cycled susceptibility.Although the experiments are performed in a reflection geometry, we probe absorptive lineshapes for the excitons since the 1L-WS2 is placed on a Ag mirror.The linear reflectivity therefore follows the imaginary part of the sample susceptibility.Thus, we take the imaginary part of the complex nonlinear susceptibility as our 2DES signal: 16,17  2 (, ,   ) = ℑ(  (, ,   )).
(S16) Here, each class of excitons (A, B) interacts with independent baths via an individual Lindblad operator.

Modeling of many-body effects
The many-body Hamiltonian introduced here (see Fig. 4a) contains two 1-quantum (1Q) states |  ⟩ and |  ⟩ and three two-quantum (2Q) states |    ⟩, |    ⟩ and |    ⟩.The mixed |    ⟩ state suppresses 2DES cross peaks and coherent oscillations via the collective ground state.This is achieved by canceling ground state bleaching (GSB) and stimulated emission (SE) pathways with ESA of opposite signs.The 2Q excitations |    ⟩ and |    ⟩ are introduced to phenomenologically account for many-body effects, 26 such as EIS, EID, or a finite bleaching of the transitions from the 1Q to the 2Q manifold.These 2Q states do not reflect bound biexciton states.In this framework, EIS is reflected by a finite energy difference between the 2Q state and the sum of the constituent 1Q excitations.EID implies an increase in dephasing when going from 1Q to 2Q states, while Pauli blocking can be accounted for by reducing the amplitude of the 1Q→2Q transition dipole moment.All three mechanisms affect the balance between 1Q and 2Q transitions and transform the Hamiltonian from a linear to nonlinear one in the applied optical field.
From the experiments, we estimate that EID dominates the nonlinearity of the sample, therefore we introduce this as the microscopic source of the nonlinearity in our simulations.Qualitatively similar results are also obtained when adding a finite amount of EIS or bleaching.The density-dependent linewidth broadening 11 and saturation of the nonlinearity seen in the fluence-dependent study of Fig. S4 are not covered by our Hamiltonian.They would require density-dependent many-body interactions and/or the inclusion of higher-lying quantum states.
To account for EID effects, we modify the system-bath interactions via  ̂, .After this modification, the altered terms in the operators read: where   have been replaced by modified rates   ′ in the first two and   ′′ in the last two cases.
These four terms represent two different EID mechanisms.The first two lines of Eq. (S17) alter the systembath interactions of the |    ⟩ and |    ⟩ state, thus reflect intra-species (A-A or B-B) manybody interactions.These terms are clearly required to reproduce the lineshape of the diagonal features of the excitons in the 2DES maps in Fig. 3.This EID mechanism neither introduces (A,B) or (B,A) cross-peaks in a 2DES map nor coherent oscillations in the dynamics.
The second case, the bottom two lines in Eq. (S17), describes inter-species interactions, where an A exciton population increases the dephasing of the B exciton transition and vice versa.This interaction directly affects the mixed |    ⟩ 2Q state responsible for suppressing signatures of A-B exciton coherences via ESA transitions.Thus, changing the amplitude, dephasing or energy of the transition from 1Q states to |    ⟩ has a similar effect as introducing a direct coupling between A and B excitons.Inter-species EID can introduce cross-peaks and coherent exciton oscillations in the nonlinear experiment without the presence of A-B exciton populations, due to a lack of explicit A-B coupling.Such a model was used in Ref. 27 to describe a possible effect of EID and EIS on cross-peaks between A and B excitons in 2DES maps of 1L-MoS2.In Ref. 27, a comparison between experiment and simulation suggested that this effect could not explain the observed cross-peaks and their dynamics.Instead, a direct coupling between A and B excitons of ~28 meV was found to reproduce the experimental observations.To simulate our experiments, we introduce only the first class of intra-species EID.Coherent couplings between A and B excitons are taken into account by a direct coupling term   , rather than by introducing indirect coupling via an inter-species EID or EIS.Such a phenomenological model reproduces our key experimental observations.It contains several free parameters to account for many-body nonlinearities.
While the chosen parameters seem realistic, based on existing microscopic calculations, 28 this does not imply that other sets of parameters may not give similar or better agreement with experiments.

Simulation parameters
We set up the system assuming the same transition dipole moments for both A and B excitons,   =   , with   ′′ =   , so that no effect of inter-species EID is introduced.We further add a fixed amount of intraspecies EID by setting   ′ = 1.1  , increasing the homogeneous linewidth by 10% (Fig. S6).All additional simulation parameters are listed in Table S1.

fs
The Hamiltonian (Fig. 3a) only contains states up to the 2Q manifold.The light-matter interaction energy  0 is kept sufficiently weak to only introduce a depletion of the ground state population on the order of 10 −7 after excitation with both pump pulses ( = 0).

Inhomogeneous broadening effects
To account for inhomogeneous broadening, the time-dependent simulations of the master equation are repeated while varying the energies of the A and B excitons by an Δ ranging from -30 to 30 meV in steps of 10 meV.A weighted averaging of these 7 pump-probe and 2DES maps then introduces inhomogeneous broadening, giving the simulations of Fig. 4. For this, we determine a Gaussian distribution with  = 14 meV standard deviation.Crosscuts of the simulated 2DES map at  = 0 fs are in Fig. S5 (dashed lines), in good agreement with experiments (solid lines).A small energy shift of ~5 meV between the diagonal profiles of experiments and simulations reflects a finite displacement of the A exciton feature from the diagonal position, not captured by the simulation.Experiments and simulations show good agreement for the crossdiagonal profiles.The deduced value for ℏ  ≈ 13.2 meV gives a homogeneous FWHM of ~26.4 meV in linear reflectivity.This is larger than the ~16 meV linewidth seen in the crossdiagonal 2DES lineshape (Fig. S5).This apparent narrowing is a consequence of EID as the dominant source of exciton nonlinearity.The pronounced narrowing effect of EID on the differential reflectivity profile is exemplarily shown in Fig. S6 for a Lorentzian without inhomogeneous broadening contribution.

Simulations based on semiconductor Bloch equations
To analyze the role of Dexter coupling microscopically, we use the semiconductor Bloch equations 29

Figure S1 :
Figure S1: Experimental 2DES setup and pulse generation.A home-built NOPA (left) is pumped with 260-fs pulses of 50 µJ energy at a center wavelength of 1030 nm and 175 kHz repetition rate.The output of the VIS-NOPA covers a spectral range from ~515-690 nm (Fig.1c).The NOPA pulses are used in a home-built 2DES setup (right).2Key components of the 2DES setup are: (i) a passively phase-stabilized TWINS interferometer, (ii) a fast line camera detector operating at 87.5 kHz and (iii) a home-built high-repetition-rate mechanical chopper.

Figure S3: a :
Figure S3: a: Raman and b: PL spectra of 1L-WS2 under 514 nm excitation at RT.

Figure S4 :
Figure S4: Fluence study of the pump-probe maps recorded on 1L-WS2 at RT. a-c: Pump-probe maps at (a) 10, (b) 50, (c) 150 µJ/cm².d: Spectral crosscuts at T = 200 fs.A broadening of the A exciton resonance with increasing pump fluence due to EID can be seen.e: Dynamics at the A exciton resonance.For high fluences >100 µJ/cm² a slight slowdown of the initial fast decay component is observed, with no signatures of exciton-exciton annihilation.f: Average pump-probe signal for T > 150 fs.The pump-probe signal saturates at ~200 µJ/cm², with a deviation from a linear fit for >40 µJ/cm².Experiments are performed at a pump fluence of ~30 µJ/cm², corresponding to an excitation density of ~1.7 ⋅ 10 12 cm -2 .

Figure S5 :
Figure S5: Estimating homogeneous and inhomogeneous FWHM by comparing diagonal and crossdiagonal cuts through the 2DES map at T = 0 fs for experiments (black and blue solid lines) and simulations (red and green dashed lines).

Figure S6 :
Figure S6: Effect of excitation-induced dephasing on reflectivity.a: Simulated reflectivity for the A exciton assuming a homogeneous linewidth that corresponds to a dephasing time of ~50 fs (black line) and a line-broadening by ~10% via EID (red line).b: Differential reflectivity as calculated from panel a.The FWHM of ~12.8 meV is much smaller than that of the absorptive linear reflectivity of ~26.4 meV.
consider optically driven interband transitions between the topmost valence band and the lowest conduction band of a 1L-TMD semiconductor.The states of electrons and holes are described by the carrier wave vector  centered around the high-symmetry point  = (, ′) and spin  = (↑, ↓).   is the bandgap energy and with the effective masses of electrons and holes,   /ℎ , the reduced mass is introduced via 1/  = 1/  e + 1/  h .The dephasing of exciton states is described via Γ.Furthermore,Σ  () = −2 ∑  | ′ +− ′ |   ′ ′ ()  ′ ,′(S19)accounts for renormalizations of the electron and hole energies.The factor two in (S18) and (S19) appears due to the degeneration of the electron and hole occupation probabilities.Moreover,Ω  () =   () + ∑  | ′ +− ′ |  ′ ,   ′ ′ () (S20)is the generalized Rabi energy.It contains the optical driving term of interband transitions via the dipole coupling   and the optical field ().For the former we employ an analytical expression from a simple two-band model derived in Ref.30.The Coulomb matrix element is given by  || =  || / || where  || is the bare Coulomb interaction and  || the dielectric screening due to the environment and the monolayer itself.For the latter an analytical expression derived by Ref. 31 is used, which describes a TMD layer of a finite width separated from a substrate below by an airgap ℎ int .In (S19) and (S20) the terms with ′ =  describe intravalley Coulomb interaction while those with ′ ≠  account for intervalley interaction which

Table S1 :
Parameters used in the simulations.