2018-05-10 11:08

1 files changed, 770 insertions, 88 deletions

diff --git a/mixed_domain/chapter.tex b/mixed_domain/chapter.tex
index 75e60ea..8128ba5 100644
--- a/mixed_domain/chapter.tex
+++ b/mixed_domain/chapter.tex
@@ -23,7 +23,7 @@ domain.  %
 

 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}  % 

+themselves. \cite{RentzepisPM1970a, MukamelShaul2000a}  %

 The ultrafast specta can be collected in the time domain or the frequency

 domain. \cite{ParkKisam1998a}  %

 

@@ -31,21 +31,21 @@ 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.  % 

+measurement.  %

 FID occurs at the frequency of the transition that has been excited by a well-defined, time-ordered

-sequence of pulses.  % 

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

+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}  % 

+  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}  % 

+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

@@ -53,7 +53,7 @@ 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.  % 

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

 

@@ -96,7 +96,7 @@ Multidimensional Spectroscopy (CMDS) as experiments transition between the two l
 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}  % 

+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}  %

@@ -115,15 +115,15 @@ 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.  % 

+broadening.  %

 In the FID limit where the coherence frequencies are defined instead by the transition, the

-diagonal character is lost.  %  

+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'}$.  % 

+$\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

@@ -160,7 +160,7 @@ from these measurement artifacts.  %
 

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

+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:

@@ -204,7 +204,7 @@ $\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{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}

@@ -227,7 +227,7 @@ The integral solution is
 \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, 

+& \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.  %

@@ -243,19 +243,19 @@ this paper.  %
 \begin{figure}

 	\includegraphics[width=\linewidth]{"mixed_domain/simulation overview"}

 	\caption[Overview of the MR-CMDS simulation.]{

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

+    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.   

+    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.  

+    help introduce our delay convention.

 	}

 	\label{mix:fig:overview}

 \end{figure}

@@ -296,10 +296,10 @@ E_L(t) = i \mu_{4}\rho_{3}(t).
 \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_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	,	

+	&\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}

@@ -330,7 +330,7 @@ In this simulation, the detection is emulated by transforming $E_{\text{tot}}(t)
 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}) = 

+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

@@ -343,39 +343,6 @@ 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.  %

 

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

-Further details of the convolution are in \autoref{mix: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 \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}  %

-

 \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

@@ -386,7 +353,7 @@ 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} 

+\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,

@@ -410,7 +377,7 @@ This is classic free induction decay (FID) evolution:  the system evolves at its
 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.  

+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 -

@@ -445,10 +412,10 @@ 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 

+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_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) \\

@@ -493,7 +460,16 @@ The build-up limit approximates well when pulses are near-resonant and arrive to
 build-up behavior is emphasized).  %

 The driven limit holds for large detunings, regardless of delay.  %

 

-\subsection{Convolution technique for inhomogeneous broadening} \label{mix:sec:convolution}  % ----

+\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.  %

+Further details of the convolution are in \autoref{mix:sec:convolution}.  %

+Dynamic broadening effects such as spectral diffusion are beyond the scope of this work.  %

 

 \begin{figure}

 	\includegraphics[width=\linewidth]{mixed_domain/convolve}

@@ -582,6 +558,587 @@ 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}

+

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

+

+\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{tab:table1} describes the relationship between our notation and the parameters that make

+up the 16 Liouville pathways.  %

+

+\begin{table*}[!hb]

+	\caption{\label{tab:table1} Parameters of each Liouville Pathway.}

+	\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}

+\end{table*}

+

+\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{GelinMaximF2009.000}.  %

+

+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{eq: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 \textbf{TODO}, 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{

+    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{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{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{section: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{

+    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{center}

+	\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}

+\end{center}

+

+$\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{eq: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:

+

+\begin{center}

+	\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}

+\end{center}

+

+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:

+

+\begin{center}

+	\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}

+\end{center}

+

+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.

+

+\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:

+

+\begin{center}

+	\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}

+\end{center}

+

+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.

+

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

 

 We now present portions of our simulated data that highlight the dependence of the spectral line

@@ -647,24 +1204,49 @@ FID character is difficult to isolate when $\Gamma_{10}\Delta_t=2.0$.  %
     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).		

+    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.  

+    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). 

+    %(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. 

+    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{

+    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"}

+	\label{fig:sqc_vs_t}

+	\caption{

+    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.

+  }

+\end{figure}

+

 \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$.

-

+Figure \textbf{TODO} 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

@@ -678,11 +1260,31 @@ 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.  %

+for $\Omega_{fx}/\Delta_{\omega}=\pm2$ manifests as a visibly disproportionate amplitude decay

+(\textbf{TODO} 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.  %

@@ -736,8 +1338,8 @@ $\Gamma_{10}\Delta_t=1$.  %
 	\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.  

+    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.

@@ -745,6 +1347,18 @@ $\Gamma_{10}\Delta_t=1$.  %
 	\label{mix:fig:pw1}

 \end{figure}

 

+\begin{figure}

+	\includegraphics[width=\textwidth]{"mixed_domain/pw1 lineshapes no mono"}

+	\label{mix:fig:pw1_no_mono}

+	\caption{

+    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.

+  }

+\end{figure}

+

 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

@@ -768,7 +1382,6 @@ 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

@@ -780,6 +1393,8 @@ The prominence of FID signal can change the resonance conditions; \autoref{mix:t
 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.  %

+\textbf{TODO} shows a simulation of \autoref{mix:fig:pw1} without monochromator frequency

+filtering.  %

 

 \begin{table}

 	\caption{\label{mix:tab:table2} Conditions for peak intensity at different pulse delays for pathway

@@ -790,7 +1405,7 @@ considered.  %
     $\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   & 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}$ & --                     &

@@ -813,7 +1428,7 @@ $\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.  % 

+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

@@ -877,10 +1492,10 @@ in unexpected ways.  %
 	\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.  

+		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.  

+    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

@@ -954,6 +1569,16 @@ vanishing signal intensities; the contour of $P=0.99$ across our systems highlig
 	\label{mix:fig:hom_2d_spectra}

 \end{figure}

 

+\begin{figure}

+	\includegraphics[width=\textwidth]{"mixed_domain/spectral evolution full"}

+	\label{mix:fig:spectral_evolution_full}

+	\caption{

+    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.

+  }

+\end{figure}

+

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

@@ -1028,11 +1653,12 @@ This narrowing, however, is unresolvable when the pulse bandwidth becomes broade
 resonance, which gives rise to a vertically elongated signal when $\Gamma_{10}\Delta_t=0.5$.  %

 

 \begin{figure}

-	\includegraphics[width=\linewidth]{"mixed_domain/wigners"}

-	\caption[Wigners.]{

-		Transient ($\omega_1$) line shapes and their dependence on $\omega_2$ frequency.  

-		The relative dephasing rate is $\Gamma_{10}\Delta_t=1$ and $\tau_{22^\prime}=0$.  

-		For each plot, the corresponding $\omega_2$ value is shown as a light gray vertical line.}

+	\includegraphics[width=\linewidth]{"mixed_domain/wigners full"}

+	\caption{

+    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}

 

@@ -1054,6 +1680,18 @@ Again, these features can resemble spectral diffusion even though our system is
 \subsection{Inhomogeneous broadening} \label{mix:sec:res_inhom}  % --------------------------------

 

 \begin{figure}

+	\includegraphics[width=\textwidth]{"mixed_domain/3PEPS"}

+	\label{mix:fig:3PEPS}

+	\caption{

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

@@ -1070,6 +1708,17 @@ Again, these features can resemble spectral diffusion even though our system is
 	\label{mix:fig:delay_inhom}

 \end{figure}

 

+\begin{figure}

+	\includegraphics[width=\textwidth]{"mixed_domain/2D delays"}

+	\label{mix:fig:2D_delays}

+	\caption{

+    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}

+

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

@@ -1099,7 +1748,7 @@ In our 2D delay plots (\autoref{mix:fig:delay_purity}, \autoref{mix:fig:delay_in
 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$.  % 

+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

@@ -1125,20 +1774,53 @@ time-ordering III is decoupled by detuning.  %
 	\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).  

+		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.  

+    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 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{

+    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_-4}

+	\caption{

+    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}

+

 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

@@ -1165,7 +1847,7 @@ Time-gating isolates different properties of the coherences and populations. Con
 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 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}

@@ -1190,7 +1872,7 @@ shape changes observed in Figures \ref{mix:fig:hom_2d_spectra} and
 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.  %  

+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).  %

@@ -1229,7 +1911,7 @@ $\omega_{2^\prime}-\omega_2$ which is not explored in our 2D frequency space.  %
 		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$.  

+		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}

@@ -1265,7 +1947,7 @@ We adopt the convention $\mathcal{E} = \left(a^2-b^2\right) / \left(a^2+b^2\righ
 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.  % 

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

@@ -1284,7 +1966,7 @@ It becomes independent of $\Delta_\text{inhom} / \Delta_\omega$ when $\Delta_\te
 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. 

+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}.  %