Influence of local fields on the dynamics of four-wave mixing signals from 2D semiconductor systems

Based on the analytic results for a single ultrashort pulse derived in equations (8) and (9), we start our discussion by calculating the two-pulse FWM signal analytically in that limit. The polarization after the second pulse depends on the time t after that pulse and the delay τ₁₂ within the pulse pair and reads

$$\begin{align}\hfill {p}_{2}\left(t,{\tau }_{12}\right)& =\left[a\enspace {\text{e}}^{\text{i}\left({\phi }_{1}+{\Phi}\right)}+b\enspace {\text{e}}^{\text{i}{\phi }_{2}}+c\enspace {\text{e}}^{\text{i}\left(2{\phi }_{2}-{\phi }_{1}-{{\Phi}}^{{\ast}}\right)}\right]{\mathrm{e}}^{-\beta t}\hfill \\ \hfill & \quad {\times}\mathrm{exp}\left(\mathrm{i}\left\{\alpha +\eta \enspace {\mathrm{e}}^{-\mathrm{I}\mathrm{m}\left({\Phi}\right)}\enspace \mathrm{cos}\left[\mathrm{R}\mathrm{e}\left({\Phi}\right)+{\phi }_{1}-{\phi }_{2}\right]\right\}\right),\hfill \end{align} \tag{ 11a }$$

with

$$\begin{equation}a\left({\tau }_{12}\right)=\frac{\mathrm{i}}{2}\enspace \mathrm{sin}\left({\theta }_{1}\right){\mathrm{cos}}^{2}\left(\frac{{\theta }_{2}}{2}\right)\mathrm{exp}\left(-\beta {\tau }_{12}\right),\end{equation} \tag{ 11b }$$

$$\begin{equation}b\left({\tau }_{12}\right)=\frac{\mathrm{i}}{2}\enspace \mathrm{sin}\left({\theta }_{2}\right)\left[1-2\enspace {\mathrm{sin}}^{2}\left(\frac{{\theta }_{1}}{2}\right){\mathrm{e}}^{-{\Gamma}{\tau }_{12}}\right],\end{equation} \tag{ 11c }$$

$$\begin{equation}c\left({\tau }_{12}\right)=-\frac{\mathrm{i}}{2}\enspace \mathrm{sin}\left({\theta }_{1}\right){\mathrm{sin}}^{2}\left(\frac{{\theta }_{2}}{2}\right)\mathrm{exp}\left(-\beta {\tau }_{12}\right),\end{equation} \tag{ 11d }$$

$$\begin{align}\hfill & \alpha \left(t,{\tau }_{12}\right)=Vt-2\frac{V}{{\Gamma}}\left[{\mathrm{sin}}^{2}\left(\frac{{\theta }_{2}}{2}\right)+\mathrm{cos}\left({\theta }_{2}\right){\mathrm{sin}}^{2}\left(\frac{{\theta }_{1}}{2}\right){\mathrm{e}}^{-{\Gamma}{\tau }_{12}}\right]\left(1-{\mathrm{e}}^{-{\Gamma}t}\right),\hfill \end{align} \tag{ 11e }$$

$$\begin{equation}{\Phi}\left({\tau }_{12}\right)=-2\frac{V}{{\Gamma}}\enspace {\mathrm{sin}}^{2}\left(\frac{{\theta }_{1}}{2}\right)\left(1-{\mathrm{e}}^{-{\Gamma}{\tau }_{12}}\right)+V{\tau }_{12},\end{equation} \tag{ 11f }$$

$$\begin{equation}\eta \left(t,{\tau }_{12}\right)=-\frac{V}{{\Gamma}}\enspace \mathrm{sin}\left({\theta }_{1}\right)\mathrm{sin}\left({\theta }_{2}\right){\text{e}}^{-\beta {\tau }_{12}}\left(1-{\mathrm{e}}^{-{\Gamma}t}\right).\end{equation} \tag{ 11g }$$

Note, that the polarization is not oscillating with the transition frequency because we are working in the rotating frame of the laser field. Directly after the second pulse at t = 0, the polarization consists of three parts depending on the phases ϕ₁, ϕ₂ and 2ϕ₂ − ϕ₁. In the field-free propagation after the second pulse, the polarization receives an additional contribution from the occupation via the local field coupling. The corresponding term is the last exponential function in equation (11a), which also contributes to the phase dependence via the difference ϕ₂ − ϕ₁ and its complex conjugate phase. This will be discussed in detail below. From equation (11) we find that in the limit Γ → 0 an imaginary part of the local field parameter only gives rise to an additional damping, depending on the pulse areas and the delay time. Thus, we expect that EID manly gives rise to a faster decay of the signals and a broadening of the corresponding spectra. Since here we are mainly interested in characteristic features in the signals and the spectra in the following we will neglect an imaginary part and assume that V is a real parameter. In the appendix B we will show results of calculations including EID, which will confirm that this expectation is indeed fulfilled.

The FWM dynamics are obtained by filtering the polarization in equation (11a) with respect to the inverse FWM phase $-{\phi }_{\mathrm{F}\mathrm{W}\mathrm{M}}^{\left(2\right)}=-\left(2{\phi }_{2}-{\phi }_{1}\right)$ and we retrieve the two-pulse FWM polarization

$$\begin{equation}{p}_{\mathrm{F}\mathrm{W}\mathrm{M}}^{\left(2\right)}\left(t,{\tau }_{12}\right)=\underset{0}{\overset{2\pi }{\int }}{p}_{2}\left(t,{\tau }_{12}\right){\text{e}}^{-\text{i}\left(2{\phi }_{2}-{\phi }_{1}\right)}\frac{\mathrm{d}{\phi }_{1}\enspace \mathrm{d}{\phi }_{2}}{{\left(2\pi \right)}^{2}}.\end{equation} \tag{ 12 }$$

Hence, the FWM dynamics consist only of that part of the polarization which carries the FWM phase as all other phase-dependencies vanish. This procedure models the typical heterodyne detection in FWM experiments [11]. Without local field interaction there is no additional phase in the free propagation. Then the FWM polarization consists only of the term ∼c in equation (11a) which reads

$$\begin{equation}{p}_{\mathrm{F}\mathrm{W}\mathrm{M}}^{\left(2\right)}\left(t,{\tau }_{12}\right)=c\enspace {\mathrm{e}}^{-\beta t}\sim {\text{e}}^{-\beta \left(t+{\tau }_{12}\right)}.\end{equation} \tag{ 13 }$$

With local field interaction, the resulting FWM dynamics are

$$\begin{equation}{p}_{\mathrm{F}\mathrm{W}\mathrm{M}}^{\left(2\right)}\left(t,{\tau }_{12}\right)={\text{e}}^{\text{i}\left(\alpha -{\Phi}\right)-\beta t}\left[c{J}_{0}\left(\eta \right)+\mathrm{i}\enspace b{J}_{1}\left(\eta \right)-a{J}_{2}\left(\eta \right)\right].\end{equation} \tag{ 14 }$$

We find that each term of the polarization contributes to the FWM dynamics with a Bessel function of a different order. The full expression, after inserting equations (11b)–(11g) into equation (14), becomes quite involved, but we will now have a more detailed look on the creation of the FWM signal to understand the origins of its different contributions. A flow chart of the pulse sequence and phase selection resulting in the FWM signal is depicted in figure 2. Initially the system is in the ground state with n = p = 0. The first pulse creates an occupation ⁰ ${n}_{1}^{+}$ and a polarization with the phase of the first ${\mathrm{p}\mathrm{u}\mathrm{l}\mathrm{s}\mathrm{e}\enspace }^{{\phi }_{1}}{p}_{1}^{+}$, marked as green wave for a laser pulse induced step. For clarity, we give the phase dependence of each quantity on the left side of the symbol and the number of the pulse on the right side, where + is directly after and − directly before the respective pulse. During the field-free propagation before the second pulse two things happen. On the one hand both quantities simply propagate in time indicated by black arrows. On the other hand the occupation contributes to the phase through the local field coupling (Φ in equation (11f)), indicated by a violet dotted arrow. The second pulse creates an occupation consisting of a part independent of the pulse phases ⁰ ${n}_{2}^{+}$ and a part depending on the pulses' phase ${\text{difference}}_{{\hspace{2pt}\phi }_{1}-{\phi }_{2}}^{{\hspace{2pt}\phi }_{2}-{\phi }_{1}}{n}_{2}^{+}$. Because the occupation is real valued, this quantity has to depend on the phase difference and its conjugate, therefore two terms appear on the left side of the symbol. As mentioned above, the polarization created by the last pulse consists of parts depending on the phases of the exciting ${\text{pulses}}^{{\hspace{2pt}\phi }_{1}}{p}_{2}^{+}{,}^{{\hspace{2pt}\phi }_{2}}{p}_{2}^{+}$ and the FWM ${\text{phase}}^{\hspace{2pt}2{\phi }_{2}-{\phi }_{1}}{p}_{2}^{+}$. After the second ${\text{pulse,}}_{{\hspace{2pt}\phi }_{1}-{\phi }_{2}}^{\hspace{2pt}{\phi }_{2}-{\phi }_{1}}{n}_{2}^{+}$ contributes to the polarization via the local field resulting in additional wave-mixing possibilities (violet dotted arrows). We illustrate this by one example:

Zoom In Zoom Out Reset image size

The simplest wave-mixing process between the ${\text{occupation}\enspace }_{{\hspace{2pt}\phi }_{1}-{\hspace{2pt}\phi }_{2}}^{\hspace{2pt}{\phi }_{2}-{\phi }_{1}}{n}_{2}^{+}$ and the ${\text{polarization}}^{{\hspace{2pt}\phi }_{2}}{p}_{2}^{+}$ that results in the FWM phase is

$$\begin{equation}{\phi }_{\mathrm{F}\mathrm{W}\mathrm{M}}=\underset{p}{\underbrace{{\phi }_{2}}}+\underset{n}{\underbrace{\left({\phi }_{2}-{\phi }_{1}\right)}},\end{equation} \tag{ 15 }$$

where the occupation enters once. It is possible to generalize this to larger numbers of admixtures ${\text{of}\enspace }_{{\hspace{2pt}\phi }_{1}-{\phi }_{2}}^{\hspace{2pt}{\phi }_{2}-{\phi }_{1}}{n}_{2}^{+}$ giving

$$\begin{equation}{\phi }_{\mathrm{F}\mathrm{W}\mathrm{M}}={\phi }_{2}+k\left({\phi }_{2}-{\phi }_{1}\right)+\left(k-1\right)\left({\phi }_{1}-{\phi }_{2}\right),\end{equation} \tag{ 16 }$$

with $k\in \mathbb{N}$, where the occupation is used in total 2k − 1 times. In the flow chart this order of the phase difference mixing is annotated on the violet dotted lines. Considering the limit of weak local field coupling |V| ≪ β, i.e., |η| ≪ 1 the mixing orders enter as η^2k−1. This directly determines the order of the corresponding Bessel function, which is one in this case. Applying this procedure also to the other polarizations we find the Bessel functions of 0th, 1st, and 2nd order.

In the case of a negative delay, the pulse alignment is inverted as depicted in figure 3. The chronologically first (second) pulse has the pulse area θ₂ (θ₁) and the phase ϕ₂ (ϕ₁). To obtain the FWM signal for negative delays, equation (11) can be re-used where the pulse areas, pulse phases and the sign of the delay have to be switched θ₁ ↔ θ₂, ϕ₁ ↔ ϕ₂, τ₁₂ → τ₂₁ = −τ₁₂. The FWM polarization then reads

$$\begin{equation}{p}_{\mathrm{F}\mathrm{W}\mathrm{M}}^{\left(2\right)}\left(t,{\tau }_{12}{< }0\right)={\text{e}}^{\text{i}\left(\alpha +2{\Phi}\right)-\beta t}\left[\mathrm{i}\enspace a{J}_{1}\left(\eta \right)-b{J}_{2}\left(\eta \right)-\mathrm{i}\enspace c{J}_{3}\left(\eta \right)\right].\end{equation} \tag{ 17 }$$

In contrast to the FWM dynamics with positive delay, the orders of the Bessel functions are different as the wave-mixing orders ${\text{of}\enspace }_{\hspace{2pt}{\phi }_{1}-{\phi }_{2}}^{\hspace{2pt}{\phi }_{2}-{\phi }_{1}}{n}_{2}^{+}$ change. In particular, as both pulses change their ordering, in the pure TLS only the ${\text{polarization}}^{\hspace{2pt}2{\phi }_{1}-{\phi }_{2}}{p}_{2}^{+}$ is created, which results in a vanishing FWM signal. Therefore, to obtain a signal with the FWM phase 2ϕ₂ − ϕ₁ at least a third order mixing via the local field is necessary. Thus, the lowest order of wave-mixing is one, which is the contribution ${\text{of}}^{\hspace{2pt}{\phi }_{2}}{p}_{2}^{+}$ mixing with the phase difference ϕ₂ − ϕ₁ of the occupation once.

Zoom In Zoom Out Reset image size

To study FWM signals with non-vanishing pulse durations Δt > 0, the Bloch equations with local field coupling in equation (2) and the phase filtering are evaluated numerically [35]. The optical field is modeled as resonant Gaussian pulses centered at t = −τ₁₂ and t = 0,

$$\begin{align}\hfill {{\Omega}}_{\mathrm{e}\mathrm{x}\mathrm{t}\left(t\right)}& =\frac{1}{\sqrt{8\pi }{\Delta}t}\left\{{\theta }_{{\;}_{2}^{1}}\enspace \mathrm{exp}\left[-\frac{{\left(t+\vert {\tau }_{12}\vert \right)}^{2}}{2{\Delta}{t}^{2}}+\mathrm{i}\enspace {\phi }_{{\;}_{2}^{1}}\right]\right.\hfill \\ \hfill & \left.\quad +{\theta }_{{\;}_{1}^{2}}\enspace \mathrm{exp}\left(-\frac{{t}^{2}}{2{\Delta}{t}^{2}}+\mathrm{i}\enspace {\phi }_{{\;}_{1}^{2}}\right)\right\},\hfill \\ \hfill & \mathrm{f}\mathrm{o}\mathrm{r}\enspace {\tau }_{12}\gtrless 0.\hfill \end{align} \tag{ 18 }$$

To account for the inverted pulse ordering, for τ₁₂ > 0 the upper indices in ${\theta }_{{\;}_{2}^{1}}$ and ${\phi }_{{\;}_{2}^{1}}$ are used, while the lower ones refer to τ₁₂ < 0. For the following discussion we keep as many system parameters as possible fixed to keep the analysis as instructive as possible. For the exciton dephasing we choose β = 0.5 ps⁻¹ which is a typical value for hBN-encapsulated TMDC monolayer [36]. The exciton decay time is usually much longer than the dephasing [37] and we choose Γ = 0. For a non-vanishing value the impact of the local field effect would additionally decay in time in a non-trivial way. To avoid this complication we disregard the decay. For the pulse duration we assume Δt = 0.1 ps as a typical value for pulse trains in FWM experiments [8]. In addition we set θ₂ = 2θ₁ which is usually applied in FWM experiments [38]. As local field coupling we consider V = 6 ps⁻¹ which is an order of magnitude larger than the dephasing rate β. In the literature similar proportions have been used to describe semiconductor quantum well systems [24, 25]. In monolayer TMDCs, the dephasing and decay rates differ between samples. Additionally, the strength of the local field is generally not known and may depend on the ambient conditions. For simplicity, it is fixed to a value such that local field effects are clearly visible. Future comparisons to experiments will help to find the correct values for the local field strength.

We begin the discussion by analyzing the influence of the local field coupling on FWM signals in the regime of small pulse areas. Exemplarily, in figure 4, the FWM dynamics for θ₁ = 0.05π of the pure TLS with V = 0 and with a local field coupling V = 6 ps⁻¹ are shown in (a) and (b), respectively.

Zoom In Zoom Out Reset image size

In the pure TLS without local field interaction, the FWM dynamics in figure 4(a) directly represent the polarization dynamics which are given by exponential decays in real time t and delay time τ₁₂ with the dephasing rate β (equation (13)). In the direction of negative delays the dynamics drop sharply when going to τ₁₂ < 0 as there is no FWM polarization created when the pulse with phase ϕ₂ arrives before the pulse with phase ϕ₁ [20]. The corresponding FWM spectrum, defined as

$$\begin{equation}{S}_{\mathrm{F}\mathrm{W}\mathrm{M}}\left(\omega ,{\tau }_{12}\right)=\left\vert \enspace \underset{-\infty }{\overset{\infty }{\int }}{p}_{\mathrm{F}\mathrm{W}\mathrm{M}}^{\left(2\right)}\left(t,{\tau }_{12}\right){\text{e}}^{\text{i}\left(\omega -{\omega }_{0}\right)t}\enspace \mathrm{d}t\enspace \right\vert \end{equation} \tag{ 19 }$$

is a Lorentzian line for positive delays, i.e., a single peak at the bare transition energy ℏω = ℏω₀ (not shown).

Upon including the local field interaction, the FWM dynamics in figure 4(b) change qualitatively. On the one hand, the signal extends to τ₁₂ < 0 and on the other hand, the maximum of the signal is shifted to larger times t > 0. Both aspects can be understood with the analytic solution for δ-pulses. For low excitation powers, the lowest order contributing to the FWM signal is called the χ⁽³⁾-regime which for t ⩾ 0 reads [23]

$$\begin{align}\hfill {p}_{\mathrm{F}\mathrm{W}\mathrm{M}}^{\left(2\right)}\left(t,{\tau }_{12}\right)& \approx {\kappa }_{2}\enspace {\mathrm{e}}^{\mathrm{i}V\left(t-{\tau }_{12}\right)-\beta \left(t+\vert {\tau }_{12}\vert \right)}\hfill \\ \hfill & \quad {\times}\left\{{\Theta}\left({\tau }_{12}\right)+2\mathrm{i}\enspace Vt\left[{\Theta}\left({\tau }_{12}\right)+{\text{e}}^{\left(-\beta +\mathrm{i}V\right)\vert {\tau }_{12}\vert }{\Theta}\left(-{\tau }_{12}\right)\right]\right\},\hfill \end{align} \tag{ 20 }$$

with ${\kappa }_{2}=-\mathrm{i}\enspace {\theta }_{1}{\theta }_{2}^{2}/8$ and the Heaviside function Θ(t). For positive delays, the dynamics consist of a part decaying as ${\text{e}}^{-\beta \left(t+{\tau }_{12}\right)}$ which reflects simply the FWM dynamics without local field interaction and a part that rises for short times $\sim \hspace{-2pt}Vt\enspace {\text{e}}^{-\beta \left(t+\vert {\tau }_{12}\vert \right)}$. This contribution results from the mixing of the ${\text{polarization}}^{\hspace{2pt}{\phi }_{2}}{p}_{2}^{+}$ with the phase-dependent ${\text{occupation}}_{\hspace{2pt}{\phi }_{1}-{\phi }_{2}}^{\hspace{2pt}{\phi }_{2}-{\phi }_{1}}{n}_{2}^{+}$, i.e., mixing order 1 in figure 2. In contrast, for negative delays, only the last term contributes, as the FWM polarization is not created in the pure TLS but the ${\text{polarization}}^{{\hspace{2pt}\phi }_{2}}{p}_{2}^{+}$ can still mix with the occupation. Here, the pulse with the phase ϕ₂ arrives first, therefore ${\text{both}}^{{\hspace{2pt}\phi }_{2}}{p}_{2}^{+}$ ${\text{and}}_{{\hspace{2pt}\phi }_{1}-{\phi }_{2}}^{{\hspace{2pt}\phi }_{2}-{\phi }_{1}}{n}_{2}^{+}$ are reduced due to dephasing before the second pulse with ϕ₁. This results in the additional decay of the signal with ${\text{e}}^{-\beta \vert {\tau }_{12}\vert }$ in equation (20). So, in total, the signal decays twice as fast for negative delays than for positive ones. The FWM spectrum also consists of a single Lorentzian-like peak that is shifted by −ℏV due to the energy renormalization from the local field coupling (not shown).

For larger pulse areas beyond the χ⁽³⁾-regime, the FWM dynamics change remarkably. In figure 5(a) the FWM dynamics for θ₁ = 0.2π are shown. As a first striking difference, they exhibit oscillations, which stem from the Bessel functions known from the δ-pulse solution in equation (14). For increasing absolute delays, the minima of this modulation shift to larger t. Their functional evolution ${\tau }_{12}^{\left(\mathrm{m}\mathrm{i}\mathrm{n}\right)}\left(t\right)$ in the plot can be estimated from the equations. It is determined by the curve for which the argument η of the Bessel functions (equation (11g)) is constant, i.e., its differential dη vanishes:

$$\begin{align}\hfill \mathrm{d}\eta & ={\text{e}}^{-\beta {\tau }_{12}}\enspace \mathrm{d}t-\beta t\enspace {\text{e}}^{-\beta {\tau }_{12}}\mathrm{d}{\tau }_{12}!{\hspace{2pt}=\hspace{2pt}}0\hfill \\ \hfill & {\Rightarrow}\enspace \frac{\mathrm{d}{\tau }_{12}}{\mathrm{d}t}=\frac{1}{\beta t}\hfill \\ \hfill & {\Rightarrow}\enspace {\tau }_{12}^{\left(\mathrm{m}\mathrm{i}\mathrm{n}\right)}=\frac{1}{\beta }\enspace \mathrm{ln}\left(t\right)+{\tau }_{0}.\hfill \end{align} \tag{ 21 }$$

This curve is included as dashed line in figure 5(a). Therefore, the exponential loss of coherence between the pulses leads to the logarithm-shaped evolution of the minima.

Zoom In Zoom Out Reset image size

In figure 5(b), the corresponding FWM spectra are shown as a function of τ₁₂. While for small absolute delays |τ₁₂| ≪ 1/β the spectrum is broad, it becomes narrow for large delays, which can be explained by considering the FWM dynamics. Because the oscillations vanish for large τ₁₂, the single maximum at t > 0 leads to a single line in the spectrum. This behavior of the FWM spectra is similar for positive and negative delays and the spectral positions of large |τ₁₂| coincide. However for negative delays the signal is weaker and decays faster, as already seen in the χ⁽³⁾-regime in equation (20). The physical meaning of the spread and position of the spectrum will be discussed below.

Also the FWM dynamics for higher pulse areas of θ₁ = 0.4π, shown in figure 5(c), exhibit similarities with the previous example like the oscillatory behavior. However, one difference is that the signal decays much faster. In addition, when looking at small delays |τ₁₂| ≈ Δt, the minima deviate from the logarithmic shape. This indicates, that there is a difference between the action of two separate and two simultaneous pulses. The reason is that the pulse area and duration have a non-trivial impact on the resulting state when the system deviates from the pure TLS [38, 39]. Therefore, the pulse area theorem [18] does not hold for extended pulses anymore [39] and the pulse area θ does not agree with the rotation angle of the Bloch vector. This makes the dynamics of the signal more involved. Details about this pulse area renormalization are given in appendix A. In addition, during the pulse overlap complex dynamics stemming from Rabi oscillations contribute to the FWM signal for large enough pulse areas [40]. Qualitatively also the FWM spectra, shown in figure 5(d), have a similar form as the ones at lower pulse area. The most striking difference is however, that the final peak energy for long delays moves from ω < ω₀ to ω > ω₀. A discussion of the differences between the simulations and the δ-pulse limit is presented in appendix B, where the respective signals in the ultrashort pulse limit are shown. As expected, some deviations between the simulations with non-vanishing pulse duration and δ-pulses appear for delays that are shorter or of the order of the pulse duration, i.e., in the regime τ₁₂ ≲ 200 fs, where the overlap has a significant influence on the generated signal. For longer delays, the δ-pulse limit turns out to be a very accurate approximation.

In the following we analyze the spread and position of the FWM spectra at small and large delays in detail. A systematic study for small delays is shown in figure 6(a), where the spectra are depicted for a continuous variation of the pulse area θ₁. The delay is fixed to twice the pulse duration τ₁₂ = 2Δt, such that the two pulses are mainly separated. At low pulse areas the spectrum is very narrow as we have discussed for the χ⁽³⁾-regime. We have already seen in the exemplary FWM spectra that the spectral width increases with growing pulse area. This spread of the spectra can be traced back to the fact, that the transition energy is modified by the occupation (see equation (3)). So, in order to develop an intuitive picture which occupations can be reached in the FWM experiment, we consider the Bloch vector representation of the systems state in figure 6(b). The Bloch vector is given by its coordinates [2 Re(p), 2 Im(p), 2n − 1] and it points to the surface of the Bloch sphere for pure states and is shorter than unity when dephasing happens. Laser pulse excitations lead to rotations of the Bloch vector around an axis that depends on the phase of the pulse and, in general, also on the laser pulse detuning. However, here we assume a resonant excitation throughout this study.

Zoom In Zoom Out Reset image size

Keeping in mind that all possible combinations of pulse phases ϕ₁ and ϕ₂ that match the FWM phase are realized during the repetition of the experiment, we search for the smallest and largest possible occupation after the two pulses. Following equation (3) the corresponding renormalized energies can be calculated to ℏω_eff which connects the spread of the FWM spectrum with the extent to which the system can be addressed by its coherences. We can interpret this reasoning as a kind of coherent control experiment [41]. In particular, such a manipulation of the occupation by varying the phases of optical pulses is performed in Ramsey fringe experiments [42].

As schematically shown in the picture the largest occupation ${N}_{+}^{\left(2\right)}$ is reached when the two pulses rotate the Bloch vector into the same direction effectively adding up their pulse areas. The smallest occupation ${N}_{-}^{\left(2\right)}$ appears if the two pulses act in opposite directions. So for δ-pulses the maximal (+) and the minimal (−) occupation are given by

$$\begin{equation}{N}_{{\pm}}^{\left(2\right)}={\mathrm{sin}}^{2}\left(\frac{\vert {\theta }_{2}{\pm}{\theta }_{1}\vert }{2}\right).\end{equation} \tag{ 22 }$$

The corresponding energies $\hslash {\omega }_{\mathrm{e}\mathrm{ff}}\left({N}_{{\pm}}^{\left(2\right)}\right)$ are marked in figure 6(a) as dashed black lines. They describe the boundaries of the spectrum very well. But we find a slight mismatch between the dashed line and the edge of the spectrum. This deviation can be traced back to a renormalization of the pulse area due to the local field interaction, as we have mentioned before. While the pulse area is defined via its rotation angle in the pure TLS, the local field interaction changes this rotation angle in a non-trivial way (see appendix A).

The FWM spectra for large delays are depicted in figure 6(c). As we have already seen in figures 5(b) and (d), the spectrum becomes a single narrow line. Its energy depends on the pulse area and is the same for positive and negative delays. These findings can again be traced back to the behavior of the Bloch vector depicted in figure 6(d). After the first Rabi rotation the polarization dephases as indicated by a dotted line, i.e., the Bloch vector moves onto the Bloch sphere's z-axis. Afterwards, the second pulse leads to another Rabi rotation, which is reduced because the length of the Bloch vector is smaller than one. The spectral position of the FWM signal is then determined by the occupation after the second pulse ${N}_{\infty }^{\left(2\right)}$. For ultrafast pulses the final occupation reads

$$\begin{align}\hfill {N}_{\infty }^{\left(2\right)}& ={\mathrm{sin}}^{2}\left(\frac{{\theta }_{2}}{2}\right)+\mathrm{cos}\left({\theta }_{2}\right){\mathrm{sin}}^{2}\left(\frac{{\theta }_{1}}{2}\right)\hfill \\ \hfill & =\frac{1}{4}\left[2-\mathrm{cos}\left({\theta }_{2}-{\theta }_{1}\right)-\mathrm{cos}\left({\theta }_{2}+{\theta }_{1}\right)\right]\hfill \end{align} \tag{ 23 }$$

and is the same for both positive and negative delays as also shown in the schematic. This is the reason why the spectra lie at the same position for positive and negative delays. Because ${N}_{\infty }^{\left(2\right)}$ does not depend on the pulse phases anymore the spectrum consists of a single sharp line. In the context of coherent control this means that for delays much longer than the dephasing time the system cannot be addressed coherently anymore. Naturally, the loss of coherence also leads to a decay of the signal strength.

With three pulses there are two delays, τ₁₂ between the first and second pulse and τ₂₃ between the second and third pulse. Like in two-pulse FWM, the pulse areas θ_1,2,3 are kept fixed while the pulse phases ϕ_1,2,3 are scanned to generate the FWM signal. Usually, the delay τ₁₂ is kept small, such that the pulses slightly overlap, and it stays fixed while the delay τ₂₃ is varied.

As explained before in the limit of ultrashort pulses all properties of the system can be calculated analytically. Following the results in equation (11) and applying equation (8) we end up with the following phase dependencies immediately after the third pulse: the occupation ${n}_{3}^{+}$ carries the phases ϕ₂ − ϕ₁, ϕ₃ − ϕ₁, ϕ₃ − 2ϕ₂ + ϕ₁ and their complex conjugates, the polarization ${p}_{3}^{+}$ carries the phases ϕ₁, ϕ₂, ϕ₃, 2ϕ₂ − ϕ₁, 2ϕ₃ − ϕ₂, 2ϕ₃ − ϕ₁, and 2ϕ₃ − 2ϕ₂ + ϕ₁. In the following free propagation additional wave mixing between the occupation and the polarization can contribute to the FWM signal. As this large number of possible phase combinations complicates the situation massively, we will only analyze the limiting case of large delays τ₂₃ ≫ 1/β analytically. So we assume that before the interaction with the third pulse the polarization has dephased entirely and we can set ${p}_{3}^{-}=0$. Choosing further τ₁₂ = 0, the polarization after the three pulses reads

$$\begin{align}\hfill & {p}_{3}\left(t,{\tau }_{23}\right)=\mathrm{exp}\left\{\mathrm{i}\left[\vartheta +\nu \enspace \mathrm{cos}\left({\phi }_{2}-{\phi }_{1}\right)+{\phi }_{3}\right]-\beta t\right\}\left[u+2w\enspace \mathrm{cos}\left({\phi }_{2}-{\phi }_{1}\right)\right],\hfill \end{align} \tag{ 24a }$$

with

$$\begin{align}\hfill & \vartheta \left(t,{\tau }_{23}\right)=Vt-\frac{2V}{{\Gamma}}\left(1-{\mathrm{e}}^{-{\Gamma}t}\right)\left\{{\mathrm{sin}}^{2}\left(\frac{{\theta }_{3}}{2}\right)+{\mathrm{e}}^{-{\Gamma}{\tau }_{23}}\enspace \mathrm{cos}\left({\theta }_{3}\right)\left[{\mathrm{sin}}^{2}\left(\frac{{\theta }_{2}}{2}\right)+\mathrm{cos}\left({\theta }_{2}\right){\mathrm{sin}}^{2}\left(\frac{{\theta }_{1}}{2}\right)\right]\right\},\hfill \end{align} \tag{ 24b }$$

$$\begin{equation}\nu \left(t,{\tau }_{23}\right)=-\frac{V}{{\Gamma}}\enspace \mathrm{sin}\left({\theta }_{1}\right)\mathrm{sin}\left({\theta }_{2}\right)\mathrm{cos}\left({\theta }_{3}\right){\mathrm{e}}^{-{\Gamma}{\tau }_{23}}\left(1-{\mathrm{e}}^{-{\Gamma}t}\right),\end{equation} \tag{ 24c }$$

$$\begin{align}\hfill & u\left({\tau }_{23}\right)=\frac{\mathrm{i}}{2}\enspace \mathrm{sin}\left({\theta }_{3}\right)\left[1-2\enspace {\mathrm{sin}}^{2}\left(\frac{{\theta }_{2}}{2}\right)-2\enspace \mathrm{cos}\left({\theta }_{2}\right){\mathrm{sin}}^{2}\left(\frac{{\theta }_{1}}{2}\right)\right]{\mathrm{e}}^{-{\Gamma}{\tau }_{23}},\hfill \end{align} \tag{ 24d }$$

$$\begin{equation}w\left({\tau }_{23}\right)=-\frac{\mathrm{i}}{4}\enspace \mathrm{sin}\left({\theta }_{1}\right)\mathrm{sin}\left({\theta }_{2}\right)\mathrm{sin}\left({\theta }_{3}\right){\mathrm{e}}^{-{\Gamma}{\tau }_{23}}.\end{equation} \tag{ 24e }$$

Note that these equations are valid for a complex local field parameter V, i.e., EID is included. However, in the following we will again concentrate on a real parameter V.

Using the FWM phase ${\phi }_{\mathrm{F}\mathrm{W}\mathrm{M}}^{\left(3\right)}={\phi }_{3}+{\phi }_{2}-{\phi }_{1}$, the three-pulse FWM dynamics are obtained by the phase integration

$$\begin{equation}{p}_{\mathrm{F}\mathrm{W}\mathrm{M}}^{\left(3\right)}\left(t,{\tau }_{23}\right)=\underset{0}{\overset{2\pi }{\int }}{p}_{3}\left(t,{\tau }_{23}\right){\text{e}}^{-\text{i}\left({\phi }_{3}+{\phi }_{2}-{\phi }_{1}\right)}\frac{\mathrm{d}{\phi }_{1}\enspace \mathrm{d}{\phi }_{2}\enspace \mathrm{d}{\phi }_{3}}{{\left(2\pi \right)}^{3}}.\end{equation} \tag{ 25 }$$

A flow chart of the single steps leading to the FWM signal is shown in figure 7. Analogous to figure 2 the first and second pulse create an occupation which consists of a part independent of the pulse ${\text{phases}}^{\hspace{2pt}0}{n}_{2}^{+}$ and a part dependent on the pulses' phase ${\text{difference}}_{\hspace{2pt}{\phi }_{1}-{\phi }_{2}}^{{\hspace{2pt}\phi }_{2}-{\phi }_{1}}{n}_{2}^{+}$. These two terms remain even for long delays because we assume that the excited state decay is much slower than the dephasing. The polarization which consists ${\text{of}\enspace }^{{\hspace{2pt}\phi }_{1}}{p}_{2}^{+}$, ${\;}^{{\phi }_{2}}{p}_{2}^{+}$ and ${\;}^{2{\phi }_{2}-{\phi }_{1}}{p}_{2}^{+}$ vanishes due to the long delay τ₂₃ which is indicated by black dotted arrows.

Zoom In Zoom Out Reset image size

As a consequence the last pulse creates polarizations exclusively from the occupations before this pulse. Therefore, the final polarization consists of a part with the phase of the third ${\text{pulse}}^{{\hspace{2pt}\phi }_{3}}{p}_{3}^{+}$, indicated by the green wave ${\text{from}}^{\hspace{2pt}0}{n}_{2}^{+}$, and a ${\text{part}\enspace }_{{\hspace{2pt}\phi }_{3}-{\phi }_{2}+{\phi }_{1}}^{{\phi }_{3}+{\hspace{2pt}\phi }_{2}-{\phi }_{1}}{p}_{3}^{+}$ which adapts the phase difference from the ${\text{occupation}}_{\hspace{2pt}{\phi }_{1}-{\phi }_{2}}^{\hspace{2pt}{\phi }_{2}-{\phi }_{1}}{n}_{2}^{+}$. There is also the possibility that the occupation does not receive an additional phase dependence and simply ${\text{becomes}}_{\hspace{2pt}{\phi }_{1}-{\phi }_{2}}^{\hspace{2pt}{\phi }_{2}-{\phi }_{1}}{n}_{3}^{+}$.

Without the local field coupling, ${\text{only}}^{\hspace{2pt}{\phi }_{3}+{\phi }_{2}-{\phi }_{1}}{p}_{3}^{+}$ contributes to the FWM signal and we see that it carries the information of the previous occupation. Through the local field coupling, the polarization mixes ${\text{with}}_{\hspace{2pt}{\phi }_{1}-{\phi }_{2}}^{\hspace{2pt}{\phi }_{2}-{\phi }_{1}}{n}_{3}^{+}$ in the following field-free propagation which is marked by violet dotted arrows. This again results in Bessel functions describing the dynamics of the FWM signal. In direct analogy to the discussion of two-pulse FWM, the minimal number of mixing processes involved determines the order of the Bessel function. The possible mixing orders are again annotated on the violet arrows.

Finally, in case of long delays τ₂₃ ≫ 1/β, the FWM dynamics read

$$\begin{equation}{p}_{\mathrm{F}\mathrm{W}\mathrm{M}}^{\left(3\right)}={\text{e}}^{\text{i}\vartheta -\beta t}\left[w{J}_{0}\left(\nu \right)+\mathrm{i}\enspace u{J}_{1}\left(\nu \right)-w{J}_{2}\left(\nu \right)\right].\end{equation} \tag{ 26 }$$

The dynamics have a similar form as the analytical result from the two-pulse FWM in equation (14). This is because the occupation ${n}_{2}^{+}$ has the same phase dependence as in the two-pulse case and lowest mixing orders are also 0, 1, and 2 ${\text{for}}^{\hspace{2pt}{\phi }_{3}+{\phi }_{2}-{\phi }_{1}}{p}_{3}^{+}$, ${\;}^{{\phi }_{3}}{p}_{3}^{+}$, ${\text{and}}^{\hspace{2pt}{\phi }_{3}-{\phi }_{2}+{\phi }_{1}}{p}_{3}^{+}$, respectively.

To discuss non-vanishing pulse durations Δt > 0, we simulate the FWM signals numerically with Gaussian laser pulses. We set all pulse areas to the same value θ₁ = θ₂ = θ₃ and fix the delay between the first two pulses to τ₁₂ = 0.2 ps. The other parameters agree with the two-pulse case, i.e., β = 0.5 ps⁻¹, Γ = 0, and Δt = 0.1 ps.

Beginning with small pulse areas, the FWM dynamics for θ₁ = 0.05π without local field interaction V = 0 and with a local field coupling of V = 6 ps⁻¹ are shown in figures 8(a) and (b), respectively. For V = 0 in (a) the FWM signal decays in real time t as the FWM polarization dephases. The signal does not change when varying positive delays τ₂₃ because we did assume a vanishing decay rate and the signal is probing the occupation dynamics. For negative delays however, the pulse with ϕ₃ arrives first. After this pulse only the polarization is phase-dependent and the next two pulses with ϕ₁ and ϕ₂ create the FWM signal. Therefore the FWM signal reflects the polarization dynamics. This quantity decays $\sim \hspace{-2pt}{\text{e}}^{-\beta \vert {\tau }_{23}\vert }$ due to the dephasing in the field-free propagation after the single pulse.

Zoom In Zoom Out Reset image size

With the local field interaction in figure 8(b), the signal changes qualitatively. Similar to two-pulse FWM, the maximum of the FWM dynamics is shifted to t > 0 due to the wave-mixing in the free propagation after all three pulses. In addition the signal is strongest for short delays τ₂₃ ≪ 1/β. These aspects can be explained by calculating the FWM dynamics in the χ⁽³⁾-regime to

$$\begin{align}\hfill {p}_{\mathrm{F}\mathrm{W}\mathrm{M}}^{{\chi }^{\left(3\right)}}& ={\kappa }_{3}\enspace {\mathrm{e}}^{\mathrm{i}Vt-\beta t}\left\{{\Theta}\left({\tau }_{23}\right)\left[1+\mathrm{i}\enspace Vt\left(1+{\mathrm{e}}^{-2\beta {\tau }_{23}}\right)\right]+{\Theta}\left(-{\tau }_{23}\right){\text{e}}^{-\beta \vert {\tau }_{23}\vert +\mathrm{i}V\vert {\tau }_{23}\vert }\left(1+2\mathrm{i}\enspace Vt\right)\right\},\hfill \end{align} \tag{ 27 }$$

with κ₃ = −i θ₁ θ₂ θ₃/4. Independent of the sign of the delay, the energy is renormalized. Hence, the FWM spectrum consists of a single line at ω = ω₀ − V (not shown).

For positive delays the first term in the curly brackets in equation (27) contributes. There, the first contribution simply stems from the pure TLS without local field coupling. The second contribution is created by the interaction with the local field and consists of two parts. The first one does not depend on the delay and stems from the mixing between the polarization after the third ${\text{pulse}}^{\hspace{2pt}{\phi }_{3}}{p}_{3}^{+}$ and the ${\text{occupation}}^{\hspace{2pt}{\phi }_{2}-{\phi }_{1}}{n}_{3}^{+}$. The second part decays with ${\mathrm{e}}^{-2\beta {\tau }_{23}}$ and stems from the ${\text{polarization}}^{\hspace{2pt}{\phi }_{2}}{p}_{3}^{+}$ which mixes ${\text{with}}^{\hspace{2pt}{\phi }_{3}-{\phi }_{1}}{n}_{3}^{+}$. As both of these quantities stem from the polarization after pulses 1 and 2, they suffer from the dephasing in the field-free propagation after the pulse pair. Therefore this contribution decays with twice the dephasing rate.

For negative delays the second term in the curly brackets in equation (27) contributes. The pulse with phase ϕ₃ arrives first and after the delay τ₂₃, the two pulses with ϕ₁, ϕ₂ excite almost simultaneously. Consequently the only contributing quantity after the first excitation ${\text{is}}^{\hspace{2pt}{\phi }_{3}}{p}_{3}^{+}$ and the entire signal decays with ${\text{e}}^{-\beta \vert {\tau }_{23}\vert }$. The FWM dynamics consist of a contribution from the pure TLS and a part from the local field interaction. The latter consists of the mixing ${\text{between}}^{\hspace{2pt}{\phi }_{3}}{p}_{3}^{+}$ ${\text{and}}^{\hspace{2pt}{\phi }_{2}-{\phi }_{1}}{n}_{3}^{+}$ and ${\text{between}}^{\hspace{2pt}{\phi }_{2}}{p}_{3}^{+}$ ${\text{and}}^{\hspace{2pt}{\phi }_{3}-{\phi }_{1}}{n}_{3}^{+}$ resulting in the term ∼Vt.

Going to larger pulse areas, the FWM spectra first change slightly. For an intermediate pulse area of θ₁ = 0.25π in figure 9(a) the FWM spectrum is initially significantly broadened and becomes narrower for larger delays but it does not evolve into a single symmetric line and remains broader than the Lorentzian in the χ⁽³⁾-regime. A detailed discussion of the linewidth at large delays will be carried out later.

Zoom In Zoom Out Reset image size

At still higher pulse areas, e.g., θ₁ = 0.45π depicted in figure 9(b), the FWM spectra change significantly. Initially, the FWM spectrum is broad and evolves into two separate lines for large τ₂₃ which is a striking difference to the previous picture. Both lines are separated by a minimum of the signal at ω = ω₀. This characteristic behavior will be explained in the context of coherent control later.

For negative delays τ₂₃ < 0 the details of the FWM spectra change significantly when comparing figures 8(a) and (b) and become quite involved. However, in all cases the intensity decays rapidly due to the dephasing as discussed for the χ⁽³⁾-regime.

While the spectra are all broadened and quite involved for small delays τ₂₃ ≈ 1 ps they become very clean for long delays τ₂₃ ≈ 5 ps. The reason for this behavior is the multitude of different possible phase combinations for short delays and the vanishing influence of the polarizations for long delays as explained before. We find pulse areas that result in a single or in two separate lines. This effect is further highlighted in figure 10(a), where the FWM spectra for large delays τ₂₃ ≫ 1/β are shown as a function of the pulse area θ. For small pulse areas in the χ⁽³⁾-regime, the spectrum appears as a single line at ω = ω₀ − V. Then the spectrum broadens for an increasing pulse area. At θ₁ ≈ π/4, the spectrum splits into two parts with a pronounced minimum at ω = ω₀ in agreement with figure 9(b). For θ₁ ≈ π/2, the two lines merge again into a single peak at ω = ω₀. Increasing the pulse area further, the single line splits into two lines again.

Zoom In Zoom Out Reset image size

For the explanation of this behavior we consider the Bloch vector picture again. Figure 10(b) shows the schematic construction of the smallest and largest possible occupation after the three pulses and a full dephasing after pulse two. The first Rabi rotations (blue and orange) can act destructively, i.e., the respective rotations of the Bloch vector compensate each other. The state of the TLS after the second pulse then points to the south pole. The third pulse rotates from there and reaches the final occupation.

$$\begin{equation}{N}_{-}^{\left(3\right)}={\mathrm{sin}}^{2}\left(\frac{{\theta }_{3}}{2}\right),\end{equation} \tag{ 28a }$$

which is the smallest possible. Note, that strictly speaking this phase combination leads to a vanishing FWM signal because both, polarization and occupation are zero after the second pulse. But at this point we are just interested in the hypothetically reached shift of the transition energy. In the other extreme case the first two pulses act constructively and the pulse areas essentially add up. Before the third pulse the state dephases entirely, which means that it is projected onto the z-axis. From there the third pulse rotates the Bloch vector on a smaller circle and reaches the final occupation

$$\begin{equation}{N}_{+}^{\left(3\right)}={\mathrm{sin}}^{2}\left(\frac{{\theta }_{3}}{2}\right)+\mathrm{cos}\left({\theta }_{3}\right){\mathrm{sin}}^{2}\left(\frac{{\theta }_{1}+{\theta }_{2}}{2}\right).\end{equation} \tag{ 28b }$$

The effective transition energies corresponding to these two limiting cases are depicted in figure 10(a) as dashed black lines. We find that they again very accurately follow the boundaries of the calculated spectra. The slight deviation again stems from the pulse area renormalization due to the local field coupling (see appendix A).

Also the reduction of the signal at ω = ω₀ (seen most clearly for pulse areas around θ₁ = 0.4π) can be explained in this picture. We consider all possible combinations of ϕ₁ and ϕ₂ that make the Bloch vector point to the equator after the second pulse. The following dephasing ends in the center of the sphere, which is the balanced mixture of ground and excited state. In this point the system is called transparent because the third pulse cannot create a polarization from there. Consequently also no FWM signal can be generated.

For the special case θ₁ = 0.5π, where the spectral peaks in figure 10(a) merge at ω = ω₀, we remember that after the full dephasing for long delays τ₂₃ all possible Bloch vectors lie on the z-axis of the Bloch sphere. From there the Bloch vector is rotated by the third pulse into the equatorial plane no matter where on the z-axis it was. Therefore the only possible occupation is ${n}_{3}^{+}=1/2$ which results in a renormalized energy of ℏω_loc = 0 (equation (3)), i.e., a single line in the spectrum at ω = ω₀. Note, that in this scenario the equator is reached after the third and not the second pulse, so the argument is not in conflict with the previous discussion of the spectral minimum at ω = ω₀.