# Sequential optical response suppression for chemical mixture
characterization

Alicia B. Magann<EMAIL_ADDRESS>Department of Chemical & Biological
Engineering, Princeton University, Princeton, New Jersey 08544, USA  Center
for Computing Research, Sandia National Laboratories, Albuquerque, New Mexico
87185, USA  Gerard McCaul<EMAIL_ADDRESS>Department of Physics, Tulane
University, New Orleans, LA 70118, USA  Herschel A. Rabitz
<EMAIL_ADDRESS>Department of Chemistry, Princeton University,
Princeton, New Jersey 08544, USA  Denys I. Bondar<EMAIL_ADDRESS>Department of Physics, Tulane University, New Orleans, LA 70118, USA

###### Abstract

The characterization of mixtures of non-interacting, spectroscopically similar
quantum components has important applications in chemistry, biology, and
materials science. We introduce an approach based on quantum tracking control
that allows for determining the relative concentrations of constituents in a
quantum mixture, using a single pulse which enhances the distinguishability of
components of the mixture and has a length that scales linearly with the
number of mixture constituents. To illustrate the method, we consider two very
distinct model systems: mixtures of diatomic molecules in the gas phase, as
well as solid-state materials composed of a mixture of components. A set of
numerical analyses are presented, showing strong performance in both settings.

##  1  Introduction

Across the sciences, methods for determining the concentrations of a mixture’s
components are in high demand. The characterization of such mixtures can be a
challenging task, especially if the components have very similar chemical and
physical properties. This can occur for example in mixtures of isomers,
isotopes, or chiral molecules, as well as in mixtures of larger biological
molecules such as amino acids or proteins. For example, spectral crosstalk
currently limits the ability to distinguish between different fluorescent
proteins in experiments [20], but if this challenge could be overcome such
that the accurate characterization of fluorescent protein mixtures were
possible, this would constitute a significant advance with implications across
the biological sciences, including in synthetic biology [27], neuroscience
[10], and cytometry [2, 42, 35].

At present, standard methods (such as chromatography [45]) often struggle to
discriminate accurately between different components, as they typically rely
on the ability to resolve specific differences (e.g. in absorption peaks)
between the species, which may not be possible for mixtures composed of very
similar constituents. Optimal dynamic discrimination (ODD) was developed to
address this challenge [17, 18], by leveraging “photonic reagents” in the form
of optimally-shaped fields, and has been demonstrated in proof-of-principle
experiments in several systems, including dye molecules [3], flavin molecules
[38, 37], amino acids [36], and fluorescent proteins [14]. In essence, ODD
makes use of the fact that as long as distinctions are present between the
Hamiltonians of the different components, they will evolve differently under
the influence of a driving field applied to the mixture. These differences in
the dynamics and the resulting optical responses can be learned, and the
premise of ODD is to optimize a field to amplify these differences over time,
thereby allowing for discrimination between the components. The procedure for
accomplishing this is based on quantum optimal control, and proceeds via an
iterative optimization loop in the laboratory.

In this article, we introduce a new procedure for characterizing mixtures of
non-interacting components based on quantum tracking control principles, which
we term _scalable, single-pulse mixture characterization_ (SSMC). Quantum
tracking control is a method for designing a field $E(t)$, $t\in[0,t_{f}]$, to
drive the evolution of an observable expectation value along a prescribed
time-dependent trajectory or “track” over a time interval $t\in[0,t_{f}]$. If
reliable models for similar systems are available, the field $E(t)$ can be
computed numerically by solving a modified dynamical equation at each time
step, such that the full field can be determined from a single forward
evolution of the system dynamics. This tracking control approach stands in
contrast to numerical field-design procedures based on quantum optimal control
that identify fields for achieving a control objective in a system via an
expensive iterative optimization procedure, where at each iteration a forward
evolution of the system dynamics is required, and hundreds or thousands of
iterations can be needed for convergence of the objective function. Numerical
studies of tracking control have been performed for numerous systems including
a qubit [21], atomic systems [4], diatomic and triatomic molecules [15, 6, 7],
molecular rotors [23], and many-body solid-state systems [26, 25].

In general, quantum tracking control can suffer from singularities appearing
in $E(t)$ due to the noninvertibility of the associated dynamical equation.
Techniques to manage singularities in this context were introduced in [47,
46], and more recently, it was found that singularities could be avoided
altogether for certain tracking control problems [4, 23, 26, 25]. In these
latter cases, quantum tracking control offers a computationally attractive
method for designing fields to control observable dynamics.

In SSMC, the task of characterizing mixtures is formulated as one such
singularity-free tracking problem, where the controlled observable is taken to
be the optical response, and the prescribed track is defined such that the
optical responses of all components in the mixture are sequentially turned
off, effectively making each species sequentially invisible for a specified
time interval. SSMC is a strategy directed towards systems where either
explicit models are available, or effective models can be estimated from
experimental data. In practice, such models will always contain some degree of
error, and in principle the SSMC fields could be used as a starting basis for
ODD iterative enhancement in the laboratory. To this end, we present a general
model-based procedure for deriving equations for fields that drive the optical
response along such tracks. When these fields are applied, the time-dependent
optical response data from the full mixture can be measured and used to
evaluate the relative concentrations of each of the components by solving a
linear system of equations.

The primary benefit of the SSMC procedure is its simplicity. Namely, the use
of tracking control enables the field amplitudes to be determined by solving a
simple analytical expression at each time step. And, while the SSMC procedure
could be implemented using other forms of quantum control (i.e. optimal
control) [39, 33], this would incur significant additional computational
costs, and should only be considered when tracking control cannot be made
singularity-free. Another key feature of this procedure is its scalability
with respect to the number of mixture constituents. That is, the field used to
accomplish the characterization has a length that scales linearly in the
number of components of the mixture under consideration, and the computational
effort required a priori and post facto to determine the relative
concentrations is also efficient with respect to the number of components (see
Fig. 3, below).

To highlight the generality of SSMC, we present proof-of-principle numerical
illustrations involving two radically different models: mixtures of diatomic
molecules in the gas phase, as well as solid-state materials composed of a
mixture of components. The SSMC procedure is shown to perform consistently
well in both settings. As a comparison we consider the use of another single-
pulse strategy, namely the use of transform-limited pulses, and show that SSMC
pulses yield superior performance in all cases considered.

The remainder of this article is organized as follows. To begin, in Section 2,
the general SSMC procedure is introduced. Then, in Section 3, the application
of SSMC to the characterization of mixtures of diatomics is considered, and a
series of numerical analyses are presented. To demonstrate the flexibility of
SSMC, this is followed in Section 4 by an application of SSMC to solid-state
systems, where associated numerical illustrations are provided. We then
conclude in Section 5 with further discussion and outlook.

##  2 Theory

Consider the optical response of a mixture of $n_{s}$ non-interacting quantum
components, whose individual optical responses at some time $t$ are denoted
$R^{(s)}(t)$, for $s=1,\cdots,n_{s}$. The optical response of the mixture is
simply given by a weighted sum over all species’ optical responses:

$R_{\text{mix}}(t)=\sum_{s=1}^{n_{s}}y_{s}R^{(s)}(t),$ (1)

where $y_{s}$ denotes the relative concentration of component $s$.

Figure 1: Diagram of naive characterization procedure using a transform-
limited pulse with a simple form. The diagram shows an input pulse $E(t)$,
which serves as a “photonic reagent” for inducing optical responses $R^{(1)}$,
$R^{(2)}$, and $R^{(3)}$ in species 1, 2, and 3, respectively. The optical
response of a mixture of these species, $R_{\text{mix}}$, is given by the
weighted sum of the species optical responses, where the weights are given by
the relative concentrations, $y_{1}$, $y_{2}$, and $y_{3}$, of the species.

A naive characterization procedure might involve the application of a
transform-limited pulse to the mixture (see Fig. 1), in order to induce an
optical response from the mixture as per Eq. (1), which could be measured and
recorded. Then, if each species optical response $R^{(s)}$ were known _a
priori_ , the measured optical response data from the mixture could be
utilized to determine the relative species concentrations by inverting Eq. (1)
with respect to the vector $(y_{1},\cdots,y_{n_{s}})^{T}$, i.e., by solving
the set of coupled equations

$\mathbf{A}\mathbf{\bar{y}}=\mathbf{R_{mix}},$ (2)

where
$\mathbf{A}=(\mathbf{R}^{(1)},\,\mathbf{R}^{(2)},\cdots,\mathbf{R}^{(n_{s})})$
is a rectangular matrix containing $\mathbf{R}^{(s)}$, $s=1,\cdots,n_{s}$,
which denote the optical responses of each of the individual species;
$\mathbf{\bar{y}}=(\bar{y}_{1},\,\bar{y}_{2},\cdots,\bar{y}_{n_{s}})^{T}$ is a
vector containing the estimated values of the relative species concentrations,
where the bars are used to distinguish between the estimated values
$\bar{y}_{s}$ and the true values $y_{s}$; and $\mathbf{R_{mix}}$ is a vector
containing the mixture optical response values. The boldface vector notation
is utilized in reference to solving the system of equations associated with
the mixture characterization problem (2), in order to denote data with time-
dependent entries. When describing the underlying physics in the material
below, we do not use boldface notation.

Eq. (2) is overdetermined, but may still be solved by (for example) using a
least squares approach. The effectiveness of inverting $\mathbf{A}$ can be
quantified using the condition number,

$\textrm{cond}(\mathbf{A})=\frac{\sigma_{max}(\mathbf{A})}{\sigma_{min}(\mathbf{A})}.$
(3)

Here, $\sigma_{max}(\mathbf{A})$ and $\sigma_{min}(\mathbf{A})$ denote the
largest and smallest singular values of $\mathbf{A}$. The larger the condition
number, then the more ill-conditioned the matrix $\mathbf{A}$ is, meaning
numerical solutions of Eq. (2) will be less accurate for estimating the
relative concentrations. In order to quantify the error $\epsilon$ associated
with the estimate of the relative concentrations, the 2-norm can be used as
per

$\epsilon=\left\lVert\mathbf{y}-\bar{\mathbf{y}}\right\rVert_{2}.$ (4)

If the constituents of a mixture possess similar chemical, physical, and
optical properties, one can expect that their optical responses collected in
$\mathbf{A}$ will also be similar, leading to large condition numbers
$\textrm{cond}(\mathbf{A})$ that result in significant errors $\epsilon$ in
the estimated relative concentrations found by solving Eq. (2).

The premise of SSMC is to design individual optical responses such that the
distinguishability between the optical responses of different species
collected in $\mathbf{A}$ is enhanced, thereby yielding a better-conditioned
system of equations in Eq. (2) that can be solved to characterize the mixture
with better accuracy. In the remainder of this section, we present the SSMC
procedure as a method for achieving more accurate and robust characterization.

The SSMC procedure begins by applying a pump pulse to the mixture, denoted by
$E_{p}(t)$, in order to induce an initial response in the mixture. We note
that as the only role of the pump pulse is to induce a response, there is
significant freedom to choose it such that other desired properties can be
imposed upon the mixture (e.g., alignment of the mixture components). In order
to promote _distinguishability_ between optical responses of the different
species, the optical responses of each of the $n_{s}$ components are then
turned off in series for some time period $T$ through the sequential
application of $n_{s}$ pulses $E_{1}(t),\cdots,E_{n_{s}}(t)$. The pulse
$E_{s}(t)$ is designed to suppress the optical response of a particular
species $s$ such that $R^{(s)}(t)=0$ for the duration of the pulse $E_{s}(t)$.

A first step in determining the $E_{1}(t),\cdots,E_{n_{s}}(t)$ that will
suppress $R^{(1)}(t),\cdots R^{(n_{s})}(t)$ is the derivation of an equation
describing the form of each $R^{(s)}(t)$ in terms of $E_{s}(t)$ and one or
multiple observable expectation values at time $t$. The derivation of such
equations relies on models for the mixture components; for example molecular
systems and solid-state systems, see Sections 3 and 4, respectively, while a
more general derivation is given in Appendix A.

Once it has been derived, the equation for
$R^{(s)}\big{(}E_{s}(t),\cdots\big{)}=0$ can simply be inverted to produce an
expression for $E_{s}(t)$. This inversion procedure is the premise of quantum
tracking control, which underlies SSMC. It is important to note that in
general, the response will be a highly nonlinear function of the driving
field. As such it is difficult to obtain an interpretation to the action of
the field $E_{s}(t)$ (however, formal mechanism analyses for SSMC applications
could be done in principle, using techniques from e.g. [29, 30, 31]).

The full SSMC field $E(t)$ is then formed by concatenating the pump pulse
$E_{p}(t)$, with each of these $n_{s}$ pulses, yielding one pulse

$E(t)=\begin{cases}E_{p}(t)&\text{if $t\in[0,T)$}\\\ E_{1}(t)&\text{if
$t\in[T,2T)$}\\\ E_{2}(t)&\text{if $t\in[2T,3T)$}\\\ \cdots&\\\
E_{n_{s}}(t)&\text{if $t\in[n_{s}T,(n_{s}+1)T)$}\end{cases}$

that achieves the goal of suppressing the optical responses of each mixture
constituent in series, thereby enhancing the distinguishability of the
elements in $\mathbf{A}$.

The design of the SSMC field $E(t)$ need only be done once for a given
combination of species, and importantly, the resultant SSMC pulse $E(t)$ is
completely independent of the values of $y_{s}$, $s=1,\cdots n_{s}$ of the
species in the mixture. Hence, regardless of the relative species
concentrations in a given mixture, the same pulse can always be used to
determine these concentrations.

After they are formed in this manner, the SSMC field $E(t)$, along with the
associated optical responses of the components $R^{(s)}(t)$,
$s=1,\cdots,n_{s}$, $t\in(0,(n_{s}+1)T]$, can be stored in a library, where a
sample entry involving three components might look like Fig. 2.111We remark
that the sample library pulse shown was formed by a simple concatenation, and
that this can lead to sharp variations in the field at the concatenation
points. A variety of approaches could be employed to improve the smoothness of
the fields at these points, such as adding a short intervening period between
segments and using this period to enforce a smooth connection between them.

Figure 2: Depiction of SSMC procedure. Following the application of the pump
pulse $E_{p}(t)$, the pulses $E_{1}(t),\cdots,E_{n_{s}}(t)$ are applied to
sequentially switch off the optical response of each species for a time period
$T$. The resultant library pulse $E(t)$, formed by a concatenation of these
pulses, can then be used on an unknown mixture to determine the relative
concentrations of the constituents.

In the laboratory, $E(t)$ serves as a photonic reagent that can be drawn from
the library, in analogy with how a chemical reagent might be drawn from a
stockroom, and applied to a mixture to determine the relative concentrations
of the different components using a single pulse. Over the course of the
pulse, the mixture’s time-dependent optical response data can be recorded via
measurements and stored as a set of discrete entries in the vector
$\mathbf{R_{mix}}$. For each of the $n_{s}$ pulses $E_{1}(t),\cdots
E_{n_{s}}(t)$, we denote the number of entries in the associated optical
response vectors by $n_{t}$, such that the total number of entries for each
species over the course of the SSMC procedure is $(1+n_{s})n_{t}$. In order to
determine the concentrations, the data associated with the pump pulse is
discarded, and only $n_{s}n_{t}$ entries are utilized. Ultimately, the
estimated species concentrations $\bar{y}_{s}$, $s=1,\cdots,n_{s}$ can be
determined by solving the set of coupled equations in Eq. (2), where in SSMC,

$\mathbf{A}=\begin{pmatrix}0&\mathbf{R}_{1}^{(2)}&\cdots&\mathbf{R}_{1}^{(n_{s})}\\\
\mathbf{R}_{2}^{(1)}&0&\cdots&\mathbf{R}_{2}^{(n_{s})}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{R}_{n_{s}}^{(1)}&\mathbf{R}_{n_{s}}^{(2)}&\cdots&0\end{pmatrix}$

is a rectangular $(n_{s}n_{t}\times n_{s})$ dimensional matrix, where
$\mathbf{R}_{1}^{(s)}$ denotes the $(n_{t}\times 1)$ vector whose elements are
the optical response values associated with species $s$ collected for
$t\in(T,2T]$ while $E_{1}(t)$ is applied, $\mathbf{R}_{2}^{(s)}$ is the
$(n_{t}\times 1)$ vector of optical response values for species $s$ collected
for $t\in(2T,3T]$, etc, and

$\mathbf{R_{mix}}=(\mathbf{R_{mix,1}},\,\mathbf{R_{mix,2}},\cdots,\mathbf{R_{mix,n_{s}}})^{T}$
(5)

is an $(n_{s}n_{t}\times 1)$ dimensional vector, where $\mathbf{R_{mix,1}}$
denotes the $(n_{t}\times 1)$ vector containing the mixture optical response
values collected for $t\in(T,2T]$ while $E_{1}(t)$ is applied,
$\mathbf{R_{mix,2}}$ is the $(n_{t}\times 1)$ vector containing mixture
optical response values collected for $t\in(2T,3T]$, etc.

Because the SSMC procedure guarantees that the optical response of each
species is zero for some period $T$, the matrix $\mathbf{A}$ will be
structured to contain zeros, which can be expected to lead to a lower
condition number and smaller error. This is in contrast to a more naive
approach, where a simple transform-limited pulse is applied to induce an
optical response, and the recorded values used to solve a similar set of
equations.

An important feature of SSMC is its scalability with respect to the number of
mixture constituents. That is, if one wishes to calculate the SSMC field for
$n_{s}$ species, one must calculate a library of $n_{s}$ different responses
of length $(n_{s}+1)T$, which yields an SSMC field $E(t)$ with a length that
scales as $O(n_{s})$. Meanwhile, the computational cost for constructing the
SSMC pulse and the associated response library is $O(n_{s}^{2})$. If one
wishes to extend this library by adding new species however, it is not
necessary to completely recalculate the SSMC pulse and species responses.
Instead, one simply appends the new species to the library, evolving it with
the already calculated pulse (at cost $n_{s}T$) before applying tracking
control to eliminate its response for a period. This length $T$ extension to
the SSMC pulse is then applied to each of the other species in the library
(provided the final state of each species has been saved) to generate the
necessary addition to the response of each species in the library. Adding new
species to the discriminating pulse and response library in this manner thus
carries a cost of only $O(n_{s})$. Both of these scalings are illustrated
schematically in Fig. 3. This scalability is of particular value in mixture
characterisation, as while the condition number (and therefore error) can
scale exponentially with $n_{s}$ (see Fig. 5), the SSMC method for reducing
these errors scales only polynomially. Additionally, if one were simply to use
a single transform-limited pulse for discrimination rather than SSMC, the
scaling would depend on whether the length of the pulse was fixed, or
lengthened by $T$. In the case of a fixed length, adding another species would
also be of cost $O(n_{s})$. This would however _increase_ the degree of
overdetermination in the system, and therefore further reduce the accuracy of
calculated concentrations. If one lengthens the overall time to counteract
this however, the transform-limited condition necessitates a modification of
the pulse. In this case, the response of each species would need to be
recalculated, leading to a cost of $O(n_{s}^{2})$.

Figure 3: Schematic illustration of the scaling of SSMC (where
$n_{+}=n_{s}+1$). Each element in this block matrix is a vector of length
$n_{t}$ representing a time period $T$. To calculate the discriminating pulse
and library of responses for $n_{s}$ species will have a cost of
$O(n^{2}_{s})$. Adding a new species to an existing library however only adds
a computational cost of  $O(n_{s})$.

For the purposes of numerical and conceptual simplicity, in what follows we
will consider only one dimensional systems. Some essential information (such
as orientation) present in a three dimensional system is inevitably lost in
these lower dimensional models, but one would expect that random orientations
in a mixture would average out these contributions to the optical response of
individual molecules. Finally, we note that in principle the SSMC procedure
can be generalised to interacting systems, but at the cost of tracing out
interacting degrees of freedom when tracking each species. Given the
considerable additional analytic and computational complexity necessary to
account for interactions, we consider only non-interacting mixtures in this
work.

##  3 Characterization of molecular mixtures

In this section, we explore an application of the SSMC procedure towards the
characterization of mixtures of $n_{s}$ non-interacting molecular species, as
could occur in e.g. atmospheric or combustion samples, where the set of
possible molecular constituents in the sample is known, whereas the relative
amounts of each are not. In this setting, we consider each species $s$ itself
as being composed of $N^{(s)}$ interacting nuclear (i.e. rotational or
vibrational), polar degrees of freedom modeled in the Born-Oppenheimer
approximation, and here assumed to be aligned along the polarization direction
of the applied field $E(t)$222In practice, the alignment of the mixture
components along the same axis as the applied field could be accomplished in a
preliminary step (e.g., using the pump pulse), using techniques such those
described in [43].. The dipole moment of a particular degree of freedom $j$
belonging to a molecular species $s$ is denoted $\mu_{j}^{(s)}$, and is a
function of the vibrational or rotational coordinate $r_{j}^{(s)}$ associated
with this degree of freedom. A typical Hamiltonian governing a single species
subsystem is given in the dipole approximation by

$\displaystyle
H^{(s)}(t)=\sum_{j=1}^{N^{(s)}}-\alpha_{j}^{(s)}\frac{\partial^{2}}{\partial
r_{j}^{{(s)}^{2}}}+V^{(s)}-\sum_{j=1}^{N^{(s)}}\mu_{j}^{(s)}E(t),$ (6)

for some scalars and $\alpha_{1}^{(s)},\cdots,\alpha_{N^{(s)}}^{(s)}$ and
potential $V^{(s)}=V(r_{1}^{(s)},\cdots,r_{N^{(s)}}^{(s)})$. Because we assume
that both $V^{(s)}$ and $\mu_{j}^{(s)}$ are functions of the coordinate(s), we
assume that $[V^{(s)},\mu_{j}^{(s)}]=0$ for all $j$.

The optical response of species $s$ is proportional to the acceleration of its
total dipole moment, given here as the sum over the dipole moments
corresponding to the different degrees of freedom, i.e.,

$R^{(s)}=\sum_{k=1}^{N^{(s)}}R_{k}^{(s)}=\sum_{k=1}^{N^{(s)}}\frac{d^{2}\langle\mu_{k}^{(s)}\rangle_{t}}{dt^{2}},$
(7)

where expressions for each $R_{k}^{(s)}$ may be obtained via the Heisenberg
equation, yielding the result,

$\displaystyle R_{k}^{(s)}$
$\displaystyle=\frac{d^{2}\langle\mu_{k}^{(s)}\rangle_{t}}{dt^{2}}=-\frac{\alpha_{k}^{(s)}}{\hbar^{2}}\big{\langle}\mathcal{B}^{(s)}\big{\rangle}_{t}+\frac{2\alpha_{k}^{(s)}}{\hbar^{2}}\bigg{\langle}\frac{\partial\mu_{k}^{(s)}}{\partial
r_{k}^{(s)}}^{2}\bigg{\rangle}_{t}E(t),$ (8)

where

$\displaystyle\mathcal{B}^{(s)}$
$\displaystyle\equiv\alpha_{k}^{(s)}\bigg{(}\frac{\partial^{4}\mu_{k}^{(s)}}{\partial
r_{k}^{(s)^{4}}}+4\frac{\partial^{3}\mu_{k}^{(s)}}{\partial
r_{k}^{(s)^{3}}}\frac{\partial}{\partial
r_{k}^{(s)}}+4\frac{\partial^{2}\mu_{k}^{(s)}}{\partial
r_{k}^{(s)^{2}}}\frac{\partial^{2}}{\partial
r_{k}^{(s)^{2}}}\bigg{)}+2\frac{\partial\mu_{k}^{(s)}}{\partial
r_{k}^{(s)}}\frac{\partial V^{(s)}}{\partial r_{k}^{(s)}},$ (9)

using $\langle\cdot\rangle_{t}\equiv\langle\psi(t)|\cdot|\psi(t)\rangle$. The
optical response of the full species $s$ containing $N$ degrees of freedom is
then given by the sum of the optical responses over all degrees of freedom, as
per

$R^{(s)}=\sum_{k=1}^{N^{(s)}}R_{k}^{(s)}=\sum_{k=1}^{N^{(s)}}\frac{\alpha_{k}^{(s)}}{\hbar^{2}}\bigg{\\{}-\big{\langle}\mathcal{B}^{(s)}\big{\rangle}_{t}+2\Bigg{\langle}\frac{\partial\mu_{k}^{(s)}}{\partial
r_{k}^{(s)}}^{2}\Bigg{\rangle}_{t}E(t)\bigg{\\}}.$ (10)

In order to obtain an expression for the laser pulse $E_{s}(t)$ that
suppresses the optical response of species $s$, we may invert the equation
$R^{(s)}(t)=0$, with $R^{(s)}$ given in Eq. (10), to yield the pulse333The
denominator of Eq. (11) cannot change sign as it contains a sum of $N^{(s)}$
non-negative expectation values. And, the likelihood of encountering a
singularity, i.e., where all of these expectation values are simultaneously
equal to zero, is very small and never occurred in our numerical analyses. The
fact that this possibility exists at all is an artifact of working with a
model in the Born-Oppenheimer approximation, and the more general and
singularity-free case is treated in the appendix.

$E_{s}(t)=\frac{\sum_{k=1}^{N}\frac{\alpha_{k}^{(s)}}{\hbar^{2}}\big{\langle}\mathcal{B}^{(s)}\big{\rangle}_{t}}{2\Big{\\{}\sum_{k=1}^{N}\frac{\alpha_{k}^{(s)}}{\hbar^{2}}\Big{\langle}\frac{\partial\mu_{k}^{(s)}}{\partial
r_{k}^{(s)}}^{2}\Big{\rangle}_{t}\Big{\\}}}.$ (11)

Then the resulting dynamical equation governing the mixture,

$\frac{d}{dt}|\psi(t)\rangle=-i(H_{0}-\mu E_{s}(\psi,t))|\psi(t)\rangle,$ (12)

is highly nonlinear due to the functional dependence of the field $E_{s}(t)$
on the state $|\psi(t)\rangle$. As such, in practice, the coupled Eqs. (11)
and (12) can be solved as follows. After evolving $|\psi(t)\rangle$ to time
$t=T$ under the influence of a pump pulse $E_{p}(t)$, begin with an initial
value for $E_{1}(t=T)$ found by evaluating Eq. (11) with the initial condition
$|\psi(T)\rangle$. The state of the mixture is then propagated forward in time
via Eq. (12) by evolving the state of each species independently under its
respective Hamiltonian forward a single time step as
$|\psi_{s}(T)\rangle\rightarrow|\psi_{s}(T+\Delta t)\rangle$ for some species
$s$. The updated states are then substituted back into Eq. (11) to determine
$E_{1}(t=T+\Delta t)$, and so forth, up until time $t=2T$. This procedure is
then repeated to determine each of the subsequent fields
$E_{2}(t),\cdots,E_{n_{s}}(t)$.

###  3.1 Numerical illustrations

Figure 4: Anharmonic Morse potentials $V(r)$ (red curve) and dipole potential
$\mu(r)$ (blue curve) as a function of the bond coordinate $r$.

As an illustration, we consider characterizing mixtures of gas-phase diatomic
molecules which differ only in their masses, as occurs when discriminating
between molecules containing different isotopes. This example is meant to
serve as a simple, proof-of-concept application study of SSMC towards
molecular systems. We model the vibrational dynamics of each of the diatomics
as nonrotating Morse oscillators, initially in their ground vibrational state
and aligned with the polarization direction of the applied field [32], with
the Hamiltonian for each species given by Eq. (6) with $N^{(s)}=1$ and
$\alpha^{(s)}=\frac{\hbar^{2}}{2m^{(s)}}$ where the reduced mass of $m^{(s)}$
is taken to be a variable parameter 444The vibrational dynamics are simulated
in coordinate space using a 1D grid with 100 points, evenly spaced between
$r=0.25$ and $r=12.25$ a.u., with differential operators represented using
finite difference formulas and with $\Delta t=2.5$ a.u. Changes to the grid
did not result in qualitative changes to the results presented in Figs. 5 and
6, or to the discussions and conclusions in Sec. 3. Codes available upon
reasonable request to the authors. . In particular, we consider mixtures of
components whose masses $m^{(1)},\cdots,m^{(n_{s})}$ take on values evenly
spaced between $(1-\Delta m)m_{ref}$ and $m_{ref}$, where we choose
$m_{ref}=1800$ a.u.555In reality, mass is not a variable parameter, and the
diatomics we consider are fictional. Nevertheless, the model parameter values
were chosen to be consistent with realistic values found in diatomics, e.g.
[19] Meanwhile, $V(r)$ is given by the anharmonic Morse potential

$V(r)=D(1-e^{-\alpha(r-r_{e})})^{2}-D.$ (13)

We consider a scenario where $r_{e}=1.3$ Å is the equilibrium bond position,
$D=37,000\textrm{ cm}^{-1}$ is the potential well depth, and
$\alpha=\frac{w_{e}}{2r_{e}\sqrt{B_{e}D}}$ defines the width of the potential,
where $w_{e}=3,000\textrm{ cm}^{-1}$ and $B_{e}=11\textrm{ cm}^{-1}$. In
addition, we assume the polarization of the field is aligned with the
molecules’ dipole moments, which are modeled in Padé-approximate form as in
[19] by the function

$\mu(r)=M_{0}\frac{(1+x)^{3}}{1+\sum_{i=1}^{4}e_{i}x^{i}},$ (14)

where $x\equiv\frac{r-r_{e}}{r_{e}}$ and the parameters are given by
$M_{0}=0.5$ Debye, $e_{1}=2$, $e_{2}=2$, $e_{3}=2$, and $e_{4}=12$ (see Fig. 4
for plots of the Morse and dipole potentials). Using this model, the SSMC
field is designed via Eq. (11), where the initial pump pulse is presumed to
have the form

$E_{p}(t)=E_{0}\sin^{2}(\pi t/T)\cos(w_{e}t),$ (15)

using $E_{0}=10^{-5}$ a.u.

Figure 5: The dependence of $\text{cond}(\mathbf{A})$ on different quantities
it shown on semilog plots for the SSMC approach (red) and a naive approach
(blue). In all cases, the SSMC procedure leads to better conditioned problems
compared with using the naive approach, as indicated by values of
$\textrm{cond}(\mathbf{A})$ that are consistently a few orders of magnitude
lower for the SSMC procedure compared with the naive procedure. In (a), the
dependence on the duration of each pulse $T$ is plotted, taking $n_{s}=10$ and
$\Delta m=0.05$. In (b), the dependence on the number of species $n_{s}$ is
plotted for $T=2,500$ a.u. and $\Delta m=0.05$. In (c), the dependence on the
relative mass difference $\Delta m$ is plotted for $T=2,500$ a.u. and
$n_{s}=10$.

Fig 5 shows how SSMC performs in comparison with a naive approach, where a
transform-limited pulse (as defined in Eq. (15)) is applied for the full time
$t\in[0,n_{s}T]$. We choose this as a benchmark to compare the performance of
SSMC against due to the fact that a transform-limited pulse is the shortest
pulse with the highest peak intensity for a given bandwidth, and thus, these
pulses will more efficiently generate non-linear effects (due to the peak
intensity) compared to other pulses of equal total intensity. That is, all
else being equal, a greater peak intensity should improve the
distinguishability of different species’ responses, as the nonlinearities will
amplify the differences in the response of each species. Despite this, it is
shown in Figs. 5(a)-(c) that as $T$, $n_{s}$, and $\Delta m$ are varied, the
SSMC procedure leads to consistently lower values of
$\textrm{cond}(\mathbf{A})$ than the naive approach by a few orders of
magnitude. Importantly, we note that even when the intensity of the transform-
limited pulses used in the naive implementation are increased by multiple
orders of magnitude (keeping the SSMC fields unchanged), a clear separation in
the performance of the two techniques persists.

Figure 6: The relation between the condition of the problem
$\textrm{cond}(\mathbf{A})$ and the error in the solution $\epsilon$ is shown.
Filled point markers correspond to the aggregated data from Fig. 5(a)-(c),
while unfilled point markers correspond to the same data when
$\mathbf{R_{mix}}$ is infected by noise, leading to higher values of
$\epsilon$. A clear, approximately linear correlation between
$\textrm{cond}(\mathbf{A})$ and $\epsilon$ is evident in both cases.

Clear trends are present in $\textrm{cond}(\mathbf{A})$ in Fig. 5(a), (b), and
(c). For example, in Fig. 5(a), $\textrm{cond}(\mathbf{A})$ decreases with $T$
for both the SSMC and naive approach. This can be understood because the
problem should become better conditioned with the addition of more data, which
is what occurs when $T$ is increased. Meanwhile, the decrease in
$\textrm{cond}(\mathbf{A})$ with $\Delta m$ in Fig. 5(c) is also to be
expected, as increasing $\Delta m$ corresponds to making the different species
in the mixture more distinct, thus making the mixture easier to characterize.
Finally, the increase in $\textrm{cond}(\mathbf{A})$ with $n_{s}$ in Fig. 5(b)
shows that as the complexity of a mixture increases, i.e. by adding more
constituents, the task of characterizing the mixture becomes increasingly
difficult. In fact, in Fig. 5(b) for SSMC it appears that the behavior of
${\rm cond}(\mathbf{A})$ is approximately linear for $n_{s}\geq 10$, while for
the naive approach it is also approximately linear for $n_{s}\geq 13$. This
linearity corresponds an exponential scaling of ${\rm cond}(\mathbf{A})$ with
respect to $n_{s}$, and if these trends continue, one can expect the relative
advantage of SSMC to remain constant as the number of mixture consituents
increases.

In Fig. 6, we report the relation between the condition of the problem,
quantified by $\textrm{cond}(\mathbf{A})$ as per Eq. (3), and the error in the
solution, quantified by $\epsilon$ as per Eq. (4). All of the data from Fig. 5
is collected in Fig. 6, such that each point marker in Fig. 6 has a
corresponding point marker in Fig. 5. The filled point markers correspond to
the errors obtained in an ideal, noise free manner, while the unfilled point
markers are associated with a noisy implementation, where each element of the
vector $\mathbf{R_{mix}}$ is infected by additive random noise, drawn from a
Gaussian distribution with zero mean and a standard deviation of
$0.001\left\lVert\mathbf{R_{mix}}\right\rVert_{\infty}$. Note that the size of
the noise variance is chosen rather arbitrarily, as the purpose of this
addition is to assess the relative performance of SSMC and the naive approach,
establish the condition number as a reliable proxy for error, and introduce a
source of noise other than machine precision in the calculation of errors.

The rationale for adding noise to $\mathbf{R_{mix}}$ is as follows. First,
within any real setting there will be errors from a number of sources. These
may be measurement errors, as well as additional unknown species which will
induce their own errors. The effect of these will be to add noises to the
linear relationship in Eq. (2):

$\left(\mathbf{A}+\mathbf{\Sigma}\right)\mathbf{\bar{y}}=\mathbf{R_{mix}}+\mathbf{\Gamma},$
(16)

where $\mathbf{\Sigma}$ represents an error in the response library matrix,
and $\mathbf{\Gamma}$ an error in the observed mixture response. Note that we
do not specify the distribution for these stochastic terms. Rearranging and
expressing this equation with explicit indices, we have

$\sum_{j}\mathbf{A}_{i,j}\mathbf{\bar{y}}_{j}=\mathbf{R_{mix}}_{,i}+\mathbf{\Gamma}_{i}+\sum_{j}\mathbf{\Sigma}_{i,j}\mathbf{\bar{y}}_{j}.$
(17)

The terms
$\mathbf{\Gamma}_{i}+\sum_{j}\mathbf{\Sigma}_{i,j}\mathbf{\bar{y}}_{j}$
represent a sum of random variables, which by the central limit theorem will
be well approximated as a single normally distributed random variable. Hence,
the set of errors present both in the individual species’ library response and
the measured total response can be collected into a single Gaussian error on
$\mathbf{R_{mix}}$. By representing this error as an additive noise to the
mixture response, we are able to capture its effects while being agnostic as
to its source. Furthermore, by assigning these errors to $\mathbf{R_{mix}}$,
we are able to retain ${\rm cond}(\mathbf{A})$ as a deterministic figure of
merit.

For both noiseless and noisy implementations depicted in Fig. 6, the true
solution $\mathbf{y}$ is taken to contain relative species concentrations
drawn at random from a uniform distribution (and subsequently normalized). We
note that there is a clear, approximately linear correlation between the
condition number $\textrm{cond}(\mathbf{A})$ and the resulting error in the
solution $\epsilon$, both when noise is added and when it is not. This
suggests that the trends displayed in Fig. 5, showing values of
$\textrm{cond}(\mathbf{A})$ that are consistently a few orders of magnitude
lower for the SSMC procedure compared with the naive procedure, can be
expected to correlate in practice to errors that are also significantly lower.

##  4 Characterization of solid-state systems

In order to demonstrate the flexibility of SSMC, we now apply it to a rather
different set of model systems. Instead of molecular mixtures, we now consider
a solid-state system, composed of non-interacting layers of different
materials. In the same way as one determines relative concentrations of
molecular species in a mixture, the relative density of each species in this
layered material can also be determined using SSMC.

For the sake of simplicity, each layer will be described as a 1D Fermi-Hubbard
model [44]. Its dynamics display a number of interesting features, but the
most relevant to the current discussion is that, like the molecular mixtures
previously considered, this system also exhibits a highly non-linear response
to driving [13, 12, 34, 40]. Furthermore, the tracking equations for the
optical response of the Hubbard model are of an entirely different form to
those derived in Sec. 3. The choice of the Hubbard model is somewhat
arbitrary, as these tracking equations for current apply to any system with an
$n$-body potential dependent only on position. This provides a good test for
the SSMC, as a way of checking that its effectiveness is not tied to a
particular type of dynamical system. A full discussion of the tracking
strategy required for many-body systems of this type (and its connection to
the tracking equations used in Sec. 3) may be found in Refs. [25, 26], but we
present a brief recapitulation here.

###  4.1 Tracking Model

We model our collection of systems using the Fermi-Hubbard model

$\displaystyle{H}^{(s)}\text{$\left(t\right)$}=$ $\displaystyle-
t_{0}\sum_{j=1\atop\sigma=\uparrow,\downarrow}^{L}\text{$\left({\rm
e}^{-i\Phi\left(t\right)}{c}_{j\sigma}^{\dagger}{c}_{j+1\sigma}+{\rm
e}^{i\Phi\left(t\right)}{c}_{j+1\sigma}^{\dagger}{c}_{j\sigma}\right)$}+U^{(s)}\sum_{j=1}^{L}{c}_{j\uparrow}^{\dagger}{c}_{j\uparrow}{c}_{j\downarrow}^{\dagger}{c}_{j\downarrow},$
(18)

where the onsite interaction potential $U^{(s)}$ is used to parametrize
different systems. The spin $\sigma$ fermionic annhilation operator for the
$j$th site is given by $c_{j\sigma}$ and $t_{0}$ is the hopping parameter.
Finally, the driving field $E(t)$ enters the model through a Peierls
substitution [11], $\Phi(t)=a\int^{t}_{0}{\rm d}t^{\prime}E(t^{\prime})$,
where $a$ is the lattice constant.

For these systems the optical response corresponds to the electronic current
defined via a continuity equation for
${\rho_{j}}=\sum_{\sigma=\uparrow,\downarrow}{c}_{j\sigma}^{\dagger}{c}_{j\sigma}$:

$\displaystyle\frac{{\rm d}{\rho_{j}}}{{\rm d}t}$
$\displaystyle=\frac{1}{a}({R}_{j}-{R}_{j-1}),$ (19) $\displaystyle{R}_{j}$
$\displaystyle=-iat_{0}\sum_{\sigma=\uparrow,\downarrow}\left({\rm
e}^{-i\Phi\left(t\right)}{c}_{j\sigma}^{\dagger}{c}_{j+1\sigma}-{\rm
h.c.}\right).$ (20)

Summing over sites, we obtain the full response $R=\sum_{j}R_{j}$,

${R}=-iat_{0}\sum_{j=1\atop\sigma=\uparrow,\downarrow}^{L}\left({\rm
e}^{-i\Phi\left(t\right)}{c}_{j\sigma}^{\dagger}{c}_{j+1\sigma}-{\rm
h.c.}\right),$ (21)

with the response expectation given by

$R^{(s)}(t)=\left\langle\psi^{(s)}(t)\left|{R}\right|\psi^{(s)}(t)\right\rangle.$
(22)

Note that this expression for the response (and therefore the tracking
equation derived from it) is identical for any system with a purely spatial
potential and a dipolar coupling to the external field. If one defines the
nearest neighbour expectation in a polar form

$\left\langle\psi^{(s)}(t)\left|\sum_{j=1\atop\sigma=\uparrow,\downarrow}^{L}{c}_{j\sigma}^{\dagger}{c}_{j+1\sigma}\right|\psi^{(s)}(t)\right\rangle=K\left(\psi^{(s)}\right){\rm
e}^{i\theta\left(\psi^{(s)}\right)},$ (23)

then an explicit relationship between the response trajectory and the control
field can be established:

$R^{(s)}(t)=-2at_{0}K(\psi^{(s)})\sin(\Phi(t)-\theta(\psi^{(s)})).$ (24)

Inverting this equation, one can find an expression for a “tracking field”
$\Phi_{T}(t,\psi^{(s)})$, which produces the desired tracking response
trajectory $R^{(s)}_{T}$ in system $\psi^{(s)}$:

$\displaystyle\Phi_{T}\left(t,\psi^{(s)}\right)$
$\displaystyle=\arcsin\left[-X(t,\psi^{(s)})\right]+\theta\left(\psi^{(s)}\right),$
(25) $\displaystyle X(t,\psi^{(s)})$
$\displaystyle=\frac{R^{(s)}_{T}\left(t\right)}{2at_{0}K\left(\psi^{(s)}\right)}.$
(26)

Substituting this tracking field into Eq. (18), we re-express the system
dynamics with a “tracking Hamiltonian”
${H}^{(s)}_{T}(R^{(s)}_{T}(t),\psi^{(s)})$ parametrised by the tracked
response trajectory $R^{(s)}_{T}(t)$:

$\displaystyle{H}^{(s)}_{T}\left(R^{(s)}_{T}(t),\psi^{(s)}\right)$
$\displaystyle=\sum_{j=1\atop\sigma=\uparrow,\downarrow}^{L}\left[P^{(s)}{\rm
e}^{-i\theta\left(\psi^{(s)}\right)}{c}_{j\sigma}^{\dagger}{c}_{j+1\sigma}+\mbox{h.c.}\right]+U^{(s)}\sum_{j=1}^{L}{c}_{j\uparrow}^{\dagger}{c}_{j\uparrow}{c}_{j\downarrow}^{\dagger}{c}_{j\downarrow},$
$\displaystyle P^{(s)}=-t_{0}$
$\displaystyle\left(\sqrt{1-X^{2}(t,\psi^{(s)})}+iX(t,\psi^{(s)})\right).$
(27)

Evolving with this tracking Hamiltonian is equivalent to generating dynamics
with Eq. (18) when the driving field is chosen such that
$\langle{R}\rangle=R^{(s)}_{T}(t)$. For SSMC, during the tracking phase for
each species, we set $R^{(s)}_{T}(t)=0$.

Finally, it is possible to show that the solutions for the tracking equations
are unique (and therefore singularity free) provided certain constraints are
satisfied. The interested reader may refer to Ref. [25] for further details.

###  4.2 Numerical illustrations

Equipped with the tracking equations for the Hubbard model, we are able to
implement SSMC and compare its performance to a transform-limited pulse. In
this case we take a pump pulse given by

$\Phi_{p}(t)=a\frac{E_{0}}{\omega_{0}}\sin^{2}\left(\frac{\pi
t}{T}\right)\sin(\omega_{0}t),$ (28)

where $E_{0}=10$MeV/cm is the electric field amplitude and $T$ is set as two
periods of angular frequency $\omega_{0}=32.9$THz. The naive discriminating
pulse will be a transform-limited pulse with the same functional form as the
pump pulse, but running for the full time $t\in[0,n_{s}T]$. For both the naive
approach and SSMC, simulations are run for $L=10$ sites at half-filling (where
the number of spin up and down electrons are equal, and the overall electron
density is one per site), with each species initially in its ground state. We
note that in general an arbitrary state (such as a Gibbs state) may be used,
but given the energy gap in the model considered is much greater than thermal
fluctuations at room temperature, the ground state is an convenient choice for
the initial preparation.

Figure 7: Performance of optical discrimination vs naive approach when
$\Delta_{U}$ is varied in a two species mixture. Here $U^{(0)}=t_{0}$ and
$U^{(1)}=U^{(0)}+\Delta_{U}$.

We first examine two systems in a mixture where $t_{0}$ is constant but the
two species have onsite repulsions $U^{(0)}$ and $U^{(1)}$. We consider the
relative accuracy of both the naive method and optical discrimination as we
vary the distance in the two species’ onsite repulsions, i.e.
$U^{(1)}=U^{(0)}+\Delta_{U}$. Fig. 7 shows that over a wide range of
$\Delta_{U}$, SSMC outperforms the naive approach. This increase in efficiency
is of particular value when $\Delta_{U}$ is small, and the systems are not
easily distinguishable from each other. While the precise dynamics of each
system will depend on both the choice of $U^{(0)}$ and $\Delta_{U}$, SSMC
retains its advantage over the naive approach regardles of the choice of
$U^{(0)}$.

The accuracy gain from using SSMC can be seen more clearly when generalising
to the $n_{s}$ species case, shown in Fig. 8. Here we find qualitatively
identitical results to those shown in Fig. 5, where there is a clear advantage
to using SSMC, and that the relative performance between the two methods
remains roughly constant as one increases $n_{s}$. Furthermore, this
consistency with the results of Sec. 3 extends to the relationship between the
condition number and error $\epsilon$.

Figure 8: Performance of optical discrimination vs naive approach when the
number of species $n_{s}$ is varied. Here, $U^{(0)}=t_{0}$ and
$U^{(n_{s})}=(1+10^{-2})t_{0}$, with each of the $n_{s}$ species being equally
spaced in $U$.

##  5 Conclusion

In this article, we have introduced SSMC as a new approach for characterizing
mixtures of non-interacting components, which uses a tailored field designed
to sequentially suppress the optical responses of the mixture components. The
duration of this SSMC field scales linearly in the number of components. We
have derived expressions for the SSMC fields needed to carry out this
procedure using quantum tracking control principles, and have numerically
studied applications of SSMC to mixtures of diatomic molecules in the gas
phase, as well as to solid-state systems. In each case, we compared the
performance of SSMC with a naive approach involving the application of a
simple, transform-limited pulse to the mixture. As illustrated in Figs. 5, 7,
and 8, we found that the SSMC procedure leads to characterization problems
that are significantly better conditioned (i.e., as quantified by condition
numbers that are consistently orders of magnitude lower, compared with those
obtained utilizing a transform limited pulse). We further found that the
condition number can reliably serve as a proxy for the error in the solution,
as demonstrated by the clear correlation between the two plotted in Fig. 6.

The numerical analyses presented in Sections 3 and 4 studied the performance
of SSMC with respect to different variables, such as the number of species in
the mixture $n_{s}$ and $\Delta m$ and $\Delta_{U}$, which represent the
similarity of species in mixtures of molecular and solid state constituents,
respectively. For each study, SSMC fields $E(t)$ were found by concatenating a
transform-limited pump pulse with a series of different pulses designed to
sequentially suppress the optical responses of the species in the mixture. For
each case presented here, the suppression order was chosen such that species
were suppressed in ascending order according to their value of $m^{(s)}$ or
$U^{(s)}$. We remark that this choice was arbitrary, and further simulations
with different orderings were performed to assess the impact of the
suppression order on the performance. For both example systems, it was found
that the suppression order did not significantly affect the values obtained
for ${\rm cond}(\mathbf{A})$ or $\epsilon$.

We also remark that although the examples in Section 3 considered molecular
systems modeled in the Born-Oppenheimer approximation, the SSMC tracking
control procedure is far more general, and can be extended in a
straightforward manner to systems outside of the Born-Oppenheimer
approximation. In addition, both of the examples considered in Sections 3 and
4 involve the application of the SSMC procedure to the characterization of
mixtures of closed quantum systems undergoing unitary dynamics. Again, the
SSMC procedure is far more general, and can be extended to address open
quantum systems with dissipative dynamics as well [4]. A full derivation
treating these more general cases is provided in Appendix A.

The SSMC procedure we have introduced stands as an alternative to the naive
approach depicted in Fig. 1, and also to other mixture characterization
methods, such as ODD. Regarding the latter, we remark that the differences
between ODD and SSMC make a comprehensive comparison of the performance of the
two approaches difficult. Namely, ODD is an experimental procedure for
iteratively identifying pulses for mixture characterization, while SSMC
involves numerical pulse design. Furthermore, there is flexibility in
selecting the objective function used in ODD to optimize the pulse in the
laboratory, and the performance of ODD will depend on this choice. For
example, if the objective function is selected to match the goals of SSMC
(i.e., to sequentially turn off the optical responses of different mixture
components), then we expect that ODD could achieve similar performance to
SSMC. However, when alternate objective functions are chosen for ODD (e.g., to
reflect measurable quantities in experiments), then it is difficult to
speculate about how ODD would perform compared with SSMC.

In the laboratory, the capabilities of laser sources (e.g. bandwidth,
intensity) continue to rapidly improve [9, 16], and new techniques for
accelerating data acquisition [1] and characterizing the high frequency non-
linear part of the response field [41] (necessary for any implementation of
tracking control) have recently been demonstrated. These developments suggest
that successful experimental implementations of SSMC for a variety of mixture
types may be feasible today, or in the near future. Furthermore, the SSMC
procedure can also be modified to better suit the capabilities of existing
laser technologies. For example, it has been shown that the tracking control
methods on which SSMC is based produce control fields which can be well
approximated by only a _few distinct frequencies_ [25]; additionally, the
method used to construct them is itself robust to noise [4]. One would
therefore expect to be able to construct an approximate library pulse for
discrimination, and provided the severity of the approximation is small
compared to other sources of error, this pulse would still provide significant
improvements in accuracy over an unshaped or transform-limited pulse.

It is also worth remarking that in practice, the computational challenge of
determining the SSMC fields for mixtures whose constituents are complex
quantum systems can be significant, given that the complexity of the
simulations is exponential in the number of degrees of freedom of the
constituents. To address this, numerous quantum dynamics approximation methods
such as time-dependent density functional theory [5], (multiconfigurational)
time-dependent Hartree(-Fock) [28, 39, 33], and tensor network methods such as
time-dependent density matrix renormalization group [8] have been developed
and shown to yield impressive performance for different classes of systems at
reduced computational cost, and could therefore serve as useful tools for use
in the numerical SSMC field design procedure. Given that any error in the
relevant tracking portion of the evolution will resemble a fluctuation around
zero response (and thus resembles an error of the type considered in Sec. 3),
we expect that errors stemming from the use of these approximation frameworks
will be tolerable, and still lead to SSMC fields that can be used to
successfully characterize mixtures. Furthermore, in the future it may even be
possible to leverage quantum computing devices in order to facilitate the
design of fields for controlling quantum systems [24].

Finally, we note that numerical tracking control is not the only method
through which one can attain the essential goal - i.e. the sequential
extinguishing of each species’ optical response. Deep learning networks have
previously been employed to experimentally determine the input field required
to generate a desired response, in a manner that is robust to noise. This
technique could easily be applied to the SSMC scheme for the calculation of
the library discrimination pulse [22].

A.B.M. acknowledges support from the U.S. Department of Energy, Office of
Science, Office of Advanced Scientific Computing Research, Department of
Energy Computational Science Graduate Fellowship (grant DE-FG02-97ER25308), as
well as helpful discussions with C. Arenz. H.A.R. acknowledges support from
the U.S. Army Research Office grant numbers (W911NF-16-1-0014 and
W911NF-19-1-0382) for theoretical foundations, and the U.S. Department of
Energy grant number (DE-FG02-02ER15344 ) for formulation of applications. G.M.
and D.I.B. are supported by Air Force Office of Scientific Research (AFOSR)
Young Investigator Research Program (grant FA9550-16-1-0254) and the Army
Research Office (ARO) (grant W911NF-19-1-0377). The views and conclusions
contained in this document are those of the authors and should not be
interpreted as representing the official policies, either expressed or
implied, of AFOSR, ARO, or the U.S. Government. The U.S. Government is
authorized to reproduce and distribute reprints for Government purposes
notwithstanding any copyright notation herein. This report was prepared as an
account of work sponsored by an agency of the United States Government.
Neither the United States Government nor any agency thereof, nor any of their
employees, makes any warranty, express or implied, or assumes any legal
liability or responsibility for the accuracy, completeness, or usefulness of
any information, apparatus, product, or process disclosed, or represents that
its use would not infringe privately owned rights. Reference herein to any
specific commercial product, process, or service by trade name, trademark,
manufacturer, or otherwise does not necessarily constitute or imply its
endorsement, recommendation, or favoring by the United States Government or
any agency thereof. The views and opinions of authors expressed herein do not
necessarily state or reflect those of the United States Government or any
agency thereof. Sandia National Laboratories is a multimission laboratory
managed and operated by National Technology & Engineering Solutions of Sandia,
LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S.
Department of Energy’s National Nuclear Security Administration under contract
DE-NA0003525. This paper describes objective technical results and analysis.
Any subjective views or opinions that might be expressed in the paper do not
necessarily represent the views of the U.S. Department of Energy or the United
States Government. SAND2022-0139 J.

## Author Contributions

A.B.M and G.M. contributed equally to this work.

## References

  * Adhikari et al. [2021] Sushovit Adhikari, Cristian L. Cortes, Xiewen Wen, Shobhana Panuganti, David J. Gosztola, Richard D. Schaller, Gary P. Wiederrecht, and Stephen K. Gray.  Accelerating ultrafast spectroscopy with compressive sensing.  _Phys. Rev. Applied_ , 15:024032, 2021.  doi: 10.1103/PhysRevApplied.15.024032.
  * Bagwell and Adams [1993] C Bruce Bagwell and Earl G Adams.  Fluorescence spectral overlap compensation for any number of flow cytometry parameters.  _Ann. N. Y. Acad. Sci._ , 677(1):167–184, 1993\.  doi: 10.1111/j.1749-6632.1993.tb38775.x.
  * Brixner et al. [2001] T Brixner, NH Damrauer, P Niklaus, and G Gerber.  Photoselective adaptive femtosecond quantum control in the liquid phase.  _Nature_ , 414(6859):57–60, 2001.  doi: 10.1038/35102037.
  * Campos et al. [2017] Andre G. Campos, Denys I. Bondar, Renan Cabrera, and Herschel A. Rabitz.  How to make distinct dynamical systems appear spectrally identical.  _Phys. Rev. Lett._ , 118:083201, 2017.  doi: 10.1103/PhysRevLett.118.083201.
  * Castro et al. [2012] A. Castro, J. Werschnik, and E. K. U. Gross.  Controlling the dynamics of many-electron systems from first principles: A combination of optimal control and time-dependent density-functional theory.  _Phys. Rev. Lett._ , 109:153603, 2012.  doi: 10.1103/PhysRevLett.109.153603.
  * Chen et al. [1995] Yu Chen, Peter Gross, Viswanath Ramakrishna, Herschel Rabitz, and Kenneth Mease.  Competitive tracking of molecular objectives described by quantum mechanics.  _J. Chem. Phys._ , 102(8001), 1995.  doi: 10.1063/1.468998.
  * Chen et al. [1997] Yu Chen, Peter Gross, Viswanath Ramakrishna, Herschel Rabitz, Kenneth Mease, and Harjinder Singh.  Control of Classical Regime Molecular Objectives -Applications of Tracking and Variations on the Theme*.  _Automatica_ , 33(9):1617–1633, 1997.  doi: 10.1016/S0005-1098(97)00077-0.
  * Doria et al. [2011] Patrick Doria, Tommaso Calarco, and Simone Montangero.  Optimal control technique for many-body quantum dynamics.  _Phys. Rev. Lett._ , 106:190501, 2011.  doi: 10.1103/PhysRevLett.106.190501.
  * Elu et al. [2020] Ugaitz Elu, Luke Maidment, Lenard Vamos, Francesco Tani, David Novoa, Michael H Frosz, Valeriy Badikov, Dmitrii Badikov, Valentin Petrov, Philip St J Russell, et al.  Seven-octave high-brightness and carrier-envelope-phase-stable light source.  _Nat. Photonics_ , pages 1–4, 2020.  doi: 10.1038/s41566-020-00735-1.
  * Feng et al. [2000] Guoping Feng, Rebecca H Mellor, Michael Bernstein, Cynthia Keller-Peck, Quyen T Nguyen, Mia Wallace, Jeanne M Nerbonne, Jeff W Lichtman, and Joshua R Sanes.  Imaging neuronal subsets in transgenic mice expressing multiple spectral variants of gfp.  _Neuron_ , 28(1):41–51, 2000.  doi: 10.1016/S0896-6273(00)00084-2.
  * Feynman [2011] Richard Feynman.  _The Feynman lectures on physics_.  Basic Books, a member of the Perseus Books Group, New York, 2011.  ISBN 0465023827.
  * Ghimire et al. [2011] Shambhu Ghimire, Anthony D. Dichiara, Emily Sistrunk, Pierre Agostini, Louis F. Dimauro, and David A. Reis.  Observation of high-order harmonic generation in a bulk crystal.  _Nat. Phys._ , 7(2):138–141, 2011.  doi: 10.1038/nphys1847.
  * Ghimire et al. [2012] Shambhu Ghimire, Anthony D. Dichiara, Emily Sistrunk, Georges Ndabashimiye, Urszula B. Szafruga, Anis Mohammad, Pierre Agostini, Louis F. Dimauro, and David A. Reis.  Generation and propagation of high-order harmonics in crystals.  _Phys. Rev. A_ , 85(4):043836, 2012.  doi: 10.1103/PhysRevA.85.043836.
  * Goun et al. [2016] Alexei Goun, Denys I. Bondar, Ali O. Er, Zachary Quine, and Herschel A. Rabitz.  Photonic reagents for concentration measurement of fluorescent proteins with overlapping spectra.  _Sci. Rep._ , 6:25827, 2016.  doi: 10.1038/srep25827.
  * Gross et al. [1993] Peter Gross, Harjinder Singh, Herschel Rabitz, Kenneth Mease, and GM Huang.  Inverse quantum-mechanical control: A means for design and a test of intuition.  _Phys. Rev. A_ , 47(6):4593, 1993.  doi: 10.1103/PhysRevA.47.4593.
  * Lesko et al. [2021] Daniel MB Lesko, Henry Timmers, Sida Xing, Abijith Kowligy, Alexander J Lind, and Scott A Diddams.  A six-octave optical frequency comb from a scalable few-cycle erbium fibre laser.  _Nat. Photonics_ , 15(4):281–286, 2021.  doi: 10.1038/s41566-021-00778-y.
  * Li et al. [2002] Baiqing Li, Gabriel Turinici, Viswanath Ramakrishna, and Herschel Rabitz.  Optimal dynamic discrimination of similar molecules through quantum learning control.  _J. Phys. Chem. B_ , 106(33):8125–8131, 2002\.  doi: 10.1021/jp0204657.
  * Li et al. [2005] Baiqing Li, Herschel Rabitz, and Jean-Pierre Wolf.  Optimal dynamic discrimination of similar quantum systems with time series data.  _J. Chem. Phys_ , 122(15):154103, 2005.  doi: 10.1063/1.1883170.
  * Li et al. [2013] Gang Li, Iouli E. Gordon, Robert J. Le Roy, Photos G. Hajigeorgiou, John A. Coxon, Peter F. Bernath, and Laurence S. Rothman.  Reference spectroscopic data for hydrogen halides. part i: Construction and validation of the ro-vibrational dipole moment functions.  _J. Quant. Spectrosc. Radiat. Transf._ , 121:78 – 90, 2013\.  doi: 10.1016/j.jqsrt.2013.02.005.
  * Lichtman and Conchello [2005] Jeff W Lichtman and José-Angel Conchello.  Fluorescence microscopy.  _Nat. Methods_ , 2(12):910–919, 2005.  doi: 10.1038/nmeth817.
  * Lidar and Schneider [2005] Daniel A Lidar and Sara Schneider.  Stabilizing qubit coherence via tracking-control.  _Quantum Inf. Comput._ , 5(4), 2005.  URL https://dl.acm.org/doi/10.5555/2011645.2011651.
  * Lohani et al. [2019] Sanjaya Lohani, Erin M. Knutson, Wenlei Zhang, and Ryan T. Glasser.  Dispersion characterization and pulse prediction with machine learning.  _OSA Continuum_ , 2(12):3438, 2019.  doi: 10.1364/OSAC.2.003438.
  * Magann et al. [2018] Alicia Magann, Tak-San Ho, and Herschel Rabitz.  Singularity-free quantum tracking control of molecular rotor orientation.  _Phys. Rev. A_ , 98:043429, 2018.  doi: 10.1103/PhysRevA.98.043429.
  * Magann et al. [2021] Alicia B. Magann, Matthew D. Grace, Herschel A. Rabitz, and Mohan Sarovar.  Digital quantum simulation of molecular dynamics and control.  _Phys. Rev. Research_ , 3:023165, Jun 2021.  doi: 10.1103/PhysRevResearch.3.023165.
  * McCaul et al. [2020a] Gerard McCaul, Christopher Orthodoxou, Kurt Jacobs, George H Booth, and Denys I Bondar.  Controlling arbitrary observables in correlated many-body systems.  _Phys. Rev. A_ , 101(5):053408, 2020a.  doi: 10.1103/PhysRevA.101.053408.
  * McCaul et al. [2020b] Gerard McCaul, Christopher Orthodoxou, Kurt Jacobs, George H Booth, and Denys I Bondar.  Driven imposters: Controlling expectations in many-body systems.  _Phys. Rev. Lett._ , 124(18):183201, 2020b.  doi: 10.1103/PhysRevLett.124.183201.
  * Medley et al. [2005] Colin D Medley, Timothy J Drake, Jeffrey M Tomasini, Richard J Rogers, and Weihong Tan.  Simultaneous monitoring of the expression of multiple genes inside of single breast carcinoma cells.  _Anal. Chem._ , 77(15):4713–4718, 2005.  doi: 10.1021/ac050881y.
  * Messina et al. [1996] Michael Messina, Kent R Wilson, and Jeffrey L Krause.  Quantum control of multidimensional systems: Implementation within the time?dependent hartree approximation.  _J. Chem. Phys._ , 104:173, 1996.  doi: 10.1063/1.470887.
  * Mitra and Rabitz [2003] Abhra Mitra and Herschel Rabitz.  Identifying mechanisms in the control of quantum dynamics through hamiltonian encoding.  _Phys. Rev. A_ , 67(3):033407, 2003.  doi: 10.1103/PhysRevA.67.033407.
  * Mitra and Rabitz [2004] Abhra Mitra and Herschel Rabitz.  Mechanistic analysis of optimal dynamic discrimination of similar quantum systems.  _J. Phys. Chem. A_ , 108(21):4778–4785, 2004\.  doi: 10.1021/jp0495390.
  * Mitra and Rabitz [2006] Abhra Mitra and Herschel Rabitz.  Quantum control mechanism analysis through field based hamiltonian encoding.  _J. Chem. Phys._ , 125(19):194107, 2006.  doi: 10.1063/1.2371079.
  * Morse [1929] Philip M. Morse.  Diatomic molecules according to the wave mechanics. ii. vibrational levels.  _Phys. Rev._ , 34:57–64, 1929.  doi: 10.1103/PhysRev.34.57.
  * Mundt and Tannor [2009] Michael Mundt and David J Tannor.  Optimal control of interacting particles: a multi-configuration time-dependent hartree–fock approach.  _New J. Phys_ , 11(10):105038, 2009.  doi: 10.1088/1367-2630/11/10/105038.
  * Murakami et al. [2018] Yuta Murakami, Martin Eckstein, and Philipp Werner.  High-Harmonic Generation in Mott Insulators.  _Phys. Rev. Lett._ , 121(5):057405, 2018.  doi: 10.1103/PhysRevLett.121.057405.
  * Perfetto et al. [2004] Stephen P Perfetto, Pratip K Chattopadhyay, and Mario Roederer.  Seventeen-colour flow cytometry: unravelling the immune system.  _Nat. Rev. Immunol._ , 4(8):648–655, 2004.  doi: 10.1038/nri1416.
  * Rondi et al. [2012] Ariana Rondi, Luigi Bonacina, A Trisorio, C Hauri, and J-P Wolf.  Coherent manipulation of free amino acids fluorescence.  _Phys. Chem. Chem. Phys._ , 14(26):9317–9322, 2012.  doi: 10.1039/C2CP23357F.
  * Roslund et al. [2011] Jonathan Roslund, Matthias Roth, Laurent Guyon, Véronique Boutou, Francois Courvoisier, Jean-Pierre Wolf, and Herschel Rabitz.  Resolution of strongly competitive product channels with optimal dynamic discrimination: Application to flavins.  _J. Chem. Phys._ , 134(3):034511, 2011.  doi: 10.1063/1.3518751.
  * Roth et al. [2009] Matthias Roth, Laurent Guyon, Jonathan Roslund, Véronique Boutou, Francois Courvoisier, Jean-Pierre Wolf, and Herschel Rabitz.  Quantum control of tightly competitive product channels.  _Phys. Rev. Lett._ , 102:253001, 2009.  doi: 10.1103/PhysRevLett.102.253001.
  * Schröder et al. [2008] Markus Schröder, José-Luis Carreón-Macedo, and Alex Brown.  Implementation of an iterative algorithm for optimal control of molecular dynamics into mctdh.  _Phys. Chem. Chem. Phys._ , 10:850, 2008.  doi: 10.1039/B714821F.
  * Silva et al. [2018] R. E.F. Silva, Igor V. Blinov, Alexey N. Rubtsov, O. Smirnova, and M. Ivanov.  High-harmonic spectroscopy of ultrafast many-body dynamics in strongly correlated systems.  _Nat. Photonics_ , 12(5):266–270, 2018.  doi: 10.1038/s41566-018-0129-0.
  * Sommer et al. [2016] A Sommer, E M Bothschafter, S A Sato, C Jakubeit, T Latka, O Razskazovskaya, H Fattahi, M Jobst, W Schweinberger, V Shirvanyan, V S Yakovlev, R Kienberger, K Yabana, N Karpowicz, M Schultze, and F Krausz.  Attosecond nonlinear polarization and light–matter energy transfer in solids.  534(7605):86–90, 2016.  doi: 10.1038/nature17650.
  * Speicher et al. [1996] Michael R Speicher, Stephen Gwyn Ballard, and David C Ward.  Karyotyping human chromosomes by combinatorial multi-fluor fish.  _Nat. Genet._ , 12(4):368–375, 1996.  doi: 10.1038/ng0496-368.
  * Stapelfeldt and Seideman [2003] Henrik Stapelfeldt and Tamar Seideman.  Colloquium: Aligning molecules with strong laser pulses.  _Rev. Mod. Phys._ , 75:543–557, Apr 2003.  doi: 10.1103/RevModPhys.75.543.
  * Tasaki [1998] Hal Tasaki.  The Hubbard model - An introduction and selected rigorous results.  _J. Phys. Cond. Mat._ , 10(20):4353–4378, 1998\.  doi: 10.1088/0953-8984/10/20/004.
  * Wilson [2000] Ian Wilson.  _Encyclopedia of separation science_.  Academic Press, San Diego, 2000.  ISBN 0122267702.
  * Zhu and Rabitz [2003] Wusheng Zhu and Herschel Rabitz.  Quantum control design via adaptive tracking.  _J. Chem. Phys._ , 119(7):3619–3625, 2003.  doi: 10.1063/1.1582847.
  * Zhu et al. [1999] Wusheng Zhu, Martina Smit, and Herschel Rabitz.  Managing singular behavior in the tracking control of quantum dynamical observables.  _J. Chem. Phys._ , 110(4):1905–1915, 1999.  doi: 10.1063/1.477857.

##  Appendix A First principles derivation

In this section, we consider a derivation for the SSMC equations for a quantum
mixture based on first principles. We begin by considering a component of a
quantum mixture which is composed of $N_{n}^{(s)}$ nuclei and $N_{e}^{(s)}$
electrons, where the subscripts distinguish between nuclei and electrons and
the superscripts label the mixture constituent they belong to. For systems
undergoing dissipative dynamics, we may model the dynamics of the state
$\rho^{(s)}(t)$ of the mixture component by the dynamical equation

$\displaystyle\frac{d}{dt}\rho^{(s)}(t)$
$\displaystyle=\mathcal{L}^{(s)}(\rho^{(s)}(t))$ (29)
$\displaystyle=-i[H^{(s)}(t),\rho^{(s)}(t)]+\mathcal{D}^{(s)}(\rho^{(s)}(t)),$

which is of Lindblad form, where the dissipator $\mathcal{D}^{(s)}(\cdot)$
describes the interactions between $\rho^{(s)}(t)$ and its environment, given
by

$\displaystyle\mathcal{D}^{(s)}(\rho^{(s)}(t))$
$\displaystyle=\sum_{j,k}\gamma_{jk}^{(s)}\Big{(}A_{jk}^{(s)}\rho^{(s)}(t)A_{jk}^{(s)\dagger}-\frac{1}{2}A_{jk}^{(s)\dagger}A_{jk}^{(s)}\rho^{(s)}(t)-\frac{1}{2}\rho^{(s)}(t)A_{jk}^{(s)\dagger}A_{jk}^{(s)}\Big{)},$
(30)

where $A_{jk}^{(s)}=|j\rangle\langle k|^{(s)}$ and $\gamma_{jk}^{(s)}$ are the
scalar coefficients weighting each of the projectors $A_{jk}^{(s)}$.
Meanwhile, the Hamiltonian $H^{(s)}(t)$ is given by

$H^{(s)}(t)=H_{0}^{(s)}+H_{c}^{(s)}(t),$ (31)

where $H_{0}^{(s)}$ describes the species time-independent drift Hamiltonian
and $H_{c}^{(s)}(t)$ describes the species time-dependent control Hamiltonian.

In the following, we consider a single species only and drop the superscript
species label. The drift Hamiltonian for any given molecular species is given
by,

$\displaystyle H_{0}$
$\displaystyle=T_{\text{nuc}}\big{(}\\{\vec{{R}}_{J}\\}\big{)}+T_{\text{el}}\big{(}\\{\vec{{r}}_{j}\\}\big{)}+T_{\text{int}}\big{(}\\{\vec{{R}}_{J}\\},\\{\vec{{r}}_{j}\\}\big{)}$
(32)
$\displaystyle=\bigg{\\{}\sum_{J=1}^{N_{n}}\frac{\vec{{P}}_{J}^{2}}{2m_{J}}+\sum_{J<K}\frac{Z_{J}Z_{K}}{\vec{R}_{JK}}\bigg{\\}}+\bigg{\\{}\sum_{j=1}^{N_{e}}\frac{\vec{{p}}_{j}^{2}}{2}+\sum_{j>k}\frac{1}{\vec{R}_{jk}}\bigg{\\}}+\bigg{\\{}-\sum_{J,j}\frac{Z_{J}}{\vec{R}_{Jj}}\bigg{\\}},$

where $\\{\vec{{R}}_{J}\\}$ and $\\{\vec{{r}}_{j}\\}$ denote the sets of
collective nuclear and electronic coordinates of a molecular species,
respectively. Vector notation is used here to denote
$\vec{{a}}=a_{x}\vec{{e}}_{x}+a_{y}\vec{{e}}_{y}+a_{z}\vec{{e}}_{z}$. We
denote the charge and mass of the $J$-th nucleus by $Z_{J}$ and $m_{J}$
respectively, we denote the momentum of the $j$-th electron by $\vec{{p}}_{j}$
and the momentum of the $J$-th nucleus by $\vec{{P}}_{J}$, and we use the
abbreviated notation $\vec{R}_{jk}\equiv|\vec{{r}}_{j}-\vec{{r}}_{k}|$,
$\vec{R}_{JK}\equiv|\vec{{R}}_{J}-\vec{{R}}_{K}|$, and
$\vec{R}_{Jj}\equiv|\vec{{R}}_{J}-\vec{{r}}_{j}|$ for the relative distances.
The control Hamiltonian describes the interaction of the system with an
applied laser field, and is modeled in the dipole approximation as

$H_{c}(t)=-\bigg{(}\sum_{j=1}^{N_{e}^{(s)}}\vec{{r}}_{j}^{(s)}-\sum_{J=1}^{N_{n}^{(s)}}Z_{J}^{(s)}\vec{{R}}_{J}^{(s)}\bigg{)}\cdot\vec{\varepsilon}(t),$
(33)

where $\vec{\varepsilon}(t)$ is the applied field amplitude at time $t$. The
optical response of a single species $s$ is proportional to the acceleration
of the dipole moment of the species,
$\vec{R}(t)\propto\frac{d^{2}}{dt^{2}}\big{\langle}\vec{\mu}\big{\rangle}_{t}$,
where the notation $\langle\vec{O}\rangle_{t}\equiv{\rm
Tr}\\{O_{x}\rho(t)\\}\vec{\mathbf{e}}_{x}+{\rm
Tr}\\{O_{y}\rho(t)\\}\vec{\mathbf{e}}_{y}+{\rm
Tr}\\{O_{z}\rho(t)\\}\vec{\mathbf{e}}_{z}$ denotes the expectation value of a
vector operator at time $t$, and we have introduced the notation

$\vec{\mu}\equiv\sum_{j=1}^{N_{e}}\vec{{r}}_{j}-\sum_{J=1}^{N_{n}}Z_{J}\vec{{R}}_{J}$
(34)

for the dipole moment of the species. In the same spirit, we also adopt the
following notation:

$\vec{p}\equiv\sum_{j=1}^{N_{e}}\vec{{p}}_{j}-\sum_{J=1}^{N_{n}}\frac{Z_{J}}{m_{J}}\vec{{P}}_{J}$
(35)
$\vec{\nabla}\vec{V}\equiv\sum_{j=1}^{N_{e}}\vec{\nabla}_{j}\frac{Z_{J}}{\vec{R}_{Jj}}-\sum_{J=1}^{N_{n}}\frac{Z_{J}}{m_{J}}\vec{\nabla}_{J}\frac{Z_{J}}{\vec{R}_{Jj}}$
(36) $Z_{tot}\equiv\sum_{J=1}^{N_{n}}Z_{J}$ (37)
$m_{tot}\equiv\sum_{J=1}^{N_{n}}m_{J},$ (38)

which leads to the following commutation relations:

$\displaystyle=-i\vec{p},\qquad\qquad[H(t),\vec{p}]$
$\displaystyle=i\vec{\nabla}\vec{V}-i\big{(}N_{e}+\tfrac{Z_{tot}}{m_{tot}}\big{)}\vec{\varepsilon}(t),$
(39)

where we remark that $\vec{\nabla}\vec{V}=\frac{\partial V_{x}}{\partial
x}\vec{\mathbf{e}}_{x}+\frac{\partial V_{y}}{\partial
y}\vec{\mathbf{e}}_{y}+\frac{\partial V_{z}}{\partial z}\vec{\mathbf{e}}_{z}$.
Our goal is to derive an expression for the species optical response when the
system is acted on by the field $\vec{\varepsilon}(t)$. We first recall that
from Eq. (29), the time evolution of a species between times $t$ and $t+dt$ is
given by

$\rho(t+dt)=\mathcal{T}e^{\int_{t}^{t+dt}\mathcal{L}(\cdot)dt^{\prime}}\rho(t),$
(40)

where $\mathcal{T}$ is the time-ordering operator. Then, using the Taylor
expansion it can be shown that eliminating $\mathcal{T}$ contributes an error
$O(dt^{3})$. A Taylor expansion can also be used to show
$\int_{t}^{t+dt}\mathcal{L}(\cdot)dt^{\prime}=\mathcal{L}_{D}(\cdot)+O(dt^{3})$
where
$\mathcal{L}_{D}(\cdot)\equiv-i\Big{[}H\big{(}t+\tfrac{dt}{2}\big{)},\cdot\Big{]}dt+\mathcal{D}(\cdot)dt$,
yielding

$\langle\vec{\mu}\rangle_{t+dt}={\rm
Tr}\Big{\\{}\vec{\mu}e^{\mathcal{L}_{D}(\cdot)dt}\rho(t)\Big{\\}}+O(dt^{3})\,.$
(41)

We now expand
$e^{\mathcal{L}_{D}(\cdot)dt}\rho(t)\approx\rho(t)+\mathcal{L}_{D}(\rho(t))dt+\mathcal{L}_{D}\big{(}\mathcal{L}_{D}(\rho(t))\big{)}\frac{dt^{2}}{2}+\cdots$.
We truncate the expansion at second order, such that the error remains
$O(dt^{3})$. We can then use the definition of the adjoint operator

${\rm Tr}\big{\\{}\vec{\mu}\mathcal{L}_{D}(\rho(t))\big{\\}}={\rm
Tr}\big{\\{}\mathcal{L}_{D}^{\dagger}(\vec{\mu})\rho(t)\big{\\}},$ (42)

where

$\mathcal{L}_{D}^{\dagger}(\vec{\mu})=i\Big{[}H\big{(}t+\tfrac{dt}{2}\big{)},\vec{\mu}\Big{]}+\mathcal{D}^{\dagger}(\vec{\mu})$
(43)

and in the above,

$\mathcal{D}^{\dagger}\big{(}\vec{\mu}\big{)}=\sum_{j,k}\gamma_{jk}\big{(}A_{jk}^{\dagger}\vec{\mu}A_{jk}-\frac{1}{2}A_{jk}^{\dagger}A_{jk}\vec{\mu}-\frac{1}{2}\vec{\mu}A_{jk}^{\dagger}A_{jk}\big{)},$
(44)

to write

$\displaystyle\langle\vec{\mu}\rangle_{t+dt}$ $\displaystyle={\rm
Tr}\big{\\{}\vec{\mu}\rho(t)\big{\\}}+{\rm
Tr}\bigg{\\{}\Big{(}\mathcal{L}_{D}^{\dagger}(\vec{\mu})dt+\mathcal{L}_{D}^{\dagger}\big{(}\mathcal{L}_{D}^{\dagger}(\vec{\mu})\big{)}\tfrac{dt^{2}}{2}\Big{)}\rho(t)\bigg{\\}}+O(dt^{3})$
(45)
$\displaystyle=\langle\vec{\mu}\rangle_{t}+\Big{\langle}\mathcal{L}_{D}^{\dagger}(\vec{\mu})dt+\mathcal{L}_{D}^{\dagger}\big{(}\mathcal{L}_{D}^{\dagger}(\vec{\mu})\big{)}\tfrac{dt^{2}}{2}\Big{\rangle}_{t}+O(dt^{3}).$

We see that

$\displaystyle\mathcal{L}_{D}^{\dagger}\big{(}\mathcal{L}_{D}^{\dagger}(\vec{\mu})\big{)}$
$\displaystyle=i\big{[}H\big{(}t+\tfrac{dt}{2}\big{)},\vec{p}\big{]}+\mathcal{D}^{\dagger}(\vec{p})+i\Big{[}H\big{(}t+\tfrac{dt}{2}\big{)},\mathcal{D}^{\dagger}(\vec{\mu})\Big{]}+\mathcal{D}^{\dagger}\big{(}\mathcal{D}^{\dagger}(\vec{\mu})\big{)}.$
(46)

Substituting the relation from Eq. (39) into Eq. (46),

$\displaystyle\mathcal{L}_{D}^{\dagger}\big{(}\mathcal{L}_{D}^{\dagger}(\vec{\mu})\big{)}$
$\displaystyle=-\vec{\nabla}\vec{V}+\big{(}N_{e}+\tfrac{Z_{tot}}{m_{tot}}\big{)}\vec{\varepsilon}(t+\tfrac{dt}{2})+\mathcal{D}^{\dagger}\big{(}\vec{p}+\mathcal{D}^{\dagger}(\vec{\mu})\big{)}+i\big{[}H\big{(}t+\tfrac{dt}{2}\big{)},\mathcal{D}^{\dagger}(\vec{\mu})\big{]}.$
(47)

From Eq. (43), we see that
$\frac{d\langle\vec{\mu}\rangle_{t}}{dt}=i\langle[H(t+\tfrac{dt}{2}),\vec{\mu}]\rangle_{t}+\langle\mathcal{D}^{\dagger}(\vec{\mu})\rangle_{t}=\langle\vec{p}\rangle_{t}+\langle\mathcal{D}^{\dagger}(\vec{\mu})\rangle_{t}$.
We substitute this and the definition

$\langle\mathcal{\vec{A}}\rangle_{t}\equiv\Big{\langle}\mathcal{D}^{\dagger}\Big{(}\vec{p}+\mathcal{D}^{\dagger}(\vec{\mu})\Big{)}\Big{\rangle}_{t}+i\Big{\langle}\Big{[}H\big{(}t+\tfrac{dt}{2}\big{)},\mathcal{D}^{\dagger}(\vec{\mu})\Big{]}\Big{\rangle}_{t}\,.$
(48)

into Eq. (47). Then, substituting Eq. (47) into Eq. (45) leads to

$\displaystyle\langle\vec{\mu}\rangle_{(t+dt)}$
$\displaystyle=\langle\vec{\mu}\rangle_{t}+dt\frac{d\langle\vec{\mu}\rangle_{t}}{dt}+\frac{dt^{2}}{2}\Big{[}\big{(}N_{e}+\tfrac{Z_{tot}}{m_{tot}}\big{)}\vec{\varepsilon}\big{(}t+\tfrac{dt}{2}\big{)}+\langle\mathcal{\vec{A}}\rangle_{t}-\langle\vec{\nabla}\vec{V}\rangle_{t}\Big{]}+O(dt^{3})\,,$
(49)

which can be rearranged to

$\displaystyle\frac{2}{dt^{2}}\big{(}\langle\vec{\mu}\rangle_{(t+dt)}-\langle\vec{\mu}\rangle_{t}\big{)}$
$\displaystyle=\frac{2}{dt}\frac{d\langle\vec{\mu}\rangle_{t}}{dt}+\big{(}N_{e}+\tfrac{Z_{tot}}{m_{tot}}\big{)}\vec{\varepsilon}\big{(}t+\tfrac{dt}{2}\big{)}+\langle\mathcal{\vec{A}}\rangle_{t}-\langle\vec{\nabla}\vec{V}\rangle_{t}+O(dt)\,.$
(50)

We now substitute in the Taylor expansion
$\langle\vec{\mu}\rangle_{(t+dt)}=\langle\vec{\mu}\rangle_{t}+\frac{d\langle\vec{\mu}\rangle_{t}}{dt}dt+\frac{d^{2}\langle\vec{\mu}\rangle_{t}}{dt^{2}}\frac{dt^{2}}{2}+O(dt^{3})$
and cancel terms to arrive at our final expression for the species optical
response $\vec{R}(t)=\frac{d^{2}\langle\vec{\mu}\rangle_{t}}{dt^{2}}$,

$\displaystyle\vec{R}(t)$
$\displaystyle=\big{(}N_{e}+\tfrac{Z_{tot}}{m_{tot}}\big{)}\vec{\varepsilon}\big{(}t+\tfrac{dt}{2}\big{)}+\langle\mathcal{\vec{A}}\rangle_{t}-\langle\vec{\nabla}\vec{V}\rangle_{t}+O(dt).$
(51)

We now return to our full system, which describes a mixture of $n_{s}$
species. The optical response of the mixture is defined as the sum of the
optical responses of the species weighted by their relative concentrations,
i.e.,

$\displaystyle\vec{R}_{\text{mix}}(t)$
$\displaystyle=\sum_{s=1}^{n_{s}}y_{s}\Big{[}\big{(}N_{e}^{(s)}+\tfrac{Z_{tot}}{m_{tot}}^{(s)}\big{)}\vec{\varepsilon}\big{(}t+\tfrac{dt}{2}\big{)}+\langle\mathcal{\vec{A}}\rangle_{t}^{(s)}-\langle\vec{\nabla}^{(s)}\vec{V}^{(s)}\rangle_{t}\Big{]}+O(dt).$
(52)

Molecular discrimination between species can be carried out by designing
$n_{s}$ pulse shapes
$\vec{\varepsilon}_{1}(t),\cdots,\vec{\varepsilon}_{n_{s}}(t)$ such that
$\vec{\varepsilon}_{s}(t)$ suppresses the optical response of species $s$,
i.e.,

$\displaystyle 0$
$\displaystyle=y_{s}\Big{[}\big{(}N_{e}^{(s)}+\tfrac{Z_{tot}}{m_{tot}}^{(s)}\big{)}\vec{\varepsilon}_{s}\big{(}t+\tfrac{dt}{2}\big{)}+\langle\mathcal{\vec{A}}\rangle_{t}^{(s)}-\langle\vec{\nabla}^{(s)}\vec{V}^{(s)}\rangle_{t}\Big{]},$
(53)

which rearranges to the final SSMC equation

$\vec{\varepsilon}_{s}(t+\tfrac{dt}{2})=-\frac{\Big{(}\langle\mathcal{\vec{A}}\rangle_{t}^{(s)}-\langle\vec{\nabla}^{(s)}\vec{V}^{(s)}\rangle_{t}\Big{)}}{\big{(}N_{e}^{(s)}+\tfrac{Z_{tot}}{m_{tot}}^{(s)}\big{)}},$
(54)

which is free of singularities due to the fact that the denominator is
restricted to be a positive number. For a closed quantum system undergoing
unitary dynamics, the expression would be the identical, but with
$\langle\mathcal{\vec{A}}\rangle_{t}^{(s)}=0$.666The molecular case discussed
in the main text in Sec. 3 can be obtained by assuming closed system dynamics,
and by assuming the Born-Oppenheimer approximation, such that the dynamics of
the nuclei and electrons are assumed to be uncoupled. Then, the illustrations
in Sec. 3.1 use a model for the vibrational dynamics of the nuclei alone.