Imaging dynamic exciton interactions and coupling in transition metal dichalcogenides

Transition metal dichalcogenides (TMDs) are regarded as a prime material platform for applications ranging from solar-energy^1,2 and lasers³ to quantum light-emitting diodes.⁴ Specifically, the rapid charge transfer and associated interlayer excitons^5–12 have received considerable attention. Interlayer excitons with nanosecond lifetimes^13–15 that are highly tunable by the twist-angle between layers¹⁰ are potential candidates for qubits.^16,17 In recent years, strain-engineering in TMD monolayers has also gained momentum with the potential for room-temperature entangled-photon sources.^18,19 Moreover, coherent coupling between excitons and trions in TMD monolayers^20,21 and intralayer excitons in TMD heterostructures⁵ has been demonstrated recently, which opens the avenue for quantum coherent control of these materials.

Depending on the sample system, different physical parameters determine the feasibility of quantum information applications. In TMD monolayers, low inhomogeneity and long intralayer exciton dephasing times are crucial. Martin et al. demonstrated picosecond dephasing times for MoSe₂ monolayers²² by measuring the homogeneous linewidth, which is inversely proportional to the dephasing time. These dephasing times are well above the previously reported values^23–26 due to the dominant contribution of inhomogeneous broadening in these samples. Meanwhile, Jakubczyk et al.^27,28 and Boule et al.²⁹ showed broader homogeneous linewidths, with the exciton optical response being in the homogeneous limit for certain areas. As established previously,^27,28 there is an expected anti-correlation between homogeneous and inhomogeneous linewidths, favoring long dephasing times in areas of large inhomogeneity. Non-radiative broadening mechanisms can weaken this anti-correlation. However, a large inhomogeneity is strongly detrimental to quantum information applications since it corresponds to several different emitters within the excitation volume. To identify promising areas for these applications, one thus needs to find uniquely suited sample areas with small inhomogeneity and large dephasing times.

For heterostructure systems, these considerations change based on the physical effects harnessed for quantum information applications. Interlayer excitons (ILEs) have been proposed as potential candidates for qubits in the past with their formation closely tied to the rapid charge transfer. Moreover, coherent coupling between excitons in the different layers of the heterostructures is another pathway toward coherent control for quantum information applications. For feasibility, identifying areas with the robustness of properties such as charge transfer, ILE lifetimes, and coherent coupling strength is thus a requirement.

However, the intricate spatial and temporal variations of excitons interacting with their environment and among themselves in TMD monolayers and heterostructures have remained mostly elusive due to limitations in the employed experimental techniques. Here, we use multi-dimensional coherent imaging spectroscopy (MDCIS) to map the distribution of dephasing times and inhomogeneity across a MoSe₂ monolayer, enabling the identification of promising and unfavorable areas for quantum information applications. MDCIS also gives insight into strain across the monolayer. Correlating it with findings from photoluminescence (PL) spectroscopy and PL-detected MDCIS, we visualize the distribution of bright and dark exciton states across the monolayer sample. We further use MDCIS to spatially map strain and conduction band changes across a MoSe₂/WSe₂ heterostructure. Despite the presence of notable, complex local strain and dielectric environment changes, we demonstrate surprising robustness of key sample properties: The rapid electron and hole transfer and coherent coupling between intralayer excitons in the MoSe₂ and WSe₂ monolayers and interlayer exciton lifetimes are robust across a majority of the heterostructure. These results have larger implications for the commonly made device application claims, strengthening the case for TMDs as the materials of choice in various applications, including quantum information.

The multi-dimensional coherent spectroscopy (MDCS) technique that is the basis for the majority of results in this paper is schematically shown in Fig. 1(a). MDCS uses a four-pulse sequence where the first two pulses (A and B) act as the pump pulse, while the third pulse (C) acts as a probe pulse and a fourth pulse (D) is used for heterodyne detection. The first pulse (A) excites a coherence between the ground state and excited state, while the second pulse (B) converts this into an excited or ground state population. The second pulse can also convert the system into a Raman-like non-radiative coherence between different excited states of the excitons, which oscillates during T. The third pulse (C) reverts the population or non-radiative coherence back into a radiative coherence between ground and excited states. Because MDCS measures the phase-resolved response, scanning the time delay τ between the first two pulses gives access to the absorption energy upon Fourier transform, while scanning the time delay t between C and D gives access to the emission energy axis. Additional dynamics, such as population decay, charge transfer, and coherent coupling, are accessible via the T time delay.

To be compatible with diffraction-limited imaging, all four pulses need to be collinear, that is, all excitations and the signal must be overlapped in the same beam. To extract the signal against a large background, each pulse is frequency tagged with an acousto-optic modulator^5,22,30 and the interference between the third-order nonlinear response (induced by the first three pulses) and the fourth pulse can subsequently be phase-sensitively detected using a lock-in amplifier. The absorption of beam D can be taken into account by measuring the reflected spectrum using the interference between C and D.

In previous implementations of imaging MDCS, spatially addressed MDCS, i.e., taking a full MDCS spectrum at specific sample points, has been used to gain a spatial understanding of the spectroscopic properties of the sample.^{27–29,31,32} In our implementation of MDCIS, instead of taking a single data point for fixed τ, T, and t positions, an entire image of the sample is taken. This yields an up to five-dimensional dataset with full temporal information along τ, T, and t for every image pixel. Subsequent Fourier transforms along the τ and t axes lead to an MDCS spectrum at every single image pixel. This process is illustrated for a hexagonal Boron Nitride (hBN) encapsulated MoSe₂ monolayer [Fig. 1(b)] in Fig. 1(c). MDCIS enables the study of linewidth evolutions across the sample and coupling dynamics by harnessing the rich spectroscopic and temporal information of MDCS and the rich spatial information of nonlinear imaging. The imaging is implemented with a laser-scanning imaging setup shown in Fig. 1(d), allowing for a diffraction-limited beam diameter (Abbe limit) of 940 nm. Because of the third-order nonlinearity of the four-wave mixing (FWM) based MDCS, the spatial resolution for the FWM measurements reaches 540 nm. The four pulses impinge on an x- and y-Galvo mirror, with the x-Galvo mirror being relayed onto the y-Galvo mirror via a 4f-setup incorporating off-axis parabolic mirrors. The deflected beams are subsequently sent through another 4f-setup containing a scan and tube lens optimized for broadband, wide field-of-view laser scanning microscopy. This arrangement images the angle of the laser beam onto the objective without changing the position on the objective aperture, allowing for aberration-free laser scanning imaging. We employ a custom lock-in amplifier³³ for rapid imaging. A Ti:sapphire laser with nearly transform-limited 35 fs pulses is used in this work. The pulses are pre-compensated with negative dispersion using an SLM-based pulse shaper to compensate for the dispersion acquired upon propagation through the setup. All experiments are performed at 6 K.

We further demonstrate the usefulness of the multi-dimensional coherent imaging spectroscopy technique, which, on a smaller scale, has been realized by spatially addressed multi-dimensional coherent spectroscopy.^{27–29,31,32}

The two samples studied in this work are an hBN encapsulated MoSe₂ monolayer and an hBN encapsulated MoSe₂/WSe₂ heterostructure. The samples were assembled using a dry-transfer technique with a stamp made of a polydimethylsiloxane cylinder with a thin film of poly(bisphenol A carbonate) on top.

The white light microscopy image of the MoSe₂ monolayer in Fig. 1(b) shows a region with cracks and bubbles toward the center but otherwise appears mostly pristine. The upper-left part of the monolayer region (shaded) is non-encapsulated, while the lower-right side of the monolayer (unshaded) is encapsulated. Given the reported instabilities and inconsistencies of non-encapsulated samples,^22,28,32 we will focus our studies on the encapsulated area of the sample. A spectrally integrated four-wave mixing (FWM) image of the MoSe₂ monolayer is shown in Fig. 1(e). For the monolayer, unless otherwise stated, data were taken at T = 50 fs. We spectrally integrate from 1630 to 1675 meV for both the absorption and emission axes. The integrated FWM will scale with the dipole moment and density of states (∝Nμ⁴), with N being the density of states and μ the dipole moment. The large area on the top shows a very strong and approximately homogeneous FWM amplitude, while the dark area toward the center of the samples coincides with cracks and bubbles visible in the white light microscopy image in Fig. 1(b). The bottom area of the sample again shows a comparatively homogeneous strength. This observation is in contrast to photoluminescence (PL)-spectroscopy measurements plotted in Fig. 1(f). Here, the integrated PL strength for the central area (yellow square) is comparable to the homogeneous area of the sample (green square). PL and PL-detected FWM (shown in the supplementary material) is only sensitive to bright, photoluminescent states, while heterodyne-detected FWM is sensitive to both bright and dark states. As discussed by Smallwood et al.,³⁴ the difference between PL strength, PL-detected FWM strength, and heterodyne-detected FWM strength suggests a highly spatially dependent distribution of bright and dark states. Furthermore, a noticeable linewidth change between the central and edge area of the monolayer can be observed in the PL spectra. This linewidth observation is corroborated by the exemplary MDCS spectra for the homogeneous bottom area of the sample [Fig. 1(g)] and the central, dark area [Fig. 1(h)]. The central area shows a more dominant inhomogeneous broadening, which manifests itself as an elongation along the diagonal (dashed line), compared to the bottom sample area. The broader linewidth in the PL spectrum in Fig. 1(f) is caused by an increase in inhomogeneity and not related to homogeneous linewidth changes. Moreover, the inhomogeneously broadened spectrum on the central area shows more weight toward lower energies, suggesting the existence of notable resonance shifts on the order of several meV.

To gain a better understanding of the underlying physics and identify promising areas for quantum information applications, a more systematic study of resonance shifts and linewidths is required. This relates to the spatial coherence of excitons, which is, in conjunction with the temporal coherence, crucial for quantum information applications and has been of interest for decades.³⁵ Spatial inhomogeneity, corresponding to local exciton energy changes due to a change in the environment, is associated with a lower spatial coherence because the energy shifts destroy the coherence. Hence, in the presence of a smaller inhomogeneity σ, the exciton coherence volume in real space is larger. A larger real space volume corresponds to a smaller k-space volume. This increases the overlap of the exciton dispersion with the light-cone, increasing the strength of the light matter interaction (and hence dipole moment μ) and decreasing the radiative lifetime T_rad. A decreased radiative lifetime decreases the overall dephasing time T₂ = T_rad + T_non-rad and via T₂ = ℏ/γ yields an increased homogeneous linewidth. By this argument, one expects larger γ for smaller σ and hence an anti-correlation between the two linewidths, as pointed out by Jakubczyk^27,28 and Boule et al.²⁹ In the idealized scenario of a radiatively limited homogeneous linewidth and no changes in the intrinsic linewidth (i.e., only changes in the radiative linewidth contribute to homogeneous linewidth changes), a perfect anti-correlation is therefore expected. However, this anti-correlation might not be perfectly linear but see a strictly monotonous decrease in γ for increasing σ. Because the Pearson correlation coefficient only captures the linear correlation, the coefficient would be close to but above −1. Furthermore, in the presence of non-radiative homogeneous broadening mechanisms, such as non-radiative electron scattering, γ can be further increased,^27,28 weakening the anti-correlation. This is reflected in the literature, where the Jakubczyk et al.²⁸ measured anti-correlations of ρ = −0.06 and ρ = −0.425 for two different MoS₂ monolayer samples.

For quantum information applications, small γ and σ values are desired. The smallest γ and σ values are hard to achieve for a given material platform because of the anti-correlation, which can only be weakened by an increase in γ due to non-radiative contributions. However, based on the chosen material platform, overall low γ and σ values can still be achieved.

Maps of γ and σ across the MoSe₂ monolayer are shown in Figs. 2(a) and 2(b), respectively. The linewidths are extracted by simultaneously fitting diagonal and cross-diagonal slices in the 2D frequency data with the analytical solutions for arbitrary amounts of γ and σ provided in Ref. 36. This is notably different from the procedure employed in Refs. 27–29, where linewidths are fitted in the temporal domain. As discussed in the supplementary material, the frequency domain is less susceptible to measurement noise because of the Fourier transform filtering out a majority of the high frequency noise. Furthermore, the frequency domain fitting procedure emphasizes the high signal areas, as elucidated upon in the supplementary material. After fitting, we further average 2 pixels × 2 pixels together to reduce pixel-to-pixel noise. The averaging leaves us with a pixel size of 700 nm, close to the spatial resolution.

The areas of lower FWM show a significantly increased σ value and a smaller γ value. Overall, γ shows drastic changes, ranging between 0.75 and 1.7 meV, more than a factor of two difference. Depending on the sample spot, σ increases by up to a factor of 10, from below 3 meV to values up to 35 meV. As pointed out in Refs. 27–29 and 37, these changes could be caused by strain, wrinkling, flake deformations, cracks, changes in the dielectric surrounding, changes in doping/trapped charges/impurities/defects, and the free carrier concentration. In this paper, we use the term “strain” to encompass all deformations of the material due to mechanical stress, including strain, wrinkling, and flake deformations. As discussed in Refs. 22 and 38, the variation in dielectric permittivity and doping is suppressed by hBN encapsulation. Furthermore, the exciton–trion ratio does not change significantly between sample points A and B from Fig. 1(b), shown in the supplementary material, suggesting no significant variation in the doping across the sample. We also see no signature of the trion in the FWM. Similar to Ref. 27, we thus exclude those changes as the main cause for the variations in the optical measurements. Similar to Refs. 27 and 28, we assign linewidth and dipole changes to strain. Strain is associated with larger exciton localization and hence larger σ, smaller γ, and smaller μ values.

The anti-correlation between γ and σ across large areas of the sample can best be visualized by plotting σ against γ, as done in Fig. 2(c). Here, sample points falling within the green and yellow rectangles [Fig. 2(a)] are plotted in green and yellow, respectively, while all other sample points are plotted in black. In the following, we use the Pearson correlation coefficient that measures the linear correlation between two sets of data A and B and is defined as ρ(A, B) = cov(A,B)/(σ_Aσ_B). Here, cov denotes the covariance of A and B and σ_A,B is their respective standard deviation. A moderately strong anti-correlation of ρ = −0.47 ± 0.1 can be observed for the two linewidths across the entirety of the sample. A closer inspection of Fig. 2(c) together with Figs. 2(a) and 2(b) suggests that there are areas of stronger anti-correlation, while other areas show a weaker correlation between the linewidths. This is further supported by plotting the two areas marked with green and yellow rectangles in Fig. 2(c). Indeed, the green sample points show a strong anti-correlation with a correlation coefficient ρ = −0.74 ± 0.13, while the yellow points show a much weaker anti-correlation with a correlation coefficient ρ = −0.33 ± 0.26. Most importantly with regard to quantum information applications, the yellow points show the values of γ around 0.9 meV with comparatively low σ values between 4.5 and 7 meV, while for the green rectangle, similar values of γ can only be found with σ above 10 meV. Hence, the yellow rectangle area toward the bottom of the sample is a more favorable area than the green rectangle for quantum information applications that rely on low inhomogeneity and slow dephasing. However, another low strain area toward the top of the sample, while showing overall lower σ, also shows higher γ, illustrating that a low strain area is not inherently a “good” area of the sample with regard to the application potential. “Good” areas, with respect to quantum information applications, are areas that (1) show linewidths in the homogeneous limit and (2) have small homogeneous linewidths (i.e., long dephasing times). With additional non-radiative contributions to γ, such as non-radiative electron scattering,^27,28 only a technique like MDCIS allows us to unambiguously identify the relevant sample areas.

Assigning the linewidth and dipole changes to strain, we can further distinguish between tensile and compressive strains via the sign of resonance energy shifts across the sample.³⁸ We plot the exciton resonance energy across the sample in Fig. 2(d). This map shows significant shifts of the resonance energy between 1645 meV and up to 1660 meV across the 10 × 15 μm² area of the sample. The areas of the strongest resonance shifts also coincide with the areas of lower integrated FWM. The observation of both red- and blueshifts across the sample points toward complex local strain dynamics involving both tensile (lower energies) and compressive (higher energies) strains, caused by the cracks and bubbles in the central area of the sample.

Having established overall strain dynamics and how intralayer excitons are affected in MDCIS across an encapsulated MoSe₂ monolayer, we can take these findings to a more complicated material system, namely, a MoSe₂/WSe₂ heterostructure.

We extend our study to a MoSe₂/WSe₂ heterostructure encapsulated in hBN. Figure 3(a) shows a white light microscopy image of the heterostructure. Three exemplary MDCS spectra taken at the three points marked in Fig. 3(a) are plotted in Fig. 3(b). For the MoSe₂/WSe₂ heterostructure, unless otherwise stated, data were taken at T = 500 fs. The two on-diagonal (dashed line) peaks are associated with the MoSe₂ and WSe₂ intralayer A-excitons. The two off-diagonal peaks are indicative of both coherent coupling and incoherent electron and hole (charge) transfer, as discussed in previous work.⁵ The spectra show significant energy shifts and varying peak strengths for both the MoSe₂ and WSe₂ resonances across the sample. We spectrally integrate over the four peaks to better visualize their strength variations across the sample. Maps of the integrated FWM for the four peaks are shown in Fig. 3(c). These figures have the same order as the peaks, with the MoSe₂ (M) peak in the lower left, the WSe₂ (W) peak in the upper right, and the MoSe₂/WSe₂ (M/W) and WSe₂/MoSe₂ (W/M) peak in the lower right and upper left, respectively. Similar to the MoSe₂ monolayer, there is a region of decreased FWM toward the upper center of the sample, which is associated with a high strain area due to wrinkles and bubbles that are caused by fabrication. Apart from this area, the peak strength is mostly homogeneous across the sample for all four peaks, except for the MoSe₂ and MoSe₂/WSe₂ peaks having a lower strength toward the bottom left of the sample. The differences in peak strengths can partially be attributed to the spatial variations of the dipole moment μ that differ for the two materials based on the local strain profile. Furthermore, the finite bandwidth of the employed laser and the reduced excitation density at sample points where the WSe₂ (MoSe₂) resonance is shifted to higher (lower) energies can contribute to the spatial peak strength variations.

The resonance energy shifts for the two resonances are plotted in Figs. 3(d) and 3(e). A large shift toward lower energies down to 1605 meV for the MoSe₂ resonance can be observed in the center-left of the sample, while the WSe₂ resonance shifts toward higher energies, up to 1680 meV in this area. Toward the bottom area of the sample, both resonances shift toward higher energies. Moreover, the MoSe₂ resonance shifts toward significantly higher energies up to 1620 meV in the upper part of the sample. To correlate the resonance shifts, we plot the MoSe₂ resonance vs the WSe₂ resonance energy in Fig. 3(f). Sample points falling within the red (blue) rectangle are plotted in red (blue), and all other data points are plotted in black. While changes in the resonance energies seem to be correlated, a Pearson coefficient of ρ = 0.07 ± 0.13 suggests otherwise. A closer examination of Fig. 3(f) together with Figs. 3(d) and 3(e) shows both strong correlation and anti-correlation, which cannot be captured using the Pearson correlation coefficient. Instead, to visualize the strong variation of correlation across the sample, we examine the select areas marked by blue and red dashed rectangles in Fig. 3(d). For the area marked by the blue rectangle (Reg D), a strong correlation with ρ = 0.54 ± 0.12 can be observed. A correlation between the two resonances is what would be expected: It is well known that both encapsulation and heterostructure formation commonly red shift the excitons,^9,22 while compressive (tensile) strain leads to red (blue) shifts.³⁷ Both monolayers should be affected in the same way by these factors. Contrary to this expectation, the area marked by the red rectangle (Reg E) shows a moderately strong anti-correlation with ρ = −0.43 ± 0.22. This points toward complex local strain dynamics, where the two monolayers experience different strains, which leads to opposite resonance shifts. One possible scenario is the compression of one monolayer, which leads to a bubble in the other monolayer, inducing tensile strain and thus the opposite sign resonance shift. The attribution to strain is further supported by the fact that the largest resonance shifts can again be observed around an area of low FWM with physical imperfections (cracks, bubbles, and wrinkles) visible in the white light microscopy image [see Fig. 3(a)].

The local strain profile also leads to large variations in the ILE PL. The integrated ILE PL plotted in Fig. 4(a) shows a strong maximum in the high-strain area of minimum FWM but appears otherwise relatively homogeneous in strength across the sample. These observations align with the observations on the MoSe₂ monolayer. The emission energy of the ILE is plotted in Fig. 4(b). The area of strong PL shows the above average energy of ILE emission around 1370 meV. The upper area of the sample also shows a higher ILE emission energy around 1355 meV, while the lower part of the sample shows lower emission energies down to 1330 meV. The strong changes in the ILE PL, including emission strength and energy, even outside the high strain area, can be further emphasized when plotting two select spectra in what appears to be a relatively homogeneous region in the bottom half of the sample. These spectra plotted in Fig. 4(c) demonstrate that even in the bottom part of the sample, ILE emission is highly spatially heterogeneous.

We further observe a strong correlation between the ILE PL-emission energy and the MoSe₂ resonance energy (measured with MDCS), as evident from Fig. 4(d) and a correlation coefficient ρ = 0.47. The correlation with the WSe₂ resonance [Fig. 4(e)] is much weaker (ρ = −0.27). This difference can be explained as follows: As Khatibi et al.³⁷ showed for TMD monolayers, along the K-point, mainly the conduction band is affected by strain, while effects on the valence band are negligible. Thus, as illustrated in Fig. 4(f), the strain-induced conduction band shift in MoSe₂ immediately affects the ILE emission energy, while the strain-induced shift in WSe₂ should have no effect. The residual anti-correlation between ILE emission energy and WSe₂ resonance energy stems from the local (anti-)correlation between MoSe₂ and WSe₂ established in the discussion above. Given these significant strain-induced changes to the heterostructure, if and how much these changes affect the coupling between excitons in the two layers is a topic of interest.

To shed some light on the coupling dynamics, we employ dynamic MDCIS, where the pump-probe delay T [see Fig. 1(a)] is varied to access coherent and incoherent coupling dynamics in the heterostructure. We plot the integrated FWM of the MoSe₂/WSe₂ (M/W) and WSe₂/MoSe₂ (W/M) coupling peaks for varying T delays in Figs. 5(a) and 5(b), respectively. Here, we have normalized both peaks by $M⋅W$⁠, with M and W being the integrated FWM of the MoSe₂ and WSe₂ peaks at T = 20 fs, respectively. This factors out the effect of local dipole moment changes; a more detailed discussion can be found in the supplementary material. All images for each peak are subsequently normalized by a single value so that the brightest pixel at T = 0 fs has a unity amplitude. Initially, the signal decays for both coupling peaks from 0 to 40 fs but recovers again for 60 fs. A smaller variation in the integrated amplitudes can be observed between 60 and 100 fs. This behavior can be assigned to coherent coupling oscillations that have been shown to occur in MoSe₂/WSe₂ heterostructures.⁵ The coherent oscillation has an amplitude of ∼50% (with respect to the peak amplitude) across the sample, the signature of strong coupling between excitons in the two layers. As discussed by Hao et al.,²¹ the amplitude is below unity due to interference between different coherent coupling contributions and exponentially decaying phase-space filling nonlinearities. In between T = 100 fs and T = 500 fs, both coupling peaks show a clear rise that we assign to the charge transfer that has been observed using MDCS in this specific heterostructure⁵ and other TMD heterostructures.⁷ While not resolved here, the charge transfer for the WSe₂/MoSe₂ peak is faster in this sample, as indicated by both previous measurements⁵ and the broader linewidths of the WSe₂ peak in Fig. 3(b). The strong amplitude change between 100 and 500 fs suggests highly efficient electron and hole transfer in the heterostructure. Most strikingly though, both coherent coupling and charge transfer appear mostly homogeneous across the sample. This is evident from the fact that the relative peak strength across the sample remains mostly unchanged along T, while the absolute peak strength changes. The homogeneity can be visualized by comparing the relative strength profiles at T = 0 fs, T = 40 fs, T = 100 fs, and T = 500 fs, as done in Fig. 6(a). Here, we plot the difference between the integrated FWM maps for the M/W peak, normalized by the M/W map at T = 0 fs. We plot the difference between T = 0 fs and T = 40 fs (leftmost map), T = 40 fs and T = 100 fs (center map), and T = 100 fs and T = 500 fs (right map). Both the initial decay and recovery, part of the coherent coupling oscillations, show a homogeneous strength profile across the entire sample. Moreover, the charge transfer rise between 100 and 500 fs shows the same homogeneous strength profile. To quantify the homogeneity of coherent coupling and charge transfer, we look at the images of the dipole-normalized coupling peaks, plotted in Fig. 5(a). With the normalization, the main physical effects contributing to a spatially inhomogeneous coupling peak map are variations in the coherent coupling strength or incoherent transfer processes. While showing some spatial inhomogeneity, 81% of the sample fall within ±20% of the average amplitude for T = 40 fs. This T-step was chosen because finite pulse effects happen at earlier T, while incoherent charge transfer dominates for later T, leading to a dominant coherent coupling contribution. Similarly, for T = 500 fs, 91% of the sample fall within ±20% of the average amplitude. This suggests that while some inhomogeneity of the coherent coupling and charge transfer remains, both are surprisingly robust toward the strain in this sample. An extensive discussion of the scheme to determine homogeneity can be found in the supplementary material. While the amplitude difference maps can be used to disentangle contributions of charge transfer from coherent coupling,⁵ noise makes quantitative results on the homogeneity inconclusive. These results are presented in the supplementary material.

Moderate resonance shifts due to strain on the order of 10–20 meV are not inherently expected to change the charge transfer dynamics significantly because the shift is comparatively small to the band offsets of hundreds of meVs.^39,40 However, the sensitivity of charge transfer to sample parameters, such as the twist angle^5,41 and lattice separation,⁴² is well documented. With the resonance shifts being an indicator of complex local strain dynamics, this strain is expected to change the interlayer spacing, among other things. These previous findings render our observation of spatially homogeneous charge transfer rather surprising. Moreover, the coherent coupling between the resonances is expected to be equally dependent on the intralayer separation. Indeed, a limited decrease in the coherent coupling strength is observable for the top part of the sample, as illustrated in Fig. 6(b). Here, we plot the integrated FWM of the MoSe₂/WSe₂ coupling peak along the three arrows drawn in Fig. 5(a). Both the green and red arrows, on the lower part of the sample, show a strong coherent coupling oscillation with only small, random changes when moving across the sample, consistent with the observations in Fig. 6(a). However, the data for the blue arrow show a weaker oscillation, which might be either caused by a reduced amplitude or a more rapid dephasing of the coherent coupling in the upper area of the sample. Given the strong resonance shifts for the MoSe₂resonance in this area and the reflectance map (see the supplementary material), we assume this to be an area of increased interlayer spacing. The increased interlayer spacing can explain the reduced coherent coupling strength/more rapid dephasing.

This specific area also shows up distinctly in the FWM decay time map plotted in Fig. 6(c). The decay time map is acquired at t = τ = 0 while moving the T delay [see Fig. 1(a)]. The FWM decay in the heterostructure for T > 50 ps is dominated by interlayer exciton decay through ground-state bleaching contributions to the signal.⁴³ The decay time can thus be taken as an indirect probe of the interlayer exciton lifetime. The bottom area of the sample again shows no notable spatial inhomogeneity in the lifetime, with values in the range between 500 and 550 ps. Given the low twist angle of this sample, these values are in good agreement with the literature.¹⁰ The upper area of the sample has on average a lower lifetime in the range of 30%–40% and larger inhomogeneity with values ranging from 300 to 550 ps. The combination of reduced interlayer exciton lifetime and reduced coherent coupling, together with the resonance shifts observed in this area, shows that the sample properties are not entirely immune to strain and defect induced changes. Nonetheless, the relatively low sensitivity toward these changes remains surprising when compared with the high inhomogeneity of the PL images.

We have examined an encapsulated MoSe₂ monolayer and an encapsulated MoSe₂/WSe₂ heterostructure using multi-dimensional coherent imaging spectroscopy. We image strain, dephasing, and inhomogeneity across the monolayer, showing a moderately strong anti-correlation between homogeneous and inhomogeneous linewidths that increases and decreases based on the area of the sample. Using MDCIS allows us to unambiguously identify promising and unfavorable areas for quantum information applications. We also visualize a spatially dependent dark vs bright state distribution using the difference between heterodyne detected FWM images, PL, and PL-detected four-wave mixing. Carrying this technique over to the heterostructure, we map the complex strain dynamics by correlating resonance energies. Our observations lead us to propose a phenomenological model for the band structure changes due to strain variations using additional interlayer exciton photoluminescence measurements. Employing dynamic multi-dimensional coherent imaging spectroscopy, we visualize coherent coupling and electron and hole transfer across the sample. While certain areas of the sample show a reduced coherent coupling and decreased interlayer exciton lifetime, we demonstrate an overall robustness of the coupling dynamics to strain and defects across the sample.

The future of these materials, especially in the realm of quantum information, is inherently coupled to the scalability and quality of fabricated devices. This work shows a reproducibility of crucial physical properties—dephasing time, coupling strength, and interlayer exciton lifetime—across large areas of the sample. The reproducibility lays the groundwork and strengthens the case for transition metal dichalcogenides as a next generation material. We further demonstrate the usefulness of the multi-dimensional coherent imaging spectroscopy technique, which, on a smaller scale, has been realized by spatially addressed multi-dimensional coherent spectroscopy.^27–29,31 However, recent technological advances in lock-in detection³³ allow for a larger scale realization of multi-dimensional coherent imaging spectroscopy and related techniques, advancing the forefront of materials and device characterization.

See the supplementary material for additional information on the experimental setup, the linewidth fitting procedure, reflection maps of the samples, and PL-detected FWM imaging data on the MoSe₂ monolayer.

We thank Blake Hipsley for the construction of the pulse shaping setup. The research at the University of Michigan was supported by NSF Grant No. 2016356. The work at the University of Washington was supported by the Department of Energy, Basic Energy Sciences, Materials Sciences and Engineering Division (Grant No. DE-SC0012509). W.G.H. was supported by the NSF Graduate Research Fellowship Program under Grant No. DGE-1762114.

E.W.M. and S.T.C. are co-founders of MONSTR Sense Technologies, LLC, which sells ultrafast spectrometers and microscopes.

The data that support the findings of this study are available from the corresponding author upon reasonable request.