\chapter{Disentangling material and instrument response}

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{Rentzepis1970,Mukamel2000}  % 
The ultrafast specta can be collected in the time domain or the frequency domain.\cite{Park1998}  %

Time-domain methods scan the pulse delays to resolve the free induction decay
(FID).\cite{Gallagher1998}   %
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{Lebedev2007, Vardeny1981, Joffre1988, Pollard1992}  % 
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{Druet1979,Oudar1980}  % 
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{Oudar1980,Wright1997a,Wright1991}  %
In this regime, both FID and driven processes are important.\cite{Pakoulev2006}  %
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{Pakoulev2007,Kohler2014,Gelin2009a}   %
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{Hamm2000} transient grating,\cite{Salcedo1978,Fourkas1992,Fourkas1992a} transient
absorption/reflection,\cite{Aubock2012,Bakker2002} photon
echo,\cite{DeBoeij1996,Patterson1984,Tokmakoff1995} two dimensional-infrared spectroscopy
(2D-IR),\cite{Hamm1999,Asplund2000,Zanni2001} 2D-electronic spectroscopy
(2D-ES),\cite{Hybl2001a,Brixner2004} and three pulse photon echo peak shift
(3PEPS)\cite{Emde1998,DeBoeij1996,DeBoeij1995,Cho1992,Passino1997} spectroscopies.  %
These methods also include fully-coherent methods involving only coherences such as Stimulated
Raman Spectroscopy (SRS),\cite{Yoon2005,McCamant2005} Doubly Vibrationally Enhanced
(DOVE),\cite{Zhao1999,Zhao1999a,Zhao2000,Meyer2003,Donaldson2007,Donaldson2008,Fournier2008} Triply
Resonant Sum Frequency (TRSF),\cite{Boyle2013a,Boyle2013,Boyle2014} Sum Frequency Generation
(SFG)\cite{Lagutchev2007}, Coherent Anti-Stokes Raman Spectroscopy
(CARS)\cite{Carlson1990b,Carlson1990a,Carlson1991} and other coherent Raman
methods\cite{Steehler1985}.  %

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{Kwac2003,Wright2016}  % 
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{Erlik2017,Smallwood2016}  %

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{Gallagher1998,Fourkas1992,Fourkas1992a} many of the experiments involve pulse widths
that are comparable to the widths of the
system.\cite{Meyer2003,Donaldson2007,Pakoulev2009,Zhao1999,Czech2015,Kohler2014}  % 
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{Wright2003} whose relative
importance depends on the relative importance of FID and driven responses.\cite{Donaldson2010}  %
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{Kohler2014,Czech2015} and have been developed on vibrational states as
well.\cite{Meyer2004}  %
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}

\begin{figure}
  \centering
	\includegraphics[width=0.5\linewidth]{"mixed_domain/WMELs"}
	\caption{
		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{fig:WMELs}
\end{figure}

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{eq: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{eq: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 Equation \ref{eq:LVN} to third order in the
electric field
interaction\cite{mukamel1995principles,Yee1978,Oudar1980,Armstrong1962,Schweigert2008} 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.  %
Fig. \ref{fig:WMELs} shows these sixteen pathways as Wave Mixing Energy Level (WMEL)
diagrams\cite{Lee1985}.  %

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{eq: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$).\footnote{$\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$).  %

Equation \ref{eq: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{eq: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.  %
Equation \ref{eq: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 \ref{eq:rho_f_int} are discussed in Appendix
\ref{sec:cw_imp}.  %

\begin{figure*}
	\includegraphics[width=\linewidth]{"mixed_domain/simulation overview"}
	\caption{
		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 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{fig:overview}
\end{figure*}

Fig. \ref{fig:overview} gives an overview of the simulations done in this work.  %
Fig. \ref{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 Fig. \ref{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.  %
Fig. \ref{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 (Fig. \ref{fig:overview}c).  %
The resulting signal is explored in 2D delay space (Fig. \ref{fig:overview}d), 2D frequency space
(Fig. \ref{fig:overview}f), and hybrid delay-frequency space (Fig. \ref{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.  %
Fig. \ref{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$.\footnote{For elucidation of the
  relationship between the generalized Liouville pathway notation and the specific parameters for
  each Liouville pathway, see Table S1 in the Supplementary Information.}  %

The electric field emitted from a Liouville pathway is proportional to the polarization created by
the third-order coherence:  %
\begin{equation}\label{eq:E_L}
E_L(t) = i \mu_{4}\rho_{3}(t).
\end{equation}
Equation \ref{eq:E_L} assumes perfect phase-matching and no pulse distortions through propagation.  Equation \ref{eq:rho_f_int} shows that the output field for this Liouville pathway is
	\begin{gather}\label{eq: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{eq:superposition}
E_{\text{tot}}= \sum_L E_L(t).
\end{equation}
For the superposition of Equation \ref{eq: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 Section
\ref{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{eq: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 Fig.
\ref{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 Fig. \ref{fig:overview}b).  %
Table S1 summarizes the arguments for each Liouville pathway.  %
Fig. \ref{fig:overview}f shows the 2D $(\omega_1, \omega_2)$ $S_{\text{tot}}$ spectrum resulting
from the pulse delay times represented in Fig. \ref{fig:overview}b.  %

\subsection{Inhomogeneity}

Inhomogeneity is isolated in CMDS through both spectral signatures, such as
line-narrowing\cite{Besemann2004,Oudar1980,Carlson1990,Riebe1988,Steehler1985}, and temporal
signatures, such as photon echoes\cite{Weiner1985,Agarwal2002}.  %
We simulate the effects of static inhomogeneous broadening by convolving the homogeneous response
with a Gaussian distribution function.  %
Further details of the convolution are in Appendix \ref{sec:convolution}.  %
Dynamic broadening effects such as spectral diffusion are beyond the scope of this work.  %

\section{Methods}  % ------------------------------------------------------------------------------

A matrix representation of differential equations of the type in Equation \ref{eq:E_L_full} was
numerically integrated for parallel computation of Liouville elements (see SI for
details).\cite{Dick1983,Gelin2005}  %
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{Oliphant2007}  %

\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{sec:evolution_SQC}

\begin{figure}
	\includegraphics[width=0.5\linewidth]{"mixed_domain/fid vs dpr"}
	\caption{
		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{fig:fid_dpr}
\end{figure}

It is illustrative to first consider the evolution of single coherences, $\rho_0 \xrightarrow{x}
\rho_1$, under various excitation conditions.  %
Fig. \ref{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 Equation \ref{eq:rho_f_int} under the various conditions in Table S1), so
the profiles in Fig. \ref{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 Appendix \ref{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 (Equation \ref{eq:sqc_rise}), but as the pulse vanishes, the
signal adopts the natural system frequency and decay rate (Equation \ref{eq:sqc_fid}).  %
For $\Gamma_{10}\Delta_t=\infty$, the signal follows the amplitude and frequency of the pulse for
all times (the driven limit, Equation \ref{eq: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 (Equation
\ref{eq: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$.  %