††thanks: Correspoding author††thanks: Current address: _Vienna Center for
Quantum Science and Technology (VCQ), Faculty of Physics, University of
Vienna, 1010 Vienna, Austria._

# Compressive quantum waveform estimation

Alex Tritt<EMAIL_ADDRESS>School of Physics & Astronomy, Monash
University, Victoria 3800, Australia.  Joshua Morris  School of Physics &
Astronomy, Monash University, Victoria 3800, Australia.  Christopher C. Bounds
School of Physics & Astronomy, Monash University, Victoria 3800, Australia.
Hamish A. M. Taylor  School of Physics & Astronomy, Monash University,
Victoria 3800, Australia.  James Saunderson  Department of Electrical and
Computer Systems Engineering, Monash University, Victoria 3800, Australia.  L.
D. Turner  School of Physics & Astronomy, Monash University, Victoria 3800,
Australia.

###### Abstract

Quantum waveform estimation, in which quantum sensors sample entire time
series, promises to revolutionize the sensing of weak and stochastic signals,
such as the biomagnetic impulses emitted by firing neurons. For long duration
signals with rapid transients, regular quantum sampling becomes prohibitively
resource intensive as it demands many measurements with distinct control and
readout. In this Manuscript, we demonstrate how careful choice of quantum
measurements, along with the modern mathematics of compressive sensing,
achieves quantum waveform estimation of sparse signals in a number of
measurements far below the Nyquist requirement. We sense synthesized neural-
like magnetic signals with radiofrequency-dressed ultracold atoms, retrieving
successful waveform estimates with as few measurements as compressive
theoretical bounds guarantee.

00footnotetext: See Supplemental Material at _[URL will be inserted by
publisher]_ for a derivation of sensitivity to sine coefficients, an
explanation of regularization parameter tuning, an explanation of how we use
the AUC statistic from signal detection theory, and recoveries from genuinely
undersampled measurements.

Quantum waveform estimation extends the applicability of precise quantum
sensors to the realm of sampling entire signals [1], at the cost of requiring
many quantum measurements. One proposed application is the non-invasive
sensing of electrical communication in a network of neurons by recording
proximal magnetic fields [2, 3, 4, 5, 6, 7, 8, 9]. Such measurements would aid
in understanding the brain on the microscale, in health and in disease [10].
Neural waveforms, composed of bipolar pulses separated by large gaps [11],
clearly contain scant information: sampling one uniformly at its Nyquist rate
[12, Sec. 7.1] will measure mostly silence. If information is so sparse, then
could one save on resources by taking far fewer quantum measurements? Naïvely
reducing a uniform sample rate, however, would ultimately lose pulses to
aliasing. In this Manuscript, we use _compressive sensing_ to demonstrate
quantum waveform estimation of a sparse, neural-like magnetic signal using
many fewer quantum measurements than the Nyquist-Shannon sampling theorem
would deem necessary.

The modern mathematics of compressive sensing [13, 14, 15] (also _compressed
sensing_ or _compressive sampling_) describes how, under certain conditions,
sparse signals can be fully recovered from an _incomplete_ set of
measurements. Compressive sensors thrive where there is benefit to taking few
measurements, e.g., cameras for wavelengths with expensive photodetectors
[16], and computed tomography scans which irradiate patients [17]. Compressive
sensing in quantum science has assisted quantum process and state tomography
[18], quantum communication [19, 20], quantum computation [21, 22, 23],
quantum annealing [24, 25, 26, 27], ghost imaging [16, 28, 29, 30], and
quantum sensor readout [31, 32]. In quantum waveform estimation, compressive
sensing has been used to denoise signals [33, 34] or disambiguate frequencies
in complete sets of quantum measurements [35]. Magesan  _et al._ [36] proposed
literal _compressive quantum sensing_ , simulating the reconstruction of a
neural-like signal from an incomplete set of an incomplete set of Walsh basis
measurements using sparse recovery. However, the proposed sensor, based on
nitrogen vacancy (NV)centers in diamond, would decohere before even a single
neural firing would complete, and no such sensor has been constructed.

Compressive quantum sensors that significantly reduce the number of required
quantum measurements would make quantum waveform estimation much more
practical. In particular, it would bring into view quantum sensor arrays
capturing entire waveforms in a single shot. To capitalize on this promise, we
need: (1) quantum sensors that remain coherent for the duration of the signal,
(2) a quantum sensing protocol to measure in a basis that collects information
uniformly, and (3) a post-processing algorithm that can robustly recover the
signal from these measurements. Compressive sensing satisfies (3); here we
demonstrate that radiofrequency-dressed cold atoms can satisfy (1) and (2).

Sensing electrical signals from a living _network_ of neurons in-vitro would
provide information about how nerve cells communicate on the network scale,
considered crucial for understanding diseases such as epilepsy [10].
Conventionally, neural pulses (action potentials) are measured using invasive
patch clamps galvanically connected to individual neurons [11]. The process
perturbs neural connections and can damage cell membranes, while not scaling
beyond microscopic regions of tissue. This motivates ongoing work to develop
non-invasive methods which infer neuronal current waveforms from their
corresponding magnetic fields. Being on the order of nanotesla at sensing
range [11, 3], quantum sensors are natural candidates for magnetometers that
could sense such neuromagnetic signals. Quantum magnetometry of living nerve
cells to date [2, 3, 4, 5, 6] has been achieved with continuous-time ensemble
sensing, averaged over many iterations of a nerve cell being stimulated,
assumed each time to evoke an identical response. In contrast, we demonstrate
a projective measurement protocol, which if extended to a quantum sensing
array, would plausibly sense a unique neural waveform, containing multiple
firings, without averaging.

Neural waveforms are an example of _sparse_ signals: those with only a small
proportion of non-zero values [see Fig. 1(a)]. Consequently, information in
neural waveforms is concentrated in short windows of time, and sampling in
time is inherently inefficient [37]. Any uniform quantum measurement of such
sparse signals will waste quantum resources (e.g., sensing duration of NV s,
or a limited number of projective measurements in a sensing array) by making
measurements that are not very informative.

We instead reduce the number of measurements by taking non-uniform samples in
a discrete sine transform (DST) [38, Eq. (37)] basis. Time sparse signals have
information spread across their broad frequency spectra. Hence, one can learn
all information about sparse signals using only an incomplete subset of
frequency measurements [39], seemingly violating the Nyquist-Shannon sampling
theorem [12, Sec. 7.1]. This is the central idea of _compressive sensing_ [14,
15]. While many possible signals are consistent with a randomly chosen
incomplete set of frequency measurements of a sparse signal, it is very
unlikely any will be as sparse as the true signal [14, Th. 6]. One then
recovers the true signal from incomplete measurements using an optimization
algorithm to find the sparsest signal consistent with the measurements and
with a linear model of the sampling process. A hurdle for applying this to
_quantum_ sensing is that, to measure a frequency coefficient, the sensor
would need to be coherent for the duration of the experiment. Zeeman coherence
times of several seconds make trapped cold atoms a promising candidate for
such a task.

Here, we outline how cold atoms can measure in a DST basis. We then describe a
recovery process that retrieves neural-like signals from an incomplete set of
such measurements. Finally, we demonstrate experimental waveform estimations
derived from incomplete data, acquired with dressed cold atoms exposed to
synthesized neural-like signals.

We model our cold ${}^{87}\mathrm{Rb}$ atom clouds as ensembles of non-
interacting spin-one systems that couple to magnetic field
$\mathbf{{\color[rgb]{.27,.545,0}B}}({\color[rgb]{.31,.58,.804}t})$ via
${\color[rgb]{.804,0,0}\hat{H}}({\color[rgb]{.31,.58,.804}t})=-{\color[rgb]{.545,.545,0}\gamma}\,\mathbf{{\color[rgb]{.27,.545,0}B}}({\color[rgb]{.31,.58,.804}t})\cdot\mathbf{{\color[rgb]{.365,.28,.545}\hat{F}}}$,
where
${\color[rgb]{.365,.28,.545}\hat{F}}_{{\color[rgb]{0.327,0.327,0.3072}X},{\color[rgb]{0.327,0.327,0.3072}Y},{\color[rgb]{0.327,0.327,0.3072}Z}}$
are hyperfine spin-1 operators and ${\color[rgb]{.545,.545,0}\gamma}$ is the
appropriate gyromagnetic ratio.

We aim to measure the Fourier sine coefficient
${\color[rgb]{.545,.545,0}y}({\color[rgb]{.536,.408,.804}f})$ of a signal
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t})$,
i.e.,

$\displaystyle{\color[rgb]{.545,.545,0}y}({\color[rgb]{.536,.408,.804}f})=$
$\displaystyle\frac{1}{{\color[rgb]{.31,.58,.804}T}}\int_{0}^{{\color[rgb]{.31,.58,.804}T}}\sin({\color[rgb]{.545,.396,.03}2\pi}{\color[rgb]{.536,.408,.804}f}{\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{0.436,0.38,0.296}\mathrm{d}}{\color[rgb]{.31,.58,.804}t},$
(1)

for different frequencies ${\color[rgb]{.536,.408,.804}f}$. To do this, we
apply three fields: a bias generating Zeeman splitting by
${\color[rgb]{.804,.408,.224}\omega}_{\color[rgb]{0,.545,.27}L}$, transverse
resonant radiofrequency dressing with Rabi frequency
${\color[rgb]{0.327,0.2808,0.2376}\Omega}$, and the signal field
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t})$
parallel to the bias. The sensor evolves under the resulting Hamiltonian

$\displaystyle{\color[rgb]{.804,0,0}\hat{H}}({\color[rgb]{.31,.58,.804}t})={\color[rgb]{.804,.408,.224}\omega}_{\color[rgb]{0,.545,.27}L}\,{\color[rgb]{.365,.28,.545}\hat{F}}_{\color[rgb]{0.327,0.327,0.3072}Z}+2{\color[rgb]{0.327,0.2808,0.2376}\Omega}\,\cos({\color[rgb]{.804,.408,.224}\omega}_{\color[rgb]{.545,0,0}\mathrm{rf}}{\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{.365,.28,.545}\hat{F}}_{\color[rgb]{0.327,0.327,0.3072}X}+{\color[rgb]{.545,.545,0}\gamma}\,{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{.365,.28,.545}\hat{F}}_{\color[rgb]{0.327,0.327,0.3072}Z},$
(2)

for duration ${\color[rgb]{.31,.58,.804}T}$, before readout of the transverse
spin projection (see Supplemental Material for details [40]). The system
approximately measures the Fourier sine coefficient of
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t})$
at frequency
${\color[rgb]{.536,.408,.804}f}={\color[rgb]{0.327,0.2808,0.2376}\Omega}/{\color[rgb]{.545,.396,.03}2\pi}$,
because the eigenvalues of its dressed frame Hamiltonian are split by
${\color[rgb]{.545,0,0}\hbar}\,{\color[rgb]{0.327,0.2808,0.2376}\Omega}$:

$\displaystyle{\color[rgb]{.804,0,0}\hat{H}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}({\color[rgb]{.31,.58,.804}t})={\color[rgb]{0.327,0.2808,0.2376}\Omega}\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{\color[rgb]{0.327,0.327,0.3072}X}+{\color[rgb]{.545,.545,0}\gamma}\,{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{.365,.28,.545}\hat{F}}_{\color[rgb]{0.327,0.327,0.3072}Z}\xrightarrow{\text{no
signal}}{\color[rgb]{0.327,0.2808,0.2376}\Omega}\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{\color[rgb]{0.327,0.327,0.3072}X}.$
(3)

The sine component is resonant with this splitting, driving transitions
between dressed states at a rate proportional to
${\color[rgb]{.545,.545,0}y}({\color[rgb]{.536,.408,.804}f})$, i.e., a shift
of the dressed-frame plane of Rabi-flopping shown in Fig. 1(b). Dressed state
populations are mapped to lab frame
${\color[rgb]{.365,.28,.545}\hat{F}}_{\color[rgb]{0.327,0.327,0.3072}Z}$
populations ${\color[rgb]{0.436,0.36,0.2656}N}_{+,0,-}$ using a $\pi/2$ pulse
[rf sequence in Fig. 1(b)], which are then measured using Stern-Gerlach
absorption imaging. As a result, we can measure different sine coefficients by
repeating this, each time selecting a different rf amplitude
${\color[rgb]{0.327,0.2808,0.2376}\Omega}$.

Figure 1:  Compressive quantum waveform estimation with rf-dressed atoms: (a)
Trapped atoms exposed to (synthesized) neural (illustrated by phreed [41] on
right) magnetic signal
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$.
Signal sampled in time-domain is illustrated as a color-coded vector
underneath plot. (b) Two example pulse sequences (upper left) differing in
their continuous dressing amplitude generate Rabi flopping at different Rabi
frequencies
${\color[rgb]{0.327,0.2808,0.2376}\Omega}_{i,{\color[rgb]{.365,.28,.545}j}}$
depicted on Bloch sphere (upper right), and as two rows of measurement matrix
${\color[rgb]{.804,.608,.112}A}_{\color[rgb]{.41,.545,.132}\Lambda}$ (lower
left). Unitary evolution with weak signal
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$
and final measurement of transverse spin projection
${\color[rgb]{.545,.545,0}y}_{\color[rgb]{.365,.28,.545}j}$ constitute sensing
a DST coefficient for each ${\color[rgb]{0.327,0.2808,0.2376}\Omega}$, shown
as linear system
${\color[rgb]{.545,.545,0}y}_{\color[rgb]{.41,.545,.132}\Lambda}={\color[rgb]{.804,.608,.112}A}_{\color[rgb]{.41,.545,.132}\Lambda}{\color[rgb]{0.327,0.327,0.3072}x}$
(lower). Protocol modeled as an (exaggerated here) underdetermined matrix
equation (lower). Equality is approximate due to measurement noise and
simplification of model. Experiment is repeated for different rf amplitudes
${\color[rgb]{0.327,0.2808,0.2376}\Omega}_{\color[rgb]{.365,.28,.545}j}$
measuring different sine coefficients
${\color[rgb]{.545,.545,0}y}_{\color[rgb]{.365,.28,.545}j}$. (c) FISTA
searches the set of signals that fit the measurements (left) for the sparsest
such signal, which it returns (right).  Figure 2:  Neural-like magnetic pulses
sensed in time domain, complete frequency domain and compressively. Solid line
is recovered signal, dashed line is signal sent to coils. (a,b) Ramsey
sampling of a (a) one- and (b) two-pulse signal. Signal is not able to be
recovered due to bias drift. (c,d) Inverse DST using the complete set of sine
coefficients recorded by the dressed atoms. (e,f) Compressive recovery using
the FISTA. Only 60 of the 99 recorded sine coefficients are used.

We discretize our model for the measurement protocol in Eq. (1) as

$\displaystyle{\color[rgb]{.545,.545,0}y}_{\color[rgb]{.365,.28,.545}j}=$
$\displaystyle\frac{1}{{\color[rgb]{.31,.58,.804}T}}\sum_{{\color[rgb]{0.327,0.27,0.1992}k}=1}^{{\color[rgb]{0.327,0.27,0.1992}K}-1}\sin[{\color[rgb]{.545,.396,.03}2\pi}({\color[rgb]{.365,.28,.545}j}\,{\color[rgb]{0.436,0.38,0.296}\Delta}{\color[rgb]{.536,.408,.804}f})({\color[rgb]{0.327,0.27,0.1992}k}\,{\color[rgb]{0.436,0.38,0.296}\Delta}{\color[rgb]{.31,.58,.804}t})]\,{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}_{\color[rgb]{0.327,0.27,0.1992}k}\,{\color[rgb]{0.436,0.38,0.296}\Delta}{\color[rgb]{.31,.58,.804}t},$
(4)

for
${\color[rgb]{.365,.28,.545}j}\in\\{1,2,\dots,{\color[rgb]{0.327,0.27,0.1992}K}-1\\}$.
Here the signal time series
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}_{\color[rgb]{0.327,0.27,0.1992}k}$
starts at ${\color[rgb]{0.327,0.27,0.1992}k}=1$ and has length
${\color[rgb]{0.327,0.27,0.1992}K}-1$. By choosing bandwidth as
${\color[rgb]{.804,.408,.224}W}\stackrel{{\scriptstyle{\color[rgb]{0.436,0.38,0.296}\text{def}}}}{{=}}{\color[rgb]{0.327,0.27,0.1992}K}\,{\color[rgb]{0.436,0.38,0.296}\Delta}{\color[rgb]{.536,.408,.804}f}=1/(2\,{\color[rgb]{0.436,0.38,0.296}\Delta}{\color[rgb]{.31,.58,.804}t})$
and sensing duration
${\color[rgb]{.31,.58,.804}T}\stackrel{{\scriptstyle{\color[rgb]{0.436,0.38,0.296}\text{def}}}}{{=}}{\color[rgb]{0.327,0.27,0.1992}K}\,{\color[rgb]{0.436,0.38,0.296}\Delta}{\color[rgb]{.31,.58,.804}t}=1/(2\,{\color[rgb]{0.436,0.38,0.296}\Delta}{\color[rgb]{.536,.408,.804}f})$,
Eq. (4) becomes the DST [38, Eq. (37)]; i.e., the matrix equation
${\color[rgb]{.545,.545,0}y}={\color[rgb]{.804,.608,.112}A}{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$,
with
${\color[rgb]{.804,.608,.112}A}_{{\color[rgb]{.365,.28,.545}j}{\color[rgb]{0.327,0.27,0.1992}k}}=\sin(\pi{\color[rgb]{.365,.28,.545}j}{\color[rgb]{0.327,0.27,0.1992}k}/{\color[rgb]{0.327,0.27,0.1992}K})/{\color[rgb]{0.327,0.27,0.1992}K}$.
Compressive sensing theory roughly says that, since
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$
is sparse, and the DCT is sensitive to all sparse signals, we can subsample
our measurements and still be able to recover the signal [14, Th. 6]: i.e.,
measure only ${\color[rgb]{.545,0,0}M}<{\color[rgb]{0.327,0.27,0.1992}K}$
components ${\color[rgb]{.545,.545,0}y}_{\color[rgb]{.365,.28,.545}j}$ of
${\color[rgb]{.545,.545,0}y}$, so that now
${\color[rgb]{.365,.28,.545}j}\in{\color[rgb]{.41,.545,.132}\Lambda}\subset\\{1,2,\dots,{\color[rgb]{0.327,0.27,0.1992}K}-1\\}$.
Equation (4) still holds, but only for
${\color[rgb]{.365,.28,.545}j}\in{\color[rgb]{.41,.545,.132}\Lambda}$; written
${\color[rgb]{.545,.545,0}y}_{\color[rgb]{.41,.545,.132}\Lambda}={\color[rgb]{.804,.608,.112}A}_{\color[rgb]{.41,.545,.132}\Lambda}{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$,
and illustrated as the matrix of sine waves of few randomly chosen frequencies
in Fig. 1(b). Therefore, we know that the signal must be one of the many
solutions ${\color[rgb]{0.327,0.327,0.3072}x}$ to the underdetermined system
${\color[rgb]{.545,.545,0}y}_{\color[rgb]{.41,.545,.132}\Lambda}={\color[rgb]{.804,.608,.112}A}_{\color[rgb]{.41,.545,.132}\Lambda}{\color[rgb]{0.327,0.327,0.3072}x}$.

This is not uniquely solvable, but because
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$
is sparse, the correct recovery,
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\star}}$,
will, with high probability, be the sparsest such
${\color[rgb]{0.327,0.327,0.3072}x}$ [14, Th. 6]. In general this is an NP-
hard problem [15, Th. 2.17], but assuming that
${\color[rgb]{.804,.608,.112}A}_{\color[rgb]{.41,.545,.132}\Lambda}$ satisfies
criteria involving the restricted isometry property (RIP), the sparest
solution of
${\color[rgb]{.545,.545,0}y}_{\color[rgb]{.41,.545,.132}\Lambda}={\color[rgb]{.804,.608,.112}A}_{\color[rgb]{.41,.545,.132}\Lambda}{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$
will also be the solution which minimizes the $\ell_{1}$-norm,
$\|{\color[rgb]{0.327,0.327,0.3072}x}\|_{1}\stackrel{{\scriptstyle{\color[rgb]{0.436,0.38,0.296}\text{def}}}}{{=}}\sum_{\color[rgb]{0.327,0.27,0.1992}k}|{\color[rgb]{0.327,0.327,0.3072}x}_{\color[rgb]{0.327,0.27,0.1992}k}|$
[14, Th. 9]. These criteria translate to requirements of measuring a system-
dependent minimum number of sine coefficients [15, Eq. (9.24)] and for the
subsampling of frequencies to have no structure (which is likely if they are
chosen randomly [15, Th. 11.23]). In practice, our linear model [Eq. (4)] will
only be approximately satisfied due to measurement noise. Hence, we relax to a
problem called the _least absolute shrinkage and selection operator (LASSO)_
[42, Eq. 3.14], where we find the
${\color[rgb]{0.327,0.327,0.3072}x}\in{\color[rgb]{.932,0,0}\mathbb{R}}^{{\color[rgb]{0.327,0.27,0.1992}K}-1}$
that minimizes
$\|{\color[rgb]{.804,.608,.112}A}_{\color[rgb]{.41,.545,.132}\Lambda}{\color[rgb]{0.327,0.327,0.3072}x}-{\color[rgb]{.545,.545,0}y}_{\color[rgb]{.41,.545,.132}\Lambda}\|_{2}^{2}+{\color[rgb]{.604,.804,.196}\lambda}\|{\color[rgb]{0.327,0.327,0.3072}x}\|_{1}$,
illustrated in Fig. 1(c). Here, the _regularization parameter_
${\color[rgb]{.604,.804,.196}\lambda}>0$ [14, Sec. 8.2.2] is chosen [40] to
produce the smallest error on simulated (using our open-source solver [43])
training data, and
$\|{\color[rgb]{.545,.545,0}y}_{\color[rgb]{.41,.545,.132}\Lambda}\|_{2}^{2}\stackrel{{\scriptstyle{\color[rgb]{0.436,0.38,0.296}\text{def}}}}{{=}}\sum_{{\color[rgb]{.365,.28,.545}j}\in{\color[rgb]{.41,.545,.132}\Lambda}}|{\color[rgb]{.545,.545,0}y}_{\color[rgb]{.365,.28,.545}j}|^{2}$
is the $\ell_{2}$-norm. To solve the LASSO, we implemented the _fast iterative
shrinkage-thresholding algorithm (FISTA)_ [44]. We chose the step-size of
FISTA using the singular values of
${\color[rgb]{.804,.608,.112}A}_{\color[rgb]{.41,.545,.132}\Lambda}$ and an an
estimate of the amplitude of
${\color[rgb]{0.327,0.327,0.3072}x}({\color[rgb]{.31,.58,.804}t})$.

To experimentally verify this sensing method, we applied the Hamiltonian Eq.
(2) to a dipole-trapped cloud of approximately $6\times 10^{5}$ laser cooled
atoms of ${}^{87}\mathrm{Rb}$ at a temperature of $1.0\ \mu\mathrm{K}$. A bias
field of $86.1\ \mu\mathrm{T}$, gave a Larmor splitting of
${\color[rgb]{.804,.408,.224}\omega}_{\color[rgb]{0,.545,.27}L}={\color[rgb]{.545,.396,.03}2\pi}\times
603\ \mathrm{kHz}$, and we tuned our rf dressing
${\color[rgb]{.804,.408,.224}\omega}_{\color[rgb]{.545,0,0}\mathrm{rf}}$ to be
resonant with this. The neural-like signal
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t})$
comprised one or two pulses, each being a single cycle sine pulse of amplitude
$143\ \mathrm{nT}$ and duration $200\ \mu\mathrm{s}$. The amplitude is much
weaker than the $50\ \mu\mathrm{T}$ pulse used by Webb  _et al._ [5], though
still stronger than sub-nanotesla fields expected in close proximity to
individual neurons. The pulse duration is much longer than the $10\
\mu\mathrm{s}$ pulse simulated by Magesan  _et al._ [36], though still shorter
than the approximately $2\ \mathrm{ms}$ duration of a typical action
potential. Starting in
${{\color[rgb]{.365,.28,.545}F}}=1,{\color[rgb]{.545,0,0}m}_{{\color[rgb]{.365,.28,.545}F}}=-1$
(a superposition of dressed states), the sensor evolves under Eq. (2), before
Zeeman populations ${\color[rgb]{0.436,0.36,0.2656}N}_{+,0,-}$ are recorded
using Stern-Gerlach absorption imaging after a readout $\pi/2$ pulse. This is
linearly sensitive to small
${\color[rgb]{.545,.545,0}y}({\color[rgb]{.536,.408,.804}f})$ via
${\color[rgb]{.545,.545,0}y}({\color[rgb]{.536,.408,.804}f})=\braket{{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}X}}}({\color[rgb]{.31,.58,.804}T})/({\color[rgb]{.545,0,0}\hbar}\,{\color[rgb]{.545,.545,0}\gamma}\,{\color[rgb]{.31,.58,.804}T})=\braket{{\color[rgb]{.365,.28,.545}\hat{F}}_{\color[rgb]{0.327,0.327,0.3072}Z}}_{\frac{\pi}{2}}/({\color[rgb]{.545,0,0}\hbar}\,{\color[rgb]{.545,.545,0}\gamma}\,{\color[rgb]{.31,.58,.804}T})=({\color[rgb]{0.436,0.36,0.2656}N}_{-}-{\color[rgb]{0.436,0.36,0.2656}N}_{+})/[{\color[rgb]{.545,.545,0}\gamma}\,{\color[rgb]{.31,.58,.804}T}\,({\color[rgb]{0.436,0.36,0.2656}N}_{-}+{\color[rgb]{0.436,0.36,0.2656}N}_{0}+{\color[rgb]{0.436,0.36,0.2656}N}_{+})]$.

The procedure was repeated, changing the value of
${\color[rgb]{0.327,0.2808,0.2376}\Omega}$ for each shot to record each
frequency component from the incomplete set
${\color[rgb]{.41,.545,.132}\Lambda}$. Our modeled neural-like pulse had
negligible sine components above ${\color[rgb]{.804,.408,.224}W}=10\
\mathrm{kHz}$ (consistent with Webb  _et al._ [3]) which fixed our highest
meaningful time resolution to be
${\color[rgb]{0.436,0.38,0.296}\Delta}{\color[rgb]{.31,.58,.804}t}=50\
\mu\mathrm{s}$ (pulses 4 samples long). We chose an experiment duration of
${\color[rgb]{.31,.58,.804}T}=5\ \mathrm{ms}$ (99 samples long), yielding a
frequency resolution of
${\color[rgb]{0.436,0.38,0.296}\Delta}{\color[rgb]{.536,.408,.804}f}=100\
\mathrm{Hz}$. We do this to sense around the crest of a single cycle of the
electrical line.

We wished to compare measurements using our compressive protocol with complete
measurements of both time and frequency samples. For the former, we used a
separate Ramsey sequence to sample the amplitude of the signal at all of the
99 points on the time-grid. Here the second $\pi/2$ pulse was in quadrature to
the first in order for the populations to be linearly sensitive to the signal
(see the coherent signal case in Mouradian  _et al._ [37]). The Ramsey
sequences approximate capturing the field at a fixed time by sensing over a
rectangular window of length $60\ \mu\mathrm{s}$ centerd on the time-grid
points. For the latter, we used the above protocol to measure a complete set
of 99 sine coefficients and recovered using a complete inverse DST. We then
obtained incomplete data
${\color[rgb]{.545,.545,0}y}_{\color[rgb]{.41,.545,.132}\Lambda}$ for the
compressive protocol by various subsamplings of these 99 sine coefficients
${\color[rgb]{.545,.545,0}y}$. This enables us to directly compare complete
inverse DST s with compressive retrievals. Of course, to benefit from
compressive sensing, one would instead take fewer measurements in the first
place. We present such experiments in Supplemental Material [40].

Figure 3:  Finding the fewest number of samples that can be used by the
compressive protocol before it fails in recovering the signal. AUC (see
Supplemental Material [40]) of $100\%$ means perfect pulse detection, $50\%$
means random categorization. AUC s from 200 random subsets of the complete
datasets were averaged to form each data point; shading shows one standard
deviation either side. One-pulse traces (light orange) draw from sine
coefficients used in Fig. 2(c), two-pulse traces (dark red) from Fig. 2(d).
Left and right insets show (average over the 200 random subsets) ROC curves
for the two-pulse signal when 20 and 60 measurements are used in the recovery
respectively. ROC s are parameterized by categorization threshold. If a
threshold perfectly distinguishes between pulses and noise, the ROC reaches
the top left and thus encloses an area (AUC plotted in main figure) of
$100\%$.

Figure 2 compares waveform measurements of repeated one- and two-pulse signals
using the Ramsey protocol, an inverse DST of the complete set of sine
coefficients, and the compressive protocol using a random set of 60 of the 99
sine coefficients. The Ramsey measurements in Figs. 2(a, b) are unable to
detect the pulses. The recovery is overwhelmed by shot-to-shot drift of the
bias magnetic field, which we later discovered was caused by suburban railway
lines $2\ \mathrm{km}$ away. Our choice to shuffle the shot order made these
slow drifts appear as white noise. Another source of deviation includes the
interference from the ac electrical lines. Furthermore, the square sampling
window of the Ramsey measurements makes the waveform estimation susceptible to
high frequency noise from surrounding instrumentation. In contrast, the
inverse DST retrieval in Figs. 2(c, d) were able to resolve the synthesized
pulses, as the DST successfully cuts out unwanted dc and high frequency
coefficients. The compressive retrievals in Figs. 2(e, f), on the other hand,
result in the cleanest signals of all the techniques, despite considering only
60 of the complete 99 sine coefficients. This is expected from the de-noising
property of the LASSO [42, Fig. 6(a)].

Compressive sensing protocols have lower bounds as to how few samples can be
used to recover a signal [15, Eq. (9.24)]. For our protocol, we wanted to
experimentally determine this limit for our protocol and dataset. Here we
compare our recovered signals to the waveform we commanded the magnetic coils
with. While the root-mean-square error is a useful metric in other contexts,
for our sparse signals we found that it had difficulty distinguishing noisy,
but otherwise _successful_ , recoveries from completely failed recoveries.
Canonical analysis of detecting signals in noise involves plotting a curve
called the receiver operating characteristic (ROC); the area under the curve
(AUC) [45, Sec. 7] is an appropriate metric for determining signal recovery
success [40]. A signal detector that never confuses signal and noise achieves
an AUC of $100\%$. In Fig. 3, the compressive protocol has AUC over $99\%$
when at least $36$ samples are used to recover the signal containing one
pulse, and at least $52$ samples for the signal containing two. For our
parameters these echo the established theoretical bounds of $34$ and $57$
samples respectively [15, Eq. (9.24)]. We thus predicted that the protocol
would work with $60$ coefficients measured, and it did, as shown in
Supplemental Material [40]. While a reduction in quantum resources by almost
half is already useful, the theoretical model predicts even larger relative
gains when measuring longer, sparser signals.

The compressive quantum sensor was able to recover neural-like signals in a
situation where Ramsey measurements could not, while using an incomplete set
of measurements. Choice of signal basis can drastically change both the
fidelity, and, when undersampling, efficiency of quantum waveform estimation.
While the implementation presented in this Letter arrays measurement in time
to estimate a repeated waveform, we envision an implementation on a quantum
sensing array of independently-controlled sensors arrayed in space, where
resource parsimony is of the essence. Such an array would be capable of
compressive waveform estimation of a never-to-be-repeated signal, in a single
shot.

We thank Anton van den Hengel, Kristian Helmerson and David Spanswick for
helpful discussions. AT, CCB, and HAMT acknowledge support through Australian
Government Research Training Program (RTP)scholarships. JS is the recipient of
an Australian Research Council Discovery Early Career Researcher Award
(project number DE210101056) funded by the Australian Government. LDT
acknowledges funding from the Australian Research Council Linkage Project
(project number LP200100082).

## References

  * Tsang  _et al._ [2011] M. Tsang, H. M. Wiseman, and C. M. Caves, Physical Review Letters 106, 090401 (2011).
  * Barry  _et al._ [2016] J. F. Barry, M. J. Turner, J. M. Schloss, D. R. Glenn, Y. Song, M. D. Lukin, H. Park, and R. L. Walsworth, Proceedings of the National Academy of Sciences 113, 14133 (2016).
  * Webb  _et al._ [2021] J. L. Webb, L. Troise, N. W. Hansen, C. Olsson, A. M. Wojciechowski, J. Achard, O. Brinza, R. Staacke, M. Kieschnick, J. Meijer, A. Thielscher, J.-F. Perrier, K. Berg-Sørensen, A. Huck, and U. L. Andersen, Scientific Reports 11, 2412 (2021).
  * Troise  _et al._ [2021] L. Troise, N. W. Hansen, C. Olsson, J. L. Webb, L. Tomasevic, J. Achard, O. Brinza, R. Staacke, M. Kieschnick, J. Meijer, A. Thielscher, H. R. Siebner, K. Berg-Sørensen, J.-F. Perrier, A. Huck, and U. L. Andersen, Laser stimulation of muscle activity with simultaneous detection using a diamond colour centre biosensor (2021), arxiv:2112.09516 [physics.bio-ph] .
  * Webb  _et al._ [2022] J. L. Webb, L. Troise, N. W. Hansen, L. F. Frellsen, C. Osterkamp, F. Jelezko, S. Jankuhn, J. Meijer, K. Berg-Sørensen, J.-F. Perrier, A. Huck, and U. L. Andersen, Physical Review Applied 17, 064051 (2022).
  * Hansen  _et al._ [2022] N. W. Hansen, J. L. Webb, L. Troise, C. Olsson, L. Tomasevic, O. Brinza, J. Achard, R. Staacke, M. Kieschnick, J. Meijer, A. Thielscher, H. R. Siebner, K. Berg-Sørensen, J.-F. Perrier, A. Huck, and U. L. Andersen, Microscopic-scale recording of brain neuronal electrical activity using a diamond quantum sensor (2022), arxiv:2208.14068 [physics, physics:quant-ph] .
  * Xu  _et al._ [2016] N. Xu, F. Jiang, Y. Tian, J. Ye, F. Shi, H. Lv, Y. Wang, J. Wrachtrup, and J. Du, Physical Review B 93, 161117 (2016).
  * Karadas  _et al._ [2018] M. Karadas, A. M. Wojciechowski, A. Huck, N. O. Dalby, U. L. Andersen, and A. Thielscher, Scientific Reports 8, 4503 (2018).
  * Parashar  _et al._ [2020] M. Parashar, K. Saha, and S. Bandyopadhyay, Communications physics 3, 174 (2020).
  * McDonald [2012] C. M. McDonald, Physical Medicine and Rehabilitation Clinics of North America 23, 495 (2012).
  * Mitterdorfer and Bean [2002] J. Mitterdorfer and B. P. Bean, Journal of Neuroscience 22, 10106 (2002).
  * Oppenheim and Willsky [1997] A. V. Oppenheim and A. S. Willsky,  _Signals & Systems_ (Pearson Education, 1997).
  * Donoho [2006] D. Donoho, IEEE Transactions on Information Theory 52, 1289 (2006).
  * Koep  _et al._ [2019] N. Koep, A. Behboodi, and R. Mathar, in  _Compressed Sensing and Its Applications: Third International MATHEON Conference 2017_, Applied and Numerical Harmonic Analysis, edited by H. Boche, G. Caire, R. Calderbank, G. Kutyniok, R. Mathar, and P. Petersen (Springer International Publishing, Cham, 2019) pp. 1–65.
  * Foucart and Rauhut [2013] S. Foucart and H. Rauhut,  _A Mathematical Introduction to Compressive Sensing_ , 2013th ed. (Birkhäuser, New York, 2013).
  * Duarte  _et al._ [2008] M. F. Duarte, M. A. Davenport, D. Takhar, J. N. Laska, T. Sun, K. F. Kelly, and R. G. Baraniuk, IEEE Signal Processing Magazine 25, 83 (2008).
  * Graff and Sidky [2015] C. G. Graff and E. Y. Sidky, Applied optics 54, C23 (2015).
  * Teo and Sánchez-Soto [2021] Y. S. Teo and L. L. Sánchez-Soto, International Journal of Quantum Information 19, 2140003 (2021).
  * Ran  _et al._ [2020] S.-J. Ran, Z.-Z. Sun, S.-M. Fei, G. Su, and M. Lewenstein, Physical Review Research 2, 033293 (2020).
  * Sherbert  _et al._ [2022] K. M. Sherbert, N. Naimipour, H. Safavi, H. C. Shaw, and M. Soltanalian, Applied Sciences 12, 7525 (2022).
  * Fontana  _et al._ [2022] E. Fontana, I. Rungger, R. Duncan, and C. Cîrstoiu, Efficient recovery of variational quantum algorithms landscapes using classical signal processing (2022), arxiv:2208.05958 [quant-ph] .
  * Seif  _et al._ [2021] A. Seif, M. Hafezi, and Y.-K. Liu, Compressed Sensing Measurement of Long-Range Correlated Noise (2021), arxiv:2105.12589 [quant-ph] .
  * Huang  _et al._ [2023] K. Huang, D. Farfurnik, A. Seif, M. Hafezi, and Y.-K. Liu, Random Pulse Sequences for Qubit Noise Spectroscopy (2023), arxiv:2303.00909 [quant-ph] .
  * Ayanzadeh  _et al._ [2020] R. Ayanzadeh, M. Halem, and T. Finin, in  _IGARSS 2020 - 2020 IEEE International Geoscience and Remote Sensing Symposium_ (2020) pp. 3517–3520.
  * Aonishi  _et al._ [2022] T. Aonishi, K. Mimura, M. Okada, and Y. Yamamoto, Quantum Science and Technology 7, 035013 (2022).
  * Ayanzadeh  _et al._ [2022] R. Ayanzadeh, J. Dorband, M. Halem, and T. Finin, in  _IGARSS 2022 - 2022 IEEE International Geoscience and Remote Sensing Symposium_ (2022) pp. 4911–4914.
  * Gunathilaka  _et al._ [2023] M. D. S. H. Gunathilaka, S. Kako, Y. Inui, K. Mimura, M. Okada, Y. Yamamoto, and T. Aonishi, Effective implementation of $l_0$-Regularised Compressed Sensing with Chaotic-Amplitude-Controlled Coherent Ising Machines, https://arxiv.org/abs/2302.12523v1 (2023).
  * Katz  _et al._ [2009] O. Katz, Y. Bromberg, and Y. Silberberg, Applied Physics Letters 95, 131110 (2009).
  * Johnson  _et al._ [2020] S. D. Johnson, P.-A. Moreau, T. Gregory, and M. J. Padgett, Applied Physics Letters 116, 260504 (2020).
  * Cao  _et al._ [2021] D.-Z. Cao, S.-H. Zhang, Y. Zhao, C. Ren, J. Zhang, B. Liang, B. Sun, and K. Wang, arXiv:2108.01037 [quant-ph] (2021), arxiv:2108.01037 [quant-ph] .
  * Arai  _et al._ [2015] K. Arai, C. Belthangady, H. Zhang, N. Bar-Gill, S. J. DeVience, P. Cappellaro, A. Yacoby, and R. L. Walsworth, Nature Nanotechnology 10, 859 (2015).
  * Holland  _et al._ [2011] D. J. Holland, M. J. Bostock, L. F. Gladden, and D. Nietlispach, Angewandte Chemie International Edition 50, 6548 (2011).
  * Müller  _et al._ [2014] C. Müller, X. Kong, J.-m. Cai, K. Melentijevic, A. Stacey, M. Markham, D. Twitchen, J. Isoya, S. Pezzagna, J. Meijer, J. F. Du, M. B. Plenio, B. Naydenov, L. P. Mcguinness, and F. Jelezko, Nature Communications 5, 4703 (2014).
  * Puentes  _et al._ [2015] G. Puentes, G. Waldherr, P. Neumann, G. Balasubramanian, and J. Wrachtrup, Scientific Reports 4, 4677 (2015).
  * Boss  _et al._ [2017] J. M. Boss, K. S. Cujia, J. Zopes, and C. L. Degen, Science 356, 837 (2017).
  * Magesan  _et al._ [2013] E. Magesan, A. Cooper, and P. Cappellaro, Physical Review A 88, 062109 (2013).
  * Mouradian  _et al._ [2021] S. L. Mouradian, N. Glikin, E. Megidish, K.-I. Ellers, and H. Haeffner, Physical Review A 103, 032419 (2021).
  * Jain [1976] A. Jain, IEEE Transactions on Communications 24, 1023 (1976).
  * Donoho and Stark [1989] D. L. Donoho and P. B. Stark, SIAM Journal on Applied Mathematics 49, 906 (1989).
  * Note [0] See Supplemental Material at _[URL will be inserted by publisher]_ for a derivation of sensitivity to sine coefficients, an explanation of regularization parameter tuning, an explanation of how we use the AUC statistic from signal detection theory, and recoveries from genuinely undersampled measurements.
  * phreed [2006] phreed, Single-neuron, openclipart (2006).
  * Santosa and Symes [1986] F. Santosa and W. W. Symes, SIAM Journal on Scientific and Statistical Computing 7, 1307 (1986).
  * Tritt  _et al._ [2023] A. Tritt, J. Morris, J. Hochstetter, R. P. Anderson, J. Saunderson, and L. D. Turner, Computer Physics Communications 287, 108701 (2023).
  * Beck and Teboulle [2009] A. Beck and M. Teboulle, SIAM Journal on Imaging Sciences 2, 183 (2009).
  * Fawcett [2006] T. Fawcett, Pattern Recognition Letters ROC Analysis in Pattern Recognition, 27, 861 (2006).
  * Blanes  _et al._ [2009] S. Blanes, F. Casas, J. A. Oteo, and J. Ros, Physics Reports 470, 151 (2009).
  * J J Sakurai (1994) [Jun John] J J Sakurai (Jun John),  _Modern Quantum Mechanics_ , rev. ed. ed. (Addison-Wesley PubCo, Reading, Mass., 1994).
  * Siminovitch  _et al._ [2003] D. Siminovitch, T. Untidt, and N. C. Nielsen, The Journal of Chemical Physics 120, 51 (2003).
  * Kay [1998] S. Kay,  _Fundamentals of Statistical Signal Processing, Volume II: Detection Theory, Volume 2_ (Pearson Education, Englewood Cliffs, N.J, 1998).

Supplemental Material

##  I Protocol for sensing Fourier sine coefficients

The lab frame Hamiltonian of our atoms is

$\displaystyle{\color[rgb]{.804,0,0}\hat{H}}({\color[rgb]{.31,.58,.804}t})={\color[rgb]{.804,.408,.224}\omega}_{\color[rgb]{0,.545,.27}L}\,{\color[rgb]{.365,.28,.545}\hat{F}}_{\color[rgb]{0.327,0.327,0.3072}Z}+2{\color[rgb]{0.327,0.2808,0.2376}\Omega}\,\cos({\color[rgb]{.804,.408,.224}\omega}_{\color[rgb]{.545,0,0}\mathrm{rf}}{\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{.365,.28,.545}\hat{F}}_{\color[rgb]{0.327,0.327,0.3072}X}+{\color[rgb]{.545,.545,0}\gamma}\,{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{.365,.28,.545}\hat{F}}_{\color[rgb]{0.327,0.327,0.3072}Z}.$
(5)

Moving into a frame ${\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}$ rotating at
${\color[rgb]{.804,.408,.224}\omega}_{\color[rgb]{.545,0,0}\mathrm{rf}}={\color[rgb]{.804,.408,.224}\omega}_{\color[rgb]{0,.545,.27}L}$
(on resonance) around
${\color[rgb]{.365,.28,.545}\hat{F}}_{\color[rgb]{0.327,0.327,0.3072}Z}$, and
taking the rotating wave approximation (RWA), the system has a Hamiltonian of

$\displaystyle{\color[rgb]{.804,0,0}\hat{H}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}({\color[rgb]{.31,.58,.804}t})={\color[rgb]{0.327,0.2808,0.2376}\Omega}\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{\color[rgb]{0.327,0.327,0.3072}X}+{\color[rgb]{.545,.545,0}\gamma}\,{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{.365,.28,.545}\hat{F}}_{\color[rgb]{0.327,0.327,0.3072}Z}.$
(6)

We can now move into a second rotating frame
${\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}$
at ${\color[rgb]{0.327,0.2808,0.2376}\Omega}$ around
${\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{\color[rgb]{0.327,0.327,0.3072}X}$,
but without taking a second RWA:

$\displaystyle{\color[rgb]{.804,0,0}\hat{H}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}({\color[rgb]{.31,.58,.804}t})={\color[rgb]{.545,.545,0}\gamma}\,{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t})\left(\sin({\color[rgb]{0.327,0.2808,0.2376}\Omega}\,{\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{\color[rgb]{0.327,0.327,0.3072}Y}-\cos({\color[rgb]{0.327,0.2808,0.2376}\Omega}\,{\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{\color[rgb]{0.327,0.327,0.3072}Z}\right).$
(7)

We assume that the signal
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$
is always small enough such that time evolution can be well approximated by a
truncated Magnus series
$\sum_{\color[rgb]{.545,0,0}m}{\color[rgb]{0.327,0.2808,0.2376}\hat{\varOmega}}_{\color[rgb]{.545,0,0}m}$
[46]. Then the time evolution operator of the spin in the second rotating
frame over the sensing duration ${\color[rgb]{.31,.58,.804}T}$ is

$\displaystyle{\color[rgb]{.804,.4,0}\hat{\mathcal{U}}}=$
$\displaystyle\exp\left({\color[rgb]{0.327,0.2808,0.2376}\hat{\varOmega}}_{1}+{\color[rgb]{0.327,0.2808,0.2376}\hat{\varOmega}}_{2}+{\color[rgb]{0.327,0.2808,0.2376}\hat{\varOmega}}_{3}+\cdots\right),$
(8) $\displaystyle{\color[rgb]{.804,.4,0}\hat{\mathcal{U}}}_{1}=$
$\displaystyle\exp\left({\color[rgb]{0.327,0.2808,0.2376}\hat{\varOmega}}_{1}\right)\quad\text{to
first order},$ (9) $\displaystyle=$
$\displaystyle\exp\left(\frac{i\,{\color[rgb]{.545,.545,0}\gamma}}{{\color[rgb]{.545,0,0}\hbar}}\int_{0}^{{\color[rgb]{.31,.58,.804}T}}{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t})\left(\sin({\color[rgb]{0.327,0.2808,0.2376}\Omega}\,{\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{\color[rgb]{0.327,0.327,0.3072}Y}-\cos({\color[rgb]{0.327,0.2808,0.2376}\Omega}\,{\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{\color[rgb]{0.327,0.327,0.3072}Z}\right){\color[rgb]{0.436,0.38,0.296}\mathrm{d}}{\color[rgb]{.31,.58,.804}t}\right)\quad\text{from
\cite[cite]{[\@@bibref{Number}{blanes_magnus_2009}{}{},
Eq.\leavevmode\nobreak\ (42)]}},$ (10) $\displaystyle=$
$\displaystyle\exp\left(\frac{i\,{\color[rgb]{.545,.545,0}\gamma}}{{\color[rgb]{.545,0,0}\hbar}}\left(\int_{0}^{{\color[rgb]{.31,.58,.804}T}}\sin({\color[rgb]{0.327,0.2808,0.2376}\Omega}\,{\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{0.436,0.38,0.296}\mathrm{d}}{\color[rgb]{.31,.58,.804}t}\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{\color[rgb]{0.327,0.327,0.3072}Y}-\int_{0}^{{\color[rgb]{.31,.58,.804}T}}\cos({\color[rgb]{0.327,0.2808,0.2376}\Omega}\,{\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{0.436,0.38,0.296}\mathrm{d}}{\color[rgb]{.31,.58,.804}t}\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{\color[rgb]{0.327,0.327,0.3072}Z}\right)\right),$
(11) $\displaystyle=$
$\displaystyle\exp\left(\frac{i}{{\color[rgb]{.545,0,0}\hbar}}\left({\color[rgb]{0.1216,0.2,0.436}S}\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{\color[rgb]{0.327,0.327,0.3072}Y}-{\color[rgb]{0,.545,.545}C}\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{\color[rgb]{0.327,0.327,0.3072}Z}\right)\right).$
(12)

Here
${\color[rgb]{0.1216,0.2,0.436}S}={\color[rgb]{.545,.545,0}\gamma}\int_{0}^{{\color[rgb]{.31,.58,.804}T}}\sin({\color[rgb]{0.327,0.2808,0.2376}\Omega}\,{\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{0.436,0.38,0.296}\mathrm{d}}{\color[rgb]{.31,.58,.804}t}={\color[rgb]{.545,.545,0}\gamma}\,{\color[rgb]{.31,.58,.804}T}\,{\color[rgb]{.545,.545,0}y}({\color[rgb]{0.327,0.2808,0.2376}\Omega}/{\color[rgb]{.545,.396,.03}2\pi})$
is proportional to the Fourier sine coefficient
${\color[rgb]{.545,.545,0}y}({\color[rgb]{0.327,0.2808,0.2376}\Omega}/{\color[rgb]{.545,.396,.03}2\pi})$
of the signal
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$
at frequency
${\color[rgb]{0.327,0.2808,0.2376}\Omega}/{\color[rgb]{.545,.396,.03}2\pi}$,
which is the quantity we want to measure. Similarly
${\color[rgb]{0,.545,.545}C}={\color[rgb]{.545,.545,0}\gamma}\int_{0}^{{\color[rgb]{.31,.58,.804}T}}\cos({\color[rgb]{0.327,0.2808,0.2376}\Omega}\,{\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t})\,{\color[rgb]{0.436,0.38,0.296}\mathrm{d}}{\color[rgb]{.31,.58,.804}t}$
is proportional to the Fourier cosine coefficient of
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$
at frequency
${\color[rgb]{.536,.408,.804}f}={\color[rgb]{0.327,0.2808,0.2376}\Omega}/{\color[rgb]{.545,.396,.03}2\pi}$.
Note that the second term in the Magnus expansion is [46, Eq. (43)]

$\displaystyle{\color[rgb]{0.327,0.2808,0.2376}\hat{\varOmega}}_{2}$
$\displaystyle=\frac{i\,{\color[rgb]{.545,.545,0}\gamma}^{2}}{2\,{\color[rgb]{.545,0,0}\hbar}}\int_{0}^{{\color[rgb]{.31,.58,.804}T}}\int_{0}^{{\color[rgb]{.31,.58,.804}t}_{1}}{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t}_{1})\,{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}({\color[rgb]{.31,.58,.804}t}_{2})\,\sin\left({\color[rgb]{0.327,0.2808,0.2376}\Omega}\,({\color[rgb]{.31,.58,.804}t}_{1}-{\color[rgb]{.31,.58,.804}t}_{2})\right)\,{\color[rgb]{0.436,0.38,0.296}\mathrm{d}}{\color[rgb]{.31,.58,.804}t}_{2}\,{\color[rgb]{0.436,0.38,0.296}\mathrm{d}}{\color[rgb]{.31,.58,.804}t}_{1}\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{\color[rgb]{0.327,0.327,0.3072}X},$
(13)

which is suppressed by being of order
$({\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}})^{2}$,
and because it drives rotations in the same axis as the dressing in the first
rotating frame. Hence we only consider the dynamics due to
${\color[rgb]{.804,.4,0}\hat{\mathcal{U}}}_{1}$.

We show later [in Eq. (18)] that we can measure
$\braket{{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}X}}({\color[rgb]{.31,.58,.804}T})}=\bra{{\color[rgb]{.545,.545,0}\psi}({\color[rgb]{.31,.58,.804}T})}{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}X}}\ket{{\color[rgb]{.545,.545,0}\psi}({\color[rgb]{.31,.58,.804}T})}=\bra{{\color[rgb]{.545,.545,0}\psi}(0)}\,{\color[rgb]{.804,.4,0}\hat{\mathcal{U}}}_{1}^{\dagger}\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}X}}\,{\color[rgb]{.804,.4,0}\hat{\mathcal{U}}}_{1}\ket{{\color[rgb]{.545,.545,0}\psi}(0)}$
for an experiment starting in some initial state
$\ket{{\color[rgb]{.545,.545,0}\psi}(0)}$. First, note that the unitary
${\color[rgb]{.804,.4,0}\hat{\mathcal{U}}}_{1}$ is a rotation on the Bloch
sphere of angle
$\sqrt{{\color[rgb]{0.1216,0.2,0.436}S}^{2}+{\color[rgb]{0,.545,.545}C}^{2}}$
around an axis lying on the
${\color[rgb]{0.327,0.327,0.3072}Y}{\color[rgb]{0.327,0.327,0.3072}Z}$ plane
at an angle
$\mathrm{arctan2}\left({\color[rgb]{0.1216,0.2,0.436}S},{\color[rgb]{0,.545,.545}C}\right)$
to the ${\color[rgb]{0.327,0.327,0.3072}Z}$ axis. Second, note that the
operator
${\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}X}}$
we are taking an expectation value of is a component of a vector operator
$\mathbf{{\color[rgb]{.365,.28,.545}\hat{F}}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}$.
Therefore, we can simply evaluate the similarity transform
${\color[rgb]{.804,.4,0}\hat{\mathcal{U}}}_{1}^{\dagger}\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}X}}\,{\color[rgb]{.804,.4,0}\hat{\mathcal{U}}}_{1}$
in this expression using standard 3-vector rotation matrices [47, Eq.
(3.10.3)]. In fact, the general result of Eq. (121) from Siminovitch  _et al._
[48] allows us to readily evaluate the similarity transform as

$\displaystyle{\color[rgb]{.804,.4,0}\hat{\mathcal{U}}}_{1}^{\dagger}\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}X}}\,{\color[rgb]{.804,.4,0}\hat{\mathcal{U}}}_{1}$
$\displaystyle=-\frac{\sin\left(\sqrt{{\color[rgb]{0.1216,0.2,0.436}S}^{2}+{\color[rgb]{0,.545,.545}C}^{2}}\right)}{\sqrt{{\color[rgb]{0.1216,0.2,0.436}S}^{2}+{\color[rgb]{0,.545,.545}C}^{2}}}\,{\color[rgb]{0.1216,0.2,0.436}S}\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}X}}-\frac{\sin\left(\sqrt{{\color[rgb]{0.1216,0.2,0.436}S}^{2}+{\color[rgb]{0,.545,.545}C}^{2}}\right)}{\sqrt{{\color[rgb]{0.1216,0.2,0.436}S}^{2}+{\color[rgb]{0,.545,.545}C}^{2}}}\,{\color[rgb]{0,.545,.545}C}\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}Y}}+\cos\left(\sqrt{{\color[rgb]{0.1216,0.2,0.436}S}^{2}+{\color[rgb]{0,.545,.545}C}^{2}}\right)\,{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}Z}}$
(14)

We can then trivially read-off the expectation value
$\braket{{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}X}}}({\color[rgb]{.31,.58,.804}T})$
when our initial state is an eigenstate of each of the spin operators. In
particular, when starting in eigenstate
$\ket{-{\color[rgb]{0.327,0.327,0.3072}Z}}$ of
${\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}Z}}$,
the expectation value is

$\displaystyle\braket{{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}X}}}({\color[rgb]{.31,.58,.804}T})/{\color[rgb]{.545,0,0}\hbar}=\frac{\sin\left(\sqrt{{\color[rgb]{0.1216,0.2,0.436}S}^{2}+{\color[rgb]{0,.545,.545}C}^{2}}\right)}{\sqrt{{\color[rgb]{0.1216,0.2,0.436}S}^{2}+{\color[rgb]{0,.545,.545}C}^{2}}}\,{\color[rgb]{0.1216,0.2,0.436}S}.$
(15)

That is (continuing with the assumption that
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$
and hence
$\sqrt{{\color[rgb]{0.1216,0.2,0.436}S}^{2}+{\color[rgb]{0,.545,.545}C}^{2}}$
are small), when starting in $\ket{-{\color[rgb]{0.327,0.327,0.3072}Z}}$, the
expectation value is linearly sensitive to the Fourier sine coefficient of
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$
that we wish to measure. Therefore we choose to start in
$\ket{-{\color[rgb]{0.327,0.327,0.3072}Z}}$. Note that, although not used in
our protocol, when starting in eigenstate
$\ket{+{\color[rgb]{0.327,0.327,0.3072}Y}}$ of
${\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}Y}}$
we are similarly linearly sensitive to the cosine coefficient of
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$
at frequency
${\color[rgb]{.536,.408,.804}f}={\color[rgb]{0.327,0.2808,0.2376}\Omega}/{\color[rgb]{.545,.396,.03}2\pi}$,
and when starting in $\ket{+{\color[rgb]{0.327,0.327,0.3072}X}}$ of
${\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}X}}$
we are similarly linearly sensitive to the spectral power (and quadratically
sensitive to the magnitude of the complex Fourier component) of
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$
at frequency
${\color[rgb]{.536,.408,.804}f}={\color[rgb]{0.327,0.2808,0.2376}\Omega}/{\color[rgb]{.545,.396,.03}2\pi}$.

At this point the analysis of the measurement protocol is entirely in the
second rotating frame. First note that a state that starts in an eigenstate of
${\color[rgb]{.365,.28,.545}\hat{F}}_{{\color[rgb]{0.327,0.327,0.3072}Z}}$ in
the lab frame will also start in an eigenstate of
${\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}Z}}$
in the second rotating frame. Thus we prepare our atoms in the lab frame
$\ket{-{\color[rgb]{0.327,0.327,0.3072}Z}}$ state, which is conveniently the
state that the magnetic catch cooling stage polarizes our atoms into. We then
let these atoms evolve such that they have an expected
$\braket{{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}X}}}({\color[rgb]{.31,.58,.804}T})$
linearly sensitive to the Fourier sine coefficient. Note that
$\braket{{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}X}}}({\color[rgb]{.31,.58,.804}T})$
in the second rotating frame is the same as
$\braket{{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}X}}}({\color[rgb]{.31,.58,.804}T})$
in the first rotating frame. We finish the pulse sequence by applying a $20\
\text{kHz}$ $\pi/2$ read-out pulse around
${\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{\color[rgb]{0.327,0.327,0.3072}Y}$
to rotate this expected value from
${\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{\color[rgb]{0.327,0.327,0.3072}X}$
to ${\color[rgb]{.365,.28,.545}\hat{F}}_{\color[rgb]{0.327,0.327,0.3072}Z}$.
The populations of the eigenstates of
${\color[rgb]{.365,.28,.545}\hat{F}}_{\color[rgb]{0.327,0.327,0.3072}Z}$ are
then imaged using a Stern-Gerlach projective measurement. Last, we record our
sine coefficients as

$\displaystyle{\color[rgb]{.545,.545,0}y}({\color[rgb]{.536,.408,.804}f})$
$\displaystyle=\frac{\braket{{\color[rgb]{.365,.28,.545}\hat{F}}^{{\color[rgb]{0.327,0.327,0.3072}\mathrm{R}}}_{{\color[rgb]{0.327,0.327,0.3072}X}}}({\color[rgb]{.31,.58,.804}T})}{{\color[rgb]{.545,0,0}\hbar}\,{\color[rgb]{.545,.545,0}\gamma}\,{\color[rgb]{.31,.58,.804}T}}$
(16)
$\displaystyle=\frac{\braket{{\color[rgb]{.365,.28,.545}\hat{F}}_{\color[rgb]{0.327,0.327,0.3072}Z}}_{\frac{\pi}{2}}}{{\color[rgb]{.545,0,0}\hbar}\,{\color[rgb]{.545,.545,0}\gamma}\,{\color[rgb]{.31,.58,.804}T}},$
(17)
$\displaystyle=\frac{1}{{\color[rgb]{.545,.545,0}\gamma}\,{\color[rgb]{.31,.58,.804}T}}\frac{{\color[rgb]{0.436,0.36,0.2656}N}_{+}-{\color[rgb]{0.436,0.36,0.2656}N}_{-}}{{\color[rgb]{0.436,0.36,0.2656}N}_{+}+{\color[rgb]{0.436,0.36,0.2656}N}_{0}+{\color[rgb]{0.436,0.36,0.2656}N}_{-}}.$
(18)

##  II Tuning of regularization parameter

The LASSO problem needs to be tuned by the regularization parameter
${\color[rgb]{.604,.804,.196}\lambda}$ [14, Sec. 8.2.2], based on the typical
noise and sparsity of measurements. To do this, we used mathematical
optimization to find a ${\color[rgb]{.604,.804,.196}\lambda}$ that produced
the smallest average error over a set of one-thousand measurements of signals.
Since collecting such data manually would require on the order of a month of
time, we instead used our open-source GPU accelerated simulation python
package _spinsim_ [43] to generate a set of simulated experiment results. Each
simulated sequence was randomly assigned 0, 1 or 2 neural-like pulses of
amplitude $1\ \mathrm{kHz}/{\color[rgb]{.545,.545,0}\gamma}$ (i.e., $143\
\mathrm{nT}$), each randomly placed in time. Imperfections were modeled in
these experiments in the form of shot-to-shot bias drift (standard deviation
of $200\ \mathrm{Hz}/{\color[rgb]{.545,.545,0}\gamma}$ consistent with lab
measurements) and Poisson-distributed atom shot noise (expected number of
atoms depreciated from order million to one thousand to also take into account
photon shot noise). The error metric we chose was the average $\ell_{1}$-norm
of the difference between true and measured signals. This is the middle-ground
between the $\ell_{0}$-norm,
$\|{\color[rgb]{0.327,0.327,0.3072}x}\|_{0}\stackrel{{\scriptstyle{\color[rgb]{0.436,0.38,0.296}\text{def}}}}{{=}}\\#\\{{\color[rgb]{0.327,0.27,0.1992}k}\text{
such that
}{\color[rgb]{0.327,0.327,0.3072}x}_{\color[rgb]{0.327,0.27,0.1992}k}=0\\}$
[14, Sec. 2.1] which favors correct supports, and the $\ell_{2}$-norm which
favors correct amplitudes. Note that we decided _against_ using area under the
receiver operating characteristic (ROC)curve (AUC)(see below) as a metric for
training since its landscape is very flat and therefore not informative for
training. We scanned ${\color[rgb]{.604,.804,.196}\lambda}$ over 200 different
values with log spacing from $100\
\mathrm{mHz}/{\color[rgb]{.545,.545,0}\gamma}$ to $10\
\mathrm{Hz}/{\color[rgb]{.545,.545,0}\gamma}$, and found that value that
returned the smallest error was ${\color[rgb]{.604,.804,.196}\lambda}=1.04\
\mathrm{Hz}/{\color[rgb]{.545,.545,0}\gamma}$. We emphasize that tuning
${\color[rgb]{.604,.804,.196}\lambda}$ this precisely is not necessary, and
using ${\color[rgb]{.604,.804,.196}\lambda}$ in the range increased the
$\ell_{1}$ norm by less than $10\%$ and yielded quality recoveries.

##  III A signal detection statistic: the area under the receiver operating
characteristic (ROC) curve (AUC)

There are three kinds of signals of interest when analyzing acquired
synthesized neural-like pulses: We have the ground-truth signal
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$,
which is a recording of the electrical signal commanded to the magnetic coils.
We have the recovered signal
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\star}}$,
which is the result of the compressive reconstructions or the Ramsey
measurements. Finally, we have the pulse template
${\color[rgb]{.804,.59,.804}p}$, which is the expected shape of each
individual pulse, but it starts at ${\color[rgb]{.31,.58,.804}t}=0$.

The highest fidelity linear method of detecting pulse locations within a
signal is by using a matched filter [49, Sec. 4.3.1]. In a matched filter, the
template ${\color[rgb]{.804,.59,.804}p}$ is translated to all the possible
starting locations in the signal, and an inner-product is taken between the
two via
${\color[rgb]{.804,.4,0}\mu}^{{\color[rgb]{.804,.4,0}\natural}}_{\color[rgb]{0.327,0.27,0.1992}k}=\sum_{\color[rgb]{.365,.28,.545}j}{\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}_{\color[rgb]{.365,.28,.545}j}\,{\color[rgb]{.804,.59,.804}p}_{{\color[rgb]{.365,.28,.545}j}+{\color[rgb]{0.327,0.27,0.1992}k}}$.
This can be written succinctly using the cross-correlation operator via
${\color[rgb]{.804,.4,0}\mu}^{{\color[rgb]{.804,.4,0}\natural}}={\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}\circledast{\color[rgb]{.804,.59,.804}p}$,
${\color[rgb]{.804,.4,0}\mu}^{{\color[rgb]{.804,.4,0}\star}}={\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\star}}\circledast{\color[rgb]{.804,.59,.804}p}$.
A high value of
${\color[rgb]{.804,.4,0}\mu}^{{\color[rgb]{.804,.4,0}\natural}}_{\color[rgb]{0.327,0.27,0.1992}k}$
or
${\color[rgb]{.804,.4,0}\mu}^{{\color[rgb]{.804,.4,0}\star}}_{\color[rgb]{0.327,0.27,0.1992}k}$
at any index corresponds to a high likelihood that a pulse is located in
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\natural}}$
or
${\color[rgb]{0.327,0.327,0.3072}x}^{{\color[rgb]{0.327,0.327,0.3072}\star}}$
respectively at time index ${\color[rgb]{0.327,0.27,0.1992}k}$. We can
classify locations as having a pulse or not by thresholding
${\color[rgb]{.804,.4,0}\mu}^{{\color[rgb]{.804,.4,0}\natural}}$ and
${\color[rgb]{.804,.4,0}\mu}^{{\color[rgb]{.804,.4,0}\star}}$. We obtain a
ground truth classification
${\color[rgb]{0,.545,.545}c}^{{\color[rgb]{0,.545,.545}\natural}}={\color[rgb]{.31,.58,.804}\mathcal{T}}_{\|{\color[rgb]{.804,.59,.804}p}\|_{2}^{2}/2}({\color[rgb]{.804,.4,0}\mu}^{{\color[rgb]{.804,.4,0}\natural}})$
by taking a threshold of the ground truth signal at a level of
$\|{\color[rgb]{.804,.59,.804}p}\|_{2}^{2}/2$. Here we define
${\color[rgb]{.31,.58,.804}\mathcal{T}}_{{\color[rgb]{0.186,0.348,0.4824}\Gamma}}({\color[rgb]{.804,.4,0}\mu}^{{\color[rgb]{.804,.4,0}\natural}})_{\color[rgb]{0.327,0.27,0.1992}k}=0$
for
${\color[rgb]{.804,.4,0}\mu}^{{\color[rgb]{.804,.4,0}\natural}}_{\color[rgb]{0.327,0.27,0.1992}k}<{\color[rgb]{0.186,0.348,0.4824}\Gamma}$,
and
${\color[rgb]{.31,.58,.804}\mathcal{T}}_{{\color[rgb]{0.186,0.348,0.4824}\Gamma}}({\color[rgb]{.804,.4,0}\mu}^{{\color[rgb]{.804,.4,0}\natural}})_{\color[rgb]{0.327,0.27,0.1992}k}=1$
otherwise.

Classifications can be made of the recovered signal
${\color[rgb]{0,.545,.545}c}^{{\color[rgb]{0,.545,.545}\star}}={\color[rgb]{.31,.58,.804}\mathcal{T}}_{{\color[rgb]{0.186,0.348,0.4824}\Gamma}}({\color[rgb]{.804,.4,0}\mu}^{{\color[rgb]{.804,.4,0}\star}})$
at various thresholds ${\color[rgb]{0.186,0.348,0.4824}\Gamma}$. For each
${\color[rgb]{0.186,0.348,0.4824}\Gamma}$, we can obtain the statistics of a
_confusion matrix_ Fawcett [45, Fig. 1] by comparing the classifications of
individual locations in
${\color[rgb]{0,.545,.545}c}^{{\color[rgb]{0,.545,.545}\star}}$ to those in
${\color[rgb]{0,.545,.545}c}^{{\color[rgb]{0,.545,.545}\natural}}$,

$\displaystyle\text{\\# True
Positives}=\mathrm{{\color[rgb]{.31,.58,.804}T}{\color[rgb]{.804,.59,.804}P}}=$
$\displaystyle\sum_{\color[rgb]{0.327,0.27,0.1992}k}{\color[rgb]{0,.545,.545}c}^{{\color[rgb]{0,.545,.545}\star}}_{\color[rgb]{0.327,0.27,0.1992}k}{\color[rgb]{0,.545,.545}c}^{{\color[rgb]{0,.545,.545}\natural}}_{\color[rgb]{0.327,0.27,0.1992}k},$
(19) $\displaystyle\text{\\# False
Positives}=\mathrm{{\color[rgb]{.365,.28,.545}F}{\color[rgb]{.804,.59,.804}P}}=$
$\displaystyle\sum_{\color[rgb]{0.327,0.27,0.1992}k}{\color[rgb]{0,.545,.545}c}^{{\color[rgb]{0,.545,.545}\star}}_{\color[rgb]{0.327,0.27,0.1992}k}(1-{\color[rgb]{0,.545,.545}c}^{{\color[rgb]{0,.545,.545}\natural}}_{\color[rgb]{0.327,0.27,0.1992}k}),$
(20) $\displaystyle\text{\\# False
Negatives}=\mathrm{{\color[rgb]{.365,.28,.545}F}{\color[rgb]{0.436,0.36,0.2656}N}}=$
$\displaystyle\sum_{\color[rgb]{0.327,0.27,0.1992}k}(1-{\color[rgb]{0,.545,.545}c}^{{\color[rgb]{0,.545,.545}\star}}_{\color[rgb]{0.327,0.27,0.1992}k}){\color[rgb]{0,.545,.545}c}^{{\color[rgb]{0,.545,.545}\natural}}_{\color[rgb]{0.327,0.27,0.1992}k},\quad\text{and}$
(21) $\displaystyle\text{\\# True
Negatives}=\mathrm{{\color[rgb]{.31,.58,.804}T}{\color[rgb]{0.436,0.36,0.2656}N}}=$
$\displaystyle\sum_{\color[rgb]{0.327,0.27,0.1992}k}(1-{\color[rgb]{0,.545,.545}c}^{{\color[rgb]{0,.545,.545}\star}}_{\color[rgb]{0.327,0.27,0.1992}k})(1-{\color[rgb]{0,.545,.545}c}^{{\color[rgb]{0,.545,.545}\natural}}_{\color[rgb]{0.327,0.27,0.1992}k}).$
(22)

Two important statistics [45, Sec. 2] can be derived from these values. The
recall (which we write ${\color[rgb]{.545,0,0}\mathrm{rc}}$, also sensitivity,
true-positive rate or hit rate in the literature) is defined as

$\displaystyle{\color[rgb]{.545,0,0}\mathrm{rc}}\stackrel{{\scriptstyle{\color[rgb]{0.436,0.38,0.296}\text{def}}}}{{=}}$
$\displaystyle\frac{\mathrm{{\color[rgb]{.31,.58,.804}T}{\color[rgb]{.804,.59,.804}P}}}{\mathrm{{\color[rgb]{.31,.58,.804}T}{\color[rgb]{.804,.59,.804}P}}+\mathrm{{\color[rgb]{.365,.28,.545}F}{\color[rgb]{0.436,0.36,0.2656}N}}},$
(23)

and is the fraction of ground-truth positives correctly detected in the
recovery.

Similarly, fallout (which we write ${\color[rgb]{.48,.215,.545}\mathrm{fo}}$,
also false-positive rate, false-alarm rate or $1-\text{specificity}$ in the
literature) is defined as

$\displaystyle{\color[rgb]{.48,.215,.545}\mathrm{fo}}\stackrel{{\scriptstyle{\color[rgb]{0.436,0.38,0.296}\text{def}}}}{{=}}$
$\displaystyle\frac{\mathrm{{\color[rgb]{.365,.28,.545}F}{\color[rgb]{.804,.59,.804}P}}}{\mathrm{{\color[rgb]{.365,.28,.545}F}{\color[rgb]{.804,.59,.804}P}}+\mathrm{{\color[rgb]{.31,.58,.804}T}{\color[rgb]{0.436,0.36,0.2656}N}}},$
(24)

and is the fraction of ground-truth negatives incorrectly detected in the
recovery. Note that these are both functions of the recovery threshold
${\color[rgb]{0.186,0.348,0.4824}\Gamma}$. The parametric curve of points
$\left({\color[rgb]{.48,.215,.545}\mathrm{fo}}({\color[rgb]{0.186,0.348,0.4824}\Gamma}),{\color[rgb]{.545,0,0}\mathrm{rc}}({\color[rgb]{0.186,0.348,0.4824}\Gamma})\right)$
that ${\color[rgb]{0.186,0.348,0.4824}\Gamma}$ traces out is called the
receiver operating characteristic (ROC)of the sensor [45, Sec. 4]. When
${\color[rgb]{0.186,0.348,0.4824}\Gamma}\to-\infty$, any noise will be
classified as a signal, meaning
$({\color[rgb]{.48,.215,.545}\mathrm{fo}}({\color[rgb]{0.186,0.348,0.4824}\Gamma}),{\color[rgb]{.545,0,0}\mathrm{rc}}({\color[rgb]{0.186,0.348,0.4824}\Gamma}))\to(1,1)$.
When ${\color[rgb]{0.186,0.348,0.4824}\Gamma}\to\infty$, none of the true
signal responses will be large enough to be classified as a signal, so
$({\color[rgb]{.48,.215,.545}\mathrm{fo}}({\color[rgb]{0.186,0.348,0.4824}\Gamma}),{\color[rgb]{.545,0,0}\mathrm{rc}}({\color[rgb]{0.186,0.348,0.4824}\Gamma}))\to(0,0)$.
If there is some threshold
${\color[rgb]{0.186,0.348,0.4824}\Gamma}^{{\color[rgb]{0.186,0.348,0.4824}\star}}$
that can completely accurately distinguish between pulses and noise, then
$({\color[rgb]{.48,.215,.545}\mathrm{fo}}({\color[rgb]{0.186,0.348,0.4824}\Gamma}),{\color[rgb]{.545,0,0}\mathrm{rc}}({\color[rgb]{0.186,0.348,0.4824}\Gamma}))=(0,1)$.

The closer the ROC curve gets to this value of $(0,1)$, the better the
classifier. This can be quantified by integrating the area under the receiver
operating characteristic (ROC)curve (AUC) [45, Sec. 7],
${\color[rgb]{.804,.608,.112}\mathrm{auc}}$. If the ROC curve crosses $(0,1)$,
then the area is a square and ${\color[rgb]{.804,.608,.112}\mathrm{auc}}=1$.
If the classifier does no better than randomly guessing, then the area is a
triangle under
${\color[rgb]{.48,.215,.545}\mathrm{fo}}={\color[rgb]{.545,0,0}\mathrm{rc}}$
and ${\color[rgb]{.804,.608,.112}\mathrm{auc}}=1/2$.

Figure 3 in the Manuscript uses AUC as a quality metric of compressive
recoveries. The insets of the figure show examples of ROC s for recoveries
made using different compression ratios, averaged for 200 different
subsamples. We found that, for our case of using data from only one signal,
the number of subsamples where the average AUC drops below $100\%$ matches
that where we would expect failed recoveries in compressive sensing
theoretical guarantees [15, Eq. (9.24)].

##  IV Genuinely undersampled compressive measurements

Figure 4:  Recoveries of two genuinely undersampled signals, each with only of
the 60 of the complete set of 99 sine coefficients measured by the cold atoms.
(a) a single pulse signal recovery. (b) a two-pulse signal recovery

While using the atoms to sample all sine coefficients is necessary to analyse
the effectiveness of quantum compressive sensing, recording a complete dataset
misses the spirit of the kind of sensor we demonstrate. In this spirit, we
show compressive recovery of signals from incomplete measurements of sine
coefficients.

We chose 60 random frequencies out of the complete DST (i.e., multiples of
$100\ \mathrm{Hz}$). We used the atoms to measure sine coefficients of a
synthesized waveform at these 60 frequencies, and only these frequencies.
FISTA then recovered the waveform based on this incomplete information.

Figure 4 shows compressive waveform estimation from genuinely incomplete data
sets. It is evident on inspection that, as per the signal detection analysis
in the previous Sec., there is a range of thresholds that will achieve correct
detection of pulses.