Phonon-assisted oscillatory exciton dynamics in monolayer MoSe2

In monolayer semiconductor transition metal dichalcogenides, the exciton-phonon interaction is expected to strongly affect the photocarrier dynamics. Here, we report on an unusual oscillatory enhancement of the neutral exciton photoluminescence with the excitation laser frequency in monolayer MoSe2. The frequency of oscillation matches that of the M-point longitudinal acoustic phonon, LA(M). Oscillatory behavior is also observed in the steady-state emission linewidth and in timeresolved photoluminescence excitation data, which reveals variation with excitation energy in the exciton lifetime. These results clearly expose the key role played by phonons in the exciton formation and relaxation dynamics of two-dimensional van der Waals semiconductors.

and relaxation of excitons remains largely unexplored. This knowledge is important for interpreting a wide range of 2D exciton phenomena and for exploring the potential of exciton-based 2D optoelectronics.
In this work, we investigate the role of exciton-phonon interaction in exciton dynamics using the model system of monolayer MoSe 2 (Fig. 1a). Performing photoluminescence excitation (PLE) spectroscopy, we observe that the neutral exciton PL intensity, as well as its linewidth, oscillates as a function of excitation energy with a period matching that of the longitudinal acoustic phonon at the M point, LA(M). Nested within the oscillations are fine structures, with linewidths one order of magnitude smaller than that of ordinary PL, originating from resonant Raman scattering.
Analysis of the emission lineshape of the neutral exciton reveals that the oscillatory behaviour also presents in the homogeneous linewidth. Moreover, time-resolved PLE shows that exciton dynamics varies with respect to excitation energy, where shorter emission lifetime is measured for off-phonon-resonance excitation. This might due to the elevated lattice temperature arising from long-wavelength (small k-vector) acoustic phonons, which enhances radiative

RESULTS AND DISCUSSIONS
In our steady-state measurement, we detect PL at 5 K while scanning the excitation energy of a continuous wave (CW) laser, i.e., PLE spectroscopy (see Method Section for experimental details). The PLE intensity plot of Fig. 1b shows excitonic emission energies as a function of laser excitation. Two luminescence peaks are identified 23 : the neutral A exciton (X 0 ), centered at 1.650 eV, and the negative trion (X -), centered at 1.624 eV. Evidently, the intensity of X 0 emission oscillates as a function of excitation laser frequency, while the behavior of Xis monotonic. Fig. 1c shows PL spectra taken at the excitation energies 1.699 eV (red) and 1.686 eV (green), which exemplify the contrasting excitation energy dependencies of X o and X -PL. Within our laser scan range, which has a high-energy limit of 1.77 eV (700 nm), five equally spaced regions of luminescence enhancement, indicated by the white arrows, can be seen in X 0 , with an average energy separation of 18.5 meV. Such oscillations of X 0 emission intensity in PLE, first reported in CdS 24 for longitudinal optical phonons, is the hallmark of resonant excitation of phonon modes. indicating neutral exciton (X 0 ) and trion (X -) emission centered at 1.650 eV and 1.624 eV, respectively. Arrows indicate regions of PL enhancement. Colour bar: counts per second. c PL spectra at two distinct excitation energies showing the variation of X 0 (but not X -) PL with excitation energy.
A closer look at Fig. 1b shows that each PLE resonance region contains several narrow peaks. Fig. 2a offers an expanded view of the spectral regime highlighted by the white square in Fig. 1b. Three narrow lines shift in parallel with the excitation laser detuning, implying a Raman scattering origin of these lines. Their sub-meV linewidths are consistent with "conventional" Raman spectra measured on a different monolayer MoSe 2 sample (Supplementary Discussion), as well as with those reported in the literature [25][26][27] . The combined spectral features of PL emission and Raman scattering give rise to the overall emission spectrum, as shown in the example of Fig. 2b. As with earlier studies in monolayer WSe 2 14 , on top of the spectrally broad features (conventional X 0 PL) sits a narrow peak arising from resonant Raman scattering. The intensity of broad X 0 PL changes gradually with excitation energy, resulting in a rising PLE background on both ends of the scan range, as apparent in Fig. 2c. This observation has been reported 28 and is most probably due to increased absorption near excitonic resonances, e.g., 1s excitonic state (1.650 eV) on the low energy side and 2s (1.830 eV) on the high energy side. According to ab initio calculations reported on monolayer MoS 2 and WS 2 29,30 , the electron-phonon interaction strength is largest for LA phonons in the vicinity of the M points. Since monolayer MoSe 2 is structurally similar to MoS 2 and WS 2 , we expect mode specific characteristics of electron-phonon coupling to qualitatively resemble those of these compounds 30 . Therefore, we assign the oscillation in PLE as overtones (harmonic series) of the LA(M). This assignment is corroborated by plotting the PL intensity at the neutral exciton resonance as a function of excess energy ( Fig. 2c), defined as the photon energy difference between the excitation laser and the X 0 resonance (1.650 eV).  As required by momentum conservation during the Raman process, the total phonon wavevector of each phonon enhanced PLE peak indicated in Fig. 2c must be zero. This requirement is easily fulfilled by and modes, but not the fundamental LA(M) mode. However, for the LA(M) overtones, the requirement can be satisfied through the following scenarios. In the case of 2LA(M), a combination of M and wavevectors conserves momentum. In 3LA(M), this requirement is met following the scheme shown in Fig. 3b, where an equilateral triangle is formed by three Mvectors, also resulting in a zero vector sum. In monolayer MoSe 2 , in addition to the K and valleys (band edges), the conduction band also has Q and valleys located close to halfway between and K/ points. The momentum carried by an LA(M) phonon therefore matches the momentum separation between the K/ and /Q valleys ( Fig. 3a and b).
Thus, following the involvement of M-point phonons, conservation of momentum stipulates that the electron be scattered between K-and -valleys; see Fig. 3c. In other words, phonon-assisted scattering occurs between the optically bright exciton X 0 with both the electron and hole in K valley, and the optically dark indirect exciton X 0 (M) with an electron in and a hole in K valley (Fig. 3c). Here, X 0 (M) can be a virtual intermediate state such that its energy is not required to match that of X 0 + LA(M). Besides, despite the estimated 0.2-eV 33 electron band energy difference between K and Q valleys, the larger effective mass of the Q-valley 34 might result in a larger exciton binding energy of X 0 (M) than that of X 0 . This binding energy difference can partially cancel the electron band energy difference between Q and K valleys, leading to a smaller energy separation between X 0 (M) and X 0 , which enhances the role of X 0 (M) as an intermediate state. A seemingly related intervalley exciton-phonon scattering is proposed to explain the strong 2LA(M) peak in excitation-dependent Raman spectroscopy of WS 2 32 , although the excitation therein is well above band edge and involves higher lying conduction bands.
We note that while oscillations due to phonon resonance feature prominently in the X 0 transition, the Xemission intensity is relatively constant, except for excitation below 1.68 eV, close to the X 0 resonance of 1.65 eV. The lack of oscillatory enhancement in Xis possibly due to its distinct radiative properties compared to X 0 , together with the availability of multiple formation pathways 35 (e.g. via the exciton-electron interaction following exciton relaxation 36 ).
For X 0 , only those inside the light cone ( ) can radiate. In contrast, Xwith a much larger range of momentum can radiate due to the electron recoil effect 23 . As a result, X 0 PL intensity depends strongly on its momentum distribution, as determined by the resonance condition of the excitation energy, while such dependency is weak for X -. Moreover, the Xformation process is largely independent of the excitation energy, because even for excitation away from the phonon resonances, optically dark excitons can still be generated at large momentum (outside the light cone), which can interact with electrons to form trions. The lack of sensitivity to the excitation energy in both trion formation and relaxation processes diminishes any oscillatory features in the Xemission. A more detailed discussion can be found in the Supplementary Discussion.
Aside from spectral information, the PLE map also offers insights into the exciton dynamics. From the fit to the spectrum taken at each excitation energy, we found that both Xand X 0 lineshapes are well-described by Voigt profiles, from which one can infer the homogeneous linewidths of the excitonic transitions, as well as the widths of the Gaussian-broadened spectral distributions of their resonances (Supplementary Discussion). The latter is associated with the inhomogeneous broadening of the excitonic transitions. Oscillatory behavior is found in the homogeneous linewidth, γ 0 , of X 0 , which is smaller for excitation resonant with phonon harmonics (Fig. 4a). Its inhomogeneous (Gaussian) width remains relatively constant, consistent with the expectation that inhomogeneous broadening should depend only weakly on excitation energy. γ 0 is associated with the coherence lifetime of the exciton, and is given by where γ is the inverse of exciton population lifetime and Γ the pure dephasing rate. Since Γ is proportional to the rate of dephasing processes such as exciton-phonon scattering 22 , it is reasonable to assume that the oscillations in γ 0 is largely due to the creation of long-wavelength phonons during the relaxation of excited (hot) excitons. Our analysis on timeresolved PLE data (more details discussed below and in the Supplementary Discussion) seems to support this interpretation. To explore the phonon-assisted dynamics directly, we measure time-resolved emission of X 0 and Xwith a streak camera. Fig. 4b presents an example of the measured spectra with the pulsed excitation centered at 1.732 eV. The time evolution of the emission is characterized by a rapid onset within a few picoseconds, followed by an exponential decay. This is shown by extracting time traces along P -P' and Q -Q' from  37,38 , leaving behind the dark excitons, which are then scattered into the light cone by phonons (Fig. 5a) and recombine at a later time, producing the observed exponential decay. When the excitation is off-resonant, excitonic relaxation results in the emission of many long-wavelength phonons (Fig. 5b), forming a phonon bath that increases the scattering rate. On the other hand, on-resonance excitation (Fig. 5c) produces only a small number of LA(M) phonons which are ineffective in the aforementioned scattering process due to momentum mismatch. This picture is consistent with the oscillations of γ 0 shown in Fig. 4a, i.e., the creation of long-wavelength phonon bath during off-resonance excitation increases exciton-phonon scattering and manifest as a broadening of the X 0 homogeneous linewidth. Procedures for PL lineshape fitting, data post-processing and analysis, and further discussions can be found in the Supplementary Discussion accompanying the paper.

DATA AVAILABILITY
Data that supports the findings of this study is available from the corresponding authors upon reasonable request.

I. PL lineshape fitting
The PL lineshapes of both neutral exciton (X 0 ) and trion (X -) are found to be inhomogeneously broadened in such a way that the broadening is on the same order of magnitude as the natural linewidth. Therefore, neither Lorentzian nor Gaussian lineshapes adequately fit the spectra, as shown in Fig. S1a. The shortcoming of these lineshapes is most evident in the "wings" of the resonances, where the best-fit Lorentzian curve overestimates the signal while the Gaussian curve underestimates it. As mentioned in the main text, we fit the spectra with Voigt profiles, or the plasma dispersion functions. Explicitly, the Voigt profile, , is given by the convolution between a Gaussian curve, , and a Lorentzian curve, : Here, A is a proportionality constant, and the average (centre) and standard of deviation of the normally (Gaussian) distributed resonance, respectively, and the half width at half maximum (HWHM) of the Lorentzian. For every PL spectrum, two separate sets of fitting parameters, A, , , and are needed for the X and X 0 resonances. Fitting is accomplished by minimizing the squared error by performing multivariable nonlinear regression. Some spectra feature prominent narrow lines attributed to resonant Raman scattering. These are fitted manually with Gaussian curves. Since the linewidths of these narrow features are much less than that of X and X 0 , their presence has minimal impact on the accuracy of the Voigt profile fits. As an example, the full-spectrum fit of Fig. 2b in the main text is shown in Fig. S1a below, along with fits using Lorentzian and Gaussian curves. The excellent fits produced with Voigt profiles for both X and X 0 resonances justify their use.  Fig. 4a in the main text, reproduced here for convenience: best-fit of the X 0 resonance. c The centre of the X 0 resonance as a function of excess energy. Shaded regions indicate regions of enhanced luminescence obtained from Fig. 1b. Error bars represent 99% confidence intervals, defined as ±2.58 × standard error of fitted parameter.  Fig. 4b. Red dashed curve is the raw streak camera measurement while blue solid curve the deconvolved data, with the IRF represented by the dotted line. b Excitation frequency-dependent rise times of X -(crosses) and X 0 (circles) emission extracted from a series of streak camera measurements in varying excitation energy. Error bars represent the standard deviations. c Excitation frequency-dependent decay rate, γ, of X -(crosses) and X 0 (circles). d Same data as in c, presented as the lifetime, τ. Reproduced from Fig. 4e.
The best-fit values of the fitting parameters are given in Table S1. Interestingly, the best-fit and for X 0 vary in an oscillatory manner with respect to the excitation photon energy, whereas those of X remain relatively constant. As shown in Fig.  S1b, of X 0 is reduced when the excitation laser is on a phonon resonance. Surprisingly, Fig. S1c reveals that also oscillates with a period corresponding to the LA(M) phonon mode. However, the oscillations are not in phase with that of because here the rising edges coincide with the phonon resonances. Artifacts from fitting can be ruled out with further analysis of the fitting error, e. g. from the lack of corresponding oscillations in the residue. The amplitude of oscillations roughly equals that of , but much smaller than the average homogeneous linewidth. This might be related to the oscillations of , in which efficient hot exciton relaxation and the subsequent smaller phonon bath is conjectured to prolong radiative lifetime (Section IV), hence reducing , with onphonon-resonance excitation. Nonetheless, the physical mechanism behind the modulation of exciton resonance seen here remains to be explored.

Table S1
Best fit parameter values for both Xand X 0 transitions. With the exception of and for the X 0 transition, all other parameters given here are found to vary only slightly with respect to the excitation energy, and without a clear oscillatory feature.

II. Exciton rise time and lifetime extraction
In time-resolved PLE measurements, due to the finite response time of the streak camera, the output is the convolution of the time-dependent emission with the instrument response function, IRF. Since the timescale of MoSe2 exciton dynamics is in the picosecond range, on the same order of magnitude as the nominal resolution (1 ps) of the streak camera, the measured time evolution of the emission will be noticeably stretched by about 1 ps. To obtain more accurate timing parameters, we perform deconvolution on the streak camera output with linear least-squares regression 1 . By scattering light from a modelocked Ti-Sapphire laser with 150 fs pulse width into the streak camera, we measured a 1.5 ps full-width-at-half-maximum (FWHM) response. This is shown in Fig. S2a (dotted line) and taken as the IRF for deconvolution. To avoid overfitting and to improve numerical stability of the deconvolution protocol, Tikhonov regularization 1 is applied. With an appropriate choice of Tikhonov factor, the noise in the reconstructed temporal profile can be suppressed without affecting the overall response, as exemplified by the time traces of the trion emission, before (dashed line) and after (solid line) deconvolution, in Fig. S2a. The same deconvolution procedure is applied to all streak camera images taken at different excitation frequencies. Expectedly, we find that the rise times of X 0 and X emission are reduced by roughly 1 ps following deconvolution, as compared with those extracted from the raw data. Nonetheless, their variations with respect to excitation frequency remain qualitatively similar with or without deconvolution.
From the deconvolved data, we obtain the rise times and lifetimes of X 0 and X with respect to excitation energy as presented in Fig. S2bd. The rise time is simply defined as the time interval for which the rising edge of the luminescence moves from 10% to 90% of peak emission. To obtain the decay rate, a linear fit to the logarithm of the falling edge is applied, where the slope gives the decay time constant, γ. The extracted excitation dependent decay rates for the first sample is shown in Fig. S2c. To facilitate comparison with the rise time, the data in the main text is presented as the lifetime, τ = 1/γ. For the conversion of the error bars from ε(γ) to ε(τ), the first order approximation is adopted, where denotes the mean value of γ. The result of the conversion is reproduced here in Fig. S2d.

III. Exciton and trion decay rates in time-resolved luminescence
The exciton dispersion curve is shown in Fig. S3a, where only excitons inside the light cone ( ) can radiatively recombine. To model the exciton decay process, we consider a simplified three-level system as shown in Fig. S3b. represents the dark exciton state outside the light cone ( ), with time dependent density , which nonradiatively decays to the vacuum state with a rate . is the bright exciton state inside the light cone ( ), with time-dependent density , which decays to , both radiatively and nonradiatively, with a total rate . The nonradiative decay rate includes relaxation from excitons to trions. Scattering with other excitons, free carriers, impurities, and long-wavelength (small k) phonons can induce interconversion from to and vice versa, with the corresponding rates given by and , respectively. Given the ~meV linewidth of excitons inside the light cone 2 , interconversion from to can be induced by both phonon emission and absorption, whose coupling matrix elements are proportional to and , respectively, where is the number of k-vector phonon. Consequentially, and increase with . and are governed by the coupled rate equations Solving the above rate equations, we find where , , and are time-independent constants determined by the initial values . The PL emission rate is given by , which has two decay time constants, and . Under a low temperature we can take 3 such that , the fast (slow) decay constant ( ) is then given by Now corresponds to the total decay rate of the bright excitons inside the light cone, and that of the dark excitons outside the light cone. : dark exciton, : bright exciton, : vacuum state. c Dispersion relation of the trion (solid black curve) and its radiative decay (red wavy arrows). Unlike the exciton, the radiative decay rate of the trion varies slowly with k, such that the scattering induced by longwavelength phonons (green curved arrows) barely changes the overall decay rate. The thickness of the wavy arrows illustrates varying decay rate with respect to wavevector.
Theory has shown that bright excitons in monolayer TMDs have a very short radiative lifetime 3 , as evidenced by the measured decay rate 2,4 of . This value is significantly larger than the inverse of PL rise time (0.2 -1 ps -1 ) and decay rate (0.2 -1 ps -1 ) shown in Fig. S2b and c. Possibly, in our time-resolved measurements, by the time the recorded luminescence intensity reaches its maximum, the bright exciton density inside the light cone is depleted. The luminescence detected at a later time predominantly comes from scatteringinduced dark-to-bright exciton conversion, whose decay rate is given by . The variation of PL decay time ( ) mainly comes from the modulation of the scattering rate, , which increases with the density of long-wavelength phonons. A related behavior due to the weak scattering rate at low temperature, termed "exciton-phonon relaxation bottleneck", has also been discussed in a recent paper 5 . In contrast to the exciton, trion has rather different radiative properties 3 . In particular, a trion with large momentum can radiatively decay, as the momentum conservation can be satisfied with the left-behind electron inherits the trion wavevector (the electron recoil effect 6 ). The radiative decay rate varies slowly with the trion wavevector 3 , such that scattering with long wavelength-phonons barely changes the overall decay rate of the trion (see Fig. S3c). Therefore, unlike the exciton, trion does not show lifetime oscillation with the excitation energy.

IV. Exciton and trion steady-state dynamics
The contrasting behaviors of exciton and trion luminescence in PLE with excitation energy can be explained by their distinct radiative properties. For the excitons, momentum conservation allows only those inside the narrow light cone ( ) to radiatively recombine. These bright excitons have an ultrafast radiative decay rate 3 of about 3 ps -1 . On the other hand, trions in a much larger momentum range can radiatively recombine, with a smoothly changing decay rate 3,6 .
The PLE spectra are determined by the steady-state population distribution of excitons/trions. We expect the steadystate total exciton/trion population to depend weakly on the excitation energy because: (i) on-or off-phonon resonance affects only the kinetic energy of the generated exciton, while the generation rate is nearly unaffected as it is determined by the LA(M) phonon emission rate; (ii) most of the trions are the outcome of subsequent relaxation from excitons 7 which is insensitive to the exciton kinetic energy (rather, trion formation rate is mostly determined by phonon emission rates of and other optical phonons 8 ). Nonetheless, the excitation energy strongly affects the exciton population distribution in k-space. When the excitation is off-phonon resonance, it generates excitons with a finite kinetic energy, resulting in excitons with a wide steady-state kspace distribution (see Fig. S4a). A large number of excitons have high energies, and cannot be efficiently scattered into the light cone for radiative recombination due to energy conservation. On the other hand, when the excitation is onphonon resonance, it generates excitons with close to zero kinetic energy. The subsequent exciton steady-state k-space distribution is narrower (Fig. S4b). When the excitation laser is tuned across phonon resonances, the variation of exciton distribution in k-space leads to the observed oscillations in PLE luminescence spectra.

Fig. S4
Exciton and trion steady-state k-space distributions for a onand b off-phonon-resonance excitation. The upper parabola represents the dispersion curve of an exciton together with a free electron (X 0 + e), while the lower that of a trion (X -). Orange arrows indicate laser generated exciton energies following phonon emission (see Fig. 5b & c in the main text), while blue curved arrows illustrate the relaxation process from a neutral exciton to a trion. Red wavy arrows represent the radiative decay process, whose thickness corresponds to the wavevector-dependent decay rate. Only excitons lying in the vicinity of valley extremum can radiatively recombine, while all trions can radiatively recombine with a slowly changing decay rate. Shaded Gaussian curves correspond to the exciton and trion k-space distributions.
Most trions are formed via the relaxation of excitons 7 . This additional step can partially reduce the difference of trion kspace distributions between exciting on-and off-phonon resonance. Meanwhile, trions in a much larger momentum range can radiatively decay with a slowly changing decay rate (Fig. S4). This further diminishes the variation of trion emission as a function of excitation energy.

V. Data from a second monolayer MoSe 2 sample.
Data from a second monolayer MoSe2 sample (Fig. S5a) exfoliated from a different bulk crystal is presented in Fig. S5b g. In this sample, the resonances for both X and X 0 are blueshifted by about 3 meV relative to the first sample. Nonetheless, the phonon resonances are identical (within the uncertainty of 1 meV) to those found in the first sample, as can be seen in the vertical line cut at the X 0 resonance shown in Fig. S5c. Here, we have included additional Raman peak assignments in addition to the LA(M) overtones where space allows. These assignments are provisional, especially for higher order modes, because multiple phonon combinations may give rise to the same energy.
The emission lineshapes of X 0 and X in the second monolayer MoSe2 is different from those in the first in the sense that they are slightly asymmetric. As a result, fitting using two Voigt profiles as is done in the previous case becomes inadequate here. Nonetheless, we found that excellent fits can be produced by adding two Gaussian curves whose peak positions and widths are constant, i.e., independent of the excitation energy, as shown in Fig. S5d. Their peak positions are each about 10 meV lower than X 0 and X resonances. Given their large linewidths of about 7.5 meV, they are unlike to be originating from defect states. It is more likely that their presence is due to contamination or inhomogeneity as a result of lattice distortion. Despite this, HWHM of X 0 emission, γ0, extracted from the lineshape analysis shows oscillations with a period close to the LA(M) phonon mode (Fig. S5e), as is the case in the first sample. However, oscillations are also noticeable in the inhomogeneous broadening width, σ, as shown in the inset of Fig S5e, although the amplitude of oscillations (about 0.15 meV) is much smaller than that of γ0. This is likely an artifact of the fit due to the weaker X 0 emission intensity and the more pronounced Raman scattering components present in the second sample. Regardless, the oscillatory behavior of the homogeneous linewidth and of the peaks position of X 0 (Fig. S5f) qualitatively agrees with results obtained from the first sample. Fig. S5 a Optical micrograph of a second monolayer MoSe2. Scale bar: 10 μm. b PLE spectra of the second sample, showing neutral exciton and trion emission centered at 1.653 eV and 1.626 eV, respectively. Colour bar: counts per second. c Vertical line cut P -P' of the PLE spectra shown in b with selected phonon peaks marked by arrows. d An example of the PL emission spectra with fit. Solid black circles are data points while the solid red line is the fit constructed from the following components: Voigt profiles for trion (solid blue line) and neutral exciton (solid greed line), and the Gaussian background (solid dark and light grey). e Best-fit homogeneous linewidth, γ0, of the X 0 resonance as a function of excess energy extracted from lineshape analysis of the PLE intensity map in b. Inset: Width of the inhomogeneous (Gaussian) broadening, σ. f Peak position of the X 0 resonance as a function of excess energy. Error bars in both e and f represent 99% confidence intervals of the fits, while shaded regions correspond to regions of enhanced luminescence seen in b. g Excitation frequency-dependent decay rates, γ's, of X -(crosses) and X 0 (circles) extracted by fitting the exponential decay time of the emission.

VI. Raman spectrum of a third monolayer MoSe 2 sample
Raman spectrum of a third monolayer MoSe2 sample is shown in Fig. S6 below. The spectrum is taken at room temperature, with an excitation laser wavelength and intensity of 514 nm and 2 mW, respectively. The resonances observed agree with the spectra obtained by other groups [9][10][11] , as well as our PLE data shown in Fig. 2c in the main text. The linewidths of the prominent peaks, e. g. and , are found to range from 0.5 to 1 meV, matching those in the PLE spectra.