Destructive Photon Echo Formation in Six‐Wave Mixing Signals of a MoSe2 Monolayer

Abstract

Monolayers of transition metal dichalcogenides display a strong excitonic optical response. Additionally encapsulating the monolayer with hexagonal boron nitride allows to reach the limit of a purely homogeneously broadened exciton system. On such a MoSe₂‐based system, ultrafast six‐wave mixing spectroscopy is performed and a novel destructive photon echo effect is found. This process manifests as a characteristic depression of the nonlinear signal dynamics when scanning the delay between the applied laser pulses. By theoretically describing the process within a local field model, an excellent agreement with the experiment is reached. An effective Bloch vector representation is developed and thereby it is demonstrated that the destructive photon echo stems from a destructive interference of successive repetitions of the heterodyning experiment.

Introduction

The spin echo^[1] is an essential effect in nuclear magnetic resonance spectroscopy and the basis for all sorts of complex radio pulse sequences^[2] that are routinely applied in medicine,^[3] chemistry,^[4] or physics.^[5–7] While spin resonances are driven by radio frequencies, optical frequencies are required to resonantly excite typical charge transitions; in this regime the analogous phenomenon is called photon echo.^[8] First demonstrations were performed on ruby crystals,^[8,9] but the photon echo spectroscopy also has a long standing history in semiconductor optics. It has been used to study different types of exciton dynamics, ranging from exciton-exciton scattering^[10–13] to exciton-phonon coupling,^[14,15] and it has been applied to 3D bulk,^[10,16,17] 2D quantum well,^[18,19] 1D nanowire,^[20] and 0D quantum dot structures.^[21–24] However, the technique is not restricted to solid state samples; it has also been applied to liquids.^[25–27] In its classical form, the photon echo is based on a four-wave mixing (FWM) process and therefore it constitutes a nonlinear process of third order (χ⁽³⁾) in the low excitation limit.^[28] Especially in single low-dimensional systems like quantum dots, it requires large effort to detect such nonlinear optical signals.^[29,30] Interestingly, due to their strong excitonic optical response, monolayers of transition metal dichalcogenides (TMDCs) show a remarkable signal strength in FWM spectroscopy.^[31–34] So far, the direct correspondence between inhomogeneous spectral broadening and photon echo duration has been used to map the inhomogeneity of TMDC monolayers.^{[31,32,34,35]}

We here exploit the strong excitonic optical response of TMDC monolayers, and explore six-wave mixing (SWM) signals from a MoSe₂ monolayer in the low excitation limit, which here represents the χ⁽⁵⁾-regime.^[36–39] We find that the signal dynamics exhibit a characteristic temporary suppression depending on the delay between the pulses. Supported by a theoretical model based on the local field effect describing exciton-exciton interaction, we explain that this suppression can be understood as the formation of a destructive photon echo.

Results Sample and Four-Wave Mixing Dynamics

We perform multi-wave mixing experiments on an hBN/MoSe₂/hBN heterostructure as schematically depicted in Figure 1. By radio frequency modulating the incoming beams, the different pulses are labeled by phases ϕ_j.^[40] After heterodyning the emitted light with a specific N-wave mixing (NWM) phase, that is a particular phase combination of the form ϕ_NWM = ∑_ja_jϕ_j ($a_j\in \mathbb{Z}$, ∑_j|a_j| = N − 1, ∑_ja_j = 1), we retrieve the corresponding nonlinear signal from the investigated monolayer.^[41] Details are given in the experimental section. We use the same sample as investigated in ref. [35] where the echo formation in two-pulse FWM signals with ϕ_FWM = 2ϕ₂ − ϕ₁ was used to study the inhomogeneity of the structure. The well known photon echo appears due to the dephasing impact of the structural inhomogeneity on the exciton's coherence as schematically depicted via Bloch vectors in Figure 2.^[42] For the sake of simplicity in the illustration we show a combination of a π/2 and a π pulse, which results in a photon echo in the full coherence. However, when considering the FWM coherence characterized by the phase ϕ_FWM the photon echo appears for any combination of pulse areas. The first laser pulse in Figure 2a having a pulse area of θ₁ = π/2 and arriving at the time t = −τ generates an exciton coherence. Because of the presence of different transition energies, originating from strain and dielectric variations, the coherences generated by the first laser pulse oscillate with different frequencies resulting in a dephasing of the total coherence. As depicted in Figure 2b, after a delay τ, that is, at the time t = 0, a second laser pulse with the pulse area θ₂ = π inverts all Bloch vectors. This is followed by a rephasing of the different coherences. The rephasing takes the same time as the dephasing and therefore the FWM signal is significantly enhanced due to constructive interferences at the time precisely given by the delay time t = τ between the two pulses.

Figure 1. Sample structure consisting of a multilayer hBN, monolayer MoSe2, multilayer hBN stack. (Left) Four-wave mixing (FWM) generated by two laser pulses with tunable delay τ. (Right) Six-wave mixing (SWM) with three excitation pulses having a tunable delay τ between pulses 2 and 3. The FWM and SWM dynamics are measured in real time t.

Figure 2. Schematic Bloch vector image of the photon echo process. a) 1st laser pulse excitation with following dissipation. b) 2nd excitation inverting all coherences and following rephasing. Laser pulse rotations are depicted in green, initial and final Bloch vectors in grey and blue, respectively, and the inhomogeneity induced coherence dynamics in red.

In Figure 3a we show the measured FWM signal from our sample as a function of the real time t after the second pulse and the delay τ. The pulse alignment is such that the pulse with phase ϕ₁ arrives at t = −τ and the pulse with ϕ₂ at t = 0. We find that the amplitude is not concentrated along the diagonal, which would represent a photon echo. Instead, starting from its maximum around (τ, t) = (0, 0) it basically decays to positive delays and times. This demonstrates that the investigated sample position is virtually homogenous because the photon echo is not present. Figure 3b shows the corresponding simulated FWM amplitude |S_FWM| within a few-level system where we take the exciton transfer into the optically uncoupled valley as well as the exciton-exciton interaction in terms of the local field model into account. It describes the optically induced dynamics of the exciton's coherence p and its occupations n and $n^{\prime }$ via^[11,43,44] [Image Omitted. See PDF][Image Omitted. See PDF][Image Omitted. See PDF]

Figure 3. Four-wave mixing dynamics as a function of real time t and delay τ. a) Experiment and b) simulation.

Here, n and $n^{\prime }$ are the occupations in the optically coupled and uncoupled valleys, respectively, Ω(t) is the time dependent Rabi frequency describing the optical excitation by co-circlarly polarized pulses, β and Γ are the dephasing and the decay rate, respectively, ℏω₀ is the exciton energy in the absence of the local field coupling, and V is the strength of the local field coupling. Compared to ref. [44] we do not find a significant impact of excitation induced dephasing as discussed in more detail in the Supporting Information. As recently studied in ref. [44] we additionally take an intervalley scattering contribution with the rate λ into account. In the special case λ = 0, the system reduces to a two-level system. We have recently investigated the local field model in the context of nonlinear optical signals focussing on FWM^[43] and pump-probe spectroscopy,^[44] showing that the handy description produces resonant optical spectroscopy signals that are consistent with experiments. In a nutshell, the local field effect leads to energy shifts of the exciton depending on the its occupation. In the limit of ultrafast optical pulses Equation (1) can be solved analytically relating the coherence p⁺ and the occupation n⁺ after the pulse to the respective values p⁻ and n⁻ immediately before the pulse [Image Omitted. See PDF][Image Omitted. See PDF]where θ is the pulse area and ϕ its phase. The approximations in Equation (2) describe the light-field induced contributions of a single pulse in first order of the pulse area, which is sufficient for reproducing the contributions of the order V² to the SWM signal, discussed in the main text of this work (other contributions, involving terms up to the second order in the pulse areas, are discussed in the Supporting Information).

Starting from the excitonic ground state characterized by $n=n_1^{-}=0$, $p=p_1^{-}=0$ the first pulse in the linear order creates [Image Omitted. See PDF][Image Omitted. See PDF]which shows that relevant exciton occupations will only be created by a second laser pulse from $p_1^{+}$. Our goal is to derive the nonlinear optical response in the lowest order of the pulse area θ because the experiments are performed with low pulse powers. In the case of the SWM signal this is the fifth order $\mathcal{O}({\theta }^5)$ (χ⁽⁵⁾-regime).

Once optical pulses have generated coherences p₀ and occupations n₀ the corresponding free propagation [Ω = 0, in Equation (1)] is governed by pure dephasing of p, exciton decay and intervalley scattering of $n^{(\prime )}$, and local field coupling between p and n. The corresponding dynamics can also be calculated analytically. Focusing on the optically addressed occupation first, its time-dependence reads [Image Omitted. See PDF]which leads to a balanced occupation between the two valleys n and $n^{\prime }$ on the timescale 1/(2λ) and a decay of both occupations with the rate Γ. As known from literature^[44–46] and as considered in this work the intervalley scattering is typically significantly faster than the decay, that is, Γ ≪ 2λ. As a simplifying approximation for the sake of interpreting the results we can therefore assume that the occupations are balanced rapidly after an optical pulse, that is, [Image Omitted. See PDF]This step might lead to slight deviations for short delays in the range τ ≈ 1/(2λ) which are however hardly visible for the parameters chosen here as shown in the Supporting Information. Further we can only optically manipulate the occupation n, while $n^{\prime }$ remains unchanged by the applied laser pulses. According to Equation (3a) these pulses moreover add phase labels ϕ. As explained at the beginning of this section and as practically applied below, we are only interested in specific phase combinations that describe the considered wave-mixing signal. Consequently, any change of the occupation n is labeled by phase factors which do not apply to the other valley $n^{\prime }$. Therefore, the latter is irrelevant for the final optical signal and we can consider $n_0^{\prime }=0$ in Equation (4b) resulting in [Image Omitted. See PDF]We want to remark that in the opposite limit of a very slow scattering rate 2λ ≪ Γ the occupation dynamics would directly be given by n(t) = n₀e^−Γt. In this case all later discussions would work in the same way. By replacing V → 2V in all the following derivations, one can even directly retrieve the corresponding equations.

Based on this approximation for the occupation dynamics we can calculate the coherence dynamics in the frame rotating with ω₀−V [Image Omitted. See PDF]In the first approximation step, we have used the approximated occuption from Equation (4c). As the exciton decay is much slower than the dephasing Γ ≪ β, in the second step we take Γ → 0. Finally, we keep the local field-induced contributions up to the second order in the exciton occupation, that is, considering Vn₀t ≪ 1. Below we will see that these are the terms which contribute to the SWM signal in the χ⁽⁵⁾-regime. In the Supporting Information, we show that only terms up to $\mathcal{O}(V^2)$ contribute to the χ⁽⁵⁾-regime of the SWM signal, while contributions with higher powers in V appear in higher orders of the optical field. The last equation tells us that the local field induced mixing of the occupation n₀ with the coherence p₀ into the coherence p(t) in the lowest order grows linearly in time with a rate given by n₀ and the local field strength V. We will use this argument later on.

We simulate the FWM signal with the phase combination 2ϕ₂ − ϕ₁ by calculating the coherence dynamics p(t) following a two-pulse sequence and filter this quantity with respect to the required phase combination as described in refs. [43, 44]. With this we find the signal dynamics $|S_{\mathrm{FWM}}{|\approx |}^{2{\varphi }_2-{\varphi }_1}{p}_2(t,\tau )|$ depicted in Figure 3b. To achieve the excellent agreement with the experiment in Figure 3a, we used a local field strength of V = 100 ps⁻¹ and pulse areas of θ₁ = θ₂/2 = θ = 0.02π, a Gaussian pulse duration of Δt = 70 fs (standard deviation), a dephasing rate of β = 3 ps⁻¹, an intervalley scattering rate λ = 4 ps⁻¹, and a decay rate of Γ = 0.6 ps⁻¹. Note that the depicted signal was calculated numerically because we considered a non-vanishing pulse duration and we did not employ the approximations mentioned in Equations (2)–(5). To set the strength of the local field coupling into context, typical values for GaAs quantum wells were reported in the range of a few meV,^[47,48] which is at least one order of magnitude smaller than considered here.

Six-Wave Mixing Dynamics

In the present study we go one step further in multi-wave mixing and consider one of the possible SWM signals generated by three laser pulses, namely ϕ_SWM = 2ϕ₃ − 2ϕ₂ + ϕ₁. In principle there are two delays in this pulse sequence that could be varied but, as schematically shown in Figure 1 (right), we set the delay between the first two pulses to τ₁₂ = 0 and only vary the second one τ₂₃ = τ. The impact of a non-vanishing τ₁₂ is discussed in the Supporting Information. In the pure two-level system this configuration probes the polarization, that is, it contains the same information as the FWM signal discussed before as shown in the Supporting Information. The measured SWM dynamics as a function of the real time t after the third pulse and the delay τ are shown in Figure 4a. Here, the two pulses with ϕ₁ and ϕ₂ arrive at t = −τ and the pulse with ϕ₃ at t = 0. The signal consists of a strong maximum at small t ≈ 0.5 ps and τ ≈ 0. Moving to negative delays τ < 0 (pulse 3 is arriving before 1 and 2), the signal is strongly damped. For positive delays τ > 0 it decays much slower on the same timescale as the FWM signal in Figure 3. We find a remarkable depression of the signal that stretches along the curved diagonal given by Equation (9) (dashed black line), as will be derived below. In correspondence with the previously described constructive signal enhancement in the photon echo, we call this pronounced signal reduction a destructive photon echo. Later, we will discuss criteria justifying the labeling of this feature as an echo effect.

Figure 4. Six-wave mixing dynamics as a function of real time t and delay τ. a) Experiment, b) simulation, with the curved dashed line depicting Equation (9).

To identify the origin of this peculiar dynamical feature we model the SWM signal within the local field model described above. We consider the same system parameters as for the FWM signal but choose equal pulse areas for all three pulses θ₁ = θ₂ = θ₃ = θ = 0.02π, in agreement with the experiment. In the Supporting Information, it is shown that the exact choice of θ in the low excitation regime does not change the SWM signal dynamics. The SWM signal is extracted via the phase combination 2ϕ₃ − 2ϕ₂ + ϕ₁ and the signal dynamics $|S_{\mathrm{SWM}}{|\approx |}^{2{\varphi }_3-2{\varphi }_2+{\varphi }_1}{p}_3(t,\tau )|$ are depicted in Figure 4b. The retrieved signal agrees very well with the respective experiment in Figure 4a showing the same characteristic suppression of the signal, that is, the destructive photon echo. We give an overview regarding the impact of the different parameters on the SWM signal in the Supporting Information. Most importantly, it is shown in the Supporting Information, that the SWM dynamics do not exhibit the suppression for small local field strengths V demonstrating that this feature is a result of exciton-exciton interaction. Other specific features in the signal's dynamics like the asymmetric decay between positive and negative delays or between τ and t were studied in detail in ref. [43] and behave similarly in SWM. The destructive photon echo effect should also be present in an inhomogeneously broadened system, where also a traditional constructive photon echo develops. As the constructive echo selects only a specific time interval of the emitted signal around t = τ, it leads to a suppression of the entire SWM signal for all other times. Consequently, the SWM signal and therefore the destructive photon echo, which bends away from the diagonal (discussed later), are only visible in the vicinity of t = τ. This aspect is studied in more detail in the Supporting Information.

Note, that the simulation shown in Figure 4b takes the non-vanishing pulse duration into account and is therefore performed numerically. In the limit of ultrafast laser pulses we can find analytical expressions for the SWM signal. Given that the experiment is carried out with weak pulse powers, we restrict the following studies on the lowest order in the optical field which is the χ⁽⁵⁾-regime. In this order we have already eight different contributions as derived in the Supporting Information. From these we will focus on the ones with the strongest local field contribution which is $\mathcal{O}(V^2)$, that is, we omit all terms $\mathcal{O}(V^1)$ and $\mathcal{O}(V^0)$. The reason for this is the absence of the destructive echo for small V. In ref. [43] we have derived a flow chart representation for the construction of nonlinear wave mixing signals. In Figure 5 we employ this procedure to disentangle the origin of the different contributions to the SWM signal. The flow chart only shows coherences p (blue) and occupations n (red) with phase combinations (green, given as left indices) relevant for the final signal; corresponding flow charts for the contributions with $\mathcal{O}(V^1)$ and $\mathcal{O}(V^0)$ are provided in the Supporting Information. Conveniently, the $\mathcal{O}(V^2)$-contributions can be derived with the approximations given in Equation (2). Note, that we introduce phase differences Δ_nm = ϕ_n − ϕ_m here. The right lower index refers to the pulse number, while the upper − (+) indicates times immediately before (after) that pulse. Starting from the excitonic ground state with n = 0, p = 0 the first pulse generates the occupation ${}^0{n}_1^{+}$ and the coherence ${}^{{\varphi }_1}{p}_1^{+}$. Note, that in the scheme we do not restrict ourselves to the lowest order contributions in the light field and thus include also the occupation which is of second order in the pulse amplitude. The second pulse arrives at the same time (τ₁₂ = 0) and creates two occupations ${}^0{n}_2^{+}$ and ${}^{|{\mathrm{\Delta }}_{21}|}{n}_2^{+}$ and the coherence ${}^{{\varphi }_2}{p}_2^{+}$. During the following propagation for the time of the delay τ two relevant things happen: On the one hand all coherences experience dephasing (blue arrows). On the other hand ${}^{{\varphi }_2}{p}_2^{+}$ is mixed with ${}^{|{\mathrm{\Delta }}_{21}|}{n}_2^{+}$ via the local field coupling resulting in the coherence ${}^{2{\varphi }_2-{\varphi }_1}{p}_3^{-}$ before the third pulse. Note, that this contribution carries the FWM phase ϕ_FWM = 2ϕ₂ − ϕ₁ which we will come back to below. After the final third pulse we have five relevant terms^[49]: The coherence ${}^{{\varphi }_3}{p}_3^{+}$ and the three occupations ${}^{|{\mathrm{\Delta }}_{32}|}{n}_3^{+}$, ${}^{|{\mathrm{\Delta }}_{32}-{\mathrm{\Delta }}_{21}|}{n}_3^{+}$, and ${}^{|{\mathrm{\Delta }}_{21}|}{n}_3^{+}$. The polarization ${}^{{\varphi }_1}{p}_1^{+}$ is not affected by the second and third pulse and just evolves into ${}^{{\varphi }_1}{p}_3^{+}$. During the remaining propagation step in real time t the three relevant SWM contributions are generated by local field mixing processes according to Equation (5).

Figure 5. Flow chart for the three main contributions to the SWM signal listing intermediate phase-filtered coherences (p, blue) and occupations (n, red). Green waved arrows show pulse induced, violet dotted ones local field induced, black ones free dynamics without, and blue ones with dephasing. The phase differences are defined as Δnm = ϕn − ϕm. The flow chart holds for any order of the optical field. However, in the linear response regime considered here, it is 0n1+=0n2+=0n3−=0.

Considering the contribution on the right first, we have to mix ${}^{{\varphi }_1}{p}_3^{+}$ with ${}^{|{\mathrm{\Delta }}_{32}|}{n}_3^{+}$ twice. This results in the phase combination ϕ₁ + 2(ϕ₃ − ϕ₂) which is the SWM phase combination. According to Equation (5), each of these local field mixing processes contributes with an amplitude of Vt resulting in the amplitude (Vt)². In addition the amplitude is damped due to the dephasing happening during the delay. Following the two paths in the diagram back to this propagation step we find that ${}^{{\varphi }_1}{p}_2^{+}\to {}^{{\varphi }_1}{p}_3^{-}$ and ${}^{{\varphi }_2}{p}_2^{+}\to {}^{{\varphi }_2}{p}_3^{-}$ contribute with a dephasing term ∼e^−βτ. The latter one is used twice due to the double local field mixing resulting in the total damping rate of e^−3βτ.

Moving on to the left contribution in Figure 5, we find that in the last propagation the coherence ${}^{{\varphi }_3}{p}_3^{+}$ is local-field mixed with ${}^{|{\mathrm{\Delta }}_{32}|}{n}_3^{+}$ and ${}^{|{\mathrm{\Delta }}_{21}|}{n}_3^{+}$ once, resulting in the SWM phase combination ϕ₃ + (ϕ₃ − ϕ₂) − (ϕ₂ − ϕ₁). Each of these processes contributes an amplitude of Vt again resulting in (Vt)². The difference to the first term is that here only ${}^{|{\mathrm{\Delta }}_{32}|}{n}_3^{+}$ stems from a coherence, namely ${}^{{\varphi }_2}{p}_2^{+}$ which experiences a dephasing during the delay. Consequently, the entire contribution is damped by e^−βτ. This already shows that this contribution is more important for the final SWM signal than the right term discussed first.

The final main contribution is the middle one in Figure 5. Here, in the last real time propagation the local field mixing happens between ${}^{{\varphi }_3}{p}_3^{+}$ and ${}^{|{\mathrm{\Delta }}_{32}-{\mathrm{\Delta }}_{21}|}{n}_3^{+}$ which contributes with an amplitude of Vt. Following the path back, we find that the occupation is created from the polarization ${}^{2{\varphi }_2-{\varphi }_1}{p}_3^{-}$ which itself is produced by a local field mixing step between ${}^{{\varphi }_2}{p}_2^{+}$ and ${}^{|{\mathrm{\Delta }}_{21}|}{n}_2^{+}$. The mixing process lasts for the delay time and therefore contributes with a factor Vτ to the final SWM amplitude. This is also the step where the only dephasing happens resulting in a final amplitude of V²τte^−βτ.

Directly comparing the three contributions in Figure 5 we find that the right one is smaller than the other two due to the stronger dephasing during the delay. We will therefore disregard this term from now on. The other two terms are of the order V² and exhibit the same dephasing with e^−βτ. From the full derivation given in the Supporting Information we find that these two terms in Figure 5 carry opposite signs. We will explain these signs later when introducing an effective Bloch vector picture to illustrate the destructive echo formation. Finally, adding the two terms we get [Image Omitted. See PDF]where we also added the dephasing rate for the propagation in t after the last pulse. We find that the two terms exactly compensate each other for t = τ. This renders our first step toward the understanding of the destructive echo. We identified two paths (signal contributions) with the same damping but with different magnitudes depending on the delay τ that act destructively in the total SWM signal. To illustrate the interplay between these two terms in Figure 6 we schematically plot the applied laser pulses and the absolute value of the final SWM signal retrieved from the previous derivations. We also include the intermediate growth of the FWM coherence lasting for the delay τ.

Figure 6. Schematic picture of the pulse sequence and signal development from the two main contributions.

The full SWM signal in the χ⁽⁵⁾-regime is derived in the Supporting Information and reads [Image Omitted. See PDF]Interestingly, we also find a suppression of the signal when only considering the terms $\mathcal{O}(V^1)$, namely for t = τ/4. We obviously do not find such a feature in our measurement, which shows that the linear order in V does not have a significant contribution.

One advantage of the spectral interferometry in the applied approach is the possibility to detect—besides the amplitude—also the phase of the SWM signal.^[50] We see that for t < τ the positive contribution ∼τt dominates while for t > τ the negative one ∼−t² is larger. Therefore we expect a phase jump when crossing the destructive echo in time t < τ → t > τ. To confirm this in Figure 7a, we plot the measured (solid blue) and calculated (dashed red) SWM signal amplitude as a function of time t for τ = 0.35 ps. We slightly adjusted the pulse duration to Δt = 80 fs to achieve the good agreement with this experiment. The finding that the simulated signal does not drop to zero shows the impact of the non-compensating contributions, like the ordinary SWM signal in Equation (7), and the influence of the considered non-vanishing pulse duration. In Figure 7b we show the corresponding SWM phase from S_SWM = |S_SWM|e^iφ. We indeed find a jump of approximately π at the destructive echo as would be expected for the two dominant terms. However, all other contributions, which are naturally present in the numerical simulation, lead to a reduction of the phase jump and a further distortion of the destructive echo dynamics. We again find that the full depression of the signal happens at t ≈ 0.4 ps which is slightly later than t = τ = 0.35 ps. This slight delay of the effect mainly stems from the non-vanishing decay rate of the exciton Γ ≠ 0.

Figure 7. Amplitude and phase dynamics of the destructive echo. a) SWM amplitude dynamics for a delay of τ = 0.35 ps. Experiment as solid blue and simulation as dashed red line. b) Respective dynamics of the phase of the SWM signal.

To include the exciton decay we have to go back to the full dynamics in Equation (5) and directly consider Vn₀(1 − e^−Γt)/Γ ≪ 1 without taking the limit Γ → 0 leading to [Image Omitted. See PDF][Image Omitted. See PDF]Now we consider the different paths in Figure 5 to determine the decay's impact on the different SWM contributions. The term ∼t² stems from the delay step ${}^{|{\mathrm{\Delta }}_{21}|}{n}_2^{+}\to {}^{|{\mathrm{\Delta }}_{21}|}{n}_3^{-}$, contributing a factor e^−Γτ, and the final double wave mixing step gives (1 − e^−Γt)². The τt-contribution has a local-field mixing step during the delay leading to (1 − e^−Γτ) and only a single local field mixing after the last pulse, that is, (1 − e^−Γt). Finally, the two contributions compensate each other for [Image Omitted. See PDF]We find that the decay additionally acts on the local field mixing processes because it dynamically reduces the occupation during the mixing. This affects the t²-contribution in a different way than the part ∼τt and the destructive echo gets delayed with respect to t = τ. Equation (9) does not depend on V as confirmed in the Supporting Information. Consequently, it also holds in both limiting cases of small and large intervalley scattering as discussed below Equation (4c) and can be expected to be correct also for intermediate λ-values. Note, that the curve describing the dynamics of the SWM signal depression in Figure 4 follows Equation (9) and is additionally shifted by the pulse width Δt to shorter times to compensate the non-vanishing pulse duration. The curve almost perfectly follows the distribution of the destructive echo.

Quasi-Bloch Vector and Phase Inhomogeneity

As discussed in the previous section, the FWM-polarization constitutes the basis for the SWM signal. Therefore, it is instructive to introduce the concept of the quasi-Bloch vector to explain the wave-mixing origin for the FWM signal in the next section. Building on this picture we will be able to explain the destructive photon echo formation in the SWM signal in the following section.

Four-Wave Mixing Quasi-Bloch Vector

The ordinary echo in FWM is most instructively visualized within the Bloch vector picture in Figure 2 where the inhomogeneity leads to different rotation speeds of the various Bloch vector realizations. In that case each exciton energy is represented by a single Bloch vector. Here, we are dealing with a slightly different situation. The wave mixing experiment is repeated numerous times with successively different phase combinations (ϕ₁, ϕ₂, ϕ₃) for each repetition of the measurement. Therefore we have to represent each run of the experiment, that is, each phase combination, by a single Bloch vector. To find the realized Bloch vectors we have a look at the coherence and occupation immediately after the second pulse^[43] [Image Omitted. See PDF][Image Omitted. See PDF]where we consider all terms up to the third order in the optical field (χ⁽³⁾-regime) because our measurements are performed with low excitation powers. Following the flow chart in Figure 5 we find that the two terms ${}^{{\varphi }_2}{p}_2^{+}$ and ${}^{|{\mathrm{\Delta }}_{21}|}{n}_2^{+}$ perform local field induced mixing during the following propagation before the third pulse (delay step). Exactly this term was identified as a contribution of the FWM signal in ref. [43]. Therefore, we will study this process in more detail. Our goal is to isolate a set of Bloch vectors that allows us to directly extract the final wave mixing signal by integrating over the entire set. To achieve this we have to already filter the coherence with the respective wave mixing phase factor. In this case we are interested in the FWM contribution, which means that we have to consider the phase factor $e^{i{\varphi }_{\mathrm{FWM}}}=e^{i(2{\varphi }_2-{\varphi }_1)}$ and consequently the filtered polarization ${}^{{\varphi }_2}{p}_2^{+}{e}^{-i{\varphi }_{\mathrm{FWM}}}$. This brings us to the two relevant terms [Image Omitted. See PDF][Image Omitted. See PDF]We use these expressions to define a quasi-Bloch vector for this FWM contribution via [Image Omitted. See PDF]with the previously introduced phase difference Δ₂₁ = ϕ₂ − ϕ₁. We call this quantity quasi-Bloch vector because it does not represent the entire density matrix, as the full Bloch vector does. When varying Δ₂₁ as it is done in experiment and simulation, we already see that this set of quasi-Bloch vectors follows a 3D Lissajous curve with a frequency relationship of 1:1:1 as depicted in Figure 8a. It appears as a circle tilted around the Re($\stackrel{\sim }{p}$)-axis and its projections on the different planes of the coordinate system form Lissajous curves with 1:1 frequency relations. This results in a distribution of occupations $\stackrel{\sim }{n}$, whose spread, represented by the tilt angle, is given by the considered pulse areas.

Figure 8. Quasi-Bloch vector of the FWM contribution. a) Lissajous curve of the initial distribution of quasi-Bloch vectors. b) Dynamics of the quasi-Bloch vectors. Initial distribution as bright open circles and for a time t > 0 after the second pulse as filled dark circles. The final integrated FWM coherence FWMp is marked as blue cross and its propagation is given by Vt.

Starting from this distribution of quasi-Bloch vectors directly after the second laser pulse their dynamics is governed by the local field mixing between $\stackrel{\sim }{p}$ and $\stackrel{\sim }{n}$ described by [see Equation (5)] [Image Omitted. See PDF]These dynamics are a rotation of the coherence $\stackrel{\sim }{p}$ depending on the occupation $\stackrel{\sim }{n}$. Particularly important is that the rotation flips its direction for opposite signs of $\stackrel{\sim }{n}$. At the same time $\stackrel{\sim }{n}$ remains unaffected. Finally, from the distribution of quasi-Bloch vectors at a given time $\stackrel{\sim }{\mathbf{v}}(t)$ we can directly retrieve the final wave mixing contribution by integrating over all possible phase combinations, that is, over Δ₂₁. We directly see that the integral over $\stackrel{\sim }{n}$ always vanishes and we are left with the final FWM coherence [Image Omitted. See PDF]

For the initial distribution of quasi-Bloch vectors in Equation (12) we find that the final FWM coherence vanishes because all vectors are equally distributed on the tilted circle in Figure 8a. However, the following dynamics lead to a non-vanishing FWM contribution. The rotation of $\stackrel{\sim }{p}$ with frequencies proportional to $\stackrel{\sim }{n}$ has two crucial consequences: i) The presence of a distribution of $\stackrel{\sim }{n}$ results in different frequencies. As already discussed in ref. [43] this means that each run of the experiment results in a slightly different emission energy of the exciton and consequently leads to a broadened FWM spectrum. Because the distribution of $\stackrel{\sim }{n}$ stems from the variation of the applied laser pulse phase combinations we call this spectral broadening local field induced phase inhomogeneity. ii) According to Equation (13) the sign of $\stackrel{\sim }{n}$ determines the rotation direction of the quasi-Bloch vectors. The initial distribution of quasi-Bloch vectors in Figure 8a is a circle that is tilted around the Re$(\stackrel{\sim }{p})$-axis and all points on this axis are not affected by the local field induced rotation ($\stackrel{\sim }{n}=0$). This axis remains a symmetry axis also for the dynamics of the quasi-Bloch vectors: one half of the circle rotates clockwise, the other counter-clockwise in exactly the same way. Consequently, when integrating over the relative phase Δ₂₁ the opposing Im$(\stackrel{\sim }{p})$ contributions compensate each other and the final FWM coherence is real, that is, $\mathrm{Re}(\stackrel{\sim }{p})\ne 0$. Because only Re($\stackrel{\sim }{p}$) will later contribute to the FWM signal, the crucial quantity determining the propagation direction of Re($\stackrel{\sim }{p}$) is $\stackrel{\sim }{n}\mathrm{Im}(\stackrel{\sim }{p})$ as can be seen in Equation (13). In Figure 8b we show the quasi-Bloch vector dynamics in the Re($\stackrel{\sim }{p}$),Im($\stackrel{\sim }{p}$)-plane after the second pulse. The bright open circles represent the initial homogeneous distribution and the filled dark circles the distribution for a time t > 0. Looking at the initial distribution in Figure 8a we find that the semicircle with Im($\stackrel{\sim }{p})<0$ also has $\stackrel{\sim }{n}<0$ and the other semicircle with Im($\stackrel{\sim }{p})>0$ has $\stackrel{\sim }{n}>0$. This means that the former rotate clockwise, while the latter rotate counter-clockwise as indicated by the curved black arrows in Figure 8b. Consequently, the weight of all quasi-Bloch vectors, which is depicted by the blue cross, moves into the positive Re($\stackrel{\sim }{p}$)-direction. The speed of this movement is given by V and the tilt angle of the initial quasi-Bloch vector distribution. This illustrates that the FWM contribution grows in time and that the local field induced phase inhomogeneity governs this process.

Six-Wave Mixing Quasi-Bloch Vector

Equipped with the quasi-Bloch vector picture we can retrieve the two most important SWM signal contributions by carrying out an analogue discussion. According to Figure 5, from all possible contributions after the third laser pulse, we only need the coherence ${}^{{\varphi }_3}{p}_3^{+}$ and the occupations ${}^{|{\mathrm{\Delta }}_{32}|}{n}_3^{+}$, ${}^{|{\mathrm{\Delta }}_{32}-{\mathrm{\Delta }}_{21}|}{n}_3^{+}$, and ${}^{|{\mathrm{\Delta }}_{21}|}{n}_3^{+}$. Omitting all other terms we get [Image Omitted. See PDF][Image Omitted. See PDF]To generate the SWM contribution ∼τt we need to local-field mix the coherence filtered with respect to the SWM phase combination 2ϕ₃ − 2ϕ₂ + ϕ₁ with the last term of the occupation in Equation (15b), that is, [Image Omitted. See PDF][Image Omitted. See PDF]From this we can directly read the respective initial quasi-Bloch vector [Image Omitted. See PDF]with α = Δ₃₂ − Δ₂₁. We see that this is the same Lissajous curve structure as in the previously discussed FWM case. The only differences are that the phase variation is now given by Δ₃₂ − Δ₂₁, which still uniformly covers all angles, a tilt around the Im($\stackrel{\sim }{p}$)-axis, and a tilt angle such that the occupation distribution $\stackrel{\sim }{n}$ has an amplitude of −θ⁴Vτ/8. The subsequent free dynamics are again governed by Equation (13).

Before discussing the shape and dynamics of this quasi-Bloch vector set we consider the other relevant SWM contribution ∼t². Therefore, we have a look at [Image Omitted. See PDF][Image Omitted. See PDF]

We are finally only interested in the coherence ^SWMp₃ describing the SWM signal. So, at this point we can already perform the integration [Image Omitted. See PDF]This leaves us with the quasi-Bloch vector for this SWM contribution reading [Image Omitted. See PDF]with α = Δ₃₂ − Δ₂₁, which is obviously a 3D Lissajous curve with the frequency ratio 2:2:1.

The two initial distributions of quasi-Bloch vectors given in Equations (17) and (20) are depicted in Figure 9 in blue and violet, respectively. We vividly see the Lissajous curves leading to the distributions in $\stackrel{\sim }{n}$-direction. While the 3D perspective of the curves is appealing, we will consider the projections to the three different planes of the coordinate system when discussing the dynamics of the quasi-Bloch vectors, which are depicted in pale colors. These projections are shown in Figure 10, where different times are sorted in rows (increasing from top to bottom) and different perspectives in columns.

Figure 9. Quasi-Bloch vector distributions after the 3rd laser pulse. The t2-contribution (violet) and τt (blue) form Lissajous curves with frequency ratios 2:2:1 and 1:1:1, respectively.

Figure 10. Projections of the quasi-Bloch vector representation of the destructive echo formation. Columns show projections of the quasi-Bloch vector defined in Equation (14) as listed on the top, the rows show different real times as given on the left. The contribution ∼tτ is shown as blue, ∼t2 as violet dots. The size of the violet dots represents the product n∼Re(p∼), positive values are larger and negative ones smaller. The situation for t = 0 is shows as bright circles for all other times. The black arrows indicate the movement of the quasi-Bloch vectors, solid (dashed) for positive (negative) n∼Re(p∼). The final SWM coherences are marked as crosses in the two right columns.

First considering the τt-contribution (blue) we find essentially the same situation as for the intermediate FWM contribution discussed in Figure 8. The only important difference is that for the SWM contribution the circle is tilted around the Im($\stackrel{\sim }{p}$)-axis. Following the same arguments as before the quasi-Bloch vector realizations move into the positive Im($\stackrel{\sim }{p}$)-direction. This movement is shown in the second column from the right in Figure 10. The initial homogeneous distribution is depicted by bright empty circles and for the considered times (written on the left) as filled dark circles. The movement is marked by the curved black arrow and the final SWM contribution after integration over all quasi-Bloch vector realizations as a dark blue cross. The speed of this movement is again constant and given by Vτ, which is proportional to the spread of $\stackrel{\sim }{n}$ and therefore given by the tilt angle of the circle.

For the t²-contribution (violet) the situation is more involved because the Lissajous curve is less trivial. In this situation the crossing point of the two loops lies at ${\stackrel{\sim }{\mathbf{v}}}_{\mathrm{SWM},t^2}=(0,\mathrm{Im}(\stackrel{\sim }{p}),0)$. Because of the vanishing occupation this point will not move in the following propagation and the Im$(\stackrel{\sim }{p})$-axis is again the symmetry axis of the dynamics. Consequently, the opposing Re$(\stackrel{\sim }{p})$ values compensate each other, such that the final SWM contribution is again purely imaginary. Therefore, we have to consider the movement into the direction governed by the product $\stackrel{\sim }{n}\mathrm{Re}(\stackrel{\sim }{p})$. When looking at the left column in Figure 10 the distribution reaches values with $\stackrel{\sim }{n}\mathrm{Re}(\stackrel{\sim }{p})>0$ (large dots) and $\stackrel{\sim }{n}\mathrm{Re}(\stackrel{\sim }{p})<0$ (small dots). Following Equation (13) the large dots move into the negative Im($\stackrel{\sim }{p}$)-direction while the small ones move into the positive Im($\stackrel{\sim }{p}$)-direction as indicated by the black horizontal arrows in the Im($\stackrel{\sim }{p}),\stackrel{\sim }{n}$-plane at t < τ. In the complex $\stackrel{\sim }{p}$-plane this leads to two counter movements marked as arrows in the right column. The small dots move up (dashed curved arrows) and the large dots move down (solid curved arrows). Initially both of these movements start with the same velocity because large and small dots have the same absolute values $|\stackrel{\sim }{n}\mathrm{Re}(\stackrel{\sim }{p})|$, that is, the Lissajous curve is symmetric. Therefore, the integrated SWM signal starts with a vanishing velocity for small t ≈ 0. However, during the following dynamics the quasi-Bloch vectors with positive Im($\stackrel{\sim }{p}$)-velocity (small dots) rotate collectively in such a way that this part of the curve obtains smaller occupations. This can be seen by focussing on the bottom-left and upper-right corners of the panels in the left column. Accordingly, those with negative Im($\stackrel{\sim }{p}$)-velocity (large dots) obtain larger $|\stackrel{\sim }{n}|$ as can be seen in the second column [Im($\stackrel{\sim }{p}),\stackrel{\sim }{n}$-plane] marked by the horizontal black arrows. In the left Lissajous curve this means that the large dots move toward the top left and bottom right corner of the plot, while the small dots move inwards. Consequently, the negative Im($\stackrel{\sim }{p}$)-velocity component gains speed because $|\stackrel{\sim }{n}\mathrm{Re}(\stackrel{\sim }{p})>0|$ (large dots) grows, while the positive component slows down as $|\stackrel{\sim }{n}\mathrm{Re}(\stackrel{\sim }{p})<0|$ (small dots) shrinks. In the complex $\stackrel{\sim }{p}$-plane (right column) this means that the upwards movement of the quasi-Bloch vectors represented by small dots (dashed arrows) slows down for increasing times t while the downwards movement of the large dots (solid arrows) speeds up. In summary, for the integrated SWM coherence, marked by the violet cross, we find an acceleration towards negative Im($\stackrel{\sim }{p}$)-values given by −V²t².

Exactly at t = τ in the third row the two integrated SWM contributions marked in the third and forth column have the same distance from the center as indicated by the dashed horizontal lines which show that they compensate each other. At this time the destructive echo appears due to the destructive quantum interference of all depicted quasi-Bloch vectors. For even larger times t < τ (bottom row) the accelerated t²-contribution is stronger than the one with constant velocity ∼τt and the SWM is non-vanishing again.

Conclusion

To conclude we again compare the traditional photon echo in FWM with the newly discovered destructive photon echo in SWM. As summarized in Table 1 we define three criteria that characterize the echo formation. Firstly, the timing of the photon echo formation is basically determined by the chosen delay τ between the two laser pulses. Slight deviations from the exact t = τ timing are used to learn about internal dynamics of the studied system.^[51] The same holds for the destructive echo in SWM. As discussed in detail the two main contributions compensate each other exactly at t = τ. Nonetheless, in the detected and simulated full signal dynamics in Figure 4 the depression happens at t > τ. The reason is that exciton decay dynamically changes the slope of the destructive echo. Secondly, the photon echo appears because all considered Bloch vectors form a constructive interference. To explain this effect only the conventional FWM signal has to be considered, which would appear as a single path in a flow chart like the one in Figure 5 (see also ref. [43]). In the case of the destructive echo we have shown that we need two paths which interfere destructively to explain the novel echo effect. Finally, the fundamental source of the photon echo is any sort of inhomogeneity, which might appear in space by a locally varying exciton energy,^[52] or in time via external noise that acts on the optical transition energy.^[30] In the case of the destructive echo we have shown that the local field together with a variation of the applied laser pulse phases results in a spectral broadening. As discussed in ref. [43] this local field induced phase inhomogeneity can already be detected in FWM spectra. This summary shows that the criteria defining the traditional photon echo can also be applied to the destructive echo effect. Therefore, we conclude that we really found a new photon echo effect.

Table 1 Comparison of criteria defining the photon echo in FWM and the destructive photon echo in SWM

	Echo (FWM)	Destructive echo (SWM)
Timing	Delay + internal dynamics	Delay + internal dynamics
Interference	Constructive (1 path)	Destructive (2 paths)
Source	Inhomogeneous broadening	Local field induced phase inhomogeneity

In addition we have developed a quasi-Bloch vector picture to illustrate the generation of the destructive echo inspired by the instructive Bloch vector image of the photon echo. We found that the SWM-relevant quasi-Bloch vectors are distributed along Lissajous curves which adds surprising aesthetics to the involved Bloch vector dynamics in SWM.

Following this proof-of-principle demonstration of a new photon echo effect the natural next tasks will be to further explore the destructive photon echo's application possibilities. One obvious step will be to measure photon echos at various spots of the sample and thereby get a measure for the spatial distribution of the local field coupling strength. We know that different pulse sequences in FWM probe different observables.^[53] Therefore we hope that novel pulse sequences in higher wave-mixing spectroscopy will be designed that should allow to access parameters governing the nonlinear light-induced dynamics in semiconductors, like the strength of interaction mechanisms, that are not easily accessible with traditional techniques. These developments are particularly attractive for materials with a large nonlinearity promoting strong wave mixing responses like TMDCs.

Experimental Section

In the experiment, we employ a mode-locked laser (Tsunami Femto provided by Spectra Physics) generating femtosecond pulses at a repetition rate of 76 MHz. The center wavelength is tuned to 752 nm, which corresponds to the A exciton transition in the investigated hBN/MoSe₂/hBN heterostructure at T = 5 K. The sample is kept in an optical He-flow cryostat and microspectroscopy is performed using an Olympus microscope objective (LCPLN50XIR, numerical aperture of 0.65), which is installed on an XYZ piezo stage from Physik Instrumente. To perform the multi-pulse, heterodyne experiments, the initial pulse train from the laser source is split into three and each of these beams is focussed into a separate acousto-optic modulator (AOM). The AOMs are driven at distinct radio-frequencies of Ω₁ = 80, Ω₂ = 79, and Ω₃ = 79.77 MHz, such that the deflected beams acquire corresponding frequency upshifts. Next, the time delays between the first two (τ₁₂) and the last two beams (τ₂₃) are introduced by a pair of mechanical delay stages. In addition, a grating-based pulse shaper is employed to correct the temporal chirp of the pulses, when passing through optical elements of the setup, especially AOMs and the microscope objective. The such prepared pulse sequence is then recombined into the same spatial mode and focussed at the heterostructure sample reaching a diffraction limited spot diameter of around 0.8 µm. The nonlinear optical signal from the sample is retrieved in the back-reflectance geometry. To isolate the desired wave-mixing response another AOM is used, operating at the heterodyne frequency generated by a home-made analogue radio-frequency mixing electronics, assembled from individual components provided by Mini-Circuits. For example, the studied SWM signal is detected at (Ω_SWM = 2Ω₃ − 2Ω₂ + Ω₁)=81.54 MHz. After being deflected from the AOM the signal is frequency downshifted by Ω_SWM and the unique wave-mixing component under consideration carries the original frequency from the laser source. The signal can now interfere with a reference beam from the same laser source, which propagates in the vicinity of the driving beams. At the same time all other optical response components present in the signal reflected from the sample are still modulated in the MHz range and thus average out completely during a single acquisition time of a few ms. Importantly, the mixing AOM operates in a Bragg configuration. Hence, simultaneously the reference beam gets deflected onto the signal, receiving a frequency upshift by Ω_SWM. This again allows for interference with the signal, which was not deflected by the AOM (operating at Ω_SWM). In total, we generate two beams in which the considered wave-mixing signal is interfering with the reference pulse, but with the opposite phase, as required by energy conservation. The spectrally-resolved interference is obtained with an imaging spectrometer (Princeton Instruments with a focal length of 750 mm) and detected with an CCD camera (PIXIS from Princeton Instruments, with an eXcelon coating). A background free, shot-noise limited detection is achieved by a balanced detection and multi-acquisition. Namely, we exploit the Bragg configuration and subtract the two mixed beams impinging at different positions on the CCD. In addition, the phase of the mixing AOM is cycled between 0 and π to overcome any classical noise from the CCD camera. The wave-mixing signal amplitude and phase are obtained from the interferogram by performing spectral interferometry.

Acknowledgements

T.H. thanks the German Academic Exchange Service (DAAD) for financial support (No. 57504619). D.W. thanks the Polish National Agency for Academic Exchange (NAWA) for financial support within the ULAM program (No. PPN/ULM/2019/1/00064). T.H., D.W., T.K., and P.M. acknowledge support from NAWA under an APM grant (No. PPI/APM/2019/1/00085). M.P., D.V., and M.B. thank the EU Graphene Flagship project. M.P., D.V., K.N., and J.K. acknowledge support from CNRS via IRP “2D materials”. M.P. and K.N. acknowledge support from TEAM programme of the Foundation for Polish Science, co-financed by the EU within the ERDFund. M.B. is supported by the ESF under the project CZ.02.2.69/0.0/0.0/20_079/0017436. K.W. and T.T. acknowledge support from the Elemental Strategy Initiative conducted by the MEXT, Japan (Grant Number JPMXP0112101001) and JSPS KAKENHI (Grant Numbers 19H05790 and JP20H00354).

Conflict of Interest

The authors declare no conflict of interest.

Data Availability Statement

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