Ultrafast spectroscopy is often collected in the mixed frequency/time domain, where pulse durations are similar to system dephasing times. In these experiments, expectations derived from the familiar driven and impulsive limits are not valid. This work simulates the mixed-domain four-wave mixing response of a model system to develop expectations for this more complex field-matter interaction. We explore frequency and delay axes. We show that these line shapes are exquisitely sensitive to excitation pulse widths and delays. Near pulse overlap, the excitation pulses induce correlations that resemble signatures of dynamic inhomogeneity. We describe these line shapes using an intuitive picture that connects to familiar field-matter expressions. We develop strategies for distinguishing pulse-induced correlations from true system inhomogeneity. These simulations provide a foundation for interpretation of ultrafast experiments in the mixed domain.

## I. INTRODUCTION

Ultrafast spectroscopy is based on using nonlinear interactions, created by multiple ultrashort (10^{−9}–10^{−15} s) pulses, to resolve spectral information on time scales as short as the pulses themselves.^{1,2} The ultrafast spectra can be collected in the time domain or in the frequency domain.^{3}

Time-domain methods scan the pulse delays to resolve the free induction decay (FID).^{4} The Fourier transform of the FID gives the ultrafast spectrum. Ideally, these experiments are performed in the impulsive limit where FID dominates the measurement. FID occurs at the frequency of the transition that has been excited by a well-defined, time-ordered sequence of pulses. Time-domain methods are compromised when the dynamics occur on faster time scales than the ultrafast excitation pulses. As the pulses temporally overlap, FID from other pulse time-orderings and emission driven by the excitation pulses both become important. These factors are responsible for the complex “coherent artifacts” that are often ignored in pump-probe and related methods.^{5–8} Dynamics faster than the pulse envelopes are best measured using line shapes in frequency domain methods.

Frequency-domain methods scan pulse frequencies to resolve the ultrafast spectrum directly.^{9,10} Ideally, these experiments are performed in the driven limit where the steady state dominates the measurement. In the driven limit, all time-orderings of the pulse interactions are equally important and FID decay is negligible. The output signal is driven at the excitation pulse frequencies during the excitation pulse width. Frequency-domain methods are compromised when the spectral line shape is narrower than the frequency bandwidth of the excitation pulses. Dynamics that are slower than the pulse envelopes can be measured in the time domain by resolving the phase oscillations of the output signal during the entire FID decay.

There is also the hybrid mixed-time/frequency-domain approach, where pulse delays and pulse frequencies are both scanned to measure the system response. This approach is uniquely suited for experiments where the dephasing time is comparable to the pulse durations so that neither frequency-domain nor time-domain approaches excel on their own.^{10–12} In this regime, both FID and driven processes are important.^{13} Their relative importance depends on pulse frequencies and delays. Extracting the correct spectrum from the measurement then requires a more complex analysis that explicitly treats the excitation pulses and the different time-orderings.^{14–16} Despite these complications, mixed-domain methods have a practical advantage: the dual frequency- and delay-scanning capabilities allow these methods to address a wide variety of dephasing rates.

The relative importance of FID and driven processes and the changing importance of different coherence pathways are important factors for understanding spectral features in all ultrafast methods. These methods include partially coherent methods involving intermediate populations such as pump-probe spectroscopy,^{17} transient grating spectroscopy,^{18–20} transient absorption/reflection spectroscopy,^{21,22} photon echo spectroscopy,^{23–25} two dimensional-infrared (2D-IR) spectroscopy,^{26–28} 2D-electronic spectroscopy (2D-ES),^{29,30} and three-pulse photon echo peak shift (3PEPS)^{23,31–34} spectroscopy. These methods also include fully coherent methods involving only coherences such as Stimulated Raman Spectroscopy (SRS),^{35,36} Doubly Vibrationally Enhanced (DOVE),^{37–43} Triply Resonant Sum Frequency (TRSF),^{44–46} Sum Frequency Generation (SFG),^{47} Coherent Anti-Stokes Raman Spectroscopy (CARS),^{48–50} and other coherent Raman methods.^{51}

This paper focuses on understanding the nature of the spectral changes that occur in Coherent Multidimensional Spectroscopy (CMDS) as experiments transition between the two limits of frequency- and time-domain methods. CMDS is a family of spectroscopies that use multiple delay and/or frequency axes to extract homogeneous and inhomogeneous broadening, as well as detailed information about spectral diffusion and chemical changes.^{52,53} For time-domain CMDS (2D-IR, 2D-ES), the complications that occur when the impulsive approximation does not strictly hold have only recently been addressed.^{54,55}

Frequency-domain CMDS methods, referred to herein as multi-resonant CMDS (MR-CMDS), have similar capabilities for measuring homogeneous and inhomogeneous broadening. Although these experiments are typically described in the driven limit,^{4,19,20} many of the experiments involve pulse widths that are comparable to the widths of the system.^{15,37,40,41,56,57} MR-CMDS then becomes a mixed-domain experiment whereby resonances are characterized with marginal resolution in both frequency and time. For example, DOVE spectroscopy involves three different pathways^{58} whose relative importance depends on the relative importance of FID and driven responses.^{59} In the driven limit, the DOVE line shape depends on the difference between the first two pulse frequencies, so the line shape has a diagonal character that mimics the effects of inhomogeneous broadening. In the FID limit where the coherence frequencies are defined instead by the transition, the diagonal character is lost. Understanding these effects is crucial for interpreting experiments, yet these effects have not been characterized for MR-CMDS.

This work considers the third-order MR-CMDS response of a 3-level model system using three ultrafast excitation beams with the commonly used four-wave mixing (FWM) phase-matching condition, $k\u2192out\u2009=\u2009k\u21921\u2212k\u21922+k\u21922\u2032$. Here, the subscripts represent the excitation pulse frequencies, $\omega 1$ and $\omega 2\u2009=\u2009\omega 2\u2032$. These experimental conditions were recently used to explore line shapes of excitonic systems^{15,57} and have been developed on vibrational states as well.^{60} Although MR-CMDS forms the context of this model, the treatment is quite general because the phase-matching condition can describe any of the spectroscopies mentioned above with the exception of SFG and TRSF, for which the model can be easily extended. We numerically simulate the MR-CMDS response with pulse durations at, above, and below the system coherence time. To highlight the role of pulse effects, we build an interpretation of the full MR-CMDS response by first showing how finite pulses affect the evolution of a coherence and then how finite pulses affect an isolated third-order pathway. When considering the full MR-CMDS response, we show that spectral features change dramatically as a function of delay, even for a homogeneous system with elementary dynamics. Importantly, the line shape can exhibit correlations that mimic inhomogeneity, and the temporal evolution of this line shape can mimic spectral diffusion. We identify key signatures that can help differentiate true inhomogeneity and spectral diffusion from these measurement artifacts.

## II. THEORY

### A. Third-order response with finite pulses

We consider a simple three-level system (states *n* = 0, 1, 2) that highlights the multidimensional line shape changes resulting from choices of the relative dephasing and detuning of the system and the temporal and spectral widths of the excitation pulses. For simplicity, we will ignore population relaxation effects: $\Gamma 11\u2009=\u2009\Gamma 00\u2009=\u20090$.

The electric field pulses, $El$, are given by

where $\omega l$ is the field carrier frequency, $k\u2192l$ is the wavevector, $\tau l$ is the pulse delay, and *c*_{l} is a slowly varying envelope. In this work, we assume normalized (real-valued) Gaussian envelopes

where $\Delta t$ is the temporal FWHM of the envelope intensity. We neglect non-linear phase effects such as chirp, so the FWHM of the frequency bandwidth is transform limited: $\Delta \omega \Delta t\u2009=\u20094ln2\u22482.77$, where $\Delta \omega $ is the spectral FWHM (intensity scale).

The Liouville-von Neumann equation propagates the density matrix, $\rho $,

Here *H*_{0} is the time-independent Hamiltonian, $\mu $ is the dipole superoperator, and $\Gamma $ contains the pure dephasing rate of the system. We perform the standard perturbative expansion of Eq. (3) to third order in the electric field interaction^{10,61–64} and restrict ourselves only to the terms that have the correct spatial wave vector $k\u2192out\u2009=\u2009k\u21921\u2212k\u21922+k\u21922\u2032$. This approximation narrows the scope to sets of three interactions, one from each field, that result in the correct spatial dependence. The set of three interactions have 3! = 6 unique time-ordered sequences, and each time-ordering produces either two or three unique system-field interactions for our system, for a total of sixteen unique system-field interaction sequences, or Liouville pathways, to consider. Figure 1 shows these sixteen pathways as Wave Mixing Energy Level (WMEL) diagrams.^{65}

We first focus on a single interaction in these sequences, where an excitation pulse, *x*, forms $\rho ij$ from $\rho ik$ or $\rho kj$. For brevity, we use $\u210f\u2009=\u20091$ and abbreviate the initial and final density matrix elements as $\rho i$ and $\rho f$, respectively. Using the natural frequency rotating frame, $\rho \u0303ij\u2009=\u2009\rho ije\u2212i\omega ijt$, the formation of $\rho f$ using pulse *x* is written as

where $\Omega fx\u2009=\u2009\kappa f\u22121\omega f\u2212\omega x(\u2009=\u2009|\omega f|\u2212\omega x)$ is the detuning, $\omega f$ is the transition frequency of the *i*th transition, $\mu f$ is the transition dipole, and $\Gamma f$ is the dephasing/relaxation rate for $\rho f$. The $\lambda f$ and $\kappa f$ parameters describe the phases of the interaction: $\lambda f\u2009=\u2009+1$ for ket-side transitions and −1 for bra-side transitions, and $\kappa f$ depends on whether $\rho f$ is formed via absorption ($\kappa f\u2009=\u2009\lambda f$) or emission ($\kappa f\u2009=\u2009\u2212\lambda f$).^{66} In the following equations, we neglect spatial dependence (*z* = 0).

Equation (4) forms the basis for our simulations. It provides a general expression for arbitrary values of the dephasing rate and excitation pulse bandwidth. The integral solution is

where $\Theta $ is the Heaviside step function. Equation (5) becomes the steady state limit expression when $\Delta t|\Gamma f+i\kappa f\Omega fx|\u2009>>\u20091$, and the impulsive limit expression results when $\Delta t|\Gamma f+i\kappa f\Omega fx|\u2009<<\u20091$. Both limits are important for understanding the multidimensional line shape changes discussed in this paper. The steady state and impulsive limits of Eq. (5) are discussed in Appendix A.

Figure 2 gives an overview of the simulations done in this work. Figure 2(a) shows an excitation pulse (gray-shaded) and examples of a coherent transient for three different dephasing rates. The color bindings to dephasing rates introduced in Fig. 2(a) will be used consistently throughout this work. Our simulations use systems with dephasing rates quantified relative to the pulse duration: $\Gamma 10\Delta t\u2009=\u20090.5,1$, or 2. The temporal axes are normalized to the pulse duration, $\Delta t$. The $\Gamma 10\Delta t\u2009=\u20092$ transient is mostly driven by the excitation pulse, while $\Gamma 10\Delta t\u2009=\u20090.5$ has a substantial free induction decay (FID) component at late times. Figure 2(b) shows a pulse sequence (pulses are shaded orange and pink) and the resulting system evolution of pathway $V\gamma $ ($00\u2192201\u21922\u203211\u2192110\u2192out00$) with $\Gamma 10\Delta t\u2009=\u20091$. The final polarization (teal) is responsible for the emitted signal, which is then passed through a frequency bandpass filter to emulate monochromator detection [Fig. 2(c)]. The resulting signal is explored in 2D delay space [Fig. 2(d)], 2D frequency space [Fig. 2(f)], and hybrid delay-frequency space [Fig. 2(e)]. The detuning frequency axes are also normalized by the pulse bandwidth, $\Delta \omega $.

We now consider the generalized Liouville pathway $L:\rho 0\u2192x\rho 1\u2192y\rho 2\u2192z\rho 3\u2192out\rho 4$, where *x*, *y*, and *z* denote properties of the first, second, and third pulses, respectively, and indices 0, 1, 2, 3, and 4 define the properties of the ground state, first, second, third, and fourth density matrix elements, respectively. Figure 2(b) demonstrates the correspondence between *x*, *y*, *z* notation and 1, 2, $2\u2032$ notation for the laser pulses with pathway $V\gamma $.^{67}

The electric field emitted from a Liouville pathway is proportional to the polarization created by the third-order coherence,

Equation (6) assumes perfect phase-matching and no pulse distortions through propagation. Equation (5) shows that the output field for this Liouville pathway is

where *R*_{L} is the third-order response function for the Liouville pathway. The total electric field will be the superposition of all the Liouville pathways,

For the superposition of Eq. (9) to be non-canceling, certain symmetries between the pathways must be broken. In general, this requires one or more of the following inequalities: $\Gamma 10\u2260\Gamma 21$, $\omega 10\u2260\omega 21$, and/or $2\mu 10\u2260\mu 21$. Our simulations use the last inequality, which is important in two-level systems ($\mu 21\u2009=\u20090$) and in systems where state-filling dominates the non-linear response, such as in semiconductor excitons. The exact ratio between $\mu 10$ and $\mu 21$ affects the absolute amplitude of the field but does not affect the multidimensional line shape. Importantly, the dipole inequality does not break the symmetry of double quantum coherence pathways (time-orderings II and IV), so such pathways are not present in our analysis.

In MR-CMDS, a monochromator resolves the driven output frequency from other nonlinear output frequencies, which in our case is $\omega m\u2009=\u2009\omega 1\u2212\omega 2+\omega 2\u2032\u2009=\u2009\omega 1$. The monochromator can also enhance spectral resolution, as we show in Sec. IV A. In this simulation, the detection is emulated by transforming $Etot(t)$ into the frequency domain, applying a narrow bandpass filter, $M(\omega )$, about $\omega 1$, and applying amplitude-scaled detection,

where $Etot(\omega )$ denotes the Fourier transform of $Etot(t)$ [see Fig. 2(c)]. For *M* we used a rectangular function of width $0.408\Delta \omega $. The arguments of $Stot$ refer to the *experimental* degrees of freedom. The signal delay dependence is parameterized with the relative delays $\tau 21$ and $\tau 22\u2032$, where $\tau nm\u2009=\u2009\tau n\u2212\tau m$ [see Fig. 2(b)]. Table S1 of the supplementary material summarizes the arguments for each Liouville pathway. Figure 2(f) shows the 2D $(\omega 1,\omega 2)$$Stot$ spectrum resulting from the pulse delay times represented in Fig. 2(b).

### B. Inhomogeneity

Inhomogeneity is isolated in CMDS through both spectral signatures, such as line-narrowing,^{10,51,68–70} and temporal signatures, such as photon echoes.^{71,72} We simulate the effects of static inhomogeneous broadening by convolving the homogeneous response with a Gaussian distribution function. Further details of the convolution are in Appendix B. Dynamic broadening effects such as spectral diffusion are beyond the scope of this work.

## III. METHODS

A matrix representation of differential equations of the type in Eq. (7) was numerically integrated for parallel computation of Liouville elements (see the supplementary material for details).^{73,74} The lower bound of integration was $2\Delta t$ before the first pulse, and the upper bound was $5\Gamma 10\u22121$ after the last pulse, with step sizes much shorter than the pulse durations. Integration was performed in the FID rotating frame; the time steps were chosen so that both the system-pulse difference frequencies and the pulse envelope were well-sampled.

The following simulations explore the four-dimensional $(\omega 1,\omega 2,\tau 21,\tau 22\u2032)$ variable space. Both frequencies are scanned about the resonance, and both delays are scanned about pulse overlap. We explored the role of the sample dephasing rate by calculating the signal for systems with dephasing rates such that $\Gamma 10\Delta t\u2009=\u20090.5,\u20091,$ and 2. Inhomogeneous broadening used a spectral FWHM, $\Delta inhom$, that satisfied $\Delta inhom/\Delta \omega \u2009=\u20090,0.5,1,$ and 2 for the three dephasing rates. For all these dimensions, both $\rho 3(t;\omega 1,\omega 2,\tau 21,\tau 22\u2032)$ and $Stot(\omega 1,\omega 2,\tau 21,\tau 22\u2032)$ are recorded for each unique Liouville pathway. Our simulations were done using the open-source SciPy library.^{75}

## IV. RESULTS

We now present portions of our simulated data that highlight the dependence of the spectral line shapes and transients on excitation pulse width, dephasing rate, detuning from resonance, pulse delay times, and inhomogeneous broadening.

### A. Evolution of single coherence

It is illustrative to first consider the evolution of single coherences, $\rho 0\u2192x\rho 1$, under various excitation conditions. Figure 3 shows the temporal evolution of $\rho 1$ with various dephasing rates under Gaussian excitation. The value of $\rho 1$ differs only by phase factors between various Liouville pathways [this can be verified by inspection of Eq. (5) under the various conditions in Table S1 of the supplementary material], so the profiles in Fig. 3 apply for the first interaction of any pathway. The pulse frequency was detuned from resonance so that frequency changes could be visualized by the color bar, but the detuning was kept slight so that it did not appreciably change the dimensionless product, $\Delta t\Gamma f+i\kappa f\Omega fx\u2248\Gamma 10\Delta t$. In this case, the evolution demonstrates the maximum impulsive character the transient can achieve. The instantaneous frequency, $d\phi /dt$, is defined as

The cases of $\Gamma 10\Delta t\u2009=\u20090(\u221e)$ agree with the impulsive (driven) expressions derived in Appendix A. For $\Gamma 10\Delta t\u2009=\u20090$, the signal rises as the integral of the pulse and has instantaneous frequency close to that of the pulse [Eq. (A3)], but as the pulse vanishes, the signal adopts the natural system frequency and decay rate [Eq. (A2)]. For $\Gamma 10\Delta t\u2009=\u2009\u221e$, the signal follows the amplitude and frequency of the pulse for all times [the driven limit, Eq. (A1)].

The other three cases show a smooth interpolation between limits. As $\Gamma 10\Delta t$ increases from the impulsive limit, the coherence within the pulse region conforms less to a pulse integral profile and more to a pulse envelope profile. Accordingly, the FID component after the pulse becomes less prominent, and the instantaneous frequency pins to the driving frequency more strongly through the course of evolution. The trends can be understood by considering the differential form of evolution [Eq. (4)] and the time-dependent balance of optical coupling and system relaxation. We note that our choices of $\Gamma 10\Delta t\u2009=\u20092.0,1.0,$ and 0.5 give coherences that have mainly driven, roughly equal driven and FID parts, and mainly FID components, respectively. The FID character is difficult to isolate when $\Gamma 10\Delta t\u2009=\u20092.0$.

Figure 4(a) shows the temporal evolution of $\rho 1$ at several values of $\Omega 1x/\Delta \omega $ with $\Gamma 10\Delta t\u2009=\u20091$ (supplementary material Fig. S3 shows the Fourier domain representation of Fig. 4(a)). As detuning increases, total amplitude decreases, FID character vanishes, and $\rho 1$ assumes a more driven character, as expected. During the excitation, $\rho 1$ maintains a phase relationship with the input field [as seen by the instantaneous frequency in Fig. 4(a)]. The coherence will persist beyond the pulse duration only if the pulse transfers energy into the system; FID evolution equates to absorption. The FID is therefore sensitive to the absorptive (imaginary) line shape of a transition, while the driven response is the composite of both absorptive and dispersive components. If the experiment isolates the latent FID response, there is consequently a narrower spectral response. This spectral narrowing can be seen in Fig. 4(a) by comparing the coherence amplitudes at *t* = 0 (driven) and at $t/\Delta t\u2009=\u20092$ (FID); amplitudes for all $\Omega fx/\Delta \omega $ values shown are comparable at *t* = 0, but the lack of FID formation for $\Omega fx/\Delta \omega \u2009=\u2009\xb12$ manifests as a visibly disproportionate amplitude decay (supplementary materials’ Fig. S4 shows explicit plots of $\rho 1(\Omega fx/\Delta \omega )$ at discrete $t/\Delta t$ values). Many ultrafast spectroscopies take advantage of the latent FID to suppress non-resonant background, improving signal to noise.^{42,47,59,76}

In driven experiments, the output frequency and line shape are fully constrained by the excitation beams. In such experiments, there is no additional information to be resolved in the output spectrum. The situation changes in the mixed domain, where $Etot$ contains the FID signal that lasts longer than the pulse duration. Figure 4(a) provides insight on how frequency-resolved detection of coherent output can enhance resolution when pulses are spectrally broad. Without frequency-resolved detection, mixed-domain resonance enhancement occurs in two ways: (1) the peak amplitude increases, and (2) the coherence duration increases due to the FID transient. Frequency-resolved detection can further discriminate against detuning by requiring that the driving frequency agrees with latent FID. The implications of discrimination are most easily seen in Fig. 4(a) with $\Omega 1x/\Delta \omega \u2009=\u2009\xb11$, where the system frequency moves from the driving frequency to the FID frequency. When the excitation pulse frequency is scanned, the resonance will be more sensitive to detuning by isolating the driven frequency (tracking the monochromator with the excitation source).

The functional form of the measured line shape can be deduced by considering the frequency domain form of Eq. (5) (assume $\rho i\u2009=\u20091$ and $\tau x\u2009=\u20090$),

where $FCx\omega $ is the Fourier transform of $cx$,

For squared-law detection of $\rho f$, the importance of the tracking monochromator is highlighted by two limits of Eq. (12) as follows:

When the transient is not frequency resolved, $sig\u2248\u222b|\rho \u0303f(\omega )|2d\omega $ and the measured line shape will be the convolution of the pulse envelope and the intrinsic (Green’s function) response [Fig. 4(b), magenta].

When the driven frequency is isolated, $sig\u2248|\rho \u0303f(\kappa f\Omega fx)|2$ and the measured line shape will give the un-broadened Green’s function [Fig. 4(b), teal].

Monochromatic detection can remove broadening effects due to the pulse bandwidth. For large $\Gamma 10\Delta t$ values, FID evolution is negligible at all $\Omega fx/\Delta \omega $ values and the monochromator is not useful. Figure 4(b) shows the various detection methods for the relative dephasing rate of $\Gamma 10\Delta t\u2009=\u20091$.

### B. Evolution of single Liouville pathway

We now consider the multidimensional response of a single Liouville pathway involving three pulse interactions. In a multi-pulse experiment, $\rho 1$ acts as a source term for $\rho 2$ (and subsequent excitations). The spectral and temporal features of $\rho 1$ that are transferred to $\rho 2$ depend on when the subsequent pulse arrives. Time-gating later in $\rho 1$ evolution will produce responses with FID behavior, while time-gating $\rho 1$ in the presence of the initial pulse will produce driven responses. An analogous relationship holds for $\rho 3$ with its source term $\rho 2$. As discussed above, the signal that time-gates FID evolution gives narrower spectra than the driven-gated signal. As a result, the spectra of even single Liouville pathways will change based on pulse delays.

The final coherence will also be frequency-gated by the monochromator. The monochromator isolates the signal at the fully driven frequency $\omega out\u2009=\u2009\omega 1$. The monochromator will induce line-narrowing to the extent that FID takes place. It effectively enforces a frequency constraint that acts as an additional resonance condition, $\omega out\u2009=\u2009\omega 1$. The driven frequency will be $\omega 1$ if *E*_{1} is the last pulse interaction (time-orderings V and VI), and the monochromator tracks the coherence frequency effectively. If *E*_{1} is not the last interaction, the output frequency may not be equal to the driven frequency, and the monochromator plays a more complex role.

We demonstrate this delay dependence using the multidimensional response of the I$\gamma $ Liouville pathway as an example (see Fig. 1). Figure 5 shows the resulting 2D delay profile of pathway I$\gamma $ signals for $\Gamma 10\Delta t\u2009=\u20091$ (left) and the corresponding $\omega 1,\omega 2$ 2D spectra at several pulse delay values (right). The spectral changes result from changes in the relative importance of driven and FID components. The prominence of the FID signal can change the resonance conditions; Table I summarizes the changing resonance conditions for each of the four delay coordinates studied. Since *E*_{1} is not the last pulse in pathway I$\gamma $, the tracking monochromator must also be considered. Supplementary materials Fig. S5 shows a simulation of Fig. 5 without monochromator frequency filtering ($M(\omega \u2212\omega 1)=1$ in Eq. 10).

Delay . | Approximate resonance conditions . | ||||
---|---|---|---|---|---|

$\tau 21/\Delta t$ . | $\tau 22\u2032/\Delta t$ . | $\rho 0\u21921\rho 1$ . | $\rho 1\u21922\rho 2$ . | $\rho 2\u21922\u2032\rho 3$ . | $\rho 3\u2192$ detection at $\omega m=\omega 1$ . |

0 | 0 | $\omega 1=\omega 10$ | $\omega 1=\omega 2$ | $\omega 1=\omega 10$ | … |

0 | −2.4 | $\omega 1=\omega 10$ | $\omega 1=\omega 2$ | $\omega 2=\omega 10$ | $\omega 1=\omega 2$ |

2.4 | 0 | $\omega 1=\omega 10$ | $\omega 2=\omega 10$ | … | $\omega 1=\omega 10$ |

2.4 | −2.4 | $\omega 1=\omega 10$ | $\omega 2=\omega 10$ | $\omega 2=\omega 10$ | $\omega 1=\omega 2$ |

Delay . | Approximate resonance conditions . | ||||
---|---|---|---|---|---|

$\tau 21/\Delta t$ . | $\tau 22\u2032/\Delta t$ . | $\rho 0\u21921\rho 1$ . | $\rho 1\u21922\rho 2$ . | $\rho 2\u21922\u2032\rho 3$ . | $\rho 3\u2192$ detection at $\omega m=\omega 1$ . |

0 | 0 | $\omega 1=\omega 10$ | $\omega 1=\omega 2$ | $\omega 1=\omega 10$ | … |

0 | −2.4 | $\omega 1=\omega 10$ | $\omega 1=\omega 2$ | $\omega 2=\omega 10$ | $\omega 1=\omega 2$ |

2.4 | 0 | $\omega 1=\omega 10$ | $\omega 2=\omega 10$ | … | $\omega 1=\omega 10$ |

2.4 | −2.4 | $\omega 1=\omega 10$ | $\omega 2=\omega 10$ | $\omega 2=\omega 10$ | $\omega 1=\omega 2$ |

When the pulses are all overlapped ($\tau 21\u2009=\u2009\tau 22\u2032\u2009=\u20090$, lower right, orange), all transitions in the Liouville pathway are simultaneously driven by the incident fields. This spectrum strongly resembles the driven limit spectrum. For this time-ordering, the first, second, and third density matrix elements have driven resonance conditions of $\omega 1\u2009=\u2009\omega 10$, $\omega 1\u2009\u2212\u2009\omega 2\u2009=\u20090$, and $\omega 1\u2009\u2212\u2009\omega 2\u2009+\u2009\omega 2\u2032\u2009=\u2009\omega 10$, respectively. The second resonance condition causes elongation along the diagonal, and since $\omega 2\u2009=\u2009\omega 2\u2032$, the first and third resonance conditions are identical, effectively making $\omega 1$ doubly resonant at $\omega 10$ and resulting in the vertical elongation along $\omega 1\u2009=\u2009\omega 10$.

The other three spectra in Fig. 5 separate the pulse sequence over time so that not all interactions are driven. At $\tau 21\u2009=\u20090$, $\tau 22\u2032\u2009=\u2009\u22122.4\Delta t$ (lower left, pink), the first two resonances remain the same as at pulse overlap (orange), but the last resonance is different. The final pulse, $E2\u2032$, is latent and probes $\rho 2$ during its FID evolution after the memory of the driven frequency is lost. There are two important consequences. First, the third driven resonance condition is now approximated by $\omega 2\u2032\u2009=\u2009\omega 10$, which makes $\omega 1$ only singly resonant at $\omega 1\u2009=\u2009\omega 10$. Second, the driven portion of the signal frequency is determined only by the latent pulse: $\omega out\u2009=\u2009\omega 2\u2032$. Since our monochromator gates $\omega 1$, we have the detection-induced correlation $\omega 1\u2009=\u2009\omega 2\u2032$. The net result is double resonance along $\omega 1\u2009=\u2009\omega 2$, and the vertical elongation of pulse overlap is strongly attenuated.

At $\tau 21\u2009=\u20092.4\Delta t,\tau 22\u2032\u2009=\u20090$ (upper right, purple), the first pulse *E*_{1} precedes the latter two, which makes the two resonance conditions for the input fields $\omega 1\u2009=\u2009\omega 10$ and $\omega 2\u2009=\u2009\omega 10$. The signal depends on the FID conversion of $\rho 1$, which gives vertical elongation at $\omega 1\u2009=\u2009\omega 10$. Furthermore, $\rho 1$ has no memory of $\omega 1$ when *E*_{2} interacts, which has two important implications. First, this means the second resonance condition is $\omega 1\u2009=\u2009\omega 2$ and the associated diagonal elongation is now absent. Second, the final output polarization frequency content is no longer functional of $\omega 1$. Coupled with the fact that *E*_{2} and $E2\u2032$ are coincident so that the final coherence can be approximated as driven by these two, we can approximate the final frequency as $\omega out\u2009=\u2009\omega 10\u2212\omega 2+\omega 2\u2032\u2009=\u2009\omega 10$. Surprisingly, the frequency content of the output is strongly independent of all pulse frequencies. The monochromator narrows the $\omega 1\u2009=\u2009\omega 10$ resonance. The $\omega 1\u2009=\u2009\omega 10$ resonance condition now depends on the monochromator slit width, the FID propagation of $\rho 1$, and the spectral bandwidth of $\rho 3$; its spectral width is not easily related to material parameters. This resonance demonstrates the importance of the detection scheme for experiments and how the optimal detection can change depending on the pulse delay time.

Finally, when all pulses are well-separated ($\tau 21\u2009=\u2009\u2212\tau 22\u2032=\u20092.4\Delta t$, upper left, cyan), each resonance condition is independent and both *E*_{1} and *E*_{2} require FID buildup to produce the final output. The resulting line shape is narrow in all directions. Again, the emitted frequency does not depend on $\omega 1$, yet the monochromator resolves the final coherence at frequency $\omega 1$. Since the driven part of the final interaction comes from $E2\u2032$, and since the monochromator track $\omega 1$, the output signal will increase when $\omega 1\u2009=\u2009\omega 2\u2032$. As a result, the line shape acquires a diagonal character.

The changes in the line shape seen in Fig. 5 have significant ramifications for the interpretations and strategies of MR-CMDS in the mixed domain. Time-gating has been used to isolate the 2D spectra of a certain time-ordering,^{13,41,60} but here we show that time-gating itself causes significant line shape changes to the isolated pathways. The phenomenon of time-gating can cause frequency and delay axes to become a functional of each other in unexpected ways.

### C. Temporal pathway discrimination

In Sec. IV B, we showed how a single pathway’s spectra can evolve with delay due to pulse effects and time gating. In real experiments, evolution with delay is further complicated by the six time-orderings/sixteen pathways present in our three-beam experiment (see Fig. 1). Each time-ordering has different resonance conditions. When the signal is collected near pulse overlap, multiple time-orderings contribute. To identify these effects, we start by considering how strongly time-orderings are isolated at each delay coordinate.

While the general idea of using time delays to enhance certain time-ordered regions is widely applied, quantitation of this discrimination is rarely explored. Because the temporal profile of the signal is dependent on both the excitation pulse profile and the decay dynamics of the coherence itself, quantitation of pathway discrimination requires simulation.

Figure 6 shows the 2D delay space with all pathways present for $\omega 1,\omega 2\u2009=\u2009\omega 10$. It illustrates the interplay of pulse width and system decay rates on the isolation of time-ordered pathways. The color bar shows the signal amplitude. The signal is symmetric about the $\tau 21\u2009=\u2009\tau 22\u2032$ line because when $\omega 1\u2009=\u2009\omega 2$, *E*_{1} and $E2\u2032$ interactions are interchangeable: $Stot(\tau 21,\tau 22\u2032)\u2009=\u2009Stot(\tau 22\u2032,\tau 21)$. The overlaid black contours represent signal “purity,” *P*, defined as the relative amount of the signal that comes from the dominant pathway at that delay value,

The dominant pathway ($maxSL\tau 21,\tau 22\u2032$) at given delays can be inferred by the time-ordered region defined in Fig. 2(d). The contours of purity generally run parallel to the time-ordering boundaries with the exception of time-ordered regions II and IV, which involve the double quantum coherences that have been neglected.

A commonly employed metric for temporal selectivity is how definitively the pulses are ordered. This metric agrees with our simulations. The purity contours have a weak dependence on $\Delta t\Gamma 10$ for $|\tau 22\u2032|/\Delta t<1$ or $|\tau 21|/\Delta t<1$, where there is significant pathway overlap and a stronger dependence at larger values where the pathways are well-isolated. Because responses decay exponentially, while pulses decay as Gaussians, there always exist delays where temporal discrimination is possible. As $\Delta t\Gamma 10\u2192\u221e$, however, such discrimination is only achieved at vanishing signal intensities; the contour of *P* = 0.99 across our systems highlights this trend.

### D. Multidimensional line shape dependence on pulse delay time

In Secs. IV A–IV C, we showed how pathway spectra and weights evolve with delay. This section ties the two concepts together by exploring the evolution of the spectral line shape over a span of $\tau 21$ delay times that include all pathways. It is a common practice to explore spectral evolution against $\tau 21$ because this delay axis shows population evolution in a manner analogous to pump-probe spectroscopies. The $k\u21922$ and $k\u21922\u2032$ interactions correspond to the pump, and the $k\u21921$ interaction corresponds to the probe. Time-orderings V and VI are the normal pump-probe time-orderings, time-ordering III is a mixed pump-probe-pump ordering (so-called pump polarization coupling), and time-ordering I is the probe-pump ordering (so-called perturbed FID). Scanning $\tau 21$ through pulse overlap complicates interpretation of the line shape due to the changing nature and balance of the contributing time-orderings. At $\tau 21>0$, time-ordering I dominates; at $\tau 21\u2009=\u20090$, all time-orderings contribute equally; at $\tau 21<0$, time-orderings V and VI dominate (Fig. 6). Conventional pump-probe techniques recognized these complications long ago,^{77,78} but the extension of these effects to MR-CMDS has not previously been done.

Figure 7 shows the MR-CMDS spectra as well as histograms of the pathway weights, while scanning $\tau 21$ through pulse overlap. The colored histogram bars and line shape contours correspond to different values of the relative dephasing rate, $\Gamma 10\Delta t$. The contour is the half-maximum of the line shape (supplementary materials Fig. S6 shows fully colored contour plots of each 2D frequency spectrum). The dependence of the line shape amplitude on $\tau 21$ can be inferred from Fig. 6.

The qualitative trend, as $\tau 21$ goes from positive to negative delays, is a change from diagonal/compressed line shapes to much broader resonances with no correlation ($\omega 1$ and $\omega 2$ interact with independent resonances). Such spectral changes could be misinterpreted as spectral diffusion, where the line shape changes from correlated to uncorrelated as population time increases due to system dynamics. The system dynamics included here, however, contain no structure that would allow for such diffusion. Rather, the spectral changes reflect the changes in the majority pathway contribution, starting with time-ordering I pathways, proceeding to an equal admixture of I, III, V, and VI, and finishing at an equal balance of V and VI when *E*_{1} arrives well after *E*_{2} and $E2\u2032$. Time-orderings I and III both exhibit a spectral correlation in $\omega 1$ and $\omega 2$ when driven, but time-orderings V and VI do not. Moreover, such a spectral correlation is forced near zero delay because the pulses time-gate the driven signals of the first two induced polarizations. The monochromator detection also plays a dynamic role because time-orderings V and VI always emit a signal at the monochromator frequency, while in time-orderings I and III, the emitted frequency is not defined by $\omega 1$, as discussed above.

When we isolate time-orderings V and VI, we can maintain the proper scaling of the FID bandwidth in the $\omega 1$ direction because our monochromator can gate the final coherence. This gating is not possible in time-orderings I and III because the final coherence frequency is determined by $\omega 2\u2032$, which is identical to $\omega 2$.

There are differences in the line shapes for the different values of the relative dephasing rate, $\Gamma 10\Delta t$. The spectral correlation at $\tau 21/\Delta t\u2009=\u20092$ decreases in strength as $\Gamma 10\Delta t$ decreases. As we illustrated in Fig. 5, this spectral correlation is a signature of the driven signal from the temporal overlap of *E*_{1} and *E*_{2}; the loss of spectral correlation reflects the increased prominence of FID in the first coherence as the field-matter interactions become more impulsive. This increased prominence of FID also reflects an increase in signal strength, as shown by $\tau 21$ traces in Fig. 6. When all pulses are completely overlapped ($\tau 21\u2009=\u20090$), each of the line shapes exhibit spectral correlation. At $\tau 21/\Delta t=\u22122$, the line shape shrinks as $\Gamma 10\Delta t$ decreases, with the elongation direction changing from horizontal to vertical. The general shrinking reflects the narrowing homogeneous linewidth of the $\omega 10$ resonance. In all cases, the horizontal line shape corresponds to the homogeneous linewidth because the narrow bandpass monochromator resolves the final $\omega 1$ resonance. The change in the elongation direction is due to the resolving power of $\omega 2$. At $\Gamma 10\Delta t\u2009=\u20090.5$, the resonance is broader than our pulse bandwidth and is fully resolved vertically. It is narrower than the $\omega 1$ resonance because time-orderings V and VI interfere to isolate only the absorptive line shape along $\omega 2$. This narrowing, however, is unresolvable when the pulse bandwidth becomes broader than that of the resonance, which gives rise to a vertically elongated signal when $\Gamma 10\Delta t\u2009=\u20090.5$.

It is also common to represent data as “Wigner plots,” where one axis is delay and the other is frequency.^{14,15,21,57} In Fig. 8, we show five $\tau 21,\omega 1$ plots for varying $\omega 2$ with $\tau 22\u2032\u2009=\u20090$. Supplementary materials’ Fig. S8 shows Wigner plots for other $\Gamma 10\Delta t$ values. The plots are the analog to the most common multidimensional experiment of transient absorption spectroscopy, where the non-linear probe spectrum is plotted as a function of the pump-probe delay. For each plot, the $\omega 2$ frequency is denoted by a vertical gray line. Each Wigner plot is scaled to its own dynamic range to emphasize the dependence on $\omega 2$. The dramatic line shape changes between positive and negative delays can be seen. This representation also highlights the asymmetric broadening of the $\omega 1$ line shape near pulse overlap when $\omega 2$ becomes non-resonant. Again, these features can resemble spectral diffusion even though our system is homogeneous.

### E. Inhomogeneous broadening

With the homogeneous system characterized, we can now consider the effect of inhomogeneity. For inhomogeneous systems, time-orderings III and V are enhanced because their final coherence will rephase to form a photon echo, whereas time-orderings I and VI will not. In delay space, this rephasing appears as a shift of the signal to time-ordered regions III and V that persists for all population times. Figure 9 shows the calculated spectra for relative dephasing rate $\Gamma 10\Delta t\u2009=\u20091$ with a frequency broadening function of width $\Delta inhom\u2009=\u20090.441\Gamma 10$. The inhomogeneity makes it easier to temporally isolate the rephasing pathways and harder to isolate the non-rephasing pathways, as shown by the purity contours.

A common metric of rephasing in delay space is the 3PEPS measurement.^{71,79–81} In 3PEPS, one measures the signal as the first coherence time, $\tau $, is scanned across both rephasing and non-rephasing pathways while keeping population time, *T*, constant. The position of the peak is measured; a peak shifted away from $\tau \u2009=\u20090$ reflects the rephasing ability of the system. An inhomogeneous system will emit a photon echo in the rephasing pathway, enhancing the signal in the rephasing time-ordering and creating the peak shift. In our 2D delay space, the $\tau $ trace can be defined if we assume that *E*_{2} and $E2\u2032$ create the population (time-orderings V and VI). The trace runs parallel to the III-V time-ordering boundary (diagonal) if $\tau 22\u2032<0$ and runs parallel to the IV-VI time-ordering boundary (horizontal) if $\tau 22\u2032>0$, and both intersect at $\tau 22\u2032\u2009=\u20090$; the $\u2212\tau 21$ value at this intersection is *T* (supplementary materials Fig. S9 illustrates how 3PEPS shifts are measured from a 2D delay plot). In our 2D delay plots (Figs. 6 and 9), the peak shift is seen as the diagonal displacement of the signal peak from the $\tau 21\u2009=\u20090$ vertical line (supplementary materials Fig. S10 shows the 3PEPS measurements of all 12 combinations of $\Gamma 10\Delta t$ and $\Delta inhom$ for every population delay surveyed). Figure 9 highlights the peak shift profile as a function of population time with the yellow trace; it is easily verified that our static inhomogeneous system exhibits a non-zero peak shift value for all population times.

The unanticipated feature of the 3PEPS analysis is the dependence on *T*. Even though our inhomogeneity is static, the peak shift is maximal at *T* = 0 and dissipates as *T* increases, mimicking spectral diffusion. This dynamic arises from signal overlap with time-ordering III, which uses *E*_{2} and *E*_{1} as the first two interactions. At *T* = 0, the $\tau $ trace gives two ways to make a rephasing pathway (time-orderings III and V) and only one way to make a non-rephasing pathway (time-ordering VI). This pathway asymmetry shifts the signal away from $\tau \u2009=\u20090$ into the rephasing direction. At large *T* (large $\tau 21$), time-ordering III is not viable and pathway asymmetry disappears. Peak shifts imply inhomogeneity only when time-orderings V and III are minimally contaminated by each other, i.e., at population times that exceed pulse overlap. This fact is easily illustrated by the dynamics of the homogeneous system (Fig. 6); even though the homogeneous system cannot rephase, there is a non-zero peak shift near *T* = 0. The contamination of the 3PEPS measurement at pulse overlap is well known and is described in some studies,^{23,72} but the dependence of pulse and system properties on the distortion has not been investigated previously. Peak-shifting due to pulse overlap is less important when $\omega 1\u2260\omega 2$ because time-ordering III is decoupled by detuning.

In frequency space, spectral elongation along the diagonal is the signature of inhomogeneous broadening. Figure 10 shows the line shape changes of a Gaussian inhomogeneous distribution (supplementary materials Fig. S7 shows all contours for Fig. 10). All systems are broadened by a distribution proportional to their dephasing bandwidth. As expected, the sequence again shows a gradual broadening along the $\omega 1$ axis, with a strong spectral correlation at early delays ($\tau 21>0$) for the more driven signals. The anti-diagonal width at early delays (e.g., Fig. 10, $\tau 21\u2009=\u20092.0\Delta t$) again depends on the pulse bandwidth and the monochromator slit width. At delay values that isolate time-orderings V and VI, however, the line shapes retain diagonal character, showing the characteristic balance of homogeneous and inhomogeneous widths.

## V. DISCUSSION

### A. An intuitive picture of pulse effects

Our chosen values of the relative dephasing time, $\Gamma 10\Delta t$, describe experiments where neither the impulsive nor the driven limit unilaterally applies. We have illustrated that in this intermediate regime, the multidimensional spectra contain attributes of both limits and that it is possible to judge when these attributes apply. In our three-pulse experiment, the second and third pulses time-gate coherences and populations produced by the previous pulse(s), and the monochromator frequency-gates the final coherence. Time-gating isolates different properties of the coherences and populations. Consequently, spectra evolve against delay. For any delay coordinate, one can develop qualitative line shape expectations by considering the following three principles:

When time-gating during the pulse, the system pins to the driving frequency with a buildup efficiency determined by resonance.

When time-gating after the pulse, the FID dominates the system response.

The emitted signal field contains both FID and driven components; the $\omega out\u2009=\u2009\omega 1$ component is isolated by the tracking monochromator.

Figure 3 illustrates principles 1 and 2 and Fig. 4 illustrates principle 2 and 3. Figure 5 provides a detailed example of the relationship between these principles and the multidimensional line shape changes for different delay times.

The principles presented above apply to a single pathway. For rapidly dephasing systems, it is difficult to achieve complete pathway discrimination, as shown in Fig. 6. In such situations, the interference between pathways must be considered to predict the line shape. The relative weight of each pathway to the interference can be approximated by the extent of pulse overlap. The pathway weights exchange when scanning across pulse overlap, which creates the dramatic line shape changes observed in Figs. 7 and 10.

### B. Conditional validity of the driven limit

We have shown that the driven limit misses details of the line shape if $\Gamma 10\Delta t\u22481$, but we have also reasoned that in certain conditions, the driven limit can approximate the response well (see principle 1). Here we examine the line shape at delay values that demonstrate this agreement. Figure 11 compares the results of our numerical simulation (third column) with the driven limit expressions for populations where $\Gamma 11\Delta t\u2009=\u20090$ (first column) or 1 (second column). The top and bottom rows compare the line shapes when $\tau 22\u2032,\tau 21=(0,0)$ and $(0,\u22124\Delta t)$, respectively. The third column demonstrates the agreement between the driven limit approximations with the simulation by comparing the diagonal and anti-diagonal cross sections of the 2D spectra.

Note the very sharp diagonal feature that appears for $(\tau 21,\tau 22\u2032)=(0,0)$ and $\Gamma 11\u2009=\u20090$; this is due to population resonance in time-orderings I and III. This expression is inaccurate: the narrow resonance is only observed when pulse durations are much longer than the coherence time. A comparison of picosecond and femtosecond studies of quantum dot exciton line shapes (Yurs *et al.*^{82} and Kohler *et al.*,^{15} respectively) demonstrates this difference well. The driven equation fails to reproduce our numerical simulations here because resonant excitation of the population is impulsive; the experiment time-gates only the rise time of the population, yet driven theory predicts the resonance to be vanishingly narrow ($\Gamma 11\u2009=\u20090$). In light of this, one can approximate this time-gating effect by substituting population lifetime with the pulse duration ($\Gamma 11\Delta t\u2009=\u20091$), which gives good agreement with the numerical simulation (third column).

When $\tau 22\u2032\u2009=\u20090$ and $\tau 21<\Delta t$, signals can also be approximated by the driven signal (Fig. 11, bottom row). Only time-orderings V and VI are relevant. The intermediate population resonance is still impulsive, but it depends on $\omega 2\u2032\u2212\omega 2$, which is not explored in our 2D frequency space.

### C. Extracting true material correlation

We have shown that pulse effects mimic the qualitative signatures of inhomogeneity. Here we address how one can extract true system inhomogeneity in light of these effects. We focus on two ubiquitous metrics of inhomogeneity: 3PEPS for the time domain and ellipticity^{83} for the frequency domain.^{52,84} In the driven (impulsive) limit, ellipticity (3PEPS) corresponds to the frequency correlation function and uniquely extracts the inhomogeneity of the models presented here. In their respective limits, the metrics give values proportional to the inhomogeneity.

Figure 12 shows the results of this characterization for all $\Delta inhom$ and $\Gamma 10\Delta t$ values explored in this work. We study how the correlations between the two metrics depend on the relative dephasing rate, $\Gamma 10\Delta t$, the absolute inhomogeneity, $\Delta inhom/\Delta \omega $, the relative inhomogeneity $\Delta inhom/\Gamma 10$, and the population time delay (supplementary materials Figs. S10–S12 show simulations for each value of the 3PEPS and ellipticity data in Fig. 12). The top row shows the correlations of the $\Delta 3PEPS/\Delta t$ 3PEPS metric that represents the normalized coherence delay time required to reach the peak intensity. The upper right graph shows the correlations for a population time delay of $T\u2009=\u20094\Delta t$ that isolates the V and VI time-orderings. For this time delay, the $\Delta 3PEPS/\Delta t$ metric works well for all dephasing times of $\Gamma 10\Delta t$ when the relative inhomogeneity is $\Delta inhom/\Delta \omega \u2009<<\u20091$. It becomes independent of $\Delta inhom/\Delta \omega $ when $\Delta inhom/\Delta \omega \u2009>\u20091$. This saturation results because the frequency bandwidth of the excitation pulses becomes smaller than the inhomogeneous width and only a portion of the inhomogeneous ensemble contributes to the 3PEPS experiment.^{71} The corresponding graph for *T* = 0 shows that a large peak shift occurs, even without inhomogeneity. In this case, the peak shift depends on pathway overlap, as discussed in Sec. IV E.

The middle row in Fig. 12 shows the ellipticity dependence on the relative dephasing rate and inhomogeneity assuming that the measurement is performed when the first two pulses are temporally overlapped ($\tau 22\u2032\u2009=\u20090$). For a $T\u2009=\u20094\Delta t$ population time, the ellipticity is proportional to the inhomogeneity until $\Delta inhom/\Delta \omega <<1$, where the excitation bandwidth is wide compared with the inhomogeneity. Unlike 3PEPS, saturation is not observed because the pulse bandwidth does not limit the frequency range scanned. The 3PEPS and ellipticity metrics are therefore complementary since 3PEPS works well for $\Delta inhom/\Delta \omega \u2009<<\u20091$ and ellipticity works well for $\Delta inhom/\Delta \omega \u2009>>\u20091$. When all pulses are temporally overlapped at *T* = 0, the ellipticity is only weakly dependent on the inhomogeneity and dephasing rate. The ellipticity is instead dominated by the dependence on the excitation pulse frequency differences of time-orderings I and III that become important at pulse overlap.

It is clear from the previous discussion that both metrics depend on the dephasing and inhomogeneity. The dephasing can be measured independently in the frequency or time domain, depending upon whether the dephasing is very fast or slow, respectively. In the mixed frequency/time domain, the measurement of the dephasing becomes more difficult. One strategy to address this challenge is to use both the 3PEPS and ellipticity metrics. The bottom row in Fig. 12 plots 3PEPS against ellipticity to show how the relationship between the metrics changes for different amounts of dephasing and inhomogeneity. The anti-diagonal contours of constant relative inhomogeneity show that these metrics are complementary and can serve to extract the system correlation parameters.

Importantly, the metrics are uniquely mapped both in the presence and absence of pulse-induced effects (demonstrated by *T* = 0 and $T\u2009=\u20094\Delta t$, respectively). The combined metrics can be used to determine correlation at *T* = 0, but the correlation-inducing pulse effects give a mapping significantly different than that at $T\u2009=\u20094\Delta t$. At *T* = 0, 3PEPS is almost nonresponsive to inhomogeneity; instead, it is an almost independent characterization of the pure dephasing. In fact, the *T* = 0 trace is equivalent to the original photon echo traces used to resolve pure dephasing rates.^{85} Both metrics are offset due to the pulse overlap effects. Accordingly, the region to the left of the homogeneous contour is non-physical because it represents observed correlations that are less than those given by pulse overlap effects. If the metrics are measured as a function of *T*, the mapping gradually changes from the left figure to the right figure in accordance with the pulse overlap. Both metrics will show a decrease, even with static inhomogeneity. If a system has spectral diffusion, the mapping at late times will disagree with the mapping at early times; both ellipticity and 3PEPS will be smaller at later times than predicted by the change in mappings alone.

## VI. CONCLUSION

This study provides a framework to describe and disentangle the influence of the excitation pulses in mixed-domain ultrafast spectroscopy. We analyzed the features of mixed-domain spectroscopy through detailed simulations of MR-CMDS signals. When pulse durations are similar to coherence times, resolution is compromised by time-bandwidth uncertainty and the complex mixture of driven and FID responses. The dimensionless quantity $\Delta t\Gamma f+i\kappa f\Omega fx$ captures the balance of driven and FID character in a single field-matter interaction. In the nonlinear experiment, with multiple field-matter interactions, this balance is also controlled by pulse delays and frequency-resolved detection. Our analysis shows how these effects can be intuitive.

The dynamic nature of pulse effects can lead to misleading changes to spectra when delays are changed. When delays separate pulses, the spectral line shapes of individual pathways qualitatively change because the delays isolate FID contributions and de-emphasize the driven response. When delays are scanned across pulse overlap, the weights of individual pathways change, further changing the line shapes. In a real system, these changes would all be present in addition to actual dynamics and spectral changes of the material.

Finally, we find that in either the frequency or time domain, pulse effects mimic signatures of ultrafast inhomogeneity. Even homogeneous systems take on these signatures. Driven character gives rise to pathway overlap peak shifting in the 2D delay response, which artificially produces rephasing near pulse overlap. Driven character also produces resonances that depend on $\omega 1\u2212\omega 2$ near pulse overlap. Determination of the homogeneous and inhomogeneous broadening at ultrashort times is only possible by performing correlation analysis in both the frequency and time domains.

## SUPPLEMENTARY MATERIAL

See supplementary material for the simulation results that are a sparse sampling of the total data simulated. The full dataset, as well as the software used to simulate, is available for download at http://dx.doi.org/10.17605/OSF.IO/EJ2XE. Additional derivations and figures can be found in the associated PDF. See Fig. S3 for a Fourier domain representation of Fig. 4(a). See Fig. S4 for explicit plots of $\rho 1(\Omega fx/\Delta \omega )$ at discrete $t/\Delta t$ values. See Fig. S5 for a representation of Fig. 5 simulated without monochromator frequency filtering [$M(\omega \u2009\u2212\u2009\omega 1)\u2009=\u20091$ for Eq. (10)]. Figure S6 shows fully colored contour plots of each 2D frequency spectrum. See Fig. S8 for Wigner plots for all $\Gamma 10\Delta t$ values. See Fig. S9 for an illustration of how 3PEPS shifts are measured from a 2D delay plot. Figure S10 shows the 3PEPS measurements of all 12 combinations of $\Gamma 10\Delta t$ and $\Delta inhom$, for every population delay surveyed. As in Fig. 7, Fig. 10 shows only the contours at the half-maximum amplitude. See Fig. S7 for all contours. The simulations for each value of the 3PEPS and ellipticity data in Fig. 12 appear in Figs. S10–S12.

## ACKNOWLEDGMENTS

This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-FG02-09ER46664.

### APPENDIX A: CHARACTERISTICS OF DRIVEN AND IMPULSIVE RESPONSE

The changes in the spectral line shapes described in this work are best understood by examining the driven/continuous wave (CW) and impulsive limits of Eqs. (5) and (7). The driven limit is achieved when pulse durations are much longer than the response function dynamics: $\Delta t|\Gamma f+i\kappa f\Omega fx|>>1$. In this limit, the system will adopt a steady state over excitation: $d\rho /dt\u22480$. Neglecting phase factors, the driven solution to Eq. (5) will be

The frequency and temporal envelopes of the excitation pulse control the coherence time evolution, and the relative amplitude and phase of the coherence are directly related to detuning from resonance.

The impulsive limit is achieved when the excitation pulses are much shorter than response function dynamics: $\Delta t|\Gamma f+i\kappa f\Omega fx|<<1$. The full description of the temporal evolution has two separate expressions: one for times when the pulse is interacting with the system and one for times after pulse interaction. Both expressions are important when describing CMDS experiments.

For times after the pulse interaction, $t\u2273\tau x+\Delta t$, the field-matter coupling is negligible. The evolution for these times, on resonance, is given by

This is classic free induction decay (FID) evolution: the system evolves at its natural frequency and decays at rate $\Gamma f$. It is important to note that while this expression is explicitly derived from the impulsive limit, FID behavior is not exclusive to impulsive excitation, as we have defined it. A latent FID will form if the pulse vanishes at a fast rate relative to the system dynamics.

For evaluating times near pulse excitation, $t\u2272\tau x+\Delta t$, we implement a Taylor expansion in the response function about zero: $e\u2212(\Gamma f+i\kappa f\Omega fx)u\u2009=\u20091\u2212(\Gamma f+i\kappa f\Omega fx)u+\cdots \u2009$. Our impulsive criterion requires that a low order expansion suffices; it is instructive to consider the result of the first order expansion of Eq. (5),

During this time, $\rho \u0303f$ builds up roughly according to the integration of the pulse envelope. The build-up is integrated because the pulse transfers energy before appreciable dephasing or detuning occurs. Contrary to the expectation of impulsive evolution, the evolution of $\rho \u0303f$ is explicitly affected by the pulse frequency, and the temporal profile evolves according to the pulse.

It is important to recognize that the impulsive limit is defined not only by having slow relaxation relative to the pulse duration but also by small detuning relative to the pulse bandwidth (as is stated in the inequality). As detuning increases, the higher orders of the response function Taylor expansion will be needed to describe the rise time, and the driven limit of Eq. (A1) will become valid. The details of this build-up time can often be neglected in impulsive approximations because build-up contributions are often negligible in the analysis; the period over which the initial excitation occurs is small in comparison to the free evolution of the system. The build-up behavior can be emphasized by the measurement, which makes Eq. (A3) important.

We now consider full Liouville pathways in the impulsive and driven limits of Eq. (7). For the driven limit, Eq. (7) can be reduced to

It is important to note that the signal depends on the multiplication of all the fields; pathway discrimination based on pulse time-ordering is not achievable because polarizations exist only when all pulses are overlapped. This limit is the basis for frequency-domain techniques. Frequency axes, however, are not independent because the system is forced to the laser frequency and influences the resonance criterion for subsequent excitations. As an example, observe that the first two resonant terms in Eq. (A4) are maximized when $\omega x\u2009=\u2009|\omega 1|$ and $\omega y\u2009=\u2009|\omega 2|$. If $\omega x$ is detuned by some value $\epsilon $, however, the occurrence of the second resonance shifts to $\omega y\u2009=\u2009|\omega 2|+\epsilon $, effectively compensating for the $\omega x$ detuning. This shifting of the resonance results in 2D line shape correlations.

If the pulses do not temporally overlap $(\tau x\u2009+\u2009\Delta t\u2272\tau y\u2009+\u2009\Delta t\u2272\tau z+\Delta t\u2272t)$, then the impulsive solution to the full Liouville pathway of Eq. (7) is

Pathway discrimination is demonstrated here because the signal is sensitive to the time-ordering of the pulses. This limit is suited for delay-scanning techniques. The emitted signal frequency is determined by the system and can be resolved by scanning a monochromator.

The driven and impulsive limits can qualitatively describe our simulated signals at certain frequency and delay combinations. Of the three expressions, the FID limit most resembles the signal when pulses are near resonance and well-separated in time (so that build-up behavior is negligible). The build-up limit is approximated well when pulses are near-resonant and arrive together (so that build-up behavior is emphasized). The driven limit holds for large detuning, regardless of delay.

### APPENDIX B: CONVOLUTION TECHNIQUE FOR INHOMOGENEOUS BROADENING

Here we describe how to transform the data of a single reference oscillator signal to that of an inhomogeneous distribution. The oscillators in the distribution are allowed to have arbitrary energies for their states, which will cause frequency shifts in the resonances. To show this, we start with a modified, but equivalent, form of Eq. (4) as follows:

We consider two oscillators with transition frequencies $\omega f$ and $\omega f\u2032\u2009=\u2009\omega f+\delta $. So long as $|\delta |\u2264\omega f$ (so that $|\omega f+\delta |\u2009=\u2009|\omega f|+\delta $ and thus the rotating wave approximation does not change), Eq. (B1) shows that the two are related by

Because both coherences are assumed to have the same initial conditions [$\rho 0(\u2212\u221e)\u2009=\u2009\rho 0\u2032(\u2212\u221e)\u2009=\u20090$], the equality also holds when both sides of the equation are integrated. The phase factor $ei\kappa f\delta \tau x$ in the substitution arises from Eq. (1), where the pulse carrier frequency maintains its phase within the pulse envelope for all delays.

The resonance translation can be extended to higher order signals as well. For a third-order signal, we compare systems with transition frequencies $\omega 10\u2032\u2009=\u2009\omega 10+a$ and $\omega 21\u2032\u2009=\u2009\omega 21+b$. The extension of Eq. (B2) to pathway $V\beta $ gives

The translation of each laser coordinate depends on which transition is made (e.g., *a* for transitions between $|0\u2009$ and $|1\u2009$ or *b* for transitions between $|1\u2009$ and $|2\u2009$), so the exact translation relation differs between pathways. We can now compute the ensemble average of the signal for pathway $V\beta $ as a convolution between the distribution function of the system, *K*(*a*, *b*), and the single oscillator response,

For this work, we restrict ourselves to a simpler ensemble where all oscillators have equally spaced levels (i.e., *a* = *b*). This makes the translation identical for all pathways and reduces the dimensionality of the convolution. Since pathways follow the same convolution, we may also perform the convolution on the total signal field,

Furthermore, since $\kappa =\u22121$ for *E*_{1} and $E2\u2032$, while $\kappa \u2009=\u20091$ for *E*_{2}, we have $eia\kappa x\tau x+\kappa y\tau y+\kappa z\tau z\u2009=\u2009e\u2212ia\tau 1\u2212\tau 2+\tau 2\u2032$ for all pathways. Equivalently, if the electric field is parameterized in terms of laser coordinates $\omega 1$ and $\omega 2$, the ensemble field can be calculated as

## REFERENCES

$\kappa f$ also has a direct relationship to the phase-matching relationship: for transitions with *E*_{2}, $\kappa f\u2009=\u20091$, and for *E*_{1} or $E2\u2032$, $\kappa f=\u22121$.

For elucidation of the relationship between the generalized Liouville pathway notation and the specific parameters for each Liouville pathway, see Table S1 in the supplementary material.

There are many ways to characterize the ellipticity of a peak shape. We adopt the convention $E\u2009=\u2009(a2\u2212b2)/(a2+b2)$, where *a* is the diagonal width and *b* is the anti-diagonal width.