Coherence and Density Dynamics of Excitons in a Single-Layer MoS2 Reaching the Homogeneous Limit.

We measure the coherent nonlinear response of excitons in a single layer of molybdenum disulfide embedded in hexagonal boron nitride, forming a h-BN/MoS2/ h-BN heterostructure. Using four-wave mixing microscopy and imaging, we correlate the exciton inhomogeneous broadening with the homogeneous one and population lifetime. We find that the exciton dynamics is governed by microscopic disorder on top of the ideal crystal properties. Analyzing the exciton ultrafast density dynamics using amplitude and phase of the response, we investigate the relaxation pathways of the resonantly driven exciton population. The surface protection via encapsulation provides stable monolayer samples with low disorder, avoiding surface contaminations and the resulting exciton broadening and modifications of the dynamics. We identify areas localized to a few microns where the optical response is totally dominated by homogeneous broadening. Across the sample of tens of micrometers, weak inhomogeneous broadening and strain effects are observed, attributed to the remaining interaction with the h-BN and imperfections in the encapsulation process.

Transition-metal dichalcogenides (TMDs) are lamellar compounds held together by van der Waals interlayer interactions. For this reason, they can be exfoliated down to a single layer (SL), similar to graphene obtained from graphite. Even though the interlayer interactions are weak, they have an important effect on the band structure of TMDs, moving them toward an indirect band gap. When thinning bulk crystals down to SLs, many TMDs can be converted to a direct band gap semiconductor, as first shown for MoS2. The discovery of efficient emission and absorption of light in SL TMDs, facilitated by creation of excitons (EXs) of high binding energy and their fast radiative recombination,[1,2] made two-dimensional TMDs candidates for next generation optoelectronics. Additionally, the symmetry and chemical composition of the atomic lattice of TMDs enable, besides their flexibility and partial transparency, a wealth of innovative application concepts.[3−7]

Among semiconducting TMDs, MoS2 has been the most studied because of its expected superior stability in atmospheric conditions. Several recent observations nevertheless challenge this expectation. The measured EX line width has been on the order of several tens of millielectronvolts, even at low temperatures.[4,6−12] These high values indicate an inhomogeneous broadening (σ) of the EX transition more than an order of magnitude above the homogenenous broadening (γ) expected in the millielectronvolt (meV) range.[13] The dominating inhomogeneities (σ ≫ γ) in SL MoS2 conceal the intrinsic properties of the EXs resulting from the underlying band structure, which is still under debate.[14−16] Presumed origins of σ observed in the experiments down to low temperature are adsorbed impurities and crystal defects such as vacancies, strongly affecting the quantum yield.[17,18] Another source of inhomogeneous broadening is the substrate onto which the two-dimensional material is typically deposited. The most common substrate used to fabricate optoelectronic devices, a thin thermal oxide (silica) layer on silicon wafers, is known to be corrugated and to contain charged impurities, generating a disordered potential landscape for EXs in the supported TMD layer.

Recently, deterministic transfer methods,[19] developed to stack two-dimensional materials in the form of so-called van der Waals heterostructures,[20] have been applied to prepare MoS2 SLs sandwiched between two thin hexagonal boron nitride (h-BN) layers.[13,21,22] In these structures, the distance between MoS2 and the charged impurities in SiO2 and possible adhered impurities on the top surface is set by the respective h-BN thicknesses and is typically several tens of nanometers (nm). Moreover, the van der Waals interaction is thought to promote a close, conformal and very flat contact at the MoS2/h-BN interface (with interfacial distances typically of a few angstroms)[23] by expelling adsorbed molecules sideways, much like a flat iron would eliminate the pleats of a clothing. In such samples, photoluminescence (PL) revealed sharp excitonic features of a few meV width, approaching the expected homogeneous limit.[13] The narrow emission allowed to deepen the understanding of the excitonic complexes involving different valleys[13] and to observe the EX’s excited states[22] in MoS2 SLs. The encapsulation strategy has also been employed to reduce σ in other TMDs.[21,21,24−26]

The question arises if σ could be suppressed with this method sufficiently to provide a dominating radiative broadening and a resulting long-range exciton-polariton formation. Toward this, can one show correlations between γ, σ, and population lifetime T1? Such fundamental issues are relevant for spectroscopists and material scientists exploring optical properties of TMD SLs and, in a broader context, for condensed matter physicists investigating two-dimensional systems. We note that the EX spectral line-shape measured in linear transmission or reflection[13,21,22] is a convolution of γ and σ. Separating them (for example, by applying line-shape fits such as a Voigt profile) requires prior knowledge of homogeneous and inhomogeneous line-shapes. For example, in the limit σ ≫ γ, the line-shape is Gaussian, such that γ cannot be reliably estimated. Interestingly, the EX emission measured in non-resonantly excited PL can show narrower line widths than those retrieved via resonant absorption. This can be understood as due to the carrier and EX relaxation selecting local potential minima at low temperatures prior to EX recombination.

To separate homogeneous from inhomogeneous broadening in the EX line-shape, nonlinear spectroscopy, specifically four-wave mixing (FWM), is particularly suited. FWM driven on an inhomogeneously broadened optical transition (for example, created by a spatially varying EX transition in TMDs) forms a photon echo.[27−30] Its temporal width is determined by σ, assuming that the excitation pulses are sufficiently short with respect to ℏ/σ. Conversely, its amplitude decay with delay time between the first two exciting pulses, i.e., τ12, is only due to the microscopic EX dephasing. In the case of a simple exponential decay, it determines the full width at half-maximum (fwhm) of the homogeneous line width γ = 2ℏ/T2, where T2 denotes the EX dephasing time. Only in the case of vanishing σ can the latter be read from the FWM transient, taking the form of a free induction decay, overcoming the necessity to scan τ12. This, however, is only known a posteriori. Because the inhomogeneous broadening is due to the spatial variation of the EX energy on the scale of its radius of a few nanometers, it can also vary on longer scales spatially across the sample surface. Furthermore, because the homogeneous broadening is also depending on this spatial variation, the measured pair (σ, γ) is a spatially varying quantity on a length scale with a lower limit given by the size of the sample region probed by the optical excitation, and an upper limit given by the size of the investigated sample. The spatial variations across the flake can be due to strain induced by the substrate or encapsulating layers,[31] the dielectric environment, the density of impurities and defects, and the free carrier concentration. These mechanisms give rise to the disorder, affecting EXs’ radiative rates[32] and, thus, γ. The disordered potential landscape results in varying EX localization lengths and produces different sets of EX energy levels,[33] determining σ and affecting a population lifetime T1.

Enhanced spatial and temporal resolution is required to experimentally investigate the above issues. It is thus instructive to study TMD SLs with FWM microspectroscopy, resolving the signal on a 100 femtosecond (fs) time and a 300 nm spatial scale.[29,30] In the employed implementation, the exciting laser pulses propagate colinearly in the same spatial mode, while the signal is discerned via optical heterodyning. This technique allows us to spatially resolve σ and γ, revealing correlations between EX’s dephasing and σ. The signal-to-noise ratio in this experiment is increased compared with more traditional two-dimensional systems, such as GaAs quantum wells, by the large oscillator strength μ of EXs in MoS2 SLs. The latter allows to generate substantial multiphoton nonlinear responses,[34,35] including FWM as its field amplitude and intensity scale as μ4 and μ8, respectively.

Results

In this work, we report FWM microscopy of two heterostructures composed of a SL MoS2 flake encapsulated by layers of high-quality h-BN.[36] Details regarding sample fabrication are provided in the Methods section. In contrast to MoSe2,[29] WSe2, and WS2[30] SLs, we find that encapsulation is essential to observe a strong, coherent, nonlinear optical response in exfoliated MoS2 SLs. By analyzing the FWM transients acquired in the first investigated sample that exhibits more disorder (sample A), we assess homogeneous and inhomogeneous contributions to the EX spectral line-shape. We find that the encapsulation leads to a global reduction of σ, down to a few meV, comparable to the homogeneous broadening. In some micron-sized areas of the sample, a decrease of σ correlated with an enhanced FWM signal can be observed. We exploit the varying inhomogeneous broadening across the sample to investigate the impact of EX disorder onto their coherence dynamics, occurring at a picosecond (ps) time scale. The retrieved correlations among σ, γ, μ, and T1 provide fundamental insights into EX localization and dynamics in a disordered two-dimensional landscape. Using the second sample, with significantly less disorder (sample B), we demonstrate the EX optical response in its homogeneous limit. Taking advantage of its high optical and structural quality, we show strong correlations between above quantities. We then assess the phonon-induced homogeneous broadening with increasing temperature and evaluate the excitation-induced dephasing with increasing EX density. Finally, the dynamics of the EX density after resonant excitation is unveiled. On a picosecond time scale, we resolve an initially dominating radiative decay,[29,37] competing with non-radiative redistribution into optically dark states. The remaining EX density is scattered back to the optically bright EX states and recombines (radiatively and non-radiatively) assisted by a three-exciton decay process visible on a nanosecond time scale. The coherent detection allows us to disentangle the EX populations in the different bright and dark states. They create a different phase of the response depending on the phase of the complex EX renormalization of the optically probed bright EXs due to Coulomb and exchange interaction with the total EX density. These contributions involve different EX populations, each corresponding to charge carriers in the inequivalent K and K’ valleys of the electronic band structure.[38] For example, broadening described by an imaginary part of the interaction is in quadrature with energy shifts described by the real part of the interaction. The different relaxation processes contribute to the FWM amplitude with different phases, imprinting unusual signatures in the measured density dynamics of the FWM amplitude due to constructive and destructive interferences.[38]

Characterization with Linear Spectroscopy

The optical pictograph of the sample A, containing an elongated MoS2 SL, displayed in Figure a, left panel, shows breaks (indicated by orange arrow) and wrinkles (red arrow) in the flake as well as air trapped in bubbles and puddles[23] (yellow arrow). In spite of these features, structurally clean areas extending across about 7 μm × 7 μm are found, such as the one enclosed with a white dashed-contour. For the initial characterization at T = 5 K, we perform hyperspectral imaging of the confocal PL (non-resonantly excited with a CW laser diode operating at 450 nm, with ∼10 μW arriving at the sample) across the spectral range around the EX emission. The experiment reveals EX center energies spanning across 30 meV, with a PL intensity varying over more than an order of magnitude and different proportions of neutral (EX) and charged (trion, TR) states, as shown in Figure S1. We also see zones in which both PL and reflectance of EX are suppressed, while a characteristic spectrally broad band[21] that was tentatively attributed to defects appears below the EX energy. A spatial map of the PL intensity across the entire flake for (636 ±  1) nm is shown in Figure a, right panel. Within this region we identify areas containing quite narrow EX emission, down to 7.1 meV fwhm, as exemplified in Figure b.

Figure 1

Optical response of the sample A, composed of a MoS2 SL embedded in between layers of h-BN. (a) Right: Spatial mapping of the confocal photoluminescence intensity (PL) at (636 ± 0.5) nm. Left: Microscope image of the sample under the white light illumination in reflection. The green dashed line indicates the position of MoS2 flake with a large extension of several tens of micrometers. Reflectance and PL measured in the area encircled with a dashed line are given in panel b. (c) Typical FWM spectral interferogram measured at the EX transition (green). The excitation intensity spectrum of the femtosecond laser is given in black. (d) FWM intensity dependence of the excitation power (in a two-beam configuration) showing the expected linear dependence with the pump power. power was fixed around 4 μW.

To inspect the linear coherent response, we performed microreflectance from the same zone, showing resonances at EX and TR as well as the B exciton, centered at ∼587 nm (not shown). In spite of the substantial improvement of the optical response with respect to previously examined free-standing MoS2 SLs,[13,39] the EX line-shape in this h-BN/MoS2/h-BN heterostructure is still affected by σ. To disentangle σ and γ, we employ three-beam FWM microscopy, inferring EX coherence and population dynamics from femtosecond to nanosecond time scales.

Coherence Dynamics via Four-Wave Mixing Microscopy

FWM spectroscopy in the k-selection regime of bare MoS2 SLs was reported in ref (28). Here, instead, we perform heterodyne FWM microscopy[30] in the configuration briefly described in the Methods section. A typical spectral interferogram of the two-beam FWM field (proportional to , where are the fields of the exciting pulses) from the neutral EX in our heterostructure, is presented in the green trace in Figure c. The measured spectrally integrated FWM intensity (amplitude squared) as a function of the power of the first arriving pulse (), given in Figure d, shows a linear dependence, consistent with the third-order regime of the FWM, up to 1 μW. The latter corresponds to an excited EX density of around 109 cm–2 per pulse, which is around 4 orders of magnitude below the EX saturation density in TMDs.

We now turn to the assessment of σ and γ by inspecting the time-resolved FWM amplitude as a function of τ12. The experimental results are given in Figure , left panels. In panel a, we see that FWM is observed for τ12 > 0 and already for τ12 > 0.2 ps takes a form of a Gaussian pulse centered close to t = τ12, i.e., the photon echo is formed. In an ideal case of a set of two-level systems and for δ pulses, the FWM signal for τ12 > 0, t > 0 can be described by a product of an exponential decay and a Gaussian shifting its maximum in time: |FWM(t, τ12)| ∝ exp[−τ12/T2] exp[−ν2(t–τ12)2/2). The echo has a constant temporal width, with a standard deviation of 1/ν and fwhm of. This quantity is linked with the fwhm of the spectral inhomogeneous broadening as . Conversely, the amplitude decay of the echo with increasing τ12 reflects the homogeneous dephasing time T2 = 2ℏ/γ. The measured signal is convoluted with the temporal duration of the applied pulses of about 150 fs, which is taken into account in the modeling presented in Figure b. The two-dimensional fit to the experimental data shown in panel a yields (γ, σ) = (2.10  ±  0.03, 12.5  ± 0.2) meV. We note that the FWM amplitude at pulse overlap (τ12 = 0), encoded in the hue level, principally reflects the EX oscillator strength μ.

Figure 2

FWM microspectroscopy carried out on sample A. (a, c, e) Time-resolved FWM amplitude for different delays τ12 at T = 5 K, showing formation of the photon echo. The disorder is decreasing from top to bottom, quantified by the increase of γ and decrease of σ. (b, d, f) Corresponding simulations.

To discuss the local character of the quantities (γ, σ, μ), let us now consider Figure e, where the FWM transient acquired within a distance of a few micrometers from the spot considered in panel a is shown. Here, the shape of the photon echo is different in several respects: (i) it is broader in real time, showing that σ is smaller; (ii) the amplitude decay along τ12 is faster, revealing a shorter T2 (and thus larger γ); and (iii) the amplitude around τ12 = 0 is larger by an order of magnitude, showing a larger μ. These changes are quantified by the fitted form of the echo, given in Figure f with the parameters (γ, σ) = (3.10  ±  0.08, 5.5  ±  0.3) meV at this position. It is worth noting that the response in panel e already deviates form the echo form, i.e., the maximum of the signal is not aligned along the diagonal, indicating a transition to a homogeneously broadened case. When further approaching this limit, the coherence dynamics displays a crossover from the photon echo toward the free induction decay, resulting in the biexponential decay, as discussed in the Figure S2. In Figure c, we present another case of the measured echo, with the intermediate line-shape parameters (γ, σ) = (2.60 ±  0.07, 7.7  ±  0.3) meV, as reflected by the simulation shown in the panel d. These examples demonstrate that the EX optical response is affected by the disorder on scales above and below the resolution of the present experiment, 300 nm. Below this resolution, the disorder leads to an effective inhomogeneous broadening of the response and a reduction of the radiative decay rate by the localization of the EXs below the optical resolution λ/2, thus leading to a mixing of dark EX states outside the radiative cone. Above this resolution, we can see the varying impact of disorder acting on EXs, as shown in the three examples.

To verify if the encapsulation can be used to virtually eliminate σ, we have processed a second heterostructure (sample B), presented in Figure S3. To exclude any aging issues, the FWM spectroscopy was performed within hours after completing the fabrication. At micron-sized areas, we measure FWM amplitude as narrow as 4.4 meV (fwhm) at T = 4.5 K, as shown by filled-blue peak in Figure a. No signatures of the photon echo can be noticed, when inspecting FWM transients versus τ12. Thus, here σ is not detectable and the broadening reaches its homogeneous limit. In that limit, time-resolved FWM takes a form of a free induction decay (FID), i.e., exponential decay from t = 0 for any τ12 > 0, with a decay constant determined by T2, as depicted by a scheme framed in panel b. A pair of explicit examples of the measured FID at 4.5 K for τ12 = 0.04 ps (yellow arrow) and τ12 = 0.6 ps (orange arrow) are given by the yellow and orange points in panel a. By deconvoluting the laser pulse duration, we retrieve T2 = (320  ±  20) fs, i.e., γ = 4.11 meV, close to the fwhm of the spectrally resolved amplitude, also in a stunning agreement with γ retrieved from time-integrated FWM as a function of τ12, shifting the focus of the discussion to panel b. Interestingly, therein at 4.5 and 10 K, we clearly detect the signal at negative delays with a FWM rise time of around 130 fs (after deconvoluting the pulse duration), close to theoretical prediction[40] of T2/2. Such contributions have previously been observed when studying homogeneously broadened EXs in GaAs quantum wells and were assigned to the local-field effect.[40]

Figure 3

Exciton coherence dynamics in a MoS2 heterostructure exhibiting low disorder (sample B), measured at different temperatures. (a) Top: FWM transients vs τ12 at T = 4.5 and 34 K showing FWM free induction decay, thus proving EX broadening in its homogeneous limit. Bottom: corresponding FWM spectra measured for τ12 = 0.5 ps directly showing the temperature-induced broadening. Red trace in the bottom-left panel is the spectral shape of the reference pulse . (b) Time-integrated FWM amplitudes as a function of τ12 measured for temperatures as indicated. The dephasing time T2 is measured from the exponential decay at τ12 > 0, and an increasing influence of the phonon dephasing is observed when increasing the temperature. At 4.5 K (blue circles), note the presence of the FWM signal for τ12 < 0. Frame: a scheme of the three-pulse FWM rephasing sequence employed to assess the coherence dynamics; τ23 is set to zero in the experiment. In the homogeneous limit (γ ≫ σ) of the EX broadening, FWM transient is a free induction decay (FID) instead of the photon echo. (c) Temperature dependence of γ retrieved from panel b (violet squares) compared to the spectral fwhm of the FWM amplitude (open green squares). The red curve is the fit to the data (see the main text).

A relevant factor influencing the EX coherent dynamics is the temperature, which determines the density of acoustic and optical phonons.[41] Increasing the temperature, and thus the phonon density, broadens γ by phonon-assisted scattering.[27,29,30,42] However, the changing phonon-scattering influences the EX relaxation dynamics.[38] To measure the impact of temperature on γ, we plot in Figure b the EX coherence dynamics, i.e., the time-integrated FWM amplitude as a function of τ12, for different temperatures. At this homogeneously broadened zone, the FWM amplitude decays as exp(−τ12/T2). With increasing temperature from 5 to 70 K, we measure shortening of the dephasing time and, thus, an increase of the homogeneous broadening from 4 to 23 meV, as marked by violet squares in panel c. In parallel, virtually the same broadening is seen in the FWM spectral amplitudes (green open squares): to directly illustrate the dominance of the homogeneous broadening mechanism through phonons at this sample position, in the bottom-right part of panel a, we present the data measured at 34 K. At higher temperatures, as exemplified for 90 K, the dephasing is too fast to be measured with our current setup, limited by the temporal resolution of around 150 fs. We fit the data with a sum of a linear and exponential activation terms: γ(T) = γ0 + aT + b/[exp(E1/kBT)–1]. For the linear term, attributed to acoustic phonons, we obtain γ0 = (3.78  ± 0.15) meV and a = (0.03  ±  0.01) meV/K. For the second term we find b = (35 ±  12) meV ad E1 = (8.3 ± 1.5) meV. In contrast to previously inspected SL TMD,[29,30] the value of activation E1 is here not consistent with the energy of optical phonons[41] of around 35 meV. We tentatively link such a more-pronounced temperature dephasing with a particularly small conduction band splitting in MoS2, favoring the population loss of bright excitons through their scattering to dark states.

While investigating dephasing mechanisms in TMD SLs, it is also instructive to determine the impact of EX–EX interactions on γ. Such interactions cause the broadening[27] of γ by 250% when increasing the EX density from 1010 to 1011 cm–2. In Figure we present dephasing curves measured for different excitaton powers P1,2,3, spanning across typically operating excitation range in our experiments. With increasing EX density from around 108 to 109cm–2, we detect a small but measurable increase of the homogeneous broadening by 15%; 0.5 meV over an order of magnitude density. Because the absorption (thus also EX density for a fixed P1,2,3 value) may vary across the flake, it is important to check that such excitation-induced dephasing is of minor importance in the applied range of P1,2,3: values of γ retrieved from the analysis of the spatially resolved FWM experiment, presented in the next section, are expected not to be significantly affected by spatially fluctuating EX–EX interactions.

Figure 4

Exciton coherence dynamics in a MoS2 heterostructure exhibiting low disorder (sample B), measured at different excitation powers. A stronger slope of the dephasing curves with increasing P1,2,3 is due to excitation-induced dephasing. Within our range of P1,2,3, we measure increase of γ by 0.5 meV when rising the EX density by an order of magnitude.

Four-Wave Mixing Mapping and Statistical Correlations

To draw a comprehensive picture of the balance between the microscopic disorder and the EX coherent dynamics, we go back to the more disordered heterostructure (sample A, Figure ). We use FWM imaging,[30] measuring photon echos, as in Figure , at the grid of spatial points of the sample. For each position, we perform two-dimensional measurements and fits as in Figure , extracting γ, σ, and the FWM amplitude, μ4. While we note that these fits are not taking into account the interaction induced nature of the signal, the resulting parameters mimic the response well. These parameters are presented as color maps in Figure a–c. Within the investigated region, γ (fwhm) spans from 1.8 to 4.3 meV. We see the spread of σ (fwhm) from 4.7 to 16 meV, with the less disordered area colored in green-blue in Figure a. This location also shows the strongest FWM amplitude, proportional to μ4, as shown in Figure b. Note that the upper limit of the measured σ is set by the spectral width of the excitation laser.

Figure 5

FWM spatial imaging performed in sample A. (a–c) Mapping of the homogeneous broadening γ, the inhomogeneous broadening σ, and of the FWM amplitude at τ12 = t = 0 proportional to μ4. The areas of the weakest disorder yield the smallest σ, the shortest T2 and largest γ, as summarized in panel e. Dashed lines in panels e and f are guides to eye. (d) Mapping of σ retrieved from fitting the spectrally resolved FWM with a Gaussian profile. (f) Correlation between the FWM amplitude measured at τ12 = t = 0 and σ. (g) Correlation between spectral FWM broadening and σ (as retrieved from the FWM delay and temporal dynamics) is visible. The diagonal is drawn with the dashed line.

A pronounced correlation between μ4 and σ is explicit in Figure f and quantified with a Pearson coefficient of . It is interpreted as follows. Smaller σ signifies weaker disorder enabling larger EX center-of-mass localization lengths,[33] thus generating large EX coherence volume in real space. In the reciprocal space, the EX wave function is thus dominated by small k-vector components and, therefore, is better contained within the light cone. This increases the light-matter interaction, and thus also μ, resulting in an enhancement of the FWM signal. Instead, at areas showing large σ, the EX spans more prominently out of the light cone. We emphasize that the smaller overlap with the light-cone results in decreasing μ (reflected by the FWM intensity) and increasing radiative lifetime: instances of such qualitative dependencies among μ, σ, and T1 are shown in Figure S4, while the statistical analysis of this effect is discussed in Figure . Close to the radiative limit (T2 = 2T1), this implies an increase of T2 and thus a decrease of the homogeneous line width γ with increasing σ. We seek for such a (σ, γ) correlation in Figure e. While some trend is visible, as marked with a dashed line, the quantitative statistical correlation is low, with . We do observe many points showing a short T2 and yet a large σ. This is attributed to other homogeneous broadening mechanisms, such as non-radiative exciton-electron scattering, which vary across the investigated area. We indeed checked that experimental points placed on the upper-right quadrant in Figure. e are measured at the areas displaying trion transitions in PL and thus exhibiting enhanced electron doping level.

Figure 6

Impact of the exciton localization range on the homogeneous line width and population decay rate measured at the low disorder limit. (a) Correlation between σ and γ, demonstrating a decrease of the homogeneous line width (increase of dephasing time T2) with increasing exciton localization range via disorder. (b) Correlation between σ and T1, demonstrating increase of the population decay time T1 (composed of radiative and non-radiative recombination channels) with increasing exciton localization range via disorder. The measurements are performed on sample B. Orange dashed lines are guides to eye.

We note that the low-disorder sample B does not show any trion absorption and emission and, thus, is characterized with a lower and more-homogeneous residual electron doping. For this reason, we have measured spatially resolved FWM, while also scanning τ12 and τ23 for each position, and performed correlation analysis as on the Sample A. In Figure a, we present the resulting (σ, γ) correlation. We unambiguously observe that with increasing σ, from values close to zero up to around 10 meV, γ decreases from 3 meV down to ≃1.5 meV. Quantitatively, the moderately strong correlation is found reflected by . The final step consists in conjugating σ with the initial EX density decay, to which we attribute population lifetime T1. The latter is inferred from τ23 dependence (see Figure S4). The result is shown in Figure b. We observe that with increasing σ, T1 increases as well, generating a strong statistical correlation with the Pearson coefficient . These experimental results obtained on sample B, confirm our assertions regarding mutual dependencies between μ, σ, γ and T1.

It is instructive to compare the parameters obtained through the temporal dynamics with the ones from the spectra. In Figure d, we present the inhomogeneous width (fwhm) of the Gaussian profiles we fit to the spectrally resolved FWM data. We find agreement between the inhomogeneous width determined by the FWM transients (panel a) and the spectral fwhm (panel d); the correlation between both is shown in panel g. In Figure , we further exploit correlations and demonstrate that σ, characterizing amount of disorder, increases with the EX center energy (obtained from the spectra). This means that the short-range (sub-resolution) disorder is dominating within the probed region and is of repulsive nature.

Figure 7

Correlations between the EX center energy, inhomogeneous σ, and homogeneous broadening γ. One can observe that the center energy varies within the line width. The center energy shifts by the amount of σ, indicating that the disorder occurs on subresolution scale and has a repulsive character. The measurement was performed on sample A.

Population Dynamics on a Nanosecond Scale

After the resonant excitation, several relaxation mechanisms play an important role.[38] They are governed by, on the one hand, the high oscillator strength[43] and thus fast radiative recombination rate and, on the other hand, by the EX conversion toward different dark states, resulting from the peculiar valley structure and available scattering channels with phonons, charge carriers, and EXs. Due to this scattering and radiative recombination, the EX phase is lost on a picosecond time scale, as shown by the results discussed in the previous sections. Nonetheless, the EX population in the dark states is evolving on a much longer time scale. Time-resolved PL performed on TMD SLs typically shows tails in a range of a few hundreds of picoseconds.[44−47] The quantitative interpretation of PL dynamics is difficult because there are many intermediate states in the scattering pathway from the initially excited electron–hole pairs to the emission of the bright EX states. Importantly, the occupation of dark EX states is not directly observable, but can be inferred indirectly, by modeling. Conversely, the experiments using phase-sensitive heterodyne detection[29,38] are sensitive not only to the bright EXs but also to the dark ones through their interaction with the optically probed bright EXs. Even more, the phase of the signal encodes the phase of the complex interaction energy, enabling the distinction between population of different dark-state reservoirs (fast direct spin-allowed, direct spin-forbidden, indirect spin-allowed, and indirect spin-forbidden).[38] Furthermore, resonant pumping generates a well-defined initial density of bright EXs (within the light cone) with a given spin-state encoded by the light helicity. Namely, the two pulses and (shifted by the radio frequencies Ω1 and Ω2, respectively), arriving in time overlap, create an EX density , oscillating at the frequency Ω2 – Ω1 = 1 MHz (see the Methods section). The FWM signal probing the density dynamics, generated by the third pulse delayed by τ23, infers the amplitude and the phase of the resulting modulation of the excitonic response. This pulse sequence is shown in Figure together with a typical dynamics measured at T = 5 K.

Figure 8

Resonantly excited EX population dynamics in a hBN/MoS2/h-BN heterostructure measured at 5 K. Co-circular polarization of , τ12 = 0.1 ps. The excitation power for each beam is 0.3 μW. The amplitude and phase retrieved by spectral interferometry are given by blue circles and green squares, respectively, along with the simultaneous fit according to the complex trial function (see discussion in the main text). Note the complex character of the measured FWM observable, inducing interference in the amplitude and phase-shifts revealing distinct densities scattering processes, while varying τ23. Top: pulse sequence employed to measure density dynamics via FWM. Data acquired on the sample A.

Qualitatively, the data for different densities and temperatures, as presented for FWM amplitudes in Figure S5, are consistently described by:

a weak signal for negative times due to previous pulses (repetition period is 13 ns), about 1 order of magnitude below the signal at τ23 = 0;

a signal amplitude rise from the negative delays τ23, given by the pulse autocorrelation, accompanied by a −π/2 phase shift;

a decay of the amplitude to less than half its value occurring for 0<τ23 < 1 ps;

subsequently, the signal amplitude rises on a time scale of 10 ps, accompanied by a phase shift of about π/2; and

for τ23 > 10 ps the amplitude decays, following a power law up to a delay of about 1 ns, while the phase shifts by about – π/2.

While the initial FWM decrease is attributed to the EX radiative recombination with the time constant given by τ1 and simultaneous scattering to the dark states, its subsequent rise (with additional features developing when increasing the temperature; see Figure S5) is surprising. To interpret the density dynamics measured with the heterodyne FWM, we recall that the technique retrieves response functions that are complex and, hence, sensitive to both amplitudes and the relative phases of the components of the signal R. Therefore the signals stemming from different densities and different interaction processes (between different EX states) interfere, which is visible when their relative contributions change along the delay τ23. In the present case, this effect is pronounced, indicating that the dynamics contains EX densities with interaction energies of significantly different phases, so that a description using a constant phase fails. In particular, to describe our data set, we introduce the following complex response function:

The meaning of the parameters is explained in Table . The response function includes a complex offset (Aof, ϕof), a two-photon absorption process (Anr, ϕnr) generating non-resonant FWM around τ23 = 0 and exponential decay processes, all multiplied by the power law decay. Note that (A, ϕ) pairs represent amplitudes and phases for the density scattering processes exhibiting exponential decays with the time constants τ. We found that a minimum of three such processes were required for a good fit to the data at 5 K, i.e., n ∈ (1, 2, 3). With increasing the temperature to 40 K, inclusion of the fourth component was necessary to model the data (not shown).

Table 1

Set of Parameters for the Fit Shown in Figure Employing the Response Function (eq )a

The values marked in blue correspond to the most relevant physical parameters of the fit.

Additionally, the third line of the above expression has been explicitly added to describe the power law visible for longer delays τ23. Introducing this term is motivated by the observed power law in the decay for different temperatures, and the fact that after some time, the EX density will be thermalized and therefore described by an overall decay process, which here appears to be EX–EX scattering leading to a power-law decay. The fitting function accounts for the temporal width of the excitation Gaussian pulses (τ0) and their repetition period (T) yielding this quite involved but analytical expression.

The simultaneously fitted dynamics to the FWM amplitude and FWM phase is shown by blue and green traces in Figure , respectively. We interpret it as follows. The initial drop is governed by the radiative recombination of EXs in the light cone competing with scattering out of the light cone leading to a fast decay with τ1 = T1 = (0.13  ± 0.04) ps. Later dynamics are characterized by τ2 = (4.2  ±  1.4) ps and comparable amplitude (A1 ≈ 0.62 × A2). Populating the indirect dark EXs results in the rise of FWM amplitude and produces a phase-shift due to modified interaction with the bright EXs. This is followed by the overall density decay via EX–EX scattering into a photon and a bright EX or nonradiative Auger recombination. To rephrase, in this model, the FWM rise (quantified by the second process with the parameters τ2, A2, and ϕ2) is due to the scattering of EXs into states with a stronger interaction with the probed bright EXs (for example, spin-forbidden direct or indirect EXs). We speculate that, within the TMD family, such scattering from bright toward dark EX configuration is the most efficient in MoS2 SLs, owing to its particularly small conduction band splitting,[48] giving a reason why it is much more pronounced here with respect to MoSe2 SLs.[38] After 10 ps, this process is completed and the subsequent decay dynamics is ruled by the power law with the power α = 0.59 ±  0.05. This value corresponds to a decay rate proportional to the EX density with the power 1 + 1/α ≈ 2.7 ≈ 3 and thus indicates triexciton scattering as dominating decay mechanism. The remaining third exponential process yields the decay constant τ3 = (56  ± 18) ps and A3 ≈ 1.4 × A1, comparable to A1. We found that, while restricting the modeling to the first two decay processes one can describe the FWM amplitude dynamics reasonably well (not shown), the third decay process was necessary to fit the phase dynamics. This exemplifies how the phase observable contains additional information, otherwise obscured when considering the FWM amplitude only. This third process is tentatively attributed to a further redistribution between the different dark exciton states, as suggested by the different phase of the resulting interaction.

Conclusions

We have shown that encapsulating MoS2 SLs in between h-BN layers drastically improves the optical quality in this material and permits us to recover giant coherent nonlinear responses of the EXs, as expected from their oscillator strength. Using two heterostructures differing in the degree of the EX disorder, we have performed three-beam FWM to infer the EX coherent and incoherent dynamics, spanning time scales from 100 fs to 1.3 ns. By measuring the coherence dynamics and time-resolved FWM, we reveal the formation of the photon echos, extracting homogeneous γ and inhomogeneous σ contributions to the EX spectral line-shape. We directly show the correlation between σ and the measured initial population loss, indicating that the latter is affected by the change in radiative lifetime on top of non-radiative channels. Importantly, at some positions of the low-disorder sample, the FWM transient appears in the form of interaction dominated free-induction decay, with no evidence of photon echo formation, showing that the EX broadening in the probed region is virtually homogeneous. FWM microscopy allowed us to reveal the impact of the local disorder on the EX’s oscillator strength and line-shape. The latter is shown to be affected by the temperature and also weakly by the excitation induced dephasing due to EX–EX interaction. The EX density dynamics measured by FWM is sensitive to both the dark and the bright EX density and can recover the complex interaction energy with the bright EXs using the amplitude and phase of the signal. Employing the modeling with complex fitting function, we identified three major EX relaxation channels, specifically: radiative decay and scattering out of the light cone (into fast EXs), scattering into dark EXs, which subsequently decay via tri-EX scattering. Further FWM transient grating investigations, thoroughly addressing a large space of parameters (the temperature, injected exciton density, valley-polarization, distinct TMD SL materials, exhibiting optically bright and dark exciton ground states, as well as the charge state, neutral EXs versus trions) are necessary to draw a comprehensive picture of the exciton relaxation dynamics. This would enable us to elaborate more accurate and predictive models, yet also more involved, with respect to the one presented in this work. To consistently describe the measured rich EX relaxation dynamics, one could for example incorporate master equations, modeling the dynamics between the different EX reservoirs.

Our methodology is well-suited to accurately assessing the impact of disorder on the EX responses in experiments with forthcoming TMD heterostructures. In particular, with the suppressed structural disorder, we could access and measure the coherent dynamics, mutual couplings and lifetime of the EX excited states. Finally, using spatially resolved FWM configuration would enable us to demonstrate long-range propagation of the coherence and of the exciton–polariton diffusion. The latter aspects are of utmost importance to reveal the exciton dispersion curve in optically active van der Waals structures.[49]

Methods

Preparation of the MoS2 van der Waals Heterostructure

In the present two samples, A and B, the natural MoS2 crystals were purchased from SPI, and the h-BN crystals were obtained from NIMS, Japan. The Si/SiO2 substrates with a 295 nm thick oxide were cleaned using acetone and isopropyl alcohol followed by nitrogen blow-drying. For the sample A (with more disorder), we used a viscoelastic stamping method for the stack preparation, while for the sample B (with less disorder), we employed a pick-up technique.

For the sample A, the h-BN crystal was placed on a scotch tape and was mechanically exfoliated onto a Si/SiO2 substrate, whereas the MoS2 crystal was exfoliated onto a PDMS layer. It has been already shown that large area monolayer MoS2 flakes could be obtained by exfoliating on PDMS.[50] A large SL MoS2 flake (15 × 80 μm2) was identified on PDMS based on optical contrast and was aligned and transferred onto the h-BN flake on Si/SiO2 by viscoelastic stamping, which is a dry method.[51] Likewise, another h-BN flake exfoliated on PDMS was stamped over the monolayer MoS2, so as to encapsulate it and to obtain a h-BN/MoS2/h-BN van der Waals heterostructure. While this process is known to yield a clean MoS2/h-BN interface at the bottom, the top interface with h-BN may encapsulate air blisters, puddles, and cracks. These may form as a result of the exerted mechanical stress, but a low amount of contaminants is expected.

For the less-disordered sample B, we use the so-called polypropylene carbonate (PPC) technique. The top and bottom h-BN and also the MoS2 are exfoliated onto Si/SiO2 substrates, respectively. Monolayer MoS2 and 15–20 nm thick h-BN flakes are identified using an optical microscope. Using a micron-thin film of PPC polymer top h-BN, monolayer MoS2 and the bottom h-BN are picked up one after another to form h-BN/MoS2/h-BN heterostructure using the procedure introduced in ref (55).

Four-Wave Mixing Microscopy

To measure coherent, resonant responses, particularly FWM spectra, we use heterodyne spectral interferometry.[52,53] We employ three laser pulses generated by the optical parametric oscillator (Radiantis Inspire) pumped by a Ti:sapphire laser (Spectra-Physics, Tsunami Femto). The pulses are chirp-corrected using a geometrical pulse-shaper, such that they arrive at the sample close to their Fourier limit, with around 150 fs duration. The beams are labeled  and are resonant with the EX transition, as displayed in Figure c (black trace). They are focused on the sample with the microscope objective (NA = 0.6) down to the diffraction limit. are frequency upshifted by distinct radio-frequencies Ω1,2,3 around 80 MHz using acousto-optic deflectors. The reflectance is collected by the same objective and spectrally dispersed by an imaging spectrometer. By exploiting optical heterodyning, we select the field component proportional to the FWM polarization proportional to , which is carried by the heterodyne beating at the Ω3 + Ω2 – Ω1 frequency, also occurring around 80 MHz. By applying the acousto-optic downshift for this mixing frequency at the detection path, the FWM spectral interference detected by a CCD camera is observed as a non-oscillating (DC) signal. Using a defined time-ordering between signal and reference pulses, the signal is resolved in amplitude and phase using spectral interferometry,[54] providing the signal in both spectral and temporal domains connected via Fourier transform. By measuring the FWM as a function of time delays τ12 (between and ) and τ23 (between and ), we gain information about the EX coherence and density dynamics, respectively.