Impact of environment on dynamics of exciton complexes in a WS$_2$ monolayer

Scientific curiosity to uncover original optical properties and functionalities of atomically thin semiconductors, stemming from unusual Coulomb interactions in the two-dimensional geometry and multi-valley band structure, drives the research on monolayers of transition metal dichalcogenides (TMDs). While recent works ascertained the exotic energetic schemes of exciton complexes in TMDs, we here employ four-wave mixing microscopy to indicate that their subpicosecond dynamics is determined by the surrounding disorder. Focusing on a monolayer WS$_2$, we observe that exciton coherence is lost primarily due to interaction with phonons and relaxation processes towards optically dark excitonic states. Notably, when temperature is low and disorder weak excitons large coherence volume results in huge oscillator strength, allowing to reach the regime of radiatively limited dephasing and we observe long valley coherence. We thus elucidate the crucial role of exciton environment in the TMDs on its dynamics and show that revealed mechanisms are ubiquitous within that family.

Introduction In spite of their illusory academic simplicity, synthetic two-dimensional (2D) materials -such as graphene, black phosphorous, and transition metal dichalcogenides (TMDs) -display stunning properties, which are also revealed in their optical responses. For instance, in monolayers (MLs) of TMDs, the reduced dielectric screening and 2D carrier confinement give rise to exotic, non-hydrogenic excitons with binding energies exceeding 0.2 eV [1], which is an asset enhancing light-matter interaction. The latter is manifested by a strong absorption and subpicosecond population lifetime, favoring formation of surface plasmon polaritons [2] and exciton-polaritons [3][4][5] with a valley degree of freedom -to name a few examples illustrating a technologydriven progress in the optics of TMDs. However, there is a need for an in-depth understanding of fundamental mechanisms governing exciton radiative and nonradiative recombination rates in various experimental settings. There is a large spread of reported values of exciton coherence and population decay [6] and little is known about their dependence on microscopic material properties and environmental factors, such as temperature, strain, dielectric surrounding and excitonic disorder on different length scales, generating inhomogeneous broadening σ.
The main obstacle to access these information, was a large size of the optically probed areas (typically, diameter of a few tens of micron), which are required to implement traditional approaches of nonlinear spectroscopy -such as angle-resolved four-wave mixing (FWM) -inferring decay times of populations T 1 and coherent polarizations T 2 in extended samples. We here overcome this difficulty, by exploiting phase-sensitive heterodyne detection. The latter permits to perform FWM spectroscopy in a microscopy configuration, attaining spatial resolution of 300 nm. Using a tungsten disulphide (WS 2 ) ML, exhibiting the strongest optical activity among all other TMDs ML [7] we observe a giant FWM response of the resonantly generated excitons and we carry out the mapping of their dephasing time (T 2 = 2 /γ), population decay time and σ. We further infer the dephasing induced by phonons, by performing FWM temperature dependence.
Additionally, two distinct types of trions (bound states of two electrons and one hole) are unambiguously identified in FWM. We show that a single electron is the ground state for optically active trions, where the additional electron and hole are within the same valley as this ground state electron (intra-valley trion), or in the opposite valley (inter-valley trion). An energetic splitting between these states due to exchange interaction was recently predicted [8], and observed [9] for WSe 2 and for WS 2 [10]. We observe the Raman quantum beats [11] resulting from this splitting, revealing coupling between both types of trions. We employ this phenomenon to measure the decay of the valley coherence T valley 2 [12][13][14][15], which appears to be signiicantly longer that previously reported.
Spectral characteristics of transitions We first perform micro-reflectance from a flake, to identify exciton (EX) and trion (TR) transitions, as shown in Fig. 1. In further experiments we employ the FWM micro-spectroscopy setup [16][17][18] (see Methods), adapted to the visible spectral range. We probe the flake with the E 1,2,3 pulsed laser beams which are spectrally centered at either EX (∼ 590 nm) or TR (∼ 600 nm), with arXiv:1709.02658v2 [cond-mat.mes-hall] 11 Oct 2017 a bandwidth about 7 nm (full width at half maximum, FWHM). The reference E R beam is focused on the surrounding SiO 2 , so that its lineshape is not affected by a strong absorption of the flake. Fig. 1 presents the resulting FWM spectral interferograms obtained on both resonances in the WS 2 flake at T=5 K. We note that the amplitude of the TR is typically an order of magnitude weaker than the EX's one. Below the TR line we further retrieve FWM of another type of valley-trion [19]. The fringe period in Fig. 1 is given by a few pico-seconds (ps) delay between the reference pulse E R and the FWM emission. Its intensity and phase are retrieved by applying spectral interferometry. The former as a function of E 1 intensity is shown in the inset, yielding the limit of the third-order χ (3) regime (where further experiments are performed) up to around 100 nW. It is worth to note that FWM can be readily detected with E 1 as low as 1 nW, generating a low carrier density of a few 10 8 cm −2 . At excitation stronger than 100 nW, FWM visibly starts to saturate, which we attribute to the absorption bleaching. In WS 2 MLs, the optically active exciton (EX) has a larger transition energy than the dark one [20], such that at low temperature the PL of EX is suppressed [10], as shown in the Supplementary Fig. S1. While this issue remains relevant in view of competing relaxation channels of the bright exciton, it is not an obstacle to drive its FWM: EX are resonantly and selectively created, generating a giant response, owing to the µ 4 scaling of the FWM, where µ is the oscillator strength.
Homogenous and inhomogeneous widths of the neutral exciton The EX spectral lineshape measured in reflectance (Fig. 1) is dominated by inhomogeneous broadening with a Gaussian distribution of around 15 meV FWHM. FWM spectroscopy has been conceived primarily to access the homogeneous broadening FWHM γ in an inhomogeneously broadened ensemble, exhibiting a spectral standard deviation of σ. The complex conju-  The non-exponential behavior is due to a constant offset, giving a noise level. c) The retrieved homogeneous broadening versus temperature γ(T) fitted with model described in the text, dashed lines represent parameters obtained from the fit. d) Exemplary measurement of the EX density dynamics for a given position. The initial decay T1 is mainly due to radiative recombination T rad and non-radiative relaxation T dark to the EX dark state, such that (T1 gate in the FWM definition, imposes phase conjugation between the first-order polarization induced by E 1 and the FWM. For τ 12 > /σ, its transient appears as a Gaussian, known as a photon echo. It is centered at t = τ 12 and has a FWHM, corrected with respect to the pulse duration, equal to 8 ln (2) /σ. Formation of such an echo is illustrated in Fig. 2 a, where time-resolved FWM amplitude of EX versus τ 12 is shown. The echo develops during the initial 0 < τ 12 < 0.5 ps. By inspecting FWM for later delays, for example τ 12 = 0.7 ps (orange trace) we retrieve inhomogeneous width of around 11 meV (FWHM). Time-integrated amplitudes of the photon echo as a function of τ 12 for different temperatures are reported in Fig. 2 b. Here, the decay reflects γ. To retrieve γ, the data are fitted with an exponential decay exp (−γτ 12 /2 ), convoluted with a Gaussian to account for a pulse duration of 0.12 ps. At T=5 K we obtain γ EX = (2.1 ± 0.1) meV, yielding dephasing time T 2 = (620 ± 20) fs. Thus EX in WS 2 shows a larger homogenous width than its counterpart in recently investigated MoSe 2 MLs [18], in line with a superior linear absorption in WS 2 with respect to MoSe 2 [7]. The temperature dependence of γ is illustrated in Fig. 2 c. The data are modeled [21](purple trace) with the following equation: The linear term [γ 0 = (2.1 ± 0.1) meV, a = (0.018 ± 0.003) meV/K] is attributed to low energy acoustic phonons. The latter term, with the b = (32 ± 6) meV, is attributed to thermal activation of optical phonons with dominant or mean energy E 1 = (37 ± 3) meV, which indeed supply a large density of states above 300 cm −1 37 meV [22]. The phonon dephasing mechanisms are therefore similar as in MoSe 2 MLs [18] and as in semiconductor quantum wells. Above T= 210 K, the FWM decay is limited by the temporal resolution, such that γ cannot be extracted, although a strong FWM is measured up to the room temperature.
Population decay A representative measurement of the population dynamics (FWM response versus τ 23 ) via three-beam FWM [17,18] is shown in Fig. 2 d. The measurement was performed at the same spot as the dephasing study, presented in Fig. 2. The population dynamics is dominated by an initial exponential decay with a constant of T 1 0.35 ps, followed by a longer dynamics described by two additional exponential decays [23] (T 1 slow 4.7 ps and T 2 slow 46 ps) that we can relate to phase space distribution via scattering processes and scattering back from the exciton dark ground state. We note that the portion of secondary excitons, decaying on a nano-second timescale, is at least an order of magnitude larger than on recently studied MoSe 2 MLs. This we associate with a dark exciton ground state in WS 2 and its bright character in MoSe 2 [20].
The obtained result (T 2 2T 1 ) indicates populationlimited dephasing. We attribute this fast decay to be due primarily to the fast radiative decay. Indeed, excitons in ML TMDs possess the radiative lifetime T rad of a few hundred femto-seconds, as recently revealed via two-colour pump-probe [24], FWM [18,25] and TMD polaritons studies [3][4][5] -all these results signify a large EX transition dipole moment and coherence volume spanning across many Bohr radii [26]. The parameter T dark describing phase space distribution via scattering processes and relaxation to the dark exciton ground state is also expected to contribute to the fast initial decay. Other nonradiative recombination processes are expected to be of minor impact, as they are not faster than the decay of secondary excitons, that is ≥46 ps (we assume that these processes have the same dynamics for both bright and dark excitons). We also note a weak role of phonons on the excitons dynamics in this low temperature range (see Fig. 2 c). To get a comprehensive view of the possible mechanisms influencing the exciton dynamics, the local insight into σ, T 2 , population decay and µ is required, and should be strengthen by imaging of these quantities across the entire flake. Crucially, such an original capability is offered by the heterodyne FWM microscopy. Thus, we now focus on the FWM mappings and analyze spatial correlations between the above parameters.
Four-wave mixing mapping In the first step, we focus on spatial variation of σ. We therefore acquire FWM spectral interferograms at τ 12 = 0.6 ps and retrieve time-resolved FWM amplitude of EX, while scanning over the flake surface. For each location, we inspect the width of the photon echo, from which we measure σ. The result is shown in Fig. 3 a. In the middle of the flake, we identify regions of a smaller inhomogeneous broadening, down to 10 meV, yet still largely dominating over γ. It is worth to note, that the largest σ, and thus most pronounced exciton localization, is measured at the borders of the flake. This is related to the strain gradients and variations of the dielectric screening by the substrate, which are expected to be strong along the edges. These locations are preferential for wrinkling, local deformations and lattice defects creating deep potential centers trapping individual emitters (see supplementary Fig. S1 d). As an origin of σ, we point toward a local strain and charges trapped on a flake. We note that suspended ML flakes displayed comparable σ (see supplementary Fig. S3), excluding the interface roughness between the SiO 2 and the flake as a principal source of inhomogeneity. Newly, it has been found that σ is reduced by encapsulating a flake in hexagonal boron nitride [27,28]: FWM performed on such heterostructure indeed has revealed reduction of σ, however its complete cancelation has not been observed (not shown). Moreover, a reduced spectral jitter was observed on non-insulating substrates [29], helping to get rid of trapped charges, further indicating decisive role of charge fluctuation on the inhomogeneous broadening.
In the next step, we focus on the area exhibiting a large variation of σ, marked with a square in Fig. 3 a. Again we perform mappings, this time varying also the delay τ 12 for each position. Such retrieved T 2 and σ are presented as color-coded maps in Fig. 3 b and c, respectively. Their spatial correlation is striking and emphasized in Fig. 3 d: the fastest dephasing is measured at the areas of smallest σ. We interpret this as follows. Center of mass of two-dimensional excitons moves within a disordered potential landscape [30], arising from the variation dielectric contrasts, strain from the substrate, uncontrollable impurities, vacancies, etc. Through the Shrödinger equation, the disorder acts on the wave-function localization in real space and thus results in its delocalization in k-space, modifying radiative rates with respect to free excitons. In other words, the disorder mixes the states inside and outside of the radiative cone, and thus creates a distribution of states with an oscillator strength reduced as compared to ones fully in the radiative cone, and spread in energy, adding up to σ. Note, that this localization is weak [31] comparing to localization resulting from deep traps [32], resulting in a distinctive emission band well below the EX emission (see Supplementary Information Fig. S1 d). As presented in Fig. 3 d, T 2 starts to decrease only for a sufficiently low σ, where the ra- diative decay time T rad becomes fast enough to compete with another channel, identified as the EX relaxation to the dark ground state. Such channel was not observed in MoSe 2 displaying a bright exciton ground state (see Supplementary Fig. S3 a) and MoS 2 (not shown). Conversely, for largest σ, the non-radiative decay dominates, as the T rad is increased through the localization. These spatial correlations demonstrate that radiative rates and dephasing of excitons in ML of TMDs are governed by exciton localization imposed by a local disorder. The slope of the transition energy owing to the strain gradient is observed in Fig. 3 e, where the center transition energy is encoded in a hue level. In Fig. 3 f we present time-integrated FWM amplitude of the EX transition (corrected by the excitation lineshape) reflecting µ. Comparing Fig. 3 f with Fig. 3 a, we note that the areas of the smallest σ yield the strongest FWM (see also Supplementary Fig. 3). This is because with decreasing σ (disorder), the spatial overlap between excitons increases, enhancing the EX interaction strength and thus resulting in a more intense time-integrated FWM. We note that the largest oscillator strength is observed in areas with the lowest transition energy.
Valley-trion structure We now turn to investigation of the trion transition (TR) [10]. The latter is formed when an additional electron occupies the lowest conduction band, as depicted in the inset of Fig. 4 a. Depending on its spin, one can form a singlet state (intra-valley trion, intra-TR) or a triplet-state (intervalley trion, inter-TR) [10,33]. To address coherence dynamics of these TR, E 1,2,3 are prepared co-circularly, selectively addressing K+ valleys. We investigate it at the same spatial location as for experiments illustrated on Fig. 2, yielding low σ and marked with a cross on Fig. 4 a, e and f. Similarly as for EX, we obtain a single exponential decay yielding the averaged TR dephasing T 2 (TR) = (440 ± 10) fs. This faster dephasing with respect to EX is attributed to fast TR relaxation into the lower lying dark states, leaving an electron with a varying momentum, and inducing additional dephasing via final state damping.
When E 1,2,3 are co-linear, K+ and K-valleys are equally excited, as linearly polarized light contains both circularly polarized components. As a result, intra-TR and inter-TR transitions are activated in both valleys. In particular, intra-TR in K+ and inter-TR in K-valley share the same ground state (corresponding to the presence of an electron in the lowest conduction band, labeled with a yellow arrow in the inset in Fig. 4 a), forming a coupled V-type system. In such a configuration, E 1 and E 2 generate valley coherence between both types of trions resulting in the Raman quantum beats [11] , as sketched in the inset in Fig. 4 a. Owing to the TR singlet-triplet splitting [10], labeled as ∆ ST , the phase of this coherence evolves when increasing τ 12 . Therefore, the measured coherence dynamics, shown in Fig. 4 a, displays beating [11] with a period T ST 0.71 ps, yielding ∆ ST = 2π /T ST = 5.8 meV.
To measure the TR density dynamics we perform the τ 23 -dependence of the FWM [18] for a fixed τ 12 = 0.2 ps. Once again, via co-circular excitation we selectively probe the intra-TR dynamics. As presented in Fig. 4 b, it displays a fast monotonous decay, owing to the relaxation into the dark exciton state, followed by a slower revival of the FWM generated by the dark excitons relaxing back into the light cone [18]. This interpretation is strengthened by comparing the TR dynamics in WS 2 with the one measured in MoSe 2 ML (see Supplementary  Fig. S4), exhibiting optically bright ground state. In the latter case, as the TR recombination leaves an electron with a varying momentum and relaxation to lower energy state is not possible, we observe a substantially longer TR radiative decay with respect to EX.
It is worth pointing out that no FWM is observed (ei-  ther at TR or EX) for opposite circular polarizations of E 1,2 and E 3 pointing towards a robust valley polarization in WS 2 MLs.
Upon co-linear excitation, the initial density dynamics displays again an oscillatory behavior: the first two pulses do not only create the intra-TR and inter-TR populations, but also induce the valley coherence between them, once more generating the Raman quantum beats [11]. All these ingredients contribute to the FWM signals, which we model with a phenomenological fitting curve (see Supporting Information) and extract the valley coherence dephasing time of T valley 2 = (1.3 ± 0.2) ps. Such a record value, much longer than previously reported [12][13][14][15] and measured on location with low σ confirms recent reports suggesting that a shallow disorder potential plays a critical role in the exciton valley coherence [34]. For longer delays the valley coherence has dephased, such that the subsequent exciton dynamics is similar for both polarization configuration. Thus both TR transitions of WS 2 are here unveiled via FWM, whereas they cannot be distinguished in reflectance owing to the spectral splitting smeared out by σ TR . Similar values of ∆ ST were retrieved when varying the position on the flake.
To conclude, we demonstrated a giant nonlinear optical response of exciton complexes in WS 2 MLs. The substantial enhancement of the FWM retrieval efficiency with respect to standard semiconductor quantum wells, was exploited to unravel the impact of a local disorder in twodimensional systems onto exciton dynamics and dephasing. The valley degree of freedom was unveiled, when considering the trion structure. Our results indicate that coherent nonlinear microscopy is suited to explore optical properties of emerging optoelectronic and optomechanical [35] devices and heterostructures [27] made of layered semiconductors. An alluring perspective is to image the exciton coherent dynamics on a nanometer areas and to conjugate it with the structural properties of TMDs, which can be revealed down to atomic scale using scanning tunneling microscopy. This could be achieved by transferring our methodology towards the nanoscopic regime [36,37], offering a spatial resolution of a few tens of nanometers. ACKNOWLEDGEMENT We acknowledge the financial support by the European Research Council (ERC) Starting Grant PICSEN (grant no. 306387), the ERC Advanced Grant MOMB (grant no. 320590), the EC Graphene Flagship project (No. 604391) and the ATOMOPTO project within the TEAM programme of the Foundation for Polish Science co-financed by the EU within the ERDFund. We also acknowledge the technical support from Nanofab facility of the Institute Néel, CNRS UGA.

METHODS
We employ an optical parametric oscillator (Inspire 50 by Radiantis pumped by Tsunami Femto by Spectra-Physics) to create a triplet of short laser pulses around 600 nm: E 1 , E 2 and E 3 , with adjustable delays τ 12 and τ 23 , as depicted in the Supplementary Fig. S1. The three beams are injected co-linearly into the microscope objective (Olympus VIS, NA=0.6), installed on a XYZ piezo stage. They are focused down to the diffraction limit of 0.6 µm, onto the sample placed in a helium-flow cryostat. E 1,2,3 are pre-chirped by using a geometrical pulse shaping [17], so as to attain close to Fourier-limited, 120 fs pulses on the sample. The WS 2 ML flake was mechanically exfoliated from a bulk crystal purchased from HQgraphene and deposited on a 90 nm thick SiO 2 substrate.
The FWM generated within the sub-wavelength (approximately half of the waist) area, diffracts in all directions. There is therefore no k-vector matching condition, on which most FWM experiments rely on. Instead, our microscopy approach imposes the signal to be selected in phase, by performing optical heterodyning. By employing acousto-optic deflectors operating at different radio-frequencies Ω 1, 2, 3 , the phases within the pulse trains E 1,2,3 are modulated by nΩ 1, 2, 3 /ν, where ν and n denote the laser repetition rate and pulse index within the train, respectively. As a result, the FWM polarization -which in the lowest, third-order is proportional to µ 4 E 1 E 2 E 3 -evolves with the phase n(Ω 3 + Ω 2 − Ω 1 )/ν. This specific phase-drift is locked onto the reference pulse Supplementary Figure S2. Photoluminescence hyperspectral imaging of a WS2 monolayer at T= 5 K. a) As the optically active exciton (EX) is an excited state, its PL is strongly quenched. TR labels the negative trion. P1 was recently interpreted as another type of valley-trion [19] coinciding with the dark EX. The origin of peaks P2, P3 and P4 remains unexplained. The reflectance is given in red, for comparison with the PL. The spatial imaging spectrally integrated across the peaks, as indicated, is shown in (b-g). L1 indicates spectrally narrow, strongly localized states, which are formed at the low energy side of P1. CW excitation at 450 nm with 5 µW at the sample.  Figure S3. Measurements at MoSe2 monolayer presented in Ref. [18] a) Correlation between homogeneous and inhomogeneous broadenings, indicating radiative rate governed by localization via disorder. Note that comparing to WS2 no saturation is observed for large inhomogeneous broadening values. We interpret the vertical spread of experimental values as resulting from activation of non-radiative recombination and fast scattering out of radiative cone, resulting in faster decoherence. Therefore, for each value of inhomogeneous broadening, only the upper point correspond to radiatively limited dephasing. b) Time-integrated FWM amplitude versus inhomogeneous broadening, the areas of the smallest σ yield the strongest FWM. This is because with decreasing σ (disorder), the spatial overlap between excitons increases, enhancing the EX interaction strength and thus resulting in a more intense time-integrated FWM response. To extract T dark we fit the following phenomenological dependence to data presented on Fig. 3 d: and obtain T dark = 0.78 ± 0.03 ps.  Crucially, the areas characterized by a fastest T1 also display short T2 and the smallest σ (compare to Fig. 3 a, b and c).

RAMAN COHERENCE FIT
We fit the FWM amplitude vs τ 23 dependence measured in colinear configuration on the TR, which is presented on Fig. 4 b, with a phenomenological model of the following form: where T valley 2 is the valley coherence decay constant, T fast and T slow describe the fast and slow population decay components, respectively, while A and B are their amplitudes, y 0 is the background signal and ∆ ST is the TR singlettriplet splitting. Table I presents the obtained fitting parameters. We performed fits on data obtained on different spots and values presented without experimental error here were typical values obtained from those fits. Therefore, only three free fitting parameters were used to fit data presented in the article (A, B and T valley 2 ).