path: root/mixed_domain/chapter.tex

\chapter{Disentangling material and instrument response}  \label{cha:mix}

\textit{This Chapter borrows extensively from \textcite{KohlerDanielDavid2017a}.}

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 which 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.  %

\section{Introduction}  % =========================================================================

Ultrafast spectroscopy is based on using nonlinear interactions, created by multiple ultrashort
($10^{-9}-10^{-15}$ s) pulses, to resolve spectral information on timescales as short as the pulses
themselves. \cite{RentzepisPM1970a, MukamelShaul2000a}  %
The ultrafast specta can be collected in the time domain or the frequency
domain. \cite{ParkKisam1998a}  %

Time-domain methods scan the pulse delays to resolve the free induction decay
(FID). \cite{GallagherSarahM1998a}   %
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. \cite{LebedevMV2007a, VardenyZ1981a, JoffreM1988a,
  PollardW1992a}  %
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. \cite{DruetSAJ1979a, OudarJL1980a}  %
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. \cite{OudarJL1980a, WrightJohnCurtis1997b, WrightJohnCurtis1991a}  %
In this regime, both FID and driven processes are important. \cite{PakoulevAndreiV2006a}  %
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. \cite{PakoulevAndreiV2007a, KohlerDanielDavid2014a, GelinMaximF2009b}   %
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 \cite{HammPeter2000a}, transient grating \cite{SalcedoJR1978a, FourkasJohnT1992a,
  FourkasJohnT1992b}, transient absorption/reflection \cite{AubockGerald2012a, BakkerHJ2002a}, photon
echo \cite{DeBoeijWimP1996a, PattersonFG1984a, TokmakoffAndrei1995a}, two dimensional-infrared
spectroscopy (2D-IR) \cite{HammPeter1999a, AsplundMC2000a, ZanniMartinT2001a}, 2D-electronic
spectroscopy (2D-ES) \cite{HyblJohnD2001b, BrixnerTobias2004a}, and three pulse photon echo peak
shift (3PEPS) \cite{EmdeMichelF1998a, DeBoeijWimP1996a, DeBoeijWimP1995a, ChoMinhaeng1992a,
  PassinoSeanA1997a} spectroscopies.  %
These methods also include fully-coherent methods involving only coherences such as Stimulated
Raman Spectroscopy (SRS) \cite{YoonSangwoon2005a, McCamantDavidW2005a}, Doubly Vibrationally
Enhanced (DOVE) \cite{ZhaoWei1999a, ZhaoWei1999b, ZhaoWei2000a, MeyerKentA2003a,
  DonaldsonPaulMurray2007b, DonaldsonPaulMurray2008a, FournierFrederic2008a},
Triply Resonant Sum Frequency (TRSF) \cite{BoyleErinSelene2013a, BoyleErinSelene2013b,
  BoyleErinSelene2014a}, Sum Frequency Generation (SFG) \cite{LagutchevAlexi2007a}, Coherent
Anti-Stokes Raman Spectroscopy (CARS) \cite{CarlsonRogerJohn1990b, CarlsonRogerJohn1990c,
  CarlsonRogerJohn1991a}, and other coherent Raman methods \cite{SteehlerJK1985a}.  %

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. \cite{KwacKijeong2003a, WrightJohnCurtis2016a}  %
For time-domain CMDS (2D-IR, 2D-ES), the complications that occur when the impulsive approximation
does not strictly hold has only recently been addressed. \cite{PerlikVaclav2017a,
  SmallwoodChristopherL2016a}  %

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, \cite{GallagherSarahM1998a, FourkasJohnT1992a, FourkasJohnT1992b} many of the experiments
involve pulse widths that are comparable to the widths of the system. \cite{MeyerKentA2003a,
  DonaldsonPaulMurray2007b,  PakoulevAndreiV2009a, ZhaoWei1999a, CzechKyleJonathan2015a,
  KohlerDanielDavid2014a}  %
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 \cite{WrightJohnCurtis2003a} whose
relative importance depends on the relative importance of FID and driven
responses. \cite{DonaldsonPaulMurray2010a}  %
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,
$\vec{k}_\text{out} = \vec{k}_1 - \vec{k}_2 + \vec{k}_{2'}$.  %
Here, the subscripts represent the excitation pulse frequencies, $\omega_1$ and $\omega_2 =
\omega_{2'}$.  %
These experimental conditions were recently used to explore line shapes of excitonic
systems, \cite{KohlerDanielDavid2014a, CzechKyleJonathan2015a} and have been developed on
vibrational states as well \cite{MeyerKentA2004a}.  %
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.  %

\section{Theory}  % ===============================================================================

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}=\Gamma_{00}=0$.  %

The electric field pulses, $\left\{E_l \right\}$, are given by:
\begin{equation} \label{mix:eqn:E_l}
E_l(t; \omega_l, \tau_l, \vec{k}_l \cdot z) = \frac{1}{2}\left[c_l(t-\tau_l)e^{i\vec{k}_l\cdot z}e^{-i\omega_l(t-\tau_l)} + c.c. \right],
\end{equation}
where $\omega_l$ is the field carrier frequency, $\vec{k}_l$ 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:  %
\begin{equation}
c_l(t) = \frac{1}{\Delta_t}\sqrt{\frac{2\ln 2}{2\pi}} \exp\left(-\ln 2 \left[\frac{t}{\Delta_t}\right]^2\right),
\end{equation}
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=4 \ln 2 \approx 2.77$, where $\Delta_{\omega}$ is the
spectral FWHM (intensity scale).  %

The Liouville-von Neumann Equation propagates the density matrix, $\bm{\rho}$:
\begin{equation} \label{mix:eqn:LVN}
\frac{d\bm{\rho}}{dt} = -\frac{i}{\hbar}\left[\bm{H_0} + \bm{\mu}\cdot \sum_{l=1,2,2^\prime} E_l(t), \bm{{\rho}}\right] + \bm{\Gamma \rho}.
\end{equation}
Here $\bm{H_0}$ is the time-independent Hamiltonian, $\bm{\mu}$ is the dipole superoperator, and
$\bm{\Gamma}$ contains the pure dephasing rate of the system.  %
We perform the standard perturbative expansion of \autoref{mix:eqn:LVN} to third order in the
electric field interaction \cite{MukamelShaul1995a, YeeTK1978a, OudarJL1980a, ArmstrongJA1962a,
  SchweigertIgorV2008a} and restrict ourselves only to the terms that have the correct spatial wave
vector $\vec{k}_{\text{out}}=\vec{k}_1-\vec{k}_2+\vec{k}_{2^\prime}$.  %
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.  %
\autoref{mix:fig:WMELs} shows these sixteen pathways as Wave Mixing Energy Level (WMEL)
diagrams \cite{LeeDuckhwan1985a}.  %

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 $\hbar=1$ and abbreviate the initial and final density matrix elements as
$\rho_i$ and $\rho_f$, respectively.  %
Using the natural frequency rotating frame, $\tilde{\rho}_{ij}=\rho_{ij} e^{-i\omega_{ij}t}$, the formation of $\rho_f$ using pulse $x$ is written as
\begin{equation} \label{mix:eqn:rho_f}
\begin{split}
\frac{d\tilde{\rho}_f}{dt} =& -\Gamma_f\tilde{\rho}_f \\
&+ \frac{i}{2} \lambda_f \mu_f c_x(t-\tau_x)e^{i\kappa_f\left(\vec{k}_x\cdot z + \omega_x \tau_x \right)}e^{i\kappa_f\Omega_{fx}t}\tilde{\rho}_i(t),
\end{split}
\end{equation}
where $\Omega_{fx}=\kappa_f^{-1}\omega_f - \omega_x (=\left|\omega_f\right| - \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=+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= \lambda_f$) or emission
($\kappa_f=-\lambda_f$).  %
$\kappa_f$ also has a direct relationship to the phase matching relationship: for transitions with
$E_2$, $\kappa_f=1$, and for $E_1$ or $E_{2^\prime}$, $\kappa_f=-1$.  %
In the following equations we neglect spatial dependence ($z=0$).  %

\autoref{mix:eqn:rho_f} 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
\begin{equation} \label{mix:eqn:rho_f_int}
\begin{split}
\tilde{\rho}_f(t) =& \frac{i}{2}\lambda_f \mu_f e^{i\kappa_f \omega_x \tau_x} e^{i\kappa_f \Omega_{fx} t} \\
&\times \int_{-\infty}^{\infty} c_x(t-u-\tau_x)\tilde{\rho}_i(t-u)\Theta(u) \\
& \qquad \quad \ \ \times e^{-\left(\Gamma_f+i\kappa_f\Omega_{fx}\right)u}du,
\end{split}
\end{equation}
where $\Theta$ is the Heaviside step function.  %
\autoref{mix:eqn:rho_f_int} becomes the steady state limit expression when $\Delta_t \left|\Gamma_f
  + i \kappa_f \Omega_{fx}\right| \gg 1$, and the impulsive limit expression results when $\Delta_t
\left|\Gamma_f + i \kappa_f \Omega_{fx}\right| \ll 1$.  %
Both limits are important for understanding the multidimensional line shape changes discussed in
this paper.  %
%The steady state and impulsive limits of Equation propagates \auotoref{mix:eqn:rho_f_int} are
% discussed in TODO
% Appendix \ref{sec:cw_imp}.  %

\autoref{mix:fig:overview} gives an overview of the simulations done in this work.  %
\autoref{mix:fig:overview}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 \autoref{mix:fig:overview}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 = 0.5, 1$, or $2$.  %
The temporal axes are normalized to the pulse duration, $\Delta_t$.  The $\Gamma_{10}\Delta_t=2$
transient is mostly driven by the excitation pulse while $\Gamma_{10} \Delta_t = 0.5$ has a
substantial free induction decay (FID) component at late times.  %
\autoref{mix:fig:overview}b shows a pulse sequence (pulses are shaded orange and pink) and the
resulting system evolution of pathway $V\gamma$ ($00 \xrightarrow{2} 01 \xrightarrow{2^\prime} 11
\xrightarrow{1} 10 \xrightarrow{\text{out}} 00$) with $\Gamma_{10}\Delta_t=1$.  %
The final polarization (teal) is responsible for the emitted signal, which is then passed through a
frequency bandpass filter to emulate monochromator detection (\autoref{mix:fig:overview}c).  %
The resulting signal is explored in 2D delay space (\autoref{mix:fig:overview}d), 2D frequency space
(\autoref{mix:fig:overview}f), and hybrid delay-frequency space (\autoref{mix:fig:overview}e).  %
The detuning frequency axes are also normalized by the pulse bandwidth, $\Delta_{\omega}$.  %

We now consider the generalized Liouville pathway $L:\rho_0 \xrightarrow{x} \rho_1 \xrightarrow{y}
\rho_2 \xrightarrow{z} \rho_3 \xrightarrow{\text{out}} \rho_4$, where $x$, $y$, and $z$ denote
properties of the first, second, and third pulse, 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.  %
\autoref{mix:fig:overview}b demonstrates the correspondence between $x$, $y$, $z$ notation and 1, 2,
$2^\prime$ notation for the laser pulses with pathway $V\gamma$.

The electric field emitted from a Liouville pathway is proportional to the polarization created by
the third-order coherence:  %
\begin{equation} \label{mix:eqn:E_L}
E_L(t) = i \mu_{4}\rho_{3}(t).
\end{equation}
\autoref{mix:eqn:E_L} assumes perfect phase-matching and no pulse distortions through propagation.
\autoref{mix:eqn:rho_f_int} shows that the output field for this Liouville pathway is
	\begin{gather} \label{mix:eqn:E_L_full}
	\begin{split}
	E_L(t) =& \frac{i}{8}\lambda_1\lambda_2\lambda_3\mu_1\mu_2\mu_3\mu_4
	e^{i\left( \kappa_1\omega_x\tau_x + \kappa_2\omega_y\tau_y + \kappa_3\omega_z\tau_z \right)}
	e^{-i\left( \kappa_3 \omega_z + \kappa_2 \omega_y + \kappa_1 \omega_x \right) t} \\
	&\times \iiint_{-\infty}^{\infty} c_z(t-u-\tau_z) c_y(t-u-v-\tau_y) c_x(t-u-v-w-\tau_x) R_L(u,v,w) dw \ dv \ du	,
	\end{split}\\
	R_L(u,v,w) = \Theta(w)e^{-\left(\Gamma_1 + i\kappa_1\Omega_{1x} \right)w}
	\Theta(v)e^{-\left(\Gamma_2 + i\left[ \kappa_1\Omega_{1x}+\kappa_2\Omega_{2y} \right] \right)v}
	\Theta(u)e^{-\left(\Gamma_3 + i\left[ \kappa_1\Omega_{1x}+\kappa_2\Omega_{2y}+\kappa_3\Omega_{3z} \right] \right)u},
	\end{gather}
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:
\begin{equation} \label{mix:eqn:superposition}
E_{\text{tot}}= \sum_L E_L(t).
\end{equation}
For the superposition of \autoref{mix:eqn:superposition} 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}\neq\Gamma_{21}$,
$\omega_{10}\neq\omega_{21}$, and/or $\sqrt{2}\mu_{10}\neq\mu_{21}$.  %
Our simulations use the last inequality, which is important in two-level systems ($\mu_{21}=0$) 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 = \omega_1 - \omega_2 + \omega_{2'} = \omega_1$.  %
The monochromator can also enhance spectral resolution, as we show in
\autoref{mix:sec:evolution_SQC}.  %
In this simulation, the detection is emulated by transforming $E_{\text{tot}}(t)$ into the
frequency domain, applying a narrow bandpass filter, $M(\omega)$, about $\omega_1$, and applying
amplitude-scaled detection:
\begin{equation} \label{mix:eqn:S_tot}
S_{\text{tot}}(\omega_1, \omega_2, \tau_{21}, \tau_{22^\prime}) =
\sqrt{ \int\left| M(\omega-\omega_1) E_{\text{tot}}(\omega) \right|^2 d\omega},
\end{equation}
where $E_{\text{tot}}(\omega)$ denotes the Fourier transform of $E_{\text{tot}}(t)$ (see
\autoref{mix:fig:overview}c).  %
For $M$ we used a rectangular function of width $0.408\Delta_{\omega}$.  %
The arguments of $S_{\text{tot}}$ refer to the \textit{experimental} degrees of freedom.  %
The signal delay dependence is parameterized with the relative delays $\tau_{21}$ and
$\tau_{22^\prime}$, where $\tau_{nm} = \tau_n - \tau_m$ (see \autoref{mix:fig:overview}b).  %
Table S1 summarizes the arguments for each Liouville pathway.  %
\autoref{mix:fig:overview}f shows the 2D $(\omega_1, \omega_2)$ $S_{\text{tot}}$ spectrum resulting
from the pulse delay times represented in \autoref{mix:fig:overview}b.  %

\begin{figure}
	\includegraphics[width=0.5\linewidth]{"mixed_domain/WMELs"}
	\caption[Sixteen triply-resonant Liouville pathways.]{
		The sixteen triply-resonant Liouville pathways for the third-order response of the system used
    here.
    Time flows from left to right.
    Each excitation is labeled by the pulse stimulating the transition; excitatons with $\omega_1$
    are yellow, excitations with $\omega_2=\omega_{2'}$ are purple, and the final emission is gray.
	}
	\label{mix:fig:WMELs}
\end{figure}

\begin{figure}
	\includegraphics[width=\linewidth]{"mixed_domain/simulation overview"}
	\caption[Overview of the MR-CMDS simulation.]{
		Overview of the MR-CMDS simulation.
		(a) The temporal profile of a coherence under pulsed excitation depends on how quickly the
    coherence dephases.  In all subsequent panes, the relative dephasing rate is kept constant at
    $\Gamma_{10}\Delta_t=1$.
		(b) Simulated evolution of the density matrix elements of a third-order Liouville pathway
    $V\gamma$ under fully resonant excitation. Pulses can be labeled both by their time of arrival
    ($x$,$y$,$z$) and by the lab lasers used to stimulate the transitions ($2$,$2^\prime$,$1$). The
    final coherence (teal) creates the output electric field.
		(c) The frequency profile of [Eccentricity ]the output electric field is filtered by a
    monochromator gating function, $M(\omega)$, and the passed components (shaded) are measured.
		(d-f) Signal is viewed  against two laser parameters, either as 2D delay (d), mixed
    delay-frequency (e), or 2D frequency plots (f).  The six time-orderings are labeled in (d) to
    help introduce our delay convention.
	}
	\label{mix:fig:overview}
\end{figure}

\subsection{Characteristics of driven and impulsive response} \label{mix:sec:cw_imp}  % -----------

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 \autoref{mix:eqn:rho_f_int} and
\ref{mix:eqn:E_L_full}.  %
The driven limit is achieved when pulse durations are much longer than the response function
dynamics: $\Delta_t \left|\Gamma_f + i \kappa_f \Omega_{fx}\right| \gg 1$.  %
In this limit, the system will adopt a steady state over excitation:  $d\rho / dt \approx 0$.  %
Neglecting phase factors, the driven solution to \autoref{mix:eqn:rho_f_int} will be
\begin{equation} \label{mix:eqn:sqc_driven}
\tilde{\rho}_f(t) = \frac{\lambda_f \mu_f}{2}
\frac{c_x(t-\tau_x)e^{i\kappa_f \Omega_{fx}t}}{\kappa_f \Omega_{fx}} \tilde{\rho}_i(t).
\end{equation}
The frequency and temporal envelope of the excitation pulse controls the coherence time evolution,
and the relative amplitude and phase of the coherence is directly related to detuning from
resonance.  %

The impulsive limit is achieved when the excitation pulses are much shorter than response function
dynamics:  $\Delta_t \left|\Gamma_f + i \kappa_f \Omega_{fx}\right| \ll 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 \gtrsim \tau_x + \Delta_t$, the field-matter coupling is
negligible.  %
The evolution for these times, on resonance, is given by
\begin{equation} \label{mix:eqn:sqc_fid}
\tilde{\rho}_f(t) =\frac{i \lambda_f\mu_f }{2} \tilde{\rho}_i(\tau_x)
\int c_x(u) du \ e^{-\Gamma_f(t-\tau_x)}.
\end{equation}
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 \lesssim \tau_x + \Delta_t$, we implement a Taylor
expansion in the response function about zero: $e^{-(\Gamma_f+i\kappa_f\Omega_{fx})u} = 1 -
(\Gamma_f+i\kappa_f\Omega_{fx})u+\cdots$.  %
Our impulsive criterion requires that a low order expansion will suffice; it is instructive to
consider the result of the first order expansion of \autoref{mix:eqn:rho_f_int}:  %
\begin{equation} \label{mix:eqn:sqc_rise}
\begin{split}
\tilde{\rho}_f(t) =& \frac{i \lambda_f\mu_f}{2} e^{-i\kappa_f\omega_x\tau_x}e^{-i\kappa_f\Omega_{fx}t} \tilde{\rho}_i(\tau_x) \\
& \times \bigg[ \left( 1-(\Gamma_f + i\kappa_f\Omega_{fx})(t-\tau_x) \right) \int_{-\infty}^{t-\tau_x} c_x(u) du \\
& \quad +(\Gamma_f + i\kappa_f\Omega_{fx}) \int_{-\infty}^{t-\tau_x} c_x(u)u \ du \bigg].
\end{split}
\end{equation}
During this time $\tilde{\rho}_f$ 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 $\tilde{\rho}_f$ 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 \autoref{mix:eqn:sqc_driven} 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 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 \autoref{mix:eqn:sqc_rise}
important.  %

We now consider full Liouville pathways in the impulsive and driven limits of
\autoref{mix:eqn:E_L_full}.  %
For the driven limit, \autoref{mix:eqn:E_L_full} can be reduced to
\begin{equation} \label{mix:eqn:E_L_driven}
\begin{split}
E_L(t) =& \frac{1}{8} \lambda_1\lambda_2\lambda_3\mu_1\mu_2\mu_3\mu_4
e^{-i(\kappa_1\omega_x\tau_x + \kappa_2\omega_y\tau_y + \kappa_3\omega_z\tau_z)} \\
& \times e^{ i(\kappa_3\omega_z + \kappa_2\omega_y + \kappa_1\omega_x)t} \\
& \times c_z(t-\tau_z)c_y(t-\tau_y)c_x(t-\tau_x) \\
& \times \frac{1}{\kappa_1\Omega_{1x}-i\Gamma_1} \frac{1}{\kappa_1\Omega_{1x} + \kappa_2\Omega_{2y} - i\Gamma_2} \\
& \times \frac{1}{\kappa_1\Omega_{1x} + \kappa_2 \Omega_{2y} + \kappa_3\Omega_{3z}-i\Gamma_3}.
\end{split}
\end{equation}
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 exists 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 \autoref{mix:eqn:E_L_driven} are
maximized when $\omega_x=\left|\omega_1\right|$ and $\omega_y=\left|\omega_2\right|$.  %
If  $\omega_x$ is detuned by some value $\varepsilon$, however, the occurrence of the second
resonance shifts to $\omega_y=\left|\omega_2\right|+\varepsilon$, 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+\Delta_t \lesssim \tau_y +\Delta_t \lesssim \tau_z
+ \Delta_t \lesssim t)$, then the impulsive solution to the full Liouville pathway of
\autoref{mix:eqn:E_L_full} is  %
\begin{equation} \label{mix:eqn:E_L_impulsive}
\begin{split}
E_L(t) =& \frac{i}{8} \lambda_1\lambda_2\lambda_3\mu_1 \mu_2 \mu_3 \mu_4 e^{i(\omega_1 + \omega_2 + \omega_3)t} \\
& \times \int c_x(w) dw \int c_y(v) dv \int c_z(u) du \\
& \times e^{-\Gamma_1(\tau_y-\tau_x)} e^{-\Gamma_2(\tau_z-\tau_y)} e^{-\Gamma(t-\tau_z)}.
\end{split}
\end{equation}
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 signal when pulses are near resonance and
well-separated in time (so that build-up behavior is negligible).  %
The build-up limit approximates well when pulses are near-resonant and arrive together (so that
build-up behavior is emphasized).  %
The driven limit holds for large detunings, regardless of delay.  %

\subsection{Inhomogeneity}  % ---------------------------------------------------------------------

Inhomogeneity is isolated in CMDS through both spectral signatures, such as
line-narrowing \cite{BesemannDanielM2004a, OudarJL1980a, CarlsonRogerJohn1990a, RiebeMichaelT1988a,
  SteehlerJK1985a}, and temporal signatures, such as photon echoes \cite{WeinerAM1985a,
  AgarwalRitesh2002a}.  %
We simulate the effects of static inhomogeneous broadening by convolving the homogeneous response
with a Gaussian distribution function.  %
Dynamic broadening effects such as spectral diffusion are beyond the scope of this work.  %

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 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 \autoref{mix:eqn:rho_f}:
\begin{equation} \label{mix:eqn:rho_f_modified}
\begin{split}
\frac{d\tilde{\rho}_f}{dt} =& -\Gamma_f\tilde{\rho}_f + \frac{i}{2}\lambda_f\mu_f c_x(t-\tau_x) \\
& \times e^{i\kappa_f\left( \vec{k}\cdot z + \omega_x \tau_x \right)} e^{-i\kappa_f\left( \omega_x-\left|\omega_f \right| \right)t}\tilde{\rho}_i(t).
\end{split}
\end{equation}

We consider two oscillators with transition frequencies $\omega_f$ and $\omega_f^\prime=\omega_f +
\delta$.  %
So long as $\left| \delta \right| \leq \omega_f$ (so that $\left| \omega_f + \delta \right| =
\left| \omega_f \right| + \delta$ and thus the rotating wave approximation does not change),
\autoref{mix:eqn:rho_f_modified} shows that the two are related by  %
\begin{equation} \label{mix:eqn:freq_translation}
\frac{d\tilde{\rho}_f^\prime}{dt}(t;\omega_x) = \frac{d\tilde{\rho}_f}{dt}(t;\omega_x-\delta)e^{i\kappa_f \delta \tau_x}.
\end{equation}

Because both coherences are assumed to have the same initial conditions
($\rho_0(-\infty)=\rho_0^\prime(-\infty)=0$), the equality also holds when both sides of the
equation are integrated.  %
The phase factor $e^{i\kappa_f\delta\tau_x}$ in the substitution arises from \autoref{mix:eqn:E_l},
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}^\prime=\omega_{10}+a$ and $\omega_{21}^\prime = \omega_{21}+b$.  %
The extension of \autoref{mix:eqn:freq_translation} to pathway $V\beta$ gives  %
\begin{equation}
\begin{split}
\tilde{\rho}_3^\prime(t;\omega_2, \omega_2^\prime, \omega_1) =& \tilde{\rho}_3(t;\omega_2-a,\omega_{2^\prime}-a,\omega_1-b) \\
&\times e^{i\kappa_2 a \tau_2} e^{i\kappa_{2^\prime} a \tau_{2^\prime}} e^{i\kappa_1 b \tau_1}.
\end{split}
\end{equation}

The translation of each laser coordinate depends on which transition is made (e.g. $a$ for
transitions between $|0\rangle$ and $|1\rangle$ or $b$ for transitions between $|1\rangle$ and
$|2\rangle$), so the exact translation relation differs between pathways.  %
We can now compute the ensemble average of signal for pathway $V\beta$ as a convolution between the
distribution function of the system, $K(a,b)$, and the single oscillator response:  %
\begin{equation}
\begin{split}
\langle \tilde{\rho}_3 (t;\omega_2,\omega_{2^\prime},\omega_1) \rangle =& \iint K(a,b)\\
& \times \tilde{\rho}_3 (t;\omega_2+a,\omega_{2^\prime}+a,\omega_1+b) \\
& \times e^{i\kappa_2 a \tau_2} e^{i\kappa_{2^\prime} a \tau_{2^\prime}} e^{i\kappa_1 b \tau_1} da \ db.
\end{split}
\end{equation}
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:
\begin{equation}
\begin{split}
\langle E_{\text{tot}}(t) \rangle =& \sum_L \mu_{4,L} \int K(a,a) \\
& \times \tilde{\rho}_{3,L}(t;\omega_x-a,\omega_y-a\omega_z-a) \\
& \times e^{ia\left( \kappa_x\tau_x + \kappa_y\tau_y + \kappa_z\tau_z \right)} da.
\end{split}
\end{equation}
Furthermore, since $\kappa=-1$ for $E_1$ and $E_{2^\prime}$, while $\kappa=1$ for $E_2$, we have
$e^{ia\left( \kappa_x\tau_x + \kappa_y\tau_y + \kappa_z\tau_z \right)} = e^{-ia\left( \tau_1 -
    \tau_2 + \tau_{2^\prime} \right)}$ 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
\begin{equation} \label{mix:eqn:convolve_final}
\begin{split}
\langle E_{\text{tot}}(t;\omega_1,\omega_2) \rangle =& \int K(a,a)E_{\text{tot}}(t;\omega_1-a,\omega_2-a) \\
&\times e^{-ia\left( \tau_1-\tau_2+\tau_{2^\prime} \right)} da.
\end{split}
\end{equation}
which is a 1D convolution along the diagonal axis in frequency space.  %
\autoref{mix:fig:convolution} demonstrates the use of \autoref{mix:eqn:convolve_final} on a
homogeneous line shape.  %

\begin{figure}
	\includegraphics[width=\linewidth]{mixed_domain/convolve}
  \caption[Convolution overview.]
  {Overview of the convolution.
    (a) The homogeneous line shape.
    (b) The distribution function, $K$, mapped onto laser coordinates.
    (c) The resulting ensemble line shape computed from the convolution.
    The thick black line represents the FWHM of the distribution function.}
	\label{mix:fig:convolution}
\end{figure}

\section{Methods}  % ==============================================================================

A matrix representation of differential equations of the type in \autoref{mix:eqn:E_L_full} was
numerically integrated for parallel computation of Liouville elements (see SI for
details). \cite{DickBernhard1983a, GelinMaximF2005a}  %
The lower bound of integration was $2\Delta_t$ before the first pulse, and the upper bound was
$5\Gamma_{10}^{-1}$ 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^\prime})$ variable space.  %
Both frequencies are scanned about the resonance, and both delays are scanned about pulse overlap.
We explored the role of sample dephasing rate by calculating signal for systems with dephasing
rates such that $\Gamma_{10}\Delta_t=0.5, 1,$ and $2$.  %
Inhomogeneous broadening used a spectral FWHM, $\Delta_{\text{inhom}}$, that satisfied
$\Delta_\text{inhom}/ \Delta_{\omega}=0,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^\prime})$ and
$S_{\text{tot}}(\omega_1, \omega_2, \tau_{21}, \tau_{22^\prime})$ are recorded for each unique
Liouville pathway.  %
Our simulations were done using the open-source SciPy library. \cite{OliphantTravisE2007a}  %

There are three details of our simulation strategy that deserve more exposition:
\begin{enumerate}[topsep=-1.5ex, itemsep=-1ex, partopsep=-1ex, parsep=1ex]
	\item Liouville pathway parameters
	\item Matrix formulation of Liouville pathways
	\item Efficient integration of the Louville Equation using the Heun method
\end{enumerate}

\subsection{Liouville pathway parameters}  % ------------------------------------------------------

Table \ref{mix:tab:table1} describes the relationship between our notation and the parameters that make
up the 16 Liouville pathways.  %

\begin{landscape}
  \begin{table}
    \begin{tabular}{l c | c c c  r r r  r r r  c c c c}
      $L$ & Liouville Pathway
      & $x$ & $y$ & $z$
      & $\lambda_1$ & $\lambda_2$ & $\lambda_3$
      & $\kappa_1$ & $\kappa_2$ & $\kappa_3$
      & $\mu_1$ & $\mu_2$ & $\mu_3$ & $\mu_4$ \\
      \hline
      I$\alpha$    & $00 \rightarrow 10 \rightarrow 11 \rightarrow 10 \rightarrow 00$
      & \multirow{3}{*}{$1$} & \multirow{3}{*}{$2$} & \multirow{3}{*}{$2^\prime$}
      & 1 & -1 & -1
      & \multirow{3}{*}{-1} & \multirow{3}{*}{1} & \multirow{3}{*}{-1}
      & $\mu_{10}$ & $\mu_{10}$ & $\mu_{10}^*$ & $\mu_{10}^*$ \\
      I$\beta$     & $00 \rightarrow 10 \rightarrow 11 \rightarrow 21 \rightarrow 11$
      & & &
      & 1 & -1 & 1
      & & &
      & $\mu_{10}$ & $\mu_{10}$ & $\mu_{21}$ & $\mu_{21}^*$ \\
      I$\gamma$    & $00 \rightarrow 10 \rightarrow 00 \rightarrow 10 \rightarrow 00$
      & & &
      & 1 & 1 & 1
      & & &
      & $\mu_{10}$ & $\mu_{10}^*$ & $\mu_{10}$ & $\mu_{10}^*$ \\
      \hline
      II$\alpha$  &  $00 \rightarrow 10 \rightarrow 20 \rightarrow 10 \rightarrow 00$
      & \multirow{2}{*}{$1$} & \multirow{2}{*}{$2^\prime$} & \multirow{2}{*}{$2$}
      & 1 & 1 & 1
      & \multirow{2}{*}{-1} & \multirow{2}{*}{-1} & \multirow{2}{*}{1}
      & $\mu_{10}$ & $\mu_{21}$ & $\mu_{21}^*$ & $\mu_{10}^*$ \\
      II$\beta$   &  $00 \rightarrow 10 \rightarrow 20 \rightarrow 21 \rightarrow 11$
      & & &
      & 1 & 1 & -1
      & & &
      & $\mu_{10}$ & $\mu_{21}$ & $\mu_{10}$ & $\mu_{21}^*$ \\
      \hline
      III$\alpha$ &  $00 \rightarrow 01 \rightarrow 11 \rightarrow 10 \rightarrow 00$
      & \multirow{3}{*}{$2$} & \multirow{3}{*}{$1$} & \multirow{3}{*}{$2^\prime$}
      & -1 & 1 & -1
      & \multirow{3}{*}{1} & \multirow{3}{*}{-1} & \multirow{3}{*}{-1}
      & $\mu_{10}$ & $\mu_{10}$ & $\mu_{10}^*$ & $\mu_{10}^*$ \\
      III$\beta$  &  $00 \rightarrow 01 \rightarrow 11 \rightarrow 21 \rightarrow 11$
      & & &
      & -1 & 1 & 1
      & & &
      & $\mu_{10}$ & $\mu_{10}$ & $\mu_{21}$ & $\mu_{21}^*$ \\
      III$\gamma$ &  $00 \rightarrow 01 \rightarrow 00 \rightarrow 10 \rightarrow 00$
      & & &
      & -1 & -1 & 1
      & & &
      & $\mu_{10}$ & $\mu_{10}^*$ & $\mu_{10}$ & $\mu_{10}^*$ \\
      \hline
      IV$\alpha$  &  $00 \rightarrow 10 \rightarrow 20 \rightarrow 10 \rightarrow 00$
      & \multirow{2}{*}{$2^\prime$} & \multirow{2}{*}{$1$} & \multirow{2}{*}{$2$}
      & 1 & 1 & 1
      & \multirow{2}{*}{-1} & \multirow{2}{*}{-1} & \multirow{2}{*}{1}
      & $\mu_{10}$ & $\mu_{21}$ & $\mu_{21}^*$ & $\mu_{10}^*$ \\
      IV$\beta$   &  $00 \rightarrow 10 \rightarrow 20 \rightarrow 21 \rightarrow 11$
      & & &
      & 1 & 1 & -1
      & & &
      & $\mu_{10}$ & $\mu_{21}$ & $\mu_{10}$ & $\mu_{21}^*$ \\
      \hline
      V$\alpha$   &  $00 \rightarrow 01 \rightarrow 00 \rightarrow 10 \rightarrow 00$
      & \multirow{3}{*}{$2$} & \multirow{3}{*}{$2^\prime$} & \multirow{3}{*}{$1$}
      & -1 & -1 & 1
      & \multirow{3}{*}{1} & \multirow{3}{*}{-1} & \multirow{3}{*}{-1}
      & $\mu_{10}$ & $\mu_{10}^*$ & $\mu_{10}$ & $\mu_{10}^*$ \\
      V$\beta$    &  $00 \rightarrow 01 \rightarrow 11 \rightarrow 21 \rightarrow 11$
      & & &
      & -1 & 1 & 1
      & & &
      & $\mu_{10}$ & $\mu_{10}$ & $\mu_{21}$ & $\mu_{21}^*$ \\
      V$\gamma$   &  $00 \rightarrow 01 \rightarrow 11 \rightarrow 10 \rightarrow 00$
      & & &
      & -1 & 1 & -1
      & & &
      & $\mu_{10}$ & $\mu_{10}$ & $\mu_{10}^*$ & $\mu_{10}^*$ \\
      \hline
      VI$\alpha$  &  $00 \rightarrow 10 \rightarrow 00 \rightarrow 10 \rightarrow 00$
      & \multirow{3}{*}{$2^\prime$} & \multirow{3}{*}{$2$} & \multirow{3}{*}{$1$}
      & 1 & 1 & 1
      & \multirow{3}{*}{-1} & \multirow{3}{*}{1} & \multirow{3}{*}{-1}
      & $\mu_{10}$ & $\mu_{10}^*$ & $\mu_{10}$ & $\mu_{10}^*$ \\
      VI$\beta$   &  $00 \rightarrow 10 \rightarrow 11 \rightarrow 21 \rightarrow 00$
      & & &
      & 1 & -1 & 1
      & & &
      & $\mu_{10}$ & $\mu_{10}$ & $\mu_{21}$ & $\mu_{21}^*$ \\
      VI$\gamma$  &  $00 \rightarrow 10 \rightarrow 11 \rightarrow 10 \rightarrow 00$
      & & &
      & 1 & -1 & -1
      & & &
      & $\mu_{10}$ & $\mu_{10}$ & $\mu_{10}^*$ & $\mu_{10}^*$ \\
    \end{tabular}
    \caption{Parameters of each Liouville Pathway.}
    \label{mix:tab:table1}
  \end{table}
\end{landscape}

\subsection{Matrix formulation}  % ----------------------------------------------------------------

\subsubsection{Generic single pathway}

In this work we explicitly incorporate our time-dependent electric fields into the Liouville
equations as done by \textcite{GelinMaximF2009a}.  %

Our third order experiment involves three successive light-matter interactions. In a generic
Liouville pathway $\rho_0 \overset{x}\rightarrow \rho_1 \overset{y}\rightarrow \rho_2
\overset{z}\rightarrow \rho_3 \overset{\mathrm{out}}\rightarrow \rho_4$, the third order
polarization can be written as three coupled differential equations:  %

\begin{eqnarray}
\frac{\mathrm{d}\tilde{\rho}_1}{\mathrm{d}t} &=& -\Gamma_1\tilde{\rho}_1 + S_1(t)\tilde{\rho}_0(t) \\
\frac{\mathrm{d}\tilde{\rho}_2}{\mathrm{d}t} &=& -\Gamma_2\tilde{\rho}_2 + S_2(t)\tilde{\rho}_1(t) \\
\frac{\mathrm{d}\tilde{\rho}_3}{\mathrm{d}t} &=& -\Gamma_3\tilde{\rho}_3 + S_3(t)\tilde{\rho}_2(t)
\end{eqnarray}

Where $S$ contains all light-matter interactions:

\begin{equation}
S_j(t) \equiv \frac{i}{2}\lambda_j\mu_je^{-i\kappa_j\omega_n\tau_n}c_n(t-\tau_n)e^{i(\kappa_j\omega_n-\omega_j)t}
\end{equation}

Here, $j$ denotes the transition and $n$ denotes the driving pulse.
See discussion of Equation 4 for identification of each term.
The emitted electric field is proportional to $\tilde{\rho}_3$ (see Equation 6), so there is no
need to explicitly consider $\rho_3 \overset{\mathrm{out}}\rightarrow \rho_4$.  %
Note that we do not include depletion due to light-matter interaction.
0 is the ground state, so $\rho_0 = 1$ and $\Gamma_0 = 0$.

Our three coupled differential equations can be cast into a matrix form:

\begin{eqnarray}
  \frac{\mathrm{d}\overline{\rho}}{\mathrm{d}t} &=&
  \overline{\overline{Q}}(t)\overline{\rho} \label{mix:eqn:generic_diff} \\
  \overline{\rho} &\equiv&
  \begin{bmatrix}
  \tilde{\rho}_0 \\
  \tilde{\rho}_1 \\
  \tilde{\rho}_2 \\
  \tilde{\rho}_3
  \end{bmatrix} \\
  \overline{\overline{Q}}(t) &\equiv&
  \begin{bmatrix}
  -\Gamma_0 & 0 & 0 & 0 \\
  S_1(t) & -\Gamma_1 & 0 & 0 \\
  0 & S_2(t) & -\Gamma_2 & 0 \\
  0 & 0 & S_3(t) & -\Gamma_3
  \end{bmatrix}
\end{eqnarray}

\autoref{mix:eqn:generic_diff} could be numerically integrated to determine $\tilde{\rho}_3$ for
this pathway.  %

\subsubsection{Full formulation}  % ----------------------------------------------------------------

As discussed throughout this chapter and shown in \autoref{mix:fig:WMELs}, there are 6 unique
time-orderings and 16 unique pathways to consider in this work.  %
One could write a separate differential equation for each pathway---this would require tabulation
of 48 density matrix elements.  %
Instead, we capitalize on the symmetry of our pathways to create a single differential equation
that is computationally cheaper.  %

\begin{figure}
	\includegraphics[width=\textwidth]{"mixed_domain/matrix flow diagram"}
	\label{mix:fig:matrix_flow_diagram}
	\caption[Liouville pathway network.]{
    Network of Liouville pathways in this work.
    Superscripts denote the field interactions that have occurred to make the density matrix
    element.
    Colors denote the pulse that is used for each transition (blue is $E_1$, red is $E_{2^\prime}$
    and green is $E_2$).
  }
\end{figure}

\autoref{mix:fig:matrix_flow_diagram} diagrams our 16 Liouville pathways as a network of
interconnected steps.  %
Some density matrix elements, such as $\tilde{\rho}_{10}$, appear multiple times because there are
several distinguishable pathways that involve that element.  %
We do not include $\tilde{\rho}_{00}^{(1-2)}$ and $\tilde{\rho}_{00}^{(2^\prime-2)}$ because they
have exactly the same amplitude as $\tilde{\rho}_{11}^{(1-2)}$ and
$\tilde{\rho}_{11}^{(2^\prime-2)}$ within our simulation. Pathways where $\rho_3$ is fed by
$\tilde{\rho}_{00}$ are sign-flipped to account for this.  %

First we define a state vector containing all nine elements in
\autoref{mix:fig:matrix_flow_diagram}:

\begin{equation}
\overline{\rho} \equiv
\begin{bmatrix}
\tilde{\rho}_{00} \\
\tilde{\rho}_{01}^{(-2)} \\
\tilde{\rho}_{10}^{(2^\prime)} \\
\tilde{\rho}_{10}^{(1)} \\
\tilde{\rho}_{20}^{(1+2^\prime)} \\
\tilde{\rho}_{11}^{(1-2)} \\
\tilde{\rho}_{11}^{(2^\prime-2)} \\
\tilde{\rho}_{10}^{(1-2+2^\prime)} \\
\tilde{\rho}_{21}^{(1-2+2^\prime)}
\end{bmatrix}
\end{equation}

To write $\overline{\overline{Q}}$, we pull most of the time dependence into a series of six variables, one for each combination of pulse (subscript) and material resonance ($A$ for $10$, $B$ for $21$).

\begin{eqnarray}
A_1 &\equiv& \frac{i}{2}\mu_{10}e^{-i\omega_1\tau_1}c_1(t-\tau_1)e^{i(\omega_1-\omega_{10})t} \\
A_2 &\equiv& \frac{i}{2}\mu_{10}e^{i\omega_2\tau_2}c_2(t-\tau_2)e^{-i(\omega_2-\omega_{10})t} \\
A_{2^\prime} &\equiv& \frac{i}{2}\mu_{10}e^{-i\omega_{2^\prime}\tau_{2^\prime}}c_{2^\prime}(t-\tau_{2^\prime})e^{i(\omega_{2^\prime}-\omega_{10})t} \\
B_1 &\equiv& \frac{i}{2}\mu_{21}e^{-i\omega_1\tau_1}c_1(t-\tau_1)e^{i(\omega_1-\omega_{21})t} \\
B_2 &\equiv& \frac{i}{2}\mu_{21}e^{i\omega_2\tau_2}c_2(t-\tau_2)e^{-i(\omega_2-\omega_{21})t} \\
B_{2^\prime} &\equiv& \frac{i}{2}\mu_{21}e^{-i\omega_{2^\prime}\tau_{2^\prime}}c_{2^\prime}(t-\tau_{2^\prime})e^{i(\omega_{2^\prime}-\omega_{21})t}
\end{eqnarray}

Each off-diagonal element in $\overline{\overline{Q}}$ has two influences determining its sign: the $\tilde{\rho}_{00}$ feeding effect discussed above and $\lambda$, $+$ for ket-side and $-$ for bra-side transitions. $\tilde{\rho}_{11}^{(1-2)} \overset{2^\prime}\rightarrow \tilde{\rho}_{10}^{(1-2+2^\prime)}$ and $\tilde{\rho}_{11}^{(2^\prime-2)} \overset{1}\rightarrow \tilde{\rho}_{10}^{(1-2+2^\prime)}$ are each doubled to account for the $\alpha$ and $\gamma$ pathways that are overlapped due to our combination of $\tilde{\rho}_{00}$ and $\tilde{\rho}_{11}$.

\begin{equation}
\overline{\overline{Q}} \equiv
\begin{bmatrix}
	0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
	-A_2 & -\Gamma_{10} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
	A_{2^\prime} & 0 & -\Gamma_{10} & 0 & 0 & 0 & 0 & 0 & 0 \\
	A_1 & 0 & 0 & -\Gamma_{10} & 0 & 0 & 0 & 0 & 0 \\
	0 & 0 & B_1 & B_{2^\prime} & -\Gamma_{20} & 0 & 0 & 0 & 0 \\
	0 & A_1 & 0 & -A_2 & 0 & -\Gamma_{11} & 0 & 0 & 0 \\
	0 & A_{2^\prime} & -A_2 & 0 & 0 & 0 & -\Gamma_{11} & 0 & 0 \\
	0 & 0 & 0 & 0 & B_2 & -2A_{2^\prime} & -2A_1 & -\Gamma_{10} & 0 \\
	0 & 0 & 0 & 0 & -A_2 & B_{2^\prime} & B_1 & 0 & -\Gamma_{21}
\end{bmatrix}
\label{eq:single_Q}
\end{equation}

Note that this approach implicitly enforces our phase matching conditions and pathways.
To simulate single time-orderings we simply remove elements from \ref{eq:single_Q}.
For example, to isolate time ordering \RomanNumeral{1}:

\begin{equation}
\overline{\overline{Q}}_\RomanNumeral{1} \equiv
\begin{bmatrix}
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & -\Gamma_{10} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & -\Gamma_{10} & 0 & 0 & 0 & 0 & 0 & 0 \\
A_1 & 0 & 0 & -\Gamma_{10} & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & -\Gamma_{20} & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & -A_2 & 0 & -\Gamma_{11} & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & -\Gamma_{11} & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & -2A_{2^\prime} & 0 & -\Gamma_{10} & 0 \\
0 & 0 & 0 & 0 & 0 & B_{2^\prime} & 0 & 0 & -\Gamma_{21}
\end{bmatrix}
\label{eq:single_Q_pw1}
\end{equation}

These equations are somewhat general to four wave mixing of systems like the one considered in this
work.  %
We have not taken all the simplifications that are possible in our specific system, such as
$\omega_{10} = \omega_{21}$ and $\Gamma_{11} = 0$.  %
See \autoref{mix:sec:scripts} for more details about our simulation, including where to find the
scripts and raw output behind this work.  %

\subsection{Heun method}

The equations defining evolution of a density matrix element under the influence of a perturbative
time-dependent coupling source can be expressed as a first order linear ordinary differential
equation with the form

\begin{equation}
\frac{dy}{dt} = -Py + Q(t)
\label{eq:ordinary_diff}
\end{equation}

For finite time steps, $\Delta$, numerical integration can be achieved by Taylor expansion about
the current time point, $t_0$:

\begin{equation}
y(t_0+\Delta) = y(t_0) + \Delta\frac{d}{dt}y(t_0) + \frac{\Delta^2}{2}\frac{d^2}{dt^2}y(t_0)+\cdots
\label{eq:taylor_expansion}
\end{equation}

The Taylor expansion must be truncated to compute.
A first order expansion (\textit{i.e.} the Euler method) will suffice for sufficiently small time
steps.
\autoref{eq:ordinary_diff} makes the first-order expansion easily evaluated:

\begin{eqnarray}
y(t_0+\Delta) &\approx& y(t_0) + \Delta\frac{d}{dt}y(t_0) \\
&=& y(t_0) + \Delta \mathpzc{f}(t_0, y(t_0))
\end{eqnarray}

We have used the operator $\mathpzc{f}(t_0, y(t_0)) \equiv -Py(t_0)+Q(t_0)$ for simplicity.

We used the Heun method (also known as the improved Euler method \cite{BlanchardP2006.000}), which
includes the second order of \autoref{eq:taylor_expansion}.
By doing this we can increase the time step while maintaining the same error tolerance.
For this case, a computationally cheap way to increase the order of the Taylor expansion without
explicitly evaluating the second derivative is to approximate the second derivative in terms of the
first derivative:

\begin{eqnarray}
\frac{d^2}{dt^2} &\approx& \frac{1}{\Delta}\left[\frac{d}{dt}y(t_0+\Delta)-\frac{d}{dt}y(t_0)\right] \\
&=& \frac{1}{\Delta}\left[\mathpzc{f}(t_0+\Delta, y(t_0+\Delta)-\mathpzc{f}(t_0, y(t_0))\right]
\end{eqnarray}

This substitution is problematic, however, because the value of $y(t_0+\Delta)$ is not known---indeed it is the desired quantity. For an approximation of the second derivative we approximate $y(t_0+\Delta)$ using the first-order expansion of $y$ about $t_0$:

\begin{equation}
y(t_0+\Delta) \approx \left[y(t_0)+\Delta\mathpzc{f}(t_0, y(t_0)\right]
\end{equation}

Note that since $y(t_0+\Delta t) = y(t_0)+\Delta\mathpzc{f}(t_0, y(t_0))+O(\Delta^2)$, the second derivative term will be

\begin{equation}
\frac{\Delta^2}{2}\frac{d^2}{dt^2}y(t_0) = \frac{\Delta}{2}\mathpzc{f}\left(t_0, y(t_0)+\Delta\mathpzc{f}(t_0,y(t_0))\right) + O(\Delta^3)
\end{equation}

which has the same error scaling as truncating the Taylor series at second order. With our second derivative substitution, the second-order Taylor expansion becomes

\begin{eqnarray}
y(t_0+\Delta) &\approx& y(t_0) + \frac{\Delta}{2}\left[\mathpzc{f}\left(t_0, y(t_0)\right) + \mathpzc{f}\left(t_0+\Delta, y(t_0)+\Delta\mathpzc{f}\left(t_0,y(t_0)\right)\right)\right] \\
&=& y(t_0) + \frac{\Delta}{2}\left[-Py(t_0)+Q(t_0)-P\left(y(t_0)+\Delta\left(-Py(t_0)+Q(t_0)\right)\right)+Q(t_0+\Delta)\right] \\
&=& y(t_0) + \frac{\Delta}{2}\left[Q(t_0) + Q(t_0+\Delta) - 2Py(t_0) + P^2\Delta y(t_0)\right]
\end{eqnarray}

Which is now easily evaluated.

It is not clear \textit{a priori} that the Heun method is more efficient than the Euler method.
Efficiency is the product of the number of time steps used and the amount of time spent at each
point.  %
The Heun method would presumably be able to use fewer time steps, but each point sampled requires
more calculation.  %
Inevitably, efficiency is also a question of external factors, such as how efficiently the code
executes each method.  %

We resorted to empirical tests to verify that the Heun methods helped our implementation of the
integration.
First we examined the convergence of the integration techniques.
\autoref{mix:fig:heun} shows the integration of a Liouville pathway using our software with Euler
and Heun methods.
The mutual convergence to the true integral value of 1 is seen with both methods.
For a given error tolerance, the Heun method requires roughly half as many points as the Euler
method.
With this in mind, we benchmarked computations where the Euler method simulations used rougly twice
as many points as the Heun method simulations.
We found the Heun method to be roughly 1.9 times faster, indicating the extra computation of the
Heun method had minimal effect on the total computation time.

We note that there are more advanced methods suited for integrating ordinary differential
equations, such as multi-stage Runge-Kutta methods.
These methods may speed up computational times even further, but we have not investigated them at
this time.

\begin{figure}
 	\includegraphics[scale=0.5]{"mixed_domain/heun"}
	\label{mix:fig:heun}
	\caption[Integration technique comparison.]{
    Comparison of numerical integration techniques for a Liouville pathway signal with different
    number of time steps.
  }
\end{figure}

\subsection{Scripts \& raw output} \label{mix:sec:scripts}  % -------------------------------------

We have made the tools and raw output of this work publicly available, including our\dots
\begin{itemize}[topsep=-1.5ex, itemsep=-1ex, partopsep=-1ex, parsep=1ex]
	\item custom simulation package
	\item raw output
	\item processing scripts
	\item figure creation scripts
\end{itemize}
This section contains details necessary to understand and use what we have shared.
All supplementary information, including this PDF, can be found on the Open Science Framework, DOI:
\href{https://dx.doi.org/10.17605/OSF.IO/EJ2XE}{10.17605/OSF.IO/EJ2XE}.

The raw and measured simulation arrays are stored within this supplementary information as a series
of Hierarchical Data Format (HDF5) files.
HDF5 is an open source format supported by many programs and programming languages, see
\href{https://www.hdfgroup.org}{https://www.hdfgroup.org} for more information.

To save space we compressed the arrays using the lossless compression format gzip, see
\href{http://www.gzip.org}{http://www.gzip.org} for more information.

If you are opening the arrays within Python we recommend using the fantastic h5py library (part of
the Scipy stack) to easily decompress and import.

\subsection{Simulation package}

We have built an open source Python package for simulating CMDS. We call it Numerical Integration of Schrödinger's Equation (NISE). It is built on top of the open source SciPy library \cite{JonesEric2011.000}, and is compatible with all versions of Python 2.7 and newer.

We have included the NISE package in this supplementary information. To use NISE first install Python and SciPy if you haven't already, see \href{https://www.scipy.org/install.html}{https://www.scipy.org/install.html}. Then place the NISE package into a directory within your \texttt{PYTHONPATH}. You should be able to \texttt{import NISE} to run and interact with simulations.

We have also built a general-purpose data analysis package called WrightTools. NISE does not require WrightTools, but many of the processing scripts and the figure creation script used in this work do. For this reason we have included the WrightTools package in the supplementary information.

Our software packages are constantly being developed. The versions kept in this supplementary information are primarily for legacy purposes, to ensure that the processing scripts are kept in their full context. Please refer to our GitHub page, \href{https://github.com/wright-group}{https://github.com/wright-group}, for the most recent version of our open source software and hardware.

\subsubsection{Raw}

The raw simulation in this work was generated using NISE and the \texttt{NISE\_iterator.py} script found in the raw folder.

In NISE, an \texttt{experiment} module is loaded to define the electric field variables and the
experimental conditions, like the phase matching.  %
In this case we used the \texttt{trive} module, which defines our normal, two-color four wave
mixing experimental conditions.  %

Within this work we have represented our data in terms of dimensionless quantities like
$\tau/\Delta_t$ and $(\omega-\omega_{10})/\Delta_\omega$.  %
The simulation within NISE was done with choices for each parameter, as tabulated below.  %
These quantities are necessary to fully understand the unprocessed arrays generated by NISE.  %

\begin{table}
  \begin{tabular}{r l}
		$\omega_{10}$ & 7000 $\mathrm{cm}^{-1}$ \\
		$\mu_{10}$ & 1 \\
		$\mu_{21}$ & 1 \\
		$\Gamma_{11}$ & 0 $\mathrm{fs}^{-1}$ \\
		$\Delta_t$ & 50 fs
	\end{tabular}
  \caption{
    Parameters used in large NISE simulation.
  }
\end{table}

$\Gamma_{10}$, $\Gamma_{21}$ and $\Gamma_{20}$ were kept equal.
Their exact value for a given run of the simulation depends on the $\Gamma_{10}\Delta_t$ quantity
as discussed in the paper.
We use the term \texttt{dpr} (dephasing pulse ratio) which is the inverse of $\Gamma_{10}\Delta_t$.

In NISE the system parameters are contained within the \texttt{hamiltonian} module, in this case
\texttt{H0}.
The \texttt{Omega} object contains all of the system attributes you would expect, including the
$\overline{\rho}(t)$ (\texttt{dm\_vector}) and $\overline{\overline{Q}}(t)$
(\texttt{gen\_matrix()}) matrices as in \autoref{mix:eqn:generic_diff}.
The matrix in the software does not account for relaxation and dephasing, that is accounted for
directly during the integration.
Within \texttt{H0} we actually define two $\overline{\overline{Q}}$ `permutations', one for
pathways in which $E_1$ arrives first (\texttt{w1\_first == True}) and one for pathways in which
$E_{2^\prime}$ arrives first (\texttt{w1\_first == False}).
This separate permutation approach is mathematically identical to the single matrix approach in
\autoref{eq:single_Q}, just slightly more computationally expensive.
\texttt{H0.Omega} allows you to define which time-orderings are included using the \texttt{TOs}
keyword argument.
This directly affects how the propagator is assembled, as discussed in \autoref{eq:single_Q_pw1}.

To generate the raw data we calculated the polarization at all coordinates within a
four-dimensional experimental array:  %

Arrays containing these points were assembled and passed into the \texttt{trive} module.
These axes and \texttt{H0} were passed into the \texttt{scan} module for numerical integration.
A single output array was saved for each scan.
To keep the output array sizes reasonable a separate simulation was done for each
$\Gamma_{10}\Delta_t$, time-ordering, and \texttt{d1} value.
For each of these simulations we saved one five-dimensional array to an HDF5 file:

The final index `timestep' contains the dependence of the output polarization on lab time.
It changes from simulation to simulation to help with computation speed.
For each simulation the timetep array is defined by a starting position (always 100 fs before the
first pulse arrives) an ending position ($5 \times \tau_{10}$ fs after the last pulse arrives), and
a step (4 fs for our longest dephasing time, 2 fs otherwise).

The output polarization is kept as a complex array in the lab frame.

\begin{table}
	\begin{tabular}{c | c | c | c}
		axis & center & full width & number of points \\ \hline
		w1 & 7000 & 3000 & 41 \\
		w2 & 7000 & 3000 & 41 \\
		d1 & 0 & 400 & 21 \\
		d2 & 0 & 400 & 21
	\end{tabular}
  \caption{
    Description of four-dimensional simulation coordinates.
  }
\end{table}

\begin{table}
	\begin{tabular}{c | c | c}
		index & name & size \\ \hline
		0 & w1 & 41 \\
		1 & w2 & 41 \\
		2 & d2 & 21 \\
		3 & permutation & 2 \\
		4 & timestep & variable \\
	\end{tabular}
  \caption{
    Simulation output format (polarization).
  }
\end{table}

\subsubsection{Measured}

As mentioned in the appendix, we introduce inhomogeneity by convolving the output with a
distribution function on the intensity level: `smearing'.
This is done in the measurement stage.
We store a separate HDF5 file for each combination of \texttt{dpr}, time ordering, and
$\Delta_{\text{inhom}}$.
Each HDF5 file contains four arrays, shown in \autoref{mix:tab:measured}

The coordinate arrays are in their native units (fs, $\mathrm{cm}^{-1}$).
The signal array is purely real, stored on the intensity level.

We measured the entire simulation space twice, one with and one without the monochromator bandpass
filter.  %

\subsubsection{Representations}

We present many representations of our simulated dataset in this work.
The script used to create these figures is \texttt{figures.py}, found in this supplementary
information.
In some cases some small additional simulation or some pre-processing step was necessary, these can
be found in the `precalculated' folder.

\begin{table}
	\begin{tabular}{c | c | c}
		alias & dimensions & shape \\ \hline
		\texttt{w1} & w1 & \texttt{(41,)} \\
		\texttt{w2} & w2 & \texttt{(41,)} \\
		\texttt{d1} & d1 & \texttt{(21,)} \\
		\texttt{d2} & w2 & \texttt{(21,)} \\
		\texttt{arr} & w1, w2, d1, d2 & \texttt{(41, 41, 21, 21)}
	\end{tabular}
  \caption{
    Simulation output format (measured).
  }
  \label{mix:tab:measured}
\end{table}

\section{Results}  % ==============================================================================

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

\subsection{Evolution of single coherence} \label{mix:sec:evolution_SQC}  % -----------------------

It is illustrative to first consider the evolution of single coherences, $\rho_0 \xrightarrow{x}
\rho_1$, under various excitation conditions.  %
\autoref{mix:fig:fid_dpr} 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 \ref{mix:eqn:rho_f_int} under the various conditions in Table S1), so
the profiles in \autoref{mix:fig:fid_dpr} 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 \left(\Gamma_f + i\kappa_f \Omega_{fx}\right)\approx \Gamma_{10}\Delta_t$.  %
In this case, the evolution demonstrates the maximum impulsive character the transient can
achieve.  %
The instantaneous frequency, $d\varphi/dt$, is defined as
\begin{equation}
\frac{d\varphi}{dt} = \frac{d}{dt} \tan^{-1}\left( \frac{\text{Im}\left(\rho_1(t)\right)}{\text{Re}\left(\rho_1(t)\right)} \right).
\end{equation}
The cases of $\Gamma_{10}\Delta_t=0 (\infty)$ agree with the impulsive (driven) expressions derived
in \autoref{mix:sec:cw_imp}.  %
For $\Gamma_{10}\Delta_t=0$, the signal rises as the integral of the pulse and has instantaneous
frequency close to that of the pulse (\autoref{mix:eqn:sqc_rise}), but as the pulse vanishes, the
signal adopts the natural system frequency and decay rate (\autoref{mix:eqn:sqc_fid}).  %
For $\Gamma_{10}\Delta_t=\infty$, the signal follows the amplitude and frequency of the pulse for
all times (the driven limit, \autoref{mix:eqn:sqc_driven}).  %

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.  %
In accordance, 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 (
\autoref{mix:eqn:rho_f}), and the time-dependent balance of optical coupling and system relaxation.  %
We note that our choices of $\Gamma_{10}\Delta_t=2.0, 1.0,$ and $0.5$ give coherences that have
mainly driven, roughly equal driven and FID parts, and mainly FID components, respectively.  %
FID character is difficult to isolate when $\Gamma_{10}\Delta_t=2.0$.  %


\autoref{mix:fig:fid_detuning}a shows the temporal evolution of $\rho_1$ at several values of
$\Omega_{1x}/\Delta_{\omega}$ with $\Gamma_{10}\Delta_t=1$.
\autoref{mix:fig:fid_vs_detuning_with_freq} shows the Fourier domain representation of
\autoref{mix:fig:fid_detuning}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 \autoref{mix:fig:fid_detuning}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 \autoref{mix:fig:fid_detuning}a by comparing the coherence
amplitudes at $t=0$ (driven)  and at $t / \Delta_t=2$ (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}=\pm2$ manifests as a visibly disproportionate amplitude decay
(\autoref{mix:fig:sqc_vs_t} shows explicit plots of $\rho_1(\Gamma_{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. \cite{LagutchevAlexi2007a, LagutchevAlexi2010a,
  DonaldsonPaulMurray2010a, DonaldsonPaulMurray2008a}  %

\autoref{mix:fig:fid_vs_detuning_with_freq} shows how a single quantum coherence evolves with
detuning in the time and frequency domain.
At large detunings (2, -2), the coherence is mostly driven, taking on the excitation pulse shape in
time and frequency.
On resonance the coherence has significant FID character, seen as narrowing in the frequency
domain.
For intermediate detuning (1, -1) the coherence is a complex mixture of driven and FID components,
resulting in a skewed frequency domain lineshape.
Time and frequency gating of this coherence results in complex multidimensional lineshapes in the
full four wave mixing experiment.

\autoref{mix:fig:sqc_vs_t} shows how the time-gated amplitude of a single quantum coherence changes
as the excitation pulse is detuned.
At larger detunings the coherence decays faster, resulting in a lineshape narrowing with increasing
time (red to yellow).
At long times, the coherence lineshape approaches the convolution of the excitation pulse (grey)
with the absorptive material response (narrow black).

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 $E_\text{tot}$ contains FID signal that lasts
longer than the pulse duration.  %
\autoref{mix:fig:fid_detuning}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 \autoref{mix:fig:fid_detuning}a with
$\Omega_{1x}/\Delta_{\omega}=\pm 1$, 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 \autoref{mix:eqn:rho_f_int} (assume $\rho_i=1$ and $\tau_x=0$):
\begin{equation} \label{mix:eqn:rho_f_int_freq}
  \tilde{\rho}_f (\omega) =
  \frac{i\lambda_f\mu_f}{2\sqrt{2\pi}} \cdot \frac{\mathcal{F}\left\{ c_x \right\}\left( \omega - \kappa_f\Omega_{fx} \right)}{\Gamma_f + i\omega},
\end{equation}
where $\mathcal{F}\left\{ c_x \right\}\left( \omega \right)$ denotes the Fourier transform of
$c_x$, which in our case gives
\begin{equation}
  \mathcal{F}\left\{ c_x \right\}\left( \omega \right) =
  \frac{1}{\sqrt{2\pi}} e^{-\frac{(\Delta_t\omega)^2}{4\ln 2}}.
\end{equation}
For squared-law detection of $\rho_f$, the importance of the tracking monochromator is highlighted
by two limits of \autoref{mix:eqn:rho_f_int_freq}:
\begin{itemize}
	\item When the transient is not frequency resolved, $\text{sig} \approx \int{\left|
        \tilde{\rho}_f(\omega) \right|^2 d\omega}$ and the measured line shape will be the
    convolution of the pulse envelope and the intrinsic (Green's function) response
    (\autoref{mix:fig:fid_detuning}b, magenta).
	\item When the driven frequency is isolated, $\text{sig} \approx \left|
      \tilde{\rho}_f(\kappa_f\Omega_{fx}) \right|^2$ and the measured line shape will give the
    un-broadened Green's function (\autoref{mix:fig:fid_detuning}b, teal).
\end{itemize}
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.  %
\autoref{mix:fig:fid_detuning}b shows the various detection methods for the relative dephasing rate of
$\Gamma_{10}\Delta_t=1$.  %

\begin{figure}
	\includegraphics[width=0.5\linewidth]{"mixed_domain/fid vs dpr"}
	\caption[Relative importance of FID and driven response for a single quantum coherence.]{
		The relative importance of FID and driven response for a single quantum coherence as a function
    of the relative dephasing rate (values of $\Gamma_{10}\Delta_t$ are shown inset).
		The black line shows the coherence amplitude profile, while the shaded color indicates the
    instantaneous frequency (see colorbar).
		For all cases, the pulsed excitation field (gray line, shown as electric field amplitude) is
    slightly detuned (relative detuning, $\Omega_{fx}/\Delta_{\omega}=0.1$).
	}
	\label{mix:fig:fid_dpr}
\end{figure}

\begin{figure}
	\includegraphics[width=0.5\linewidth]{"mixed_domain/fid vs detuning"}
	\caption[Pulsed excitation of a single quantum coherence and its dependance on pulse detuning.]{
		Pulsed excitation of a single quantum coherence and its dependence on the pulse detuning.  In
    all cases, relative dephasing is kept at $\Gamma_{10}\Delta_t = 1$.
		(a) The relative importance of FID and driven response for a single quantum coherence as a
    function of the detuning (values of relative detuning, $\Omega_{fx}/\Delta_{\omega}$, are shown
    inset).
		The color indicates the instantaneous frequency (scale bar on right), while the black line
    shows the amplitude profile.  The gray line is the electric field amplitude.
		%Comparison of the temporal evolution of single quantum coherences at different detunings
    %(labeled inset).
		(b)  The time-integrated coherence amplitude as a function of the detuning.  The integrated
    amplitude is collected both with (teal) and without (magenta) a tracking monochromator that
    isolates the driven frequency components.  %$\omega_m=\omega_\text{laser}$.
		For comparison, the Green's function of the single quantum coherence is also shown (amplitude
    is black, hashed; imaginary is black, solid).
		In all plots, the gray line is the electric field amplitude.
	}
	\label{mix:fig:fid_detuning}
\end{figure}

\begin{figure}
	\includegraphics[width=\textwidth]{"mixed_domain/fid vs detuning with freq"}
	\label{mix:fig:fid_vs_detuning_with_freq}
	\caption[Frequency domain representation of a single quantum coherence vs pulse detuning.]{
    Numerical simulation of a single quantum coherence under pulsed excitation
    ($\Gamma_{10}\Delta_t=1$) at different detunings (labelled inset).
    The coherence is shown in both the time (left column) and frequency (right column) domain.
    The color indicates the frequency (scale on right), while the black line shows the amplitude
    profile. The excitation profile is shown as a grey line.
  }
\end{figure}

\begin{figure}
	\includegraphics[width=\textwidth]{"mixed_domain/SQC lineshapes against t"}
	\caption[Time-gated amplitude of a single quantum coherence vs pulse detuning.]{
    Amplitude of a single quantum coherence under pulsed excitation as a function of detuning (x
    axis) and delay after excitation (line color, scale on right) for the three
    $\Gamma_{10}\Delta_t$ values considered in this work.
    For comparison, the excitation pulse lineshape (grey), absorptive material response (black,
    narrow) and magnitude material response (black, wide) are also shown. Each peak is normalized
    to its own maximum amplitude.
  }
	\label{mix:fig:sqc_vs_t}
\end{figure}

\subsection{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, signal that time-gates FID evolution gives narrower spectra than 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 signal at the fully driven frequency $\omega_\text{out} = \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_\text{out}=\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 \autoref{mix:fig:WMELs}).  %
\autoref{mix:fig:pw1} shows the resulting 2D delay profile of pathway I$\gamma$ signals for
$\Gamma_{10}\Delta_t=1$ (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 FID signal can change the resonance conditions; \autoref{mix:tab:table2} 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.  %
\autoref{mix:fig:pw1_no_mono} shows a simulation of \autoref{mix:fig:pw1} without monochromator
frequency filtering.  %

When the pulses are all overlapped ($\tau_{21}=\tau_{22^\prime}=0$, 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=\omega_{10}$, $\omega_1-\omega_2=0$, and
$\omega_1-\omega_2+\omega_{2^\prime}=\omega_{10}$, respectively.  %
The second resonance condition causes elongation along the diagonal, and since
$\omega_2=\omega_{2^\prime}$, 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=\omega_{10}$.  %

The other three spectra of \autoref{mix:fig:pw1} separate the pulse sequence over time so that not all
interactions are driven.  %
At $\tau_{21}=0$, $\tau_{22^\prime}=-2.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, $E_{2^\prime}$, is latent and probes $\rho_2$ during its FID evolution after
memory of the driven frequency is lost.  %
There are two important consequences.  %
Firstly, the third driven resonance condition is now approximated by
$\omega_{2^\prime}=\omega_{10}$, which makes $\omega_1$ only singly resonant at
$\omega_1=\omega_{10}$.  %
Secondly, the driven portion of the signal frequency is determined only by the latent pulse:
$\omega_{\text{out}}=\omega_{2^\prime}$.  %
Since our monochromator gates $\omega_1$, we have the detection-induced correlation
$\omega_1=\omega_{2^\prime}$.  %
The net result is double resonance along $\omega_1=\omega_2$, and the vertical elongation of pulse
overlap is strongly attenuated.  %

At $\tau_{21}=2.4\Delta_t,\tau_{22^\prime}=0$ (upper right, purple), the first pulse $E_1$ precedes
the latter two, which makes the two resonance conditions for the input fields
$\omega_1=\omega_{10}$ and $\omega_2=\omega_{10}$.  %
The signal depends on the FID conversion of $\rho_1$, which gives vertical elongation at
$\omega_1=\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 $\omega_1=\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 $E_{2^\prime}$ are coincident, so that the final coherence can
be approximated as driven by these two, we can approximate the final frequency as
$\omega_{\text{out}} = \omega_{10}-\omega_2+\omega_{2^\prime} = \omega_{10}$.  %
Surprisingly, the frequency content of the output is strongly independent of all pulse
frequencies.  %
The monochromator narrows the $\omega_1=\omega_{10}$ resonance.  %
The $\omega_1=\omega_{10}$ resonance condition now depends on the monochromator slit width, the FID
propagation of $\rho_1$, 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}=-\tau_{22^\prime}=2.4\Delta_t$, upper left,
cyan), each resonance condition is independent and both $E_1$ and $E_2$ require FID buildup to
produce 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 $E_{2^\prime}$, and since the
monochromator track $\omega_1$, the output signal will increase when
$\omega_1=\omega_{2^\prime}$.  %
As a result, the line shape acquires a diagonal character.  %

The changes in line shape seen in \autoref{mix:fig:pw1} 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 \cite{MeyerKentA2004a, PakoulevAndreiV2006a, DonaldsonPaulMurray2007b}, 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 functional of each other
in unexpected ways.  %

\begin{figure}
	\includegraphics[width=\linewidth]{"mixed_domain/pw1 lineshapes"}
	\caption[2D frequency response of a single Liouville pathway at different delay values.]{
		Changes to the 2D frequency response of a single Liouville pathway (I$\gamma$) at different
    delay values.
    The normalized dephasing rate is $\Gamma_{10}\Delta_t = 1$.
		Left:  The 2D delay response of pathway I$\gamma$ at triple resonance.
		Right: The 2D frequency response of pathway I$\gamma$ at different delay values.
    The delays at which the 2D frequency plots are collected are indicated on the delay plot;
    compare 2D spectrum frame color with dot color on 2D delay plot.
		}
	\label{mix:fig:pw1}
\end{figure}

\begin{figure}
	\includegraphics[width=\textwidth]{"mixed_domain/pw1 lineshapes no mono"}
	\caption[2D frequency response of a single Liouville pathway without a tracking monochromator.]{
    Pathway \RomanNumeral{1}$\gamma$ temporal response in the 2D pulse delay space at triple
    resonance (left) and the corresponding 2D frequency plots at different delay values.
    The delays at which the 2D frequency plots are collected are indicated on the delay plot;
    compare 2D spectrum frame color with dot color on 2D delay plot.
    Unlike elsewhere in this work, signal here was not filtered by a tracking monochromator.
  }
 	\label{mix:fig:pw1_no_mono}
\end{figure}

\begin{table}
	\begin{tabular}{c c | c c c c}
		\multicolumn{2}{c}{Delay} & \multicolumn{4}{|c}{Approximate Resonance Conditions} \\
		$\tau_{21}/\Delta_t$ & $\tau_{22^\prime}/\Delta_t$ & $\rho_0\xrightarrow{1}\rho_1$ &
    $\rho_1\xrightarrow{2}\rho_2$ & $\rho_2\xrightarrow{2^\prime}\rho_3$ & $\rho_3\rightarrow$
    detection at $\omega_m=\omega_1$ \\
		\hline\hline
		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$ \\
	\end{tabular}
  \caption{
    Conditions for peak intensity at different pulse delays for pathway I $\gamma$.
  }
  \label{mix:tab:table2}
\end{table}

\subsection{Temporal pathway discrimination}  % ---------------------------------------------------

In the last section 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 \autoref{mix:fig:WMELs}).  %
Each time-ordering has different resonance conditions.  %
When 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.  %

\autoref{mix:fig:delay_purity} shows the 2D delay space with all pathways present for
$\omega_1,\omega_2=\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.  %
Signal is symmetric about the $\tau_{21}=\tau_{22^\prime}$ line because when $\omega_1=\omega_2$,
$E_1$ and $E_{2^\prime}$ interactions are interchangeable:
$S_\text{tot}(\tau_{21},\tau_{22^\prime})=S_\text{tot}(\tau_{22^\prime}, \tau_{21})$.  %
The overlaid black contours represent signal ``purity,'' $P$, defined as the relative amount of
signal that comes from the dominant pathway at that delay value:
\begin{equation} \label{mix:eqn:P}
P(\tau_{21},\tau_{22^\prime})=\frac{\max \left\{S_L\left( \tau_{21},\tau_{22^\prime} \right)\right\}}
{\sum_L S_L\left( \tau_{21},\tau_{22^\prime} \right)}.
\end{equation}
The dominant pathway ($\max{\left\{ S_L \left( \tau_{21},\tau_{22^\prime} \right) \right\}}$) at
given delays can be inferred by the time-ordered region defined in \autoref{mix:fig:overview}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
$\left|\tau_{22^\prime}\right|/\Delta_t < 1$ or $\left|\tau_{21}\right|/\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}\rightarrow \infty$, however, such discrimination is only achieved at
vanishing signal intensities; the contour of $P=0.99$ across our systems highlights this trend.  %

\begin{figure}
	\includegraphics[width=\linewidth]{"mixed_domain/delay space ratios"}
	\caption[2D delay response for different relative dephasing rates.]{
		Comparison of the 2D delay response for different relative dephasing rates (labeled atop each
    column).
		All pulses are tuned to exact resonance.
		In each 2D delay plot, the signal amplitude is depicted by the colors.
		The black contour lines show signal purity, $P$ (see \autoref{mix:eqn:P}), with purity values
    denoted on each contour.
		The small plots above each 2D delay plot examine a $\tau_{22^\prime}$ slice of the delay
    response ($\tau_{21}=0$).
    The plot shows the total signal (black), as well as the component time-orderings VI (orange), V
    (purple), and III or I (teal).
	}
	\label{mix:fig:delay_purity}
\end{figure}

\subsection{Multidimensional line shape dependence on pulse delay time}  % ------------------------

In the previous sections 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 $\vec{k}_2$ and $\vec{k}_{2^\prime}$ interactions correspond to the pump, and the $\vec{k}_1$
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}=0$, all time-orderings contribute
equally; at $\tau_{21}<0$ time-orderings V and VI dominate (\autoref{mix:fig:delay_purity}).  %
Conventional pump-probe techniques recognized these complications long ago, \cite{BritoCruzCH1988a,
  PalfreySL1985a} but the extension of these effects to MR-CMDS has not previously been done.  %

\autoref{mix:fig:hom_2d_spectra} 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.  %
The dependence of the line shape amplitude on $\tau_{21}$ can be inferred from
\autoref{mix:fig:delay_purity}.  %

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 $E_{2^\prime}$.  %
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 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 always 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 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^\prime}$ 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=2$ decreases in strength as $\Gamma_{10}\Delta_t$
decreases.  %
As we illustrated in \autoref{mix:fig:pw1}, this spectral correlation is a signature of driven signal
from 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 \autoref{mix:fig:delay_purity}.  %
When all pulses are completely overlapped, ($\tau_{21}=0$), each of the line shapes exhibit
spectral correlation.  %
At $\tau_{21}/\Delta_t=-2$, 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 elongation direction is due to the resolving power of $\omega_2$.  %
At $\Gamma_{10}\Delta_t=0.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=0.5$.  %

It is also common to represent data as ``Wigner plots,'' where one axis is delay and the other is
frequency. \cite{KohlerDanielDavid2014a, AubockGerald2012a, CzechKyleJonathan2015a,
  PakoulevAndreiV2007a}  %
In \autoref{mix:fig:wigners} we show five $\tau_{21},\omega_1$ plots for varying $\omega_2$ with
$\tau_{22^\prime}=0$.  %
The plots are the analogue 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.  %

\begin{figure}
	\includegraphics[width=\linewidth]{"mixed_domain/spectral evolution"}
	\caption[Evolution of the 2D frequency response.]{
		Evolution of the 2D frequency response as a function of $\tau_{21}$ (labeled inset) and the
    influence of the relative dephasing rate ($\Gamma_{10}\Delta_t=0.5$ (red), $1.0$ (green), and
    $2.0$ (blue)).
		In all plots $\tau_{22^\prime}=0$.  To ease comparison between different dephasing rates, the
    colored line contours (showing the half-maximum) for all three relative dephasing rates are
    overlaid.
		The colored histograms below each 2D frequency plot show the relative weights of each
    time-ordering for each relative dephasing rate.
		Contributions from V and VI are grouped together because they have equal weights at
    $\tau_{22^\prime}=0$.
	}
	\label{mix:fig:hom_2d_spectra}
\end{figure}

\begin{figure}
	\includegraphics[width=\textwidth]{"mixed_domain/spectral evolution full"}
	\caption[Evolution of the 2D frequency response, with all contours shown.]{
    Spectral evolution of the homogeneous exciton resonance as a function of $\tau_{21}$, with
    $\tau_{22^\prime}=0$.
    The 50\% contour is darkened to ease comparison with Figure 7.
  }
	\label{mix:fig:spectral_evolution_full}
\end{figure}

\begin{figure}
	\includegraphics[width=\linewidth]{"mixed_domain/wigners full"}
	\caption[Simulated Wigner spectra.]{
    Mixed $\tau_{21}$, $\omega_1$ plots for each $\Gamma_{10}$ value simulated in this work.
    For each plot, the corresponding $\omega_2$ value is shown as a gray vertical line.
    Each plot is separately normalized.
  }
	\label{mix:fig:wigners}
\end{figure}

\subsection{Inhomogeneous broadening} \label{mix:sec:res_inhom}  % --------------------------------

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 signal to time-ordered regions III and V that
persists for all population times.  %
\autoref{mix:fig:delay_inhom} shows the calculated spectra for relative dephasing rate
$\Gamma_{10}\Delta_t=1$ with a frequency broadening function of width
$\Delta_{\text{inhom}}=0.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. \cite{WeinerAM1985a, FlemingGrahmR1998a, DeBoeijWimP1998a, SalvadorMayroseR2008a}  %
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=0$ reflects the rephasing
ability of the system.  %
An inhomogeneous system will emit a photon echo in the rephasing pathway, enhancing 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 $E_2$ and $E_{2^\prime}$ create
the population (time-orderings V and VI).  %
The trace runs parallel to the III-V time-ordering boundary (diagonal) if $\tau_{22^\prime}<0$ and
runs parallel to the IV-VI time-ordering boundary (horizontal) if $\tau_{22^\prime}>0$, and both
intersect at $\tau_{22^\prime}=0$; the $-\tau_{21}$ value at this intersection is $T$.  %
In our 2D delay plots (\autoref{mix:fig:delay_purity}, \autoref{mix:fig:delay_inhom}), the peak
shift is seen as the diagonal displacement of the signal peak from the $\tau_{21}=0$ vertical
line.  %
\autoref{mix:fig:delay_inhom} highlights the peak shift profile as a function of population time with
the yellow trace; it is easily verified that our static inhomogneous 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 ,and merely reflects $E_1$ and $E_2$ temporal overlap.  %
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 signal away from $\tau=0$ 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 homogeneous system (Fig.
\autoref{mix:fig:delay_purity}); 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, \cite{DeBoeijWimP1996a, AgarwalRitesh2002a} 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\neq\omega_2$ because
time-ordering III is decoupled by detuning.  %

In frequency space, spectral elongation along the diagonal is the signature of inhomogeneous
broadening.  %
\autoref{mix:fig:inhom_2d_spectra} shows the line shape changes of a Gaussian inhomogeneous
distribution.  %
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. \autoref{mix:fig:inhom_2d_spectra},
$\tau_{21}=2.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 width.  %

\begin{figure}
	\includegraphics[width=\textwidth]{"mixed_domain/3PEPS"}
	\label{mix:fig:3PEPS}
	\caption[3PEPS tutorial.]{
    Extraction of 3PEPS peak shifts from MR-CMDS delay space. Left-hand plot: thick colored lines
    denote contours of constant $\tau$ for $T=0, 1, 2, 3$.
    Dots indicate the fitted peak shift for each $\tau$ contour.
    Right-hand plot: numerically simulated amplitude traces (solid), Gaussian fits (transparent)
    and fit centers (vertical lines) for each $T$ (colors matched).
  }
\end{figure}

\begin{figure}
	\includegraphics[width=0.5\linewidth]{"mixed_domain/inhom delay space ratios"}
	\caption[2D delay response with inhomogeneity.]{
		2D delay response for $\Gamma_{10}\Delta_t=1$ with sample inhomogeneity.  %
		All pulses are tuned to exact resonance.  %
		The colors depict the signal amplitude.  %
		The black contour lines show signal purity, $P$ (see \autoref{mix:eqn:P}), with purity values
    denoted on each contour.  %
		The thick yellow line denotes the peak amplitude position that is used for 3PEPS analysis.  %
		The small plot above each 2D delay plot examines a $\tau_{22^\prime}$ slice of the delay
    response ($\tau_{21}=0$).  %
		The plot shows the total signal (black), as well as the component time-orderings VI (orange), V
    (purple), III (teal, dashed), and I (teal, solid).  %
	}
	\label{mix:fig:delay_inhom}
\end{figure}

\begin{figure}
	\includegraphics[width=\textwidth]{"mixed_domain/2D delays"}
	\label{mix:fig:2D_delays}
	\caption[2D delay response for all combinations of inhomogeneity, dephasing rate.]{
    2D delay scans at $\omega_1=\omega_2=\omega_{10}$ for all 12 combinations of $\Gamma_{10}$
    (rows) and $\Delta_{inhom}$ (columns) simulated in this work.
    The 3PEPS shift trace is plotted in yellow, annotated to indicate the magnitude of the $\tau$
    shift at $T=0$ and $T=4\Delta_t$.
  }
\end{figure}

\begin{figure}
	\includegraphics[width=\linewidth]{"mixed_domain/inhom spectral evolution"}
	\caption[Spectral evolution of an inhomogenious system.]{
		Same as \autoref{mix:fig:hom_2d_spectra}, but each system has inhomogeneity
    ($\Delta_{\text{inhom}}=0.441\Gamma_{10}$).
		Relative dephasing rates are $\Gamma_{10}\Delta_t=0.5$ (red), $1.0$ (green), and $2.0$ (blue).
		In all plots $\tau_{22^\prime}=0$.
		To ease comparison between different dephasing rates, the colored line shapes of all three
    systems are overlaid.
		Each 2D plot shows a single representative contour (half-maximum) for each
    $\Gamma_{10}\Delta_t$ value.
		The colored histograms below each 2D frequency plot show the relative weights of each
    time-ordering for each 2D frequency plot.
		In contrast to \autoref{mix:fig:hom_2d_spectra}, inhomogeneity makes the relative contributions of
    time-orderings V and VI unequal.
	}
	\label{mix:fig:inhom_2d_spectra}
\end{figure}

\begin{figure}
	\includegraphics[width=\textwidth]{"mixed_domain/inhom spectral evolution full"}
	\label{mix:fig:inhom_spectral_evolution_full}
	\caption[Spectral evolution of an inhomogeneous system, with all contours shown.]{
    Spectral evolution of the exciton resonance as a function of $\tau_{21}$, with
    $\tau_{22^\prime}=0$.
    For each system $\Delta_{inhom}=0.441\Gamma_{10}$.
    The 50\% contour is darkened to ease comparison with Figure 10.
  }
\end{figure}

\begin{figure}
	\includegraphics[width=\textwidth]{"mixed_domain/2D frequences at zero"}
	\label{mix:fig:2D_frequencies_at_zero}
	\caption[Eccentricity at zero delay.]{
    2D frequency scans at $\tau_{21}=\tau_{22^\prime}=0$ for all 12 combinations of $\Gamma_{10}$
    (columns) and $\Delta_{inhom}$ (rows) simulated in this work.
    The eccentricity of each spectrum is inset and represented by the yellow ellipse (50\%
    contour).
  }
\end{figure}

\begin{figure}
	\includegraphics[width=\textwidth]{"mixed_domain/2D frequences at -4"}
	\label{mix:fig:2D_frequencies_at large population time.at_-4}
	\caption[Eccentricity at large population time.]{
    2D frequency scans at large $T$ ($\tau_{22^\prime}=0$, $\tau_{21}=-4\Delta_t$) for all 12
    combinations of $\Gamma_{10}$ (columns) and $\Delta_{inhom}$ (rows) simulated in this work.
    The eccentricity of each spectrum is inset and represented by the yellow ellipse (50\%
    contour).
  }
\end{figure}

\section{Discussion}  % ---------------------------------------------------------------------------

\subsection{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 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:
\begin{enumerate}
	\item When time-gating during the pulse, the system pins to the driving frequency with a buildup efficiency determined by resonance.
	\item When time-gating after the pulse, the FID dominates the system response.
	\item The emitted signal field contains both FID and driven components; the $\omega_{\text{out}} = \omega_1$ component is isolated by the tracking monochromator.
\end{enumerate}
\autoref{mix:fig:fid_dpr} illustrates principles 1 and 2 and \autoref{mix:fig:fid_detuning}
illustrates principle 2 and 3.  %
\autoref{mix:fig:pw1} 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 \autoref{mix:fig:delay_purity}.  %
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 Figures \ref{mix:fig:hom_2d_spectra} and
\ref{mix:fig:inhom_2d_spectra}.  %

\subsection{Conditional validity of the driven limit}

We have shown that the driven limit misses details of the line shape if $\Gamma_{10} \Delta_t
\approx 1$, 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.  %
Fig. \ref{mix:fig:steady_state} compares the results of our numerical simulation (third column) with
the driven limit expressions for populations where $\Gamma_{11}\Delta_t=0$ (first column) or $1$
(second column).  %
The top and bottom rows compare the line shapes when $\left(\tau_{22^\prime},
  \tau_{21}=(0,0)\right)$ and $(0,-4\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.  %

% TODO:  [ ] population resonance is not clear
Note the very sharp diagonal feature that appears for $(\tau_{21},\tau_{22^\prime}) = (0,0)$ and $
\Gamma_{11}=0$; 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
(\textcite{YursLenaA2011a} and \textcite{KohlerDanielDavid2014a}, 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}=0$).  %
In light of this, one can approximate this time-gating effect by substituting population lifetime
with the pulse duration ($\Gamma_{11}\Delta_t=1$), which gives good agreement with the numerical
simulation (third column).  %

When $\tau_{22^\prime}=0$ and $\tau_{21}<\Delta_t$, signals can also be approximated by driven
signal (\autoref{mix:fig:steady_state} bottom row).  %
Only time-orderings V and VI are relevant.  %
The intermediate population resonance is still impulsive but it depends on
$\omega_{2^\prime}-\omega_2$ which is not explored in our 2D frequency space.  %

\begin{figure}
	\includegraphics[width=\linewidth]{"mixed_domain/steady state"}
	\caption[Conditional validity of the driven limit.]{
		Comparing approximate expressions of the 2D frequency response with the directly integrated
    response.  %
		$\Gamma_{10}\Delta_t=1$.  %
		The top row compares the 2D response of all time-orderings ($\tau_{21}=0$) and the bottom row
    compares the response of time-orderings V and VI ($\tau_{21}=-4\Delta_t$).  %
		First column:  The driven limit response.  Note the narrow diagonal resonance for $\tau_{21}=0$.
		Second column:  Same as the first column, but with ad hoc substitution $\Gamma_{11}=\Delta_t$.
		Third column:  The directly integrated response.  %
	}
	\label{mix:fig:steady_state}
\end{figure}

\subsection{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 for the frequency domain \cite{KwacKijeong2003a, OkumuraKo1999a}.  %
There are many ways to characterize the ellipticity of a peak shape.  %
We adopt the convention $\mathcal{E} = \left(a^2-b^2\right) / \left(a^2+b^2\right)$, where $a$ is
the diagonal width and $b$ is the antidiagonal width.  %
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.  %

\autoref{mix:fig:metrics} shows the results of this characterization for all $\Delta_\text{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_\text{inhom} / \Delta_\omega$, the
relative inhomogeneity $\Delta_\text{inhom} / \Gamma_{10}$, and the population time delay.  %
The top row shows the correlations of the $\Delta_\text{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 = 4\Delta_t$ that
isolates the V and VI time-orderings.  %
For this time delay, the $\Delta_\text{3PEPS} / \Delta_t$ metric works well for all dephasing times
of $\Gamma_{10}\Delta_t$ when the relative inhomogeneity is $\Delta_\text{inhom} / \Delta_\omega
\ll 1$.  %
It becomes independent of $\Delta_\text{inhom} / \Delta_\omega$ when $\Delta_\text{inhom} /
\Delta_\omega > 1$.  %
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. \cite{WeinerAM1985a}  %
The corresponding graph for $T = 0$ shows a large peak shift occurs, even without inhomogeneity.
In this case, the peak shift depends on pathway overlap, as discussed in
\autoref{mix:sec:res_inhom}.  %

The middle row in \autoref{mix:fig:metrics} shows the ellipticity dependence on the relative dephasing
rate and inhomogeneity assuming the measurement is performed when the first two pulses are
temporally overlapped ($\tau_{22^\prime}=0$).  %
For a $T=4\Delta_t$ population time, the ellipticity is proportional to the inhomogeneity until
$\Delta_\text{inhom} / \Delta_\omega \ll 1$ where the excitation bandwidth is wide compared with
the inhomogeneity.  %
Unlike 3PEPS, saturation is not observed because pulse bandwidth does not limit the frequency range
scanned.  %
The 3PEPS and ellipticity metrics are therefore complementary since 3PEPS works well for
$\Delta_\text{inhom} / \Delta_\omega \ll 1$ and ellipticity works well for $\Delta_\text{inhom} /
\Delta_\omega \gg 1$.  %
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, 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 \autoref{mix:fig:metrics} 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 = 4\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 at $T = 4\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. \cite{AartsmaThijsJ1976a}  %
Both metrics are offset due to the pulse overlap effects.  %
Accordingly, the region to the left of homogeneous contour is non-physical, because it represents
observed correlations that are less than that 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.  %

\begin{figure}
	\includegraphics[width=0.5\linewidth]{"mixed_domain/metrics"}
	\caption[Metrics of correlation.]{
		Temporal (3PEPS) and spectral (ellipticity) metrics of correlation and their relation to the
    true system inhomogeneity.  %
		The left column plots the relationship at pulse overlap ($T=0$) and the right column plots the
    relationship at a delay where driven correlations are removed ($T=4\Delta_t$).  %
    For the ellipticity measurements, $\tau_{22^\prime}=0$.  %
		In each case, the two metrics are plotted directly against system inhomogeneity (top and middle
    row) and against each other (bottom row).  %
		Colored lines guide the eyes for systems with equal relative dephasing rates
    ($\Gamma_{10}\Delta_t$, see upper legend), while the area of the data point marker indicates
    the relative inhomogeneity ($\Delta_{\text{inhom}}/\Gamma_{10}$, see lower legend).  %
		Gray lines indicate contours of constant relative inhomogeneity (scatter points with the same
    area are connected).  %
}
	\label{mix:fig:metrics}
\end{figure}

\section{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 response.  %
The dimensionless quantity $\Delta_t \left(\Gamma_f + i\kappa_f \Omega_{fx}\right)$ 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 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 frequency or time domain, pulse effects mimic signatures of
ultrafast inhomogeneity.  %
Even homogeneous systems take on these signatures.  %
For mixed domain experiments, pulse effects induce spectral ellipticity and photon echo signatures,
even in homogeneous systems.  %
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-\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 domain.  %