# Efficient measurement of the time-dependent cavity field through compressed
sensing

Fang Zhao  Center for Quantum Technology Research and Key Laboratory of
Advanced Optoelectronic Quantum Architecture and Measurements (MOE), School of
Physics, Beijing Institute of Technology, Beijing 100081, China  China Academy
of Engineering Physics, Beijing 100088, People’s Republic of China  Qing Zhao
Center for Quantum Technology Research and Key Laboratory of Advanced
Optoelectronic Quantum Architecture and Measurements (MOE), School of Physics,
Beijing Institute of Technology, Beijing 100081, China  Dazhi Xu
<EMAIL_ADDRESS>Center for Quantum Technology Research and Key Laboratory of
Advanced Optoelectronic Quantum Architecture and Measurements (MOE), School of
Physics, Beijing Institute of Technology, Beijing 100081, China

###### Abstract

We propose a method based on compressed sensing (CS) to measure the evolution
processes of the states of a driven cavity quantum electrodynamics system. In
precisely reconstructing the coherent cavity field amplitudes, we have to
prepare the same states repetitively and each time perform one measurement
with short sampling intervals considering the quantum nature of measurement
and the Nyquist–Shannon sampling theorem. However, with the help of CS, the
number of measurements can be exponentially reduced without loss of the
recovery accuracy. We use largely detuned atoms and control their interactions
with the cavity field to modulate coherent state amplitudes according to the
scheme encoded in the sensing matrix. The simulation results show that the CS
method efficiently recovers the amplitudes of the coherent cavity field even
in the presence of noise.

##  I Introduction

A measurement usually has destructive consequences for the quantum state,
which prevents repetitive probing of the quantum system. Even though the
quantum non-demolition measurement allows the measurement of an observable
without destroying the system (Braginsky  _et al._ , 1980), this method has
limitations. It generally fails in the case of frequent measurements owing to
the quantum Zeno effect (Misra and Sudarshan, 1977; Itano  _et al._ , 1990). A
typical method of recording the time evolution of a quantum system is to
prepare an ensemble of identical quantum systems and perform only one
measurement for each copy at a different sampling instant. However, this
approach requires enormous discrete samplings to reconstruct the original
continuous-time dynamic signals. The minimal sampling rate of this approach is
given by the Nyquist–Shannon sampling theorem (Nyquist, 1928; Shannon, 1949).

Fortunately, the theory of compressed sensing (CS) developed by mathematicians
nearly two decades ago (Candès  _et al._ , 2006; Candès and Tao, 2006; Donoho,
2006) has proven to be able to exactly recover sparse signals with only a
small set of linear random measurements. As sparsity is a general property of
almost all meaningful signals, CS techniques have been widely applied in
various areas, such as medical imaging (Lustig  _et al._ , 2007), radar
(Potter  _et al._ , 2010; Liu  _et al._ , 2019), and geophysics (Lin and
Herrmann, 2007). Recently, accompanying the fast development of quantum
sensing, many attempts have been made to use CS techniques to improve the
speed of quantum measurements. Most of these quantum CS schemes concern time-
independent signals; e.g., quantum imaging (Katz  _et al._ , 2009; Bromberg
_et al._ , 2009; Zerom  _et al._ , 2011), quantum tomography (Gross  _et al._
, 2010; Cramer  _et al._ , 2010; Shabani  _et al._ , 2011; Rudinger and Joynt,
2015; Riofrío  _et al._ , 2017), and two-dimensional spectroscopy (Dunbar  _et
al._ , 2013; Roeding  _et al._ , 2017).

To measure time-dependent signals in a compressive way, one needs to control
the dynamics of the quantum system effectively. In an interesting work
(Magesan  _et al._ , 2013), the authors proposed a method of compressively
measuring the time-varying magnetic field via a single spin, where a sequence
of $\pi$-pulses flips the spin according to the requirements of the CS theory.
However, the spin in Ref. (Magesan  _et al._ , 2013) is along the magnetic
field; thus, the Hamiltonians at different times commute with each other
($\left[H(t),H(t^{\prime})\right]=0$), which only covers a particular type of
quantum dynamics. On the basis of a cavity quantum electrodynamics (QED)
system, we show that the idea of the compressive measurement of the time-
dependent signals can be generalized to a quantum system with the property
$\left[H(t),H(t^{\prime})\right]\neq 0$ ($t\neq t^{\prime}$) in this work.

The cavity QED system is an essential platform for quantum computation (Ladd
_et al._ , 2010) and quantum sensing (Malnou  _et al._ , 2019; Lewis-Swan  _et
al._ , 2020). We employ CS techniques to continuously measure the state of the
cavity field driven by a time-varying classical source. Our purpose is to
reconstruct the coherent cavity amplitude with high accuracy while
exponentially reducing the number of measurements. To this end, we send
sequences of largely detuned atoms through the cavity at randomly chosen
instants during the evolution processes, such that the cavity field can be
controlled according to the assigned sensing matrix. Each measurement is
performed at the end of one complete time evolution process. With only a small
set of under-sampled measurement results, the convex optimization algorithm
can successfully recover the coherent amplitude even in the presence of white
noise. Our approach expands the application range of CS techniques in the
measurements of quantum systems.

The paper is arranged as follows. Section II revisits the theory of CS.
Section III briefly introduces the driven cavity QED system and the
Nyquist–Shannon sampling method for the cavity field. Section IV shows how to
apply CS techniques to the cavity QED system. Section V presents the recovery
results and the error analysis. Section VI summarizes our results and
discusses the possible future work.

##  II Compressed Sensing

The Nyquist–Shannon theorem states that to exactly recover a continuous-time
signal, the sampling frequency should be at least twice the signal’s highest
frequency (Nyquist, 1928; Shannon, 1949). If the signal has a sparse
representation in a set of known bases, the CS theory allows a much lower
sampling frequency than that required by the Nyquist–Shannon theorem. Before
applying the CS techniques to the cavity QED system, we briefly introduce the
basic knowledge of CS theory in this section. In the linear measurement
models, the signal to be measured can be modeled using an $N$-dimensional
vector $\vec{x}\in\mathbb{R}^{N}$, and the measurement is formulated as an
$M\times N$ matrix $\mathcal{A}$ acting on $\vec{x}$:

$\vec{y}=\mathcal{A}\cdot\vec{x},$ (1)

where the $M$-dimensional vector $\vec{y}\in\mathbb{R}^{M}$ is the measurement
results and $M$ is the number of measurements. If $M=N$ and $\mathcal{A}$ is a
full-rank matrix, the signal can be uniquely and exactly reconstructed
according to $\vec{x}=\mathcal{A}^{-1}\cdot\vec{y}$. However, we suppose the
number of measurements $M$ is smaller than the dimension $N$ of the signal. In
that case, there are infinite numbers of $\vec{x}$ satisfying Eq. (1), and it
is thus difficult to tell which one is the actual signal without further
information. Fortunately, real-world signals are generally structured, and
sparsity is one of the essential prior characteristics of the signals. A
signal is called $S$-sparse if it has no more than $S$ ($S\ll N$) non-zero
components. The CS technique takes advantage of sparsity to recover the signal
from an appreciable reduced number of samplings.

Mathematically, theory says if two conditions are satisfied, (1) $\vec{x}$ is
$S$-sparse and (2) $\mathcal{A}$ satisfies the restricted isometry property
(RIP) of order $2S$, then the unique $S$-sparse $\vec{x}$ can be exactly
recovered by solving the convex optimization problem:

$\vec{x}=\arg\min\|\vec{x}^{\prime}\|_{1}\leavevmode\nobreak\ \text{suject
to}\leavevmode\nobreak\ \vec{y}=\mathcal{A}\cdot\vec{x}^{\prime},$ (2)

where $\vec{x}^{\prime}$ is the recovered signal and
$\|\vec{x}^{\prime}\|_{1}\coloneqq\sum_{i}|x_{i}^{\prime}|$ is the
$l_{1}$-norm of vector $\vec{x}^{\prime}$. The RIP of order $S$ requires
$\mathcal{A}$ satisfying

$(1-\delta_{S})\|\vec{x}\|_{2}^{2}\leq\|\mathcal{A}\cdot\vec{x}\|_{2}^{2}\leq(1+\delta_{S})\|\vec{x}\|_{2}^{2}$
(3)

for all $S$-sparse vectors $\vec{x}$. Here, $\|\vec{x}\|_{2}^{2}$ is the
$l_{2}$-norm or Euclidean distance of $\vec{x}$, and $\delta_{S}\in(0,1)$ is
an $S$-dependent number (Candès  _et al._ , 2006; Baraniuk  _et al._ , 2008).
The RIP of order $2S$ is a sufficient condition of the sensing matrix
$\mathcal{A}$ for successfully recovering the unique $S$-sparse signal, which
can be intuitively understood as a restriction on $\mathcal{A}$ such that it
approximately preserves the distance between any two $S$-sparse vectors. A
smaller $\delta_{S}$ corresponds to $\mathcal{A}$ being more similar to an
isometry transformation. It is usually difficult to construct a matrix that
strictly satisfies the RIP. However, it has been proved that a random matrix
can satisfy the RIP with very high probability (Candès  _et al._ , 2006),
which makes the application of CS techniques extremely convenient.

There are various recovery algorithms suitable for solving the convex
optimization problem of Eq. (2). We use the orthogonal matching pursuit (OMP)
algorithm (Pati  _et al._ , 1993), which is accurate and efficient for the
goals of this work. If the signal is not sparse in the bases in which we
perform the measurements, we first find an orthogonal transformation $\Phi$
such that $\vec{x}_{\Phi}=\Phi\cdot\vec{x}$ is sparse. The transformation
$\Phi$ is generally known for a class of signals in prior. For instance, the
images are usually sparse after wavelet transformation. As a result, we can
reformulate the constraint in the optimization problem of Eq. (2) as
$\vec{y}=\mathcal{A}_{\Phi}\cdot\vec{x}_{\Phi}^{\prime}$, where the
transformed sensing matrix $\mathcal{A}_{\Phi}=\mathcal{A}\cdot\Phi^{-1}$ is a
random matrix satisfying the RIP as long as $\mathcal{A}$ is random.

Notwithstanding that the goal of CS is to reduce the number of linear
measurements $M$ as much as possible, it is evident that $M$ cannot be made
arbitrarily small to extract sufficient information of the signal. If the
sensing matrix $\mathcal{A}$ satisfies the RIP of order $2S$ with
$\delta_{2S}\in(0,\frac{1}{2}]$, then the lower bound of $M$ is given by

$M\geq C\cdot S\log_{2}(\frac{N}{S}),$ (4)

where $C$ is a constant (Baraniuk  _et al._ , 2008). This theorem shows that
CS can reduce the number of linear measurements from the order of
$\mathcal{O}(N)$ to the order of $\mathcal{O}(\log N)$.

Furthermore, as the idea of CS is to extract the information of only the most
significant elements of the signal, CS techniques naturally resist low noise
in the signal. In the presence of low noise with strength $\epsilon$ in the
signal, we can generalize the convex optimization problem of Eq. (2) as

$\vec{x}=\arg\min\|\vec{x}^{\prime}\|_{1}\leavevmode\nobreak\ \text{suject
to}\leavevmode\nobreak\
\|\vec{y}-\mathcal{A}\cdot\vec{x}^{\prime}\|_{2}^{2}\leq\epsilon.$ (5)

The noise resistance property is another advantage of the CS techniques for
our purpose of reconstructing the time evolution of the cavity field.

##  III Nyquist–Shannon sampling of the driven cavity field

Techniques used for the instantaneous measurement of the coherent field or the
average photon number inside the microwave cavity have been well developed;
e.g., homodyne detection (Lvovsky and Raymer, 2009) and quantum non-demolition
measurement (Brune  _et al._ , 1992; Nogues  _et al._ , 1999). Continuous
probing of the varying cavity field not only provides more information about
the dynamics of the quantum system but also is helpful to the quantum
precision measurement of external parameters. However, an efficient method of
continuously interrogating the cavity QED system is still lacking owing to the
fragility of the quantum system. To this end, we adopt CS techniques to
accelerate the sampling processes. For preparation, this section introduces
the driven cavity system with the usual continuous measurement procedure
without CS techniques.

The system that we are interested in is a microwave cavity driven by a
classical time-dependent radio-frequency field, the Hamiltonian of which in
the rotating frame is

$H(t)=if(t)(e^{-i\Delta t}a^{\dagger}-e^{i\Delta t}a).$ (6)

Here, $a(a^{\dagger})$ is the annihilation (creation) operator of the cavity
field, and $\Delta=\omega_{\text{d}}-\omega_{0}$ is the detuning between the
driving field frequency $\omega_{\text{d}}$ and the cavity mode $\omega_{0}$.
The amplitude of the driving field $f(t)$ is assumed to be a real function of
time. The evolution operator corresponding to $H(t)$ is given by

$U(t,t_{0})=e^{i\phi(t,t_{0})}D[\alpha(t,t_{0})],$ (7)

where $\exp[i\phi(t,t_{0})]$ is a global phase factor with

$\phi(t,t_{0})=\int_{t_{0}}^{t}ds\int_{t_{0}}^{s}ds^{\prime}f(s)f(s^{\prime})\sin\left[\Delta(s-s^{\prime})\right]$
(8)

and $D[\alpha]\equiv\exp[\alpha a^{\dagger}-\alpha^{\ast}a]$ is the
displacement operator with the time-dependent parameter

$\alpha(t_{2},t_{1})=\int_{t_{1}}^{t_{2}}dsf(s)e^{-i\Delta s}.$ (9)

It is straightforward to find that $\alpha(t_{2},t_{1})$ satisfies the
relation $\alpha(t_{3},t_{1})=\alpha(t_{3},t_{2})+\alpha(t_{2},t_{1})$ with
$t_{3}\geq t_{2}\geq t_{1}$. The cavity is assumed to be initially prepared in
the vacuum state $|\psi(t_{0})\rangle=|0\rangle$; the average photon number
$\bar{n}(t)=\langle\psi(t)|a^{\dagger}a|\psi(t)\rangle=|\alpha(t,t_{0})|^{2}$
is then a function of time, where
$|\psi(t)\rangle=U(t,t_{0})|\psi(t_{0})\rangle$.

The coherent state amplitude $\alpha$ of a steady field can be obtained from a
balanced homodyne measurement (Hansen  _et al._ , 2001; Lvovsky and Raymer,
2009). However, the continuous measurement of the amplitude of a varying
photon field is not easy because the measurements inevitably perturb the
quantum state and alter its photon statistics even with the help of the
quantum non-demolition measurement (Ueda  _et al._ , 1992). Moreover, the
repeated measurement of the cavity field can freeze the evolution of the
cavity state owing to the quantum Zeno effect (Bernu  _et al._ , 2008; Xu  _et
al._ , 2011).

Therefore, it is impossible to make all the measurements within one unitary
evolution in recovering the continuous signal $\alpha(t,t_{0})$ within an
interval $t\in[t_{0},t_{N}]$. Instead, one straightforward approach is to
prepare $N$ identical initial states repetitively and each time let the system
evolve according to Eq. (7) until a different moment $t_{n}\equiv
n\tau_{\text{B}}$ $(n=1,2,\dots,N)$ when the measurement is carried out. Here,
$\tau_{\text{B}}$ is the uniform sampling interval, which is chosen as
$\tau_{\text{B}}^{-1}\geq 2B$ according to the Nyquist–Shannon sampling
theorem. The sampling sequence is illustrated in Fig. 1(a). The total number
of measurements is then $N=t_{N}/\tau_{\text{B}}$. The resource consumed in
the experiments linearly grows with the number of measurements $N$, which is
proportional to the bandwidth of the signal.

In the language of signal processing, the time-domain discretized coherent
state amplitude $\alpha(t,t_{0})$ is represented by an $N$-dimensional column
vector $\vec{\alpha}=(\alpha_{1},\dots,\alpha_{N})^{\text{T}}$ with components
$\alpha_{n}\equiv\alpha(t_{n},t_{0})$, which is the signal that we aim to
recover. The above Nyquist–Shannon sampling procedure is a trivial version of
linear measurement described by Eq. (1), where $\vec{\alpha}$ corresponds to
the vector $\vec{x}$, $\vec{y}$ is an $N$-dimensional column vector, each
component of which is measured in one experiment. We therefore have
$y_{n}=\alpha_{n}$, and $\mathcal{A}$ is a trivial $N\times N$ unit matrix
without considering experimental error.

Figure 1: Schematic illustration of the sampling processes of the state of the
cavity field driven by a time-dependent signal: (a) the Nyquist–Shannon
sampling sequence with uniform sampling interval $\tau_{\text{B}}^{-1}=2B$;
(b) the CS sampling scheme requires $M$ samplings each at the end of one
complete evolution process, during which the phase of the field is altered
through interaction with $K$ largely detuned atoms passing through the cavity
at randomly chosen moments.

##  IV Compressed sampling of the driven cavity field

The advantage of adopting CS techniques relies on exponentially reducing the
sampling rate in reconstructing the dynamic parameters, such as
$\alpha(t,t_{0})$ and $f(t)$. Then, only $M$ ($M\ll N$) samplings are needed
to construct a low-dimensional column vector $\vec{y}$. The critical step
toward recovering $\vec{\alpha}$ is the design and realization of the sensing
matrix $\mathcal{A}$ with the experimental feasible quantum control methods.

Mathematically, each row of the sensing matrix $\vec{\mathcal{A}}^{m}$
$(m=1,2,\dots,M)$ acting on the vector $\vec{\alpha}$ yields a linear
combination of the components of the coherent amplitude $\vec{\alpha}$.
Specifically, in the cavity QED setup, such linear manipulation of the cavity
field can be realized by sending an elaborately designed sequence of largely
detuned two-level atoms through the cavity. The effective Hamiltonian
describing the interaction between the cavity field and atom is given by

$H_{\text{A}}=\xi a^{\dagger}a(|e\rangle\langle e|-|g\rangle\langle g|).$ (10)

Here, $|e\rangle$ and $|g\rangle$ are the atomic excited and ground states,
respectively. The effective coupling strength
$\xi=g^{2}/(\omega_{0}-\omega_{\text{a}})$ is defined by the Rabi frequency
$g$, the cavity mode frequency $\omega_{0}$ and the frequency splitting
$\omega_{\text{a}}$ between the atomic levels.

Each atom is assumed to interact with the cavity field for a very short time
$\tau_{\text{A}}$ before leaving the cavity. The corresponding evolution
operator is $\exp(-iH_{\text{A}}\tau_{\text{A}})$. We also assume that $\xi$
is much larger than the average strength of the driven field $f(t)$ during
$\tau_{\text{A}}$; thus, the unitary evolution $U(t+\tau_{\text{A}},t)$ can be
neglected when the atom is interacting with the cavity field. If the atom is
prepared in the initial state $|e\rangle$ and the cavity field in the coherent
state $|\alpha\rangle$ at the instant that the atom enters the cavity, the
state of the atom–cavity system becomes

$e^{-iH_{\text{A}}\tau_{\text{A}}}|\alpha\rangle\otimes|e\rangle=|\alpha
e^{-i\xi\tau_{\text{A}}}\rangle\otimes|e\rangle$ (11)

once the atom leaves the cavity. It is seen that the atom does not alter the
average photon number $\bar{n}=|\alpha|^{2}$ because the dispersive coupling
in Eq. (10) conserves the photon number. However, the cavity field acquires a
phase factor $e^{-i\xi\tau_{\text{A}}}$ associated with the atomic state.

Once the atoms are prepared in the excited state $|e\rangle$, they will not
evolve because of the dispersive coupling. As a result, one can only consider
the cavity field state, whose evolution operator is given by a projection on
$|e\rangle$ as

$R(\tau_{\text{A}})=\langle
e|e^{-iH_{\text{A}}\tau_{\text{A}}}|e\rangle=e^{-i\xi\tau_{\text{A}}a^{\dagger}a}.$
(12)

We use this phase-modulation operator $R(\tau_{\text{A}})$ as a tool with
which to manipulate the cavity field according to the arrangement set by the
sensing matrix.

### Procedures of CS measurements in the cavity QED system

All the requirements for the CS measurements have been prepared. We list the
main procedures of CS measurements of the cavity field as follows.

(i) Initially prepare the cavity field in the vacuum state $|0\rangle$.

(ii) Drive the cavity field using the classical external source $f(t)$ from
time $t_{0}$ to $t_{N}$. During the evolution, largely detuned two-level atoms
prepared in the excited state pass through the cavity field separately at $K$
randomly chosen moments $t_{1},t_{2},\dots,t_{K}$. All the $K$ moments should
be positive integer multiples of $\tau_{\text{B}}$. The atom–photon
interaction time can be tuned using the velocity of the atom and is set at
$\tau_{\text{A}}^{\ast}=\pi/\xi$, during which the free evolution of the
cavity field can be neglected.

(iii) At time $t_{N}$, the coherent field amplitude $\Lambda_{\text{f}}$ is
measured via homodyne detection, which completes one compressed sampling of a
single data point.

(iv) Return to the beginning and repeat steps (i) to (iii) $M$ times. The $M$
measured results $\Lambda_{\text{f}}^{m}$ ($m=1,2,\dots,M$ labels the $m$th
measurement) form a column vector
$\vec{\Lambda}_{\text{f}}=(\Lambda_{\text{f}}^{1},\Lambda_{\text{f}}^{2},\dots,\Lambda_{\text{f}}^{M})^{\text{T}}$,
which is the signal obtained by compressive sampling.

(v) The time-dependent cavity field amplitude is recovered by the convex
optimization algorithm with the input $\vec{\Lambda}_{\text{f}}$ and the
sensing matrix constructed in step (ii).

A schematic illustration of the CS approach is shown in Fig. 1 (b).

### Sensing matrix

We now take the $m$th measurement as an example to illustrate how to construct
the sensing matrix $\mathcal{A}$ based on the dynamic evolution of the cavity
state. The initial state of the cavity field is the vacuum state
$|\psi(t_{0})\rangle=|0\rangle$, and the $K$ injected atoms are always in the
product state $|\mathcal{E}\rangle\equiv\otimes_{j=1}^{K}|e\rangle_{j}$.
Before the first atom enters the cavity at the randomly chosen moment
$t_{1}^{m}$, the cavity field evolves into a coherent state

$|\psi^{m}(t_{1}^{m})\rangle=U(t_{1}^{m},t_{0})|\psi(t_{0})\rangle=|\alpha(t_{1}^{m},t_{0})\rangle.$
(13)

At $t_{1}^{m}$, the first atom starts interacting with the cavity field for a
short interval $\tau_{\text{A}}^{\ast}$. When the interaction ends, the sign
of the cavity field amplitude is flipped as

$R(\tau_{\text{A}}^{\ast})|\psi^{m}(t_{1}^{m})\rangle=|-\alpha(t_{1}^{m},t_{0})\rangle.$
(14)

Here, the role of the excited state atom is akin to the $\pi$-pulse that flips
the spin in the NMR system. As in Ref. (Magesan  _et al._ , 2013), the varying
magnetic field is detected using CS techniques, where the spin is flipped by
the $\pi$-pulses serving as a time-domain modulation of the magnetic field
signal. In our work, the excited state atoms alter the sign of the field
amplitude, which is the modulation of the integrated signal of the external
parameter $f(t)$ [Eq. (9)]. We note that $\tau_{\text{A}}^{\ast}=\pi/\xi$ is
not a unique option; any value except $\tau_{A}$ equaling integer multiplies
of $2\pi/\xi$ will work, and we here consider only the simplest case.

$U(t_{j+1}^{m},t_{j}^{m})$ and $R(\tau_{\text{A}}^{\ast})$ alternately govern
the subsequent time evolution according to the randomly generated time
sequence $t_{1}^{m},t_{2}^{m},\dots,t_{K}^{m}$. We can use the method of
induction to derive the expression of the cavity state. Generally, the state
of the cavity field when the $j$th atom leaves the cavity can be written as

$|\psi^{m}(t_{j}^{m})\rangle=e^{i\Theta_{j}^{m}}|\Lambda_{j}^{m}\rangle,$ (15)

where $\Lambda_{j}^{m}$ is the coherent state amplitude at $t_{j}^{m}$ and
$\Theta_{j}^{m}$ is an accumulated global phase. During the subsequent
evolution from $t_{j}^{m}$ to $t_{j+1}^{m}$, the cavity is first driven by
$U(t_{j+1}^{m},t_{j}^{m})$ and then interacts with the $(j+1)$th atom.
Therefore, when the $(j+1)$th atom leaves the cavity, the state of the cavity
field reads

$\displaystyle|\psi^{m}(t_{j+1}^{m})\rangle$
$\displaystyle=R(\tau_{\text{A}}^{\ast})U(t_{j+1}^{m},t_{j}^{m})|\psi^{m}(t_{j}^{m})\rangle$
$\displaystyle=e^{i\Theta_{j+1}^{m}}|-[\Lambda_{j}^{m}+\alpha(t_{j+1}^{m},t_{j}^{m})]\rangle,$
(16)

from which it is straightforward to find the recursive relation
$\Lambda_{j+1}^{m}=-\Lambda_{j}^{m}-\alpha(t_{j+1}^{m},t_{j}^{m})$.With the
initial value $\Lambda_{1}^{m}=-\alpha(t_{1}^{m},t_{0}^{m})$, we find

$\Lambda_{j}^{m}=\sum_{k=1}^{j}(-1)^{j-k+1}\alpha(t_{k}^{m},t_{k-1}^{m}).$
(17)

The accumulated global phase can be derived similarly as
$\Theta_{j}^{m}=\sum_{k=1}^{j}\text{Im}\left[\alpha(t_{k}^{m},t_{k-1}^{m})\Lambda_{k-1}^{m\ast}\right]$,
though it is not related to this work.

Finally, the cavity freely evolves from $t_{K}^{m}$ to the end time $t_{N}$
(where $t_{N}$ is the same for all $M$ measurement processes), and the final
state reads

$|\psi^{m}(t_{N})\rangle=U(t_{N},t_{K}^{m})|\psi^{m}(t_{K}^{m})\rangle=e^{i\Theta_{K+1}^{m}}|\Lambda_{\text{f}}^{m}\rangle.$
(18)

Here, we denote $t_{K+1}^{m}\equiv t_{N}$ and

$\Lambda_{\text{f}}^{m}=\sum_{k=1}^{K+1}(-1)^{K-k+1}\alpha(t_{k}^{m},t_{k-1}^{m}).$
(19)

It is seen from Eq. (17) that the CS signal $\Lambda_{\text{f}}^{m}$ is a
linear combination of the $K+1$ randomly split segments of the original signal
$\alpha(t,t_{0})$. To express $\Lambda_{\text{f}}^{m}$ with the components of
the discretized signal
$\vec{\alpha}=(\alpha_{1},\dots,\alpha_{N})^{\text{T}}$, we define the
increase in $\alpha(t,t_{0})$ during a step of evolution $\tau_{\text{B}}$ by
$\beta_{n}\equiv\alpha(n\tau_{\text{B}},(n-1)\tau_{\text{B}})$, and thus
$\alpha_{n}=\sum_{l=1}^{n}\beta_{l}$. In this manner, Eq. (17) can be
rewritten in the form

$\Lambda_{\text{f}}^{m}=\vec{\mathcal{A}}^{m}\cdot\vec{\beta},$ (20)

Here, $\vec{\beta}=(\beta_{1},\beta_{2},\dots,\beta_{N})^{\text{T}}$ is an
$N$-dimensional column vector and
$\vec{\mathcal{A}}^{m}=(\mathcal{A}_{1}^{m},\mathcal{A}_{2}^{m},\dots,\mathcal{A}_{N}^{m})$
is a row vector whose elements are

$\mathcal{A}_{n}^{m}=\sum_{k=1}^{K+1}(-1)^{K-k+1}H_{k}^{m}(n\tau_{\text{B}}),$
(21)

where $n=1,\dots,N,m=1,\dots,M$ and

$H_{k}^{m}(t)=\begin{cases}1,&t_{k-1}^{m}<t\leq t_{k}^{m},\\\
0,&\text{otherwise}.\end{cases}$ (22)

When all $M$ measurements are complete, the information that we obtained about
$\vec{\alpha}$ is compressed into the vector
$\vec{\Lambda}_{\text{f}}=(\Lambda_{\text{f}}^{1},\Lambda_{\text{f}}^{2},\dots,\Lambda_{\text{f}}^{M})^{\text{T}}$.
By arranging the $M$ row vectors $\vec{\mathcal{A}}^{m}$ into an $M\times N$
dimensional matrix
$\mathcal{A}=[\vec{\mathcal{A}}^{1},\vec{\mathcal{A}}^{2},\dots,\vec{\mathcal{A}}^{M}]$,
we can express the CS measurement of the cavity field in the form

$\vec{\Lambda}_{\text{f}}=\mathcal{A}\cdot\vec{\beta}.$ (23)

A comparison with Eq. (1) reveals that $\vec{\Lambda}_{\text{f}}$ plays the
role of the measurement results $\vec{y}$, $\mathcal{A}$ the sensing matrix,
and $\vec{\beta}$ the signal $\vec{x}$. Once we successfully obtain the
recovered result $\vec{\beta}^{\prime}$, which should coincide with
$\vec{\beta}$ with high accuracy, the coherent field amplitude can be
straightforwardly reconstructed using
$\alpha_{n}=\sum_{l=1}^{n}\beta_{l}^{\prime}$.

So far, we have systematically introduced the CS approach to measure the time-
dependent cavity field, especially describing how to construct the sensing
matrix with coherent operations. It is noted that the methods used in this
cavity QED setup can be generalized to other systems of circuit QED,
optomechanics, and nitrogen-vacancy centers for monitoring specific dynamic
variables.

Figure 2: (color online) Recovered amplitudes of the coherent cavity field
driven by (a) periodical squared and (b) randomly generated protocols. The
time-dependent driving fields $f(t)$ are shown in the upper row. The middle
and bottom rows are the real and imaginary parts of the coherent amplitudes;
the red solid lines correspond to the signals $\vec{\alpha}^{\prime}$
recovered taking the CS approach in the presence of noise, and the blue dashed
lines correspond to the original signals calculated using Eq. (9). The
recovered signals coincide with the original signals with mean-squared error
(MSE) less than $5\times 10^{-4}$. We choose the Nyquist–Shannon sampling
interval to be $\omega_{0}^{-1}$ and set the total sampling times as $N=1000$.
The numbers of CS measurements are $M=220$ and $M=200$ for driving protocols
(a) and (b), respectively. The detuning between the driving field and cavity
mode frequencies is set as $\Delta=0.02\omega_{0}$.

##  V Recovery results

In this section, we take two representative driving protocols as examples to
demonstrate the numerical simulation of the signal recovery via CS techniques
in cavity QED setups. In the first protocol, we design the driving field
$f(t)$ as periodical square pulses with period $\tau=(t_{N}-t_{0})/5$ and a
20% duty cycle [upper row in Fig. 2(a)]. In the second protocol, the driving
field $f(t)$ is a randomly generated smooth function of time [upper row in
Fig. 2(b)]. The original values of the real and imaginary parts of
$\alpha(t,t_{0})$ are calculated using Eq. (9) with respect to the two cases
and shown by blue dashed lines in the middle and bottom rows of Fig. 2,
respectively.

As mentioned in Sec. II, besides dramatically reducing the sampling number,
another advantage of CS is its robustness against noise. In consideration of
practical cases, we add white noise to the magnitude of the driving strength
as $f_{\xi}(t)=f(t)+0.05\omega_{0}\xi$, where $\xi$ is a stochastic variable
uniformly distributed in $[-1,1]$. In the simulation, we randomly generate an
$M\times N$ sensing matrix $\mathcal{A}$ according to Eq. (21) and then apply
$\mathcal{A}$ to $\vec{\beta}_{\xi}$ in obtaining the sample vector
$\vec{\Lambda}_{\text{f}}$. Here, $\vec{\beta}_{\xi}$ is the signal mixed with
noise.

To solve the optimization problems defined in Eq. (5), $\vec{\beta}_{\xi}$
should be transformed into a sparse representation in advance. We assume the
signals $\vec{\beta}_{\xi}$ are sparse in the discrete cosine bases, and the
corresponding transformation $\Phi$ is expressed by

$\Phi_{ij}=\sqrt{\frac{2-\delta_{1,i}}{N}}\cos\left[\frac{\pi}{2N}(i-1)(2j-1)\right],i,j=1,2,\dots,N.$

Equation (23) can then be rewritten as
$\vec{\Lambda}_{\text{f}}=(\mathcal{A}\Phi^{-1})\cdot(\Phi\vec{\beta}_{\xi})$,
where $\Phi\vec{\beta}_{\xi}$ is a sparse vector and $\mathcal{A}\Phi^{-1}$
remains a random matrix. The choice of sparse representation is usually not
unique, and the use of different sparse representations (e.g., discrete
Fourier bases and Walsh bases) will result in different sparsity of the signal
$\Phi\vec{\beta}_{\xi}$. However, finding the representation with minimum
sparsity is not the goal of this work, and the discrete cosine bases work well
for the signals driven by the two different protocols here.

The total evolution time is fixed and discretized into $N=1000$ uniform
intervals, which is the number of measurements required by the Nyquist–Shannon
sampling theorem. In the discrete cosine bases, the sparsity of
$\vec{\beta}_{\xi}$ is approximately $S=50$ for the periodic square wave
driving protocol and $S=40$ for the randomly driving protocol. Therefore, the
numbers of measurements needed for the CS approach are estimated as $\sim
S\log_{2}(N/S)$ as $M=220$ ($M=200$) in the first (second) protocol. Each row
of the sensing matrix $\mathcal{A}$ comprises $N$ elements of $\pm 1$ with a
total of $K$ randomly flips of $1\leftrightarrow-1$. In the simulation, we set
$K=30$ ($K=20$) for the first (second) protocol.

The real and imaginary parts of $\vec{\beta}_{\xi}$ are separately recovered
from the real and imaginary parts of $\vec{\Lambda}_{\text{f}}$, respectively,
via the OMP algorithm. We denote the recovered signal as
$\vec{\beta}^{\prime}$, and the coherent field amplitudes $\vec{\alpha}$ are
then obtained straightforwardly according to
$\alpha_{n}=\sum_{l=1}^{n}\beta_{l}^{\prime}$ and are shown by red solid lines
in the middle and bottom rows of Fig. 2, corresponding to the two driving
protocols. We see that the recovered field amplitudes coincide well with the
theoretically calculations even in the presence of noise. The MSE of both the
real and imaginary parts of $\vec{\alpha}$ are less than $5\times 10^{-4}$.

Figure 3: (color online) Probability of successful recovery for the randomly
driving protocol as a function of $M$. Different numbers of flips $K$ during a
single measurement are presented with different colored markers.

In Fig. 2, the number of measurements $M$ is set around the theoretical lower
bound, which is large enough to guarantee the effectiveness of CS.
Nonetheless, as the sensing matrix $\mathcal{A}$ is randomly generated, we
cannot guarantee the recovery is successful every time, even though the
failure case is rare. We want $M$ to be as small as possible meanwhile
maintaining good recovery quality. To this end, we need to check the
probability of successful recovery with respect to different values of $M$ and
$K$. For every parameter combination of $\\{M,K\\}$, we carry out $1000$
simulations and sum the MSEs of the real and imaginary parts of each
$\vec{\beta}^{\prime}$. We judge the recovery to be successful if the MSE is
less than $2\times 10^{-3}$; otherwise, the recovery has failed. The
probabilities of successful recovery for the randomly driving protocol are
shown in Fig. 3. When $M\geq 220$, the success probability quickly converges
to above 99.5%, and if we lower the success criterion of the MSE to $5\times
10^{-3}$ the success probability stabilizes at 1 for $M\geq 180$.

Another important result obtained from Fig. 3 is that $K$ has an approximate
optimal value for the success probability at a fixed $M$. This result can be
roughly explained in that each column of $\Phi^{-1}$ is an orthogonal vector
of the discrete cosine bases, and then if the row vectors of $\mathcal{A}$ are
almost homogeneously filled with $1$ or $-1$ (small $K$ case), or almost
alternately filled with $1$ and $-1$ (large $K$ case), then the magnitudes of
many elements of $\mathcal{A}\Phi^{-1}$ will be greatly canceled, which yields
a poor quality sensing matrix $\mathcal{A}\Phi^{-1}$. In the randomly driving
case, the success probability roughly increases with $K$ when $K$ is less than
10 and then decreases rapidly when $K$ is larger than 20. This favorable
result implies that we only need a small number of atoms to modulate the
cavity field in each measurement, and the error introduced by the imperfect
control of the atom–photon interaction thus does not have an appreciable
effect.

##  VI CONCLUSIONS AND OUTLOOK

We proposed a measurement approach based on CS techniques to monitor the
dynamic evolution of a driven cavity QED system. Taking this approach, the
time-varying coherent state of the cavity field can be probed and
reconstructed with high accuracy but with an exponentially reduced number of
measurements compared with the number required by the Nyquist–Shannon sampling
theorem. The key issue is the realization of the random sensing matrix via
experimental realizable quantum operations. We suggest that sending a sequence
of largely detuned two-level atoms through the cavity at time points indicated
by the sensing matrix can modulate the cavity field in the designed CS method.

The simulation results show that the reconstructed amplitudes of the coherent
cavity field coincide well with the original amplitudes. Even in the presence
of weak white noise, the error can still be suppressed with MSE no more than
$5\times 10^{-4}$. Additionally, the dependence of the probability of
successful recovery on the numbers of measurements and the number of atoms in
a single measurement was investigated, further strengthening the fact that CS
techniques can greatly reduce the volume of resources consumed in experiments
without affecting the recovery accuracy.

In summary, the quantum measurement method developed in this work reveals the
great potential of CS techniques in quantum sensing and metrology, and more
efforts should made to incorporate the ingenious concept of CS with the
precise quantum control and measurement techniques.

###### Acknowledgements.

This study is supported by the National Natural Science Foundation of China
(Grant No. 12075025).

## References

  * Braginsky  _et al._ (1980) V. B. Braginsky, Y. I. Vorontsov, and K. S. Thorne, Science 209, 547 (1980).
  * Misra and Sudarshan (1977) B. Misra and E. C. G. Sudarshan, J. Math. Phys. 18, 756 (1977).
  * Itano  _et al._ (1990) W. M. Itano, D. J. Heinzen, J. J. Bollinger, and D. J. Wineland, Phys. Rev. A 41, 2295 (1990).
  * Nyquist (1928) H. Nyquist, Transactions of the American Institute of Electrical Engineers 47, 617 (1928).
  * Shannon (1949) C. Shannon, Proceedings of the IRE 37, 10 (1949).
  * Candès  _et al._ (2006) E. J. Candès, J. Romberg, and T. Tao, IEEE Trans. Inf. Theory 52, 489 (2006).
  * Candès and Tao (2006) E. J. Candès and T. Tao, IEEE Trans. Inf. Theory 52, 5406 (2006).
  * Donoho (2006) D. L. Donoho, IEEE Trans. Inf. Theory 52, 1289 (2006).
  * Lustig  _et al._ (2007) M. Lustig, D. Donoho, and J. M. Pauly, Magn. Reson. Med. 58, 1182 (2007).
  * Potter  _et al._ (2010) L. C. Potter, E. Ertin, J. T. Parker, and M. Cetin, Proc. IEEE 98, 1006 (2010).
  * Liu  _et al._ (2019) S. Liu, X.-R. Yao, X.-F. Liu, D.-Z. Xu, X.-D. Wang, B. Liu, C. Wang, G.-J. Zhai, and Q. Zhao, Opt. Express 27, 22138 (2019).
  * Lin and Herrmann (2007) T. T. Y. Lin and F. J. Herrmann, Geophysics 72, SM77 (2007).
  * Katz  _et al._ (2009) O. Katz, Y. Bromberg, and Y. Silberberg, Appl. Phys. Lett. 95, 131110 (2009).
  * Bromberg  _et al._ (2009) Y. Bromberg, O. Katz, and Y. Silberberg, Phys. Rev. A 79, 053840 (2009).
  * Zerom  _et al._ (2011) P. Zerom, K. W. C. Chan, J. C. Howell, and R. W. Boyd, Phys. Rev. A 84, 061804 (2011).
  * Gross  _et al._ (2010) D. Gross, Y.-K. Liu, S. T. Flammia, S. Becker, and J. Eisert, Phys. Rev. Lett. 105, 150401 (2010).
  * Cramer  _et al._ (2010) M. Cramer, M. B. Plenio, S. T. Flammia, R. Somma, D. Gross, S. D. Bartlett, O. Landon-Cardinal, D. Poulin, and Y.-K. Liu, Nat. Commun. 1, 149 (2010).
  * Shabani  _et al._ (2011) A. Shabani, R. L. Kosut, M. Mohseni, H. Rabitz, M. A. Broome, M. P. Almeida, A. Fedrizzi, and A. G. White, Phys. Rev. Lett. 106, 100401 (2011).
  * Rudinger and Joynt (2015) K. Rudinger and R. Joynt, Phys. Rev. A 92, 052322 (2015).
  * Riofrío  _et al._ (2017) C. A. Riofrío, D. Gross, S. T. Flammia, T. Monz, D. Nigg, R. Blatt, and J. Eisert, Nat. Commun. 8, 15305 (2017).
  * Dunbar  _et al._ (2013) J. A. Dunbar, D. G. Osborne, J. M. Anna, and K. J. Kubarych, J. Phys. Chem. Lett. 4, 2489 (2013).
  * Roeding  _et al._ (2017) S. Roeding, N. Klimovich, and T. Brixner, J. Chem. Phys. 146, 084201 (2017).
  * Magesan  _et al._ (2013) E. Magesan, A. Cooper, and P. Cappellaro, Phys. Rev. A 88, 062109 (2013).
  * Ladd  _et al._ (2010) T. D. Ladd, F. Jelezko, R. Laflamme, Y. Nakamura, C. Monroe, and J. L. O’Brien, Nature 464, 45 (2010).
  * Malnou  _et al._ (2019) M. Malnou, D. A. Palken, B. M. Brubaker, L. R. Vale, G. C. Hilton, and K. W. Lehnert, Phys. Rev. X 9, 021023 (2019).
  * Lewis-Swan  _et al._ (2020) R. J. Lewis-Swan, D. Barberena, J. A. Muniz, J. R. K. Cline, D. Young, J. K. Thompson, and A. M. Rey, Phys. Rev. Lett. 124, 193602 (2020).
  * Baraniuk  _et al._ (2008) R. Baraniuk, M. Davenport, R. DeVore, and M. Wakin, Constr. Approx. 28, 253 (2008).
  * Pati  _et al._ (1993) Y. C. Pati, R. Rezaiifar, and P. S. Krishnaprasad, In Proc. of the 27th Asilomar Conference on Signals, Systems and Computers 1, 40 (1993).
  * Lvovsky and Raymer (2009) A. I. Lvovsky and M. G. Raymer, Rev. Mod. Phys. 81, 299 (2009).
  * Brune  _et al._ (1992) M. Brune, S. Haroche, J. M. Raimond, L. Davidovich, and N. Zagury, Phys. Rev. Lett. 45, 5193 (1992).
  * Nogues  _et al._ (1999) G. Nogues, A. Rauschenbeutel, S. Osnaghi, M. Brune, J. M. Raimond, and S. Haroche, Nature 400, 239 (1999).
  * Hansen  _et al._ (2001) H. Hansen, T. Aichele, C. Hettich, P. Lodahl, A. I. Lvovsky, J. Mlynek, and S. Schiller, Opt. Lett. 26, 1714 (2001).
  * Ueda  _et al._ (1992) M. Ueda, N. Imoto, H. Nagaoka, and T. Ogawa, Phys. Rev. A 46, 2859 (1992).
  * Bernu  _et al._ (2008) J. Bernu, S. Deléglise, C. Sayrin, S. Kuhr, I. Dotsenko, M. Brune, J. M. Raimond, and S. Haroche, Phys. Rev. Lett. 101, 180402 (2008).
  * Xu  _et al._ (2011) D. Z. Xu, Q. Ai, and C. P. Sun, Phys. Rev. A 83, 022107 (2011).


# Comments on “A Gibbs sampler for a class of random convex polytopes”

Kentaro Hoffman, Jan Hannig and Kai Zhang
University of North Carolina at Chapel Hill

###### Abstract

In this comment we discuss relative strengths and weaknesses of simplex and
Dirichlet Dempster-Shafer inference as applied to multi-resolution tests of
independence.

##  1 Introduction

We would like to congratulate the authors for their contribution to this
longstanding open problem in mathematical statistics. Their clever
implementation of MCMC to obtain simplex based Dempster-Shafer (DS) samples
for parameters of multinomial distribution and its connection to graph theory
is extremely thought provoking. We expect that the contributions seen in this
paper will have impacts for years to come.

In our comment, we would like to contribute our thoughts particularly on the
applications of the various DS approaches to nonparametric tests of
independence. Nonparametric dependence detection is a classical statistical
problem but recently gains great interest from both statisticians and computer
scientists due to its applications in machine learning. See Hoeffding (1948);
Székely et al. (2007); Gretton et al. (2007); Reshef et al. (2011); Heller et
al. (2013); Chatterjee (2020); Shi et al. (2020), and references therein.

One general approach in nonparametric tests of independence is the multi-
resolution approach. See some recent work by Ma and Mao (2019); Zhang (2019);
Lee et al. (2019); Gorsky and Ma (2018); Zhang et al. (2021). Some advantages
of this approach include uniform consistency, minimax optimal power, clear
interpretability and efficient computation.

The multi-resolution approach reduces the test of independence to the test of
discrete uniformity over multinomial distributions. In this paper, we study
the performance of such tests, particularly when the sample size is small
compared to the resolution. We show potential gains in power and computational
scalability from using an alternative Dirichlet DS method based on unpublished
manuscript of Lawrence et al. (2009).

##  2 Simplex and Dirichlet DS for Multinomial Data

A crucial part to multi-resolution tests of independence is being able to
conduct inference for parameters
$\mathbf{p}=(p_{1},...p_{k})\in\mathcal{P}=\\{\mathbf{p}:\ 0\leq p_{i},\
\sum_{i=1}^{k}p_{i}=1\\}$ using observation $\mathbf{z}=(z_{1},...z_{k})$
following multinomial$(n,\mathbf{p})$ distribution. In particular, we will
compare two different DS approaches to inference parameters $\mathbf{p}$.

The first, which we shall refer to as simplex DS, refers to the method used in
the discussed paper. This version of the DS is based on inverting a very
natural data generating algorithm proposed by Dempster et al. (1966) and
described in detail in Section 2 of the discussed paper. In particular, a
sample from the multinomial distribution is generated by sampling ancillary
variables uniformly at random on a simplex and selecting the multinomial
category depending on which section of the simplex the ancillary variable
falls into, see Figure 1a in the discussed paper. To generate DS samples one
needs to invert this data generating algorithm, i.e., sample new ancillary
variables uniformly and find partitions of the simplex, determining the
parameters $\mathbf{p}$, that would reproduces the observed data. The main
challenge is that such partitions of the simplex exist only if the ancillary
random variables fall into a rather complicated small polytope. The discussed
paper proposes an ingenious Gibbs sampler that does just that. Once such
ancillary variable $u^{(t)}$ is sampled, a linear programming problems need to
be solved to obtain the particle $\mathcal{F}(u^{(t)})$ of parameters
$\mathbf{p}$ that are needed for inference.

While this Gibbs sampler offers a theoretically interesting representation
using shortest path algorithm on graphs, we find that this algorithm has some
drawbacks in terms of scalability. For example, generating 1000 MCMC samples a
7-Dimensional multinomial takes approximately 1 minute. However, converting
the Gibbs sampler’s results $u^{(t)}$ into a convex polytope
$\mathcal{F}(u^{(t)})$ takes nearly 20 times as long as running the Gibbs
sampler. This is due to the large number of vertexes in the polytopes that
result from the Simplex method. The simulation resulted in polygons with an
average of 62 and a max over 250 edges. This places a large computational
burden on the simplex method.

Time Running the Gibbs Sampler | Time converting MCMC to convex polytope
---|---
2.6 seconds | 53.6 seconds

Figure 1: (Top): Runtime comparison of Gibbs sampler vs convex polytope
computation. (Bottom): The number of edges in 1000 Convex polytopes from the
Simplex method.

The other DS method was proposed in an unpublished manuscript of Lawrence et
al. (2009) and we will be referring to as the Dirichlet DS. While Lawrence et
al. (2009) obtain Dirchlet DS using a clever use of Dempster’s rule of
combination starting from $k$ binomial DSs. In particular, let
$\mathbf{W}^{\star}$ follow Dirichlet$(1,z_{1},\ldots,z_{k})$ distribution and
the Dirchlet DS polytope

$\mathcal{F}(\mathbf{W}^{\star})=\\{\mathbf{p}:W_{i}^{\star}\leq p_{i},\
i=1,\ldots,k\\}$ (2.1)

is a simplex with vertexes
$(W_{1}^{\star},W_{2}^{\star},\ldots,W_{k}^{\star})$,
$(W_{1}^{\star}+W_{0}^{\star},W_{2}^{\star},\ldots,W_{k}^{\star})$,
$(W_{1}^{\star},W_{2}^{\star}+W_{0}^{\star},\ldots,W_{k}^{\star}),\ldots,$
$(W_{1}^{\star},W_{2}^{\star},\ldots,W_{k}^{\star}+W_{0}^{\star})$ (Hannig et
al., 2016).

Dirichlet DS has a number of clear computational advantages. First, a
independent sampling from Dirichlet$(1,z_{1},\ldots,z_{k})$ can be done
directly using methods built into most software packages. Additionally the
Dirichlet DS polytope is always a simplex, so for $k=7$ it yields a convex
polytope with 8 vertexes, instead of a maximum of 250. This makes the
Dirichlet DS to scale to much larger problems, e.g. multinomial distributions
with thousands of categories, which is important for testing independence at
high dimensions and/or high resolutions.

Next, many estimators of the multinomial proportions have the following
invariance property: If we merge two categories, the estimate of the merged
proportion is the sum of the proportions that are being merged. Investigation
of (2.1) reveals that Dirichlet DS has this invariance property. Consequently,
inference on proportions for categories in which we have observations is not
influenced by addition or deletion of empty categories. As seen in Section 4.1
of the discussed paper, simplex DS does not have this invariance property, but
whether this invariance property is desirable could be a matter of opinion.

Finally, let us assume that there were originally $n+r$ of the multinoulli
experiments but we were only given results of $n$ and the other $r$ are
missing. Notice that we are not assuming that these $r$ multinoullis are
missing at random, in fact they could be deleted deliberately. The Dirichlet
DS naturally accommodates this additional information by changing the
distribution of $\mathbf{W}^{\star}$ to Dirichlet$(1+r,z_{1},\ldots,z_{k})$
consequently enlarging the polytope in (2.1). Such an enlargement of the DS
polytopes is sometimes called weakening. It is not clear to us if simplex DS
could be weakened to accommodate missing observations.

##  3 DS Test of Independence

For a pair of continuous random variables $(X,Y)$ distributed on $[0,1]^{2}$
we would like to test using $n$ i.i.d. samples $(X_{1},Y_{1})...(X_{n},Y_{n})$
if $X$ is independent of $Y$. An important step towards this goal is to be
able to test if $(X,Y)$ is uniform, or alternatively at what resolutions we
have sufficient power to make a rejection. We adopt the notation of Gorsky and
Ma (2018). For a joint sample space $\Omega=\Omega_{x}\times\Omega_{y}$ it is
possible to create a coarse-to-fine discretization of $\Omega$. At resolution
$k$, we partition $\Omega$ into sets as follows:

$\Omega=\bigcup_{i=1}^{k}\bigcup_{j=1}^{k}I_{i}^{k}\times I_{j}^{k}$

Where $i\in\\{1,k\\}$ $I_{i}^{k}=[\frac{i-1}{k},\frac{i}{k})$. At the coarsest
level, resolution 2, $\Omega$ is discretized into four different pieces:
$\Omega_{11}=[0,\frac{1}{2}]\times[0,\frac{1}{2}]$,
$\Omega_{12}=[0,\frac{1}{2}]\times(\frac{1}{2},1]$,
$\Omega_{21}=(\frac{1}{2},1]\times[0,\frac{1}{2}]$,
$\Omega_{22}=(\frac{1}{2},1]\times(\frac{1}{2},1]$ with

$\Omega=\Omega_{11}\cup\Omega_{21}\cup\Omega_{12}\cup\Omega_{22}$

These sets $\\{I_{i}^{k}\times I_{j}^{k}\\}$ each define statistics

$Z_{i,j}=|(X_{i},Y_{j})\in\Omega_{i,j}|$

which naturally describe a $k\times k$ contingency table. Under the classical
multinomial sampling scheme for $k\times k$ contingency table, the statistics
$Z_{i,j}$ are distributed:

$(Z_{11}...Z_{kk}|\sum_{i,j}Z_{i,j})\sim multinomial(n,(p_{11},...p_{kk}))$

To test if $(X,Y)$ are uniform at resolution k, it suffices to discretize
$(X_{1},Y_{1})...(X_{n},Y_{n})$ into a $k\times k$ contingency table and used
the counts $z_{ij}$ in each bin to obtain the distribution of random polytope
$\Delta\subset\mathcal{P}=\\{\mathbf{p}:\ 0\leq p_{ij},\
\sum_{i=1}^{k}\sum_{j=1}^{k}p_{ij}=1\\}$ following either simplex or Dirichlet
DS distribution. The null hypothesis of uniformity should be rejected at this
resolution if the null-hypothesis point $(k^{-2},\ldots,k^{-2})$ lies well
outside of the bulk of the distribution of $\Delta$. To quantify this we will
find the upper and lower probability of the complement of the smallest ball
centered on the point estimator
$\mathbf{\widehat{p}}=\left(\frac{z_{11}+k^{-2}}{n+1},\ldots,\frac{z_{kk}+k^{-2}}{n+1}\right)$
and containing the null-hypothesis point $(k^{-2},\ldots,k^{-2})$.

To this end we generate $m$ samples $\\{\Delta_{1},...\Delta_{m}\\}$ from
either the simplex or Dirichlet DS. Then, we compute the upper and lower
distances $\mathcal{U}=\\{u_{1},...u_{m}\\}$,
$\mathcal{L}=\\{l_{1},...l_{m}\\}$, where

$l_{i}=\inf\\{\|\mathbf{y}-\mathbf{\widehat{p}}\|_{2}:\
\mathbf{y}\in\Delta_{i}\\},\quad
u_{i}=\sup\\{\|\mathbf{y}-\mathbf{\widehat{p}}\|_{2}:\
\mathbf{y}\in\Delta_{i}\\}.$

Leveraging the convexity of $\Delta_{i}$, $u_{i}$ can be computed as
$u_{i}=\max\\{\|\mathbf{e}-\mathbf{\widehat{p}}\|_{2}:\ \mathbf{e}\in
V(\Delta_{i})\\},$ where $V(\Delta_{i})$ is the set of vertexes of the convex
polytope $\Delta_{i}$, and $l_{i}$ can be found by solving a convex
optimization problem. We then use $u_{i}$ and $l_{i}$ and the distance
$r_{center}=\|\mathbf{\widehat{p}}-(\frac{1}{k^{2}},...\frac{1}{k^{2}})\|_{2}$
to define upper and lower $p$-values as

$\widehat{p}_{upper}(\\{u_{1},...u_{m}\\})=\frac{|\mathcal{U}\geq
r_{center}|}{m},\quad\widehat{p}_{lower}(\\{l_{1},...l_{m}\\})=\frac{|\mathcal{L}\geq
r_{center}|}{m}.$

The upper $p$-value can be then used as usual $p$-value for testing the null
hypothesis of uniformity at resolution $k$. The gap between the upper and
lower $p$-value is specific to DS and measures lack of knowledge. In
particular if the gap between upper and lower $p$-value indicates inability to
make decision about independence at this resolution.

##  4 Simulation Results

To compare the performance of the Simplex and Dirichlet methods, we consider a
simple test of independence. First, 100 data sets of sample size $n=30$ are
generated under either the following null or alternate hypotheses:

$H_{0}\sim\text{Beta}(1,1)^{2},\quad H_{1}\sim\text{Beta}(1,2)^{2}.$ (4.1)

Each of these data sets is then discretized into 2$\times$2, 3$\times$3, and
$6\times 6$ contingency tables and each table is tested for independence using
the simplex DS, the Dirichlet DS, and the classical $\chi^{2}$ tests. To
generate the $p$-values for the test of independence, both the simplex and
Dirichlet DS generate 200 polytopes with a burn-in of 300 for the former. The
purpose of the low sample size ($n=30$) in this simulation is to demonstrate
that as the resolution $k$ and the number of multinomial categories $k^{2}$
increases with sample size held constant, the uncertainty indicated between
the gap between the upper and lower $p$-values increases.

Figure 2: Empirical CDFs of the upper and lower $p$-values for $H_{0}$
analyzed using three tests: (Top Left): Simplex DS, (Top Right) Dirichlet DS,
and (Bottom): $\chi^{2}$. The $x$-axis is the nominal $p$-value, the $y$-axis
is proportion of $p$-values below $p$-value cutoff.

In Figure 2 we present plots of the Empirical CDFs of the upper and lower
$p$-values under the assumption that $H_{0}$ is true. Well calibrated
$p$-values follow a uniform distribution which CDF is represented by the
$45^{\text{o}}$ line. As expected, we see that $p$-values empirical CDFs from
the $\chi^{2}$ test closely follow this dotted line. The upper $p$-values for
both DS tests are below the dotted line, showing that these $p$-values are
conservative, i.e., sub-uniform. Next, we see that while the upper $p$-values
for the Dirichlet and simplex method behave similarly, the lower $p$-values of
the simplex method are more skewed towards rejecting. Consequently, Dirichlet
DS has a much smaller gap between the upper and lower $p$-values than the
simplex DS. Finally, we remark that there are no $p$-values for the 6x6
simplex method as the computation timed out after 2hr without producing a
simplex.

Figure 3:  Empirical CDFs of the upper and lower $p$-values for $H_{0}$
analyzed using three tests: (Top Left): Simplex DS, (Top Right) Dirichlet DS,
and (Bottom): $\chi^{2}$. The $x$-axis is the nominal $p$-value, the $y$-axis
is proportion of $p$-values below the cutoff.

In Figure 3 we show empirical CDFs based on data generated under the alternate
hypothesis in (4.1). All three tests correctly lean towards rejecting the null
hypothesis. In terms of the power of the lower $p$-value, the simplex method
performs similarly to the Dirichlet method. However, the empirical CDFs of
upper $p$-values for the simplex method is lower than their corresponding
empirical CDFs for the Dirichlet DS empirical Cr. This indicates that for the
Dirichlet method has more power to reject $H_{0}$. In addition, we can see the
gaps between lower and upper $p$-value plots increase as the resolution
increases in both the simplex and Dirichlet DS.

As for runtime comparisons, the difference is substantial. Generating one
polytope under Dirichlet DS at the $3\times 3$ level takes approximately 2
seconds while a similar polytope takes nearly 30 seconds under Simplex DS. The
difference in runtime comparison gets larger with the $6\times 6$ level, where
the Dirichlet DS is still under 5 seconds while the Simplex DS timed out after
at least 2 hours.

Method | Runtime to generate one polytope
---|---
Dirichlet Method for $2\times 2$ | 1.76 seconds
Dirichlet Method for $3\times 3$ | 2.15 seconds
Dirichlet Method for $6\times 6$ | 4.15 seconds
Simplex Method $2\times 2$ | 5.31 seconds
Simplex Method $3\times 3$ | 29.60 seconds
Simplex Method $6\times 6$ | Timed out at $>2$hr

##  5 Discussion

Based on this small simulation example, it appears that the Dirichlet DS might
be better suited for performing the test for independence than the simplex DS
both in terms of speed and frequentist performance. Dirichlet DS also has a
relatively easy way to adapt to situations with missing data with a potential
to be applied to adversarial attack scenarios. The one negative is that unlike
simplex DS, the Dirichlet DS does not appear to be fiducial, i.e., based on an
inverse of a data generating algorithm. On balance, it is not clear to us in
what practical situations the simplex DS would be the preferable choice over
the Dirichlet DS when conducting statistical inference.

## Acknowledgements

Jan Hannig’s research was supported in part by the National Science Foundation
under Grant No. IIS-1633074 and DMS-1916115. Kai Zhang’s research was
supported in part by the National Science Foundation under Grant No.
DMS-1613112, IIS-1633212, and DMS-1916237.

## References

  * Chatterjee (2020) Chatterjee, S. (2020).  A new coefficient of correlation.  Journal of the American Statistical Association, in press.
  * Dempster et al. (1966) Dempster, A. P. et al. (1966).  New methods for reasoning towards posterior distributions based on sample data.  The Annals of Mathematical Statistics 37(2), 355–374.
  * Gorsky and Ma (2018) Gorsky, S. and L. Ma (2018).  Multiscale Fisher’s independence test for multivariate dependence.  arXiv preprint arXiv:1806.06777.
  * Gretton et al. (2007) Gretton, A., K. Fukumizu, C. Teo, L. Song, B. Schölkopf, and A. Smola (2007).  A kernel statistical test of independence.  Nips 20, 585–592.
  * Hannig et al. (2016) Hannig, J., H. Iyer, R. C. S. Lai, and T. C. M. Lee (2016).  Generalized fiducial inference: A review and new results.  Journal of the American Statistical Association 111(515), 1346–1361.
  * Heller et al. (2013) Heller, R., Y. Heller, and M. Gorfine (2013).  A consistent multivariate test of association based on ranks of distances.  Biometrika 100(2), 503–510.
  * Hoeffding (1948) Hoeffding, W. (1948).  A non-parametric test of independence.  The Annals of Mathematical Statistics, 546–557.
  * Lawrence et al. (2009) Lawrence, E. C., S. Vander Wiel, C. Liu, and J. Zhang (2009).  A new method for multinomial inference using Dempster-Shafer theory.  Technical report, Los Alamos National Lab.(LANL), Los Alamos, NM (United States).
  * Lee et al. (2019) Lee, D., K. Zhang, and M. R. Kosorok (2019).  Testing independence with the binary expansion randomized ensemble test.  arXiv preprint arXiv:1912.03662.
  * Ma and Mao (2019) Ma, L. and J. Mao (2019).  Fisher exact scanning for dependency.  Journal of the American Statistical Association 114(525), 245–258.
  * Reshef et al. (2011) Reshef, D. N., Y. A. Reshef, H. K. Finucane, S. R. Grossman, G. McVean, P. J. Turnbaugh, E. S. Lander, M. Mitzenmacher, and P. C. Sabeti (2011).  Detecting novel associations in large data sets.  Science 334(6062), 1518–1524.
  * Shi et al. (2020) Shi, H., M. Hallin, M. Drton, and F. Han (2020).  Rate-optimality of consistent distribution-free tests of independence based on center-outward ranks and signs.  arXiv preprint arXiv:2007.02186.
  * Székely et al. (2007) Székely, G. J., M. L. Rizzo, and N. K. Bakirov (2007, 12).  Measuring and testing dependence by correlation of distances.  The Annals of Statistics 35(6), 2769–2794.
  * Zhang (2019) Zhang, K. (2019).  BET on independence.  Journal of the American Statistical Association 114(528), 1620–1637.
  * Zhang et al. (2021) Zhang, K., Z. Zhao, and W. Zhou (2021).  BEAUTY powered BEAST.  arXiv preprint arXiv:2103.00674.


Efficient Hotspot Switching in Plasmonic Nanoantennas using Phase-shaped Laser
Pulses controlled by Neural Networks

Alberto Comin∗ and Achim Hartschuh

Department of Chemistry and Center for NanoScience (CeNS), LMU Munich, 81377,
Munich, Germany

<EMAIL_ADDRESS>

Abstract:  We present a novel procedure for manipulating the near-field of
plasmonic nanoantennas using neural network-controlled laser pulse-shaping. As
model systems we numerically studied the spatial distribution of the second
harmonic response of L-shaped nanoantennas illuminated by broadband laser
pulses. We first show that a trained neural network can be used to predict the
relative intensity of the second-harmonic hotspots of the nanoantenna for a
given spectral phase and that it can be employed to deterministically switch
individual hotspots on and off on sub-diffraction length scale by shaping the
spectral phase of the laser pulse. We then demonstrate that a neural network
trained on a $90~{}\mathrm{nm}\times 150~{}\mathrm{nm}$ nano-L can in addition
efficiently predict the hotspot intensities in an antenna with different
aspect ratio after minimal further training for varying spectral phases. These
results could lead to novel applications of machine-learning and optical
control to nanoantennas and nanophotonics components.

OCIS codes: 190.7110 Ultrafast nonlinear optics, 250.5403 Plasmonics

## References and links

  * [1] L. Novotny and N. van Hulst, “Antennas for light,” Nature Photonics 5, 83–90 (2011).
  * [2] L. Piatkowski, N. Accanto, and N. F. van Hulst, “Ultrafast Meets Ultrasmall: Controlling Nanoantennas and Molecules,” ACS Photonics 3, 1401–1414 (2016).
  * [3] M. I. Stockman, “Nanoplasmonics: past, present, and glimpse into future,” Optics Express 19, 22029–22106 (2011).
  * [4] M. Gu, X. Li, and Y. Cao, “Optical storage arrays: A perspective for future big data storage,” Light: Science and Applications 3 (2014).
  * [5] B. Ashall, J. F. López-Barberá, É. McClean-Ilten, and D. Zerulla., “Highly efficient broadband ultrafast plasmonics,” Optics Express 21, 27383 (2013).
  * [6] M. I. Stockman, K. Kneipp, S. I. Bozhevolnyi, S. Saha, A. Dutta, J. Ndukaife, N. Kinsey, H. Reddy, U. Guler, V. M. Shalaev, A. Boltasseva, B. Gholipour, H. N. S. Krishnamoorthy, K. F. MacDonald, C. Soci, N. I. Zheludev, V. Savinov, R. Singh, P. Groß, C. Lienau, M. Vadai, M. L. Solomon, D. R. Barton, M. Lawrence, J. A. Dionne, S. V. Boriskina, R. Esteban, J. Aizpurua, X. Zhang, S. Yang, D. Wang, W. Wang, T. W. Odom, N. Accanto, P. M. de Roque, I. M. Hancu, L. Piatkowski, N. F. van Hulst, and M. F. Kling, “Roadmap on plasmonics,” Journal of Optics 20, 043001 (2018).
  * [7] T. Brixner, F. J. García de Abajo, J. Schneider, C. Spindler, and W. Pfeiffer, “Ultrafast adaptive optical near-field control,” Physical Review B 73, 125437 (2006).
  * [8] G. Piredda, C. Gollub, R. de Vivie-Riedle, and A. Hartschuh, “Controlling near-field optical intensities in metal nanoparticle systems by polarization pulse shaping,” Applied Physics B 100, 195–206 (2010).
  * [9] M. Aeschlimann, M. Bauer, D. Bayer, T. Brixner, S. Cunovic, F. Dimler, A. Fischer, W. Pfeiffer, M. Rohmer, C. Schneider, F. Steeb, C. Strüber, and D. V. Voronine, “Spatiotemporal control of nanooptical excitations.” Proceedings of the National Academy of Sciences of the United States of America 107, 5329–33 (2010).
  * [10] D. Brinks, M. Castro-Lopez, R. Hildner, and N. F. van Hulst, “Plasmonic antennas as design elements for coherent ultrafast nanophotonics,” Proceedings of the National Academy of Sciences 110, 18386–18390 (2013).
  * [11] Y. Kojima, Y. Masaki, and F. Kannari, “Control of ultrafast plasmon pulses by spatiotemporally phase-shaped laser pulses,” Journal of the Optical Society of America B 33, 2437 (2016).
  * [12] M. Aeschlimann, M. Bauer, D. Bayer, T. Brixner, S. Cunovic, A. Fischer, P. Melchior, W. Pfeiffer, M. Rohmer, C. Schneider, C. Strüber, P. Tuchscherer, and D. V. Voronine, “Optimal open-loop near-field control of plasmonic nanostructures,” New Journal of Physics 14, 033030 (2012).
  * [13] J. S. Huang, D. V. Voronine, P. Tuchscherer, T. Brixner, and B. Hecht, “Deterministic spatiotemporal control of optical fields in nanoantennas and plasmonic circuits,” Physical Review B 79, 195441 (2009).
  * [14] T. Brixner, F. J. García de Abajo, J. Schneider, C. Spindler, and W. Pfeiffer, “Ultrafast adaptive optical near-field control,” Physical Review B 73, 125437 (2006).
  * [15] P. Tuchscherer, C. Rewitz, D. V. Voronine, F. J. García de Abajo, W. Pfeiffer, and T. Brixner, “Analytic coherent control of plasmon propagation in nanostructures,” Optics Express 17, 14235–14259 (2009).
  * [16] T. Baumert, T. Brixner, V. Seyfried, M. Strehle, and G. Gerber, “Femtosecond pulse shaping by an evolutionary algorithm with feedback,” Applied Physics B: Lasers and Optics 65, 779–782 (1997).
  * [17] A. Meyer-Baese and V. Schmid, “Genetic Algorithms,” in “Pattern Recognition and Signal Analysis in Medical Imaging,” (Elsevier, 2014), pp. 135–149.
  * [18] I. Goodfellow, Y. Bengio, and A. Courville, _Deep Learning_ (MIT Press, 2016). http://www.deeplearningbook.org.
  * [19] V. Myroshnychenko, E. Carbó-Argibay, I. Pastoriza-Santos, J. Pérez-Juste, L. M. Liz-Marzán, and F. J. García de Abajo, “Modeling the Optical Response of Highly Faceted Metal Nanoparticles with a Fully 3D Boundary Element Method,” Advanced Materials 20, 4288–4293 (2008).
  * [20] U. Hohenester and A. Trügler, “{MNPBEM} – a matlab toolbox for the simulation of plasmonic nanoparticles,” Computer Physics Communications 183, 370 – 381 (2012).
  * [21] Y. Lecun, Y. Bengio, and G. Hinton, “Deep learning,” Nature 521, 436–444 (2015).
  * [22] R. Esteban, R. Vogelgesang, J. Dorfmüller, A. Dmitriev, C. Rockstuhl, C. Etrich, and K. Kern, “Direct near-field optical imaging of higher order plasmonic resonances,” Nano Letters 8, 3155–3159 (2008).
  * [23] S. Simoncelli, Y. Li, E. Cortés, and S. A. Maier, “Nanoscale control of molecular self-assembly induced by plasmonic hot-electron dynamics,” ACS Nano 12, 2184–2192 (2018). PMID: 29346720.
  * [24] F. J. García de Abajo, “Optical excitations in electron microscopy,” Rev. Mod. Phys. 82, 209–275 (2010).
  * [25] G. Bachelier, J. Butet, I. Russier-Antoine, C. Jonin, E. Benichou, and P.-F. Brevet, “Origin of optical second-harmonic generation in spherical gold nanoparticles: Local surface and nonlocal bulk contributions,” Physical Review B 82, 235403 (2010).
  * [26] R. Czaplicki, J. Mäkitalo, R. Siikanen, H. Husu, J. Lehtolahti, M. Kuittinen, and M. Kauranen, “Second-harmonic generation from metal nanoparticles: Resonance enhancement versus particle geometry,” Nano Letters 15, 530–534 (2015).

##  1 Introduction

Plasmonic nanoantennas can efficiently focus broadband optical fields to form
nanometer sized ultrafast hotspots [1, 2]. They could thus become a key
component of future nanophotonic devices, combining high storage density and
fast processing of information [3, 4, 5, 6]. The ability to deterministically
control the brightness of individual hotspots within a nanoantenna by varying
the input optical fields would add enormous flexibility to antenna schemes.
Towards this goal, several groups have explored the possibility of controlling
the near-field in plasmonic structures by tuning the spectrum, polarization
and spectral phase profile of the input laser excitation [7, 8, 3, 9, 10, 11].

Spectral amplitude shaping could be applied in the simple case in which the
nanoantenna features spectrally distinct plasmonic resonances connected to
hotspots at different positions. For such nanoantennas, it is possible to lit
each hotspot individually by tuning the color of the excitation light, albeit
at the expense of reduced spectral bandwidth and thus temporal resolution. In
more general cases, spatio-temporal control of plasmonic near-fields is
understood to have two main control mechanisms [7]. Efficient tuning can be
achieved by polarization pulse shaping which exploits the interference of
plasmonic near-field modes with different polarization responses. Polarization
pulse shaping has wide applicability and was successfully used to
experimentally demonstrate sub-wavelength hotspot switching on different
antenna systems [12, 7]. The second control mechanism is based on spectral
phase shaping without the need of polarization control. In this scheme, the
spatial control of non-linear responses is achieved by imprinting a spectral
phase profile on the incident laser pulses which is set to compensate the
phase response of the nanoantenna at a particular position [10, 13, 14].

Whereas for simple systems the parameters for the coherent control of the
hotspot position can mostly be derived analytically, the optimum pulse
characteristics for more complex multi-modal systems cannot be predicted [15].
Although evolutionary algorithms, such as genetic algorithms (GA), could be
employed as versatile optimization tools for such problem sets, they typically
require a large number of experimental iterations and are thus limited in
efficiency [16]. Furthermore, the results obtained by GA are hard to
generalize and replicate, while the GA must be re-optimized for each sample
and experimental configuration. For example, when an external perturbation
shifts the plasmonic resonance of the nanoantenna, the GA can only adapt by
crossover and random mutation, without taking into account previously learned
information about the sample [17].

We propose a combined use of GA and neural-network (NN) as a more efficient
and deterministic way to control the near-field in plasmonic nanoantennas. NN
have a layered structured which can encode more general information in the
first layers and more specific information in the last layers [18]. This is
one of the main reasons why NNs became very popular: It is possible to re-use
a pre-trained NN on a different domain after just minimal training of the
final layers [18].

In the following, we will show that a NN consisting of only four fully
connected layers can accurately encode the dependence of the near-field on the
spectral phase of the incident laser pulses. [10, 13, 14]. Moreover, a NN
trained on a specific nanoantenna can also be used with minimal further
training on nanoantennas with different size and aspect ratio.

As an example of the efficacy of this approach, we apply the GA-NN combination
to achieve second harmonic generation (SH) hotspot switching in L-shaped
plasmonic nanoantennas by means of spectral-phase shaping. L-shaped
nanoantennas were selected as simple model systems supporting multiple
plasmonic resonances within the spectrum of ultrafast Ti:Sapphire lasers.
Whereas polarization pulse shaping could be exploited as an additional
powerful degree of freedom, here we wanted to limit the complexity of the
optimization scheme. Although sub-wavelength and second-harmonic hotspot
switching have been separately demonstrated [12, 10], the control of SHG in
sub-diffraction nanoantennas was not shown before.

##  2 Setup

In order to use NNs for optimal control on real nanoantennas, a reasonable
feasibility step is to train them using realistic simulated data. Producing
high-quality nanostructures and measuring them is still a time-consuming task,
but a vast simulated dataset can be produced in much shorter time. Plasmonic
nanostructures are particularly convenient to simulate: Their optical
properties are dominated by the surface density of charge and current, and it
is possible to accurately model them using a boundary element method (BEM)
[19]. For this purpose we applied a customized version of the Matlab MNPBEM
toolbox [20], which was extended to perform non-linear optics simulations:
additional details are provided in the appendix in section 5.2.

Fig. 1: a) Flow diagram of the training phase of the neural network using the
population generated by a genetic algorithm. A random initial population of 20
spectral phases is chosen. At each iteration, the fitness is calculated
according to a pre-specified goal: The fittest individuals are kept unchanged
and the rest undergo mutation or cross-over. The least fit individuals are
optimized by the NN using back-propagation. b) Training history of the GA
without assistance from the NN. The graph reports the fitness of the best and
worst individuals, and the average fitness. The fitness function (vertical
axis) was chosen to maximize the relative SHG flux at a target hotspot of a
$90~{}\mathrm{nm}\times 250~{}\mathrm{nm}$ gold L-shaped nanoantenna,
indicated by an arrow in the inset. c) Similar to b) but with the GA assisted
by an un-trained NN. d) similar to b) but with the GA assisted by a pre-
trained NN: The fitness of the best individuals converges after just one
iteration, the fitness of the worst individual remains lower due to the random
mutations introduced by the GA.

We trained the NN using the populations produced by a genetic algorithm, as
illustrated in Figure 1a. Using a GA has a double advantage: First, it allows
to have a fast feedback on the nanoantenna design, e.g. if it is likely or not
to be controllable withing the given experimental parameters. Second, it
generates a varied dataset of both near-optimal and pseudo-random solutions.
At each iteration the population was used to train the NN and then sub-divided
according to the relative fitness: a few of the fittest individuals were kept
unchanged, 80% of the remaining fittest were used to create crossover children
and the rest underwent random mutation. The least fit individuals were
optimized by the NN using a back-propagation algorithm [21]. The choice of
optimizing the least fit individuals allows the possibility to start with an
untrained NN. As the NN was trained, its accuracy increased until it started
to accelerate the training of the GA. When the NN was sufficiently trained, it
could be directly used to produce optimal solutions.

We note that an alternative and simpler way to generate a set of training
spectral phases would be to use a random numbers generator. However, random
phases tend to have fast oscillations that increase the pulse duration and
reduce any nonlinearity.

For the results in this paper, the NN was trained to predict the second-
harmonic (SH) flux at the hotspots on the surface of the nanoantenna utilizing
the spectral phase of the laser excitation pulse as input. For efficiency
reasons, the spectral phase was specified at 6 nodal points, evenly
distributed within 2 standard deviations from the central frequency of the
laser excitation pulse. The phases were then interpolated on a finer mesh,
using piecewise cubic Hermite interpolation (pchip). An advantage of this
approach is that it produces reasonably smooth phase profiles, without fast
oscillations. The spectral phase at the nodal points was bounded between $\pm
100~{}\mathrm{rad}$, an interval chosen to match the capability of a standard
$128~{}\mathrm{pixel}$ pulse-shaper. The spectral amplitude was taken to be
Gaussian with central frequency of $375~{}\mathrm{THz}$ ($800~{}\mathrm{nm}$)
and standard deviation of $26.5~{}\mathrm{THz}$, corresponding to a full-
width-half-max (FWHM) temporal intensity duration of $10~{}\mathrm{fs}$. These
constrains will make the results easier to test in an experimental setup
equipped with a pulse-shaper and a femtosecond laser source.

##  3 Results and Discussion

The efficacy of the NN to accelerate the training of a GA is demonstrated in
Figure 1. Panel b) shows the convergence of the GA without assistance of the
NN for a $90~{}\mathrm{nm}\times 250~{}\mathrm{nm}$ L-shaped gold nanoantenna.
The bars indicate the fitness of the best and worst individuals and the
average fitness. The fitness function was chosen to maximize the relative SH
flux at a specific hotspot, indicated in the inset by an arrow:
$R_{\phi}=h_{\phi,i}/\underset{i\neq j}{\max}(h_{\phi,j})$, where $R_{\phi}$
is the fitness, $h_{\phi,i}$ is the SH flux at the target hotspot ’$i$’
resulting from a specific spectral phase profile $\phi$ and $\underset{i\neq
j}{\max}(h_{\phi,j})$ is the maximum SH flux using the same phase $\phi$ over
all the other hotspots. The inset shows the generated distribution of the SH
field at the outer surface of the nanoantenna for a flat phase pulse. Figure
1c,d shows that a un-trained neural network can reduce the number of needed
iterations substantially and that, with a pre-trained NN, the fitness of the
best individual converges to the optimal value in just one iteration.

Optimizing the relative local intensity of the SH allows to switch the
position of the brightest hotspot.

Fig. 2: Hotspot switching in a $90~{}\mathrm{nm}\times 250~{}\mathrm{nm}$ gold
nanoantenna. a) Spectral phases which maximize the relative flux intensity of
two hotspots (labeled ’1’ and ’2’) as found by the GA+NN algorithm. The
hotspot positions are indicated by arrows in panels b,c). The laser spectrum
(orange shaded area) and the nanoantenna absorption spectrum (green shaded
area) are also shown. b,c) SHG flux intensity at the outer surface of the
nano-antenna corresponding to the optimal phases shown in panel a). Panels
d,e,f) report similar information as panels a,b,c) but with a different
optimization goal: to maximize the absolute, rather than the relative, SH flux
intensity at a specific hotspot

The idea of hotspot switching is illustrated in Figure 2. Panel a) shows the
laser spectrum, the nanoantenna absorption spectrum and the spectral phases
which maximize the relative SH flux $R_{\phi}$ as defined above at the
hotspots labeled ’1’ and ’2’ in b) and c), respectively. Panel b) shows the SH
surface field: It can be recognized that the location $x$ of the brightest
hotspot changes according to the optimization goal $F_{\phi}(x)$. The color
scale was normalized to the maximum SH surface flux for a flat phase laser
pulse. It can be recognized that maximizing the relative intensity of a
specific hotspot comes with a reduction of the overall SHG flux intensity of a
factor of five.

Whilst optimizing the relative hotspot intensity is interesting, it result in
a overall decrease of the SH intensity, which might be detrimental for actual
applications. A different optimization goal is to maximize the absolute value
of the SH flux for a given hotspot. The fitness is now defined by
$A_{\phi}=h_{\phi}/h_{0}$, where $h_{\phi}$ is again the SH flux at the target
hotspot resulting from a specific spectral phase $\phi$ and $h_{0}$ is the SH
flux over all the nanoantenna surface for a flat spectral phase profile. In
this case we are not guaranteed that the target hotspot will be the brightest
one, but the resulting local fields will be larger.

We present the results of this kind of optimization in Figure 2d,e,f. For the
a 90 nm250 nm gold nanoantenna the intensity gain was about a factor of 5 as
compared to the optimization in Figure 3b, c as can be seen from the
respective color bars.

Fig. 3: a) Relative SH flux intensity for a target hotspot of a gold L-shaped
nanoantenna with size $90~{}\mathrm{nm}\times 250~{}\mathrm{nm}$ for different
spectral phase profiles of the incident laser pulse. The blue line refers to
the simulated value, the red dashed line refers to the value predicted by the
NN. The inset shows the nanoantenna with target hotspots marked. Panel b,c)
shows the performance of same NN but used to predict the SHG hotspot
intensities for a nanoantenna with different size: $90~{}\mathrm{nm}\times
150~{}\mathrm{nm}$ with no further training, and after training only the
output layer of the NN for 10 epochs. d,e,f) report similar results as a,b,c)
but for a different hotspot. The figure illustrates the flexibility of the NN
with respect to the sample size and aspect ratio.

The results shown in Figure 3 demonstrate the possibility of switching the
position of the main hotspot between different corners of an L-shaped
nanoantennas. The displacement distance, for the illustrated case, was
$100~{}\mathrm{nm}$, well below the diffraction limit for the laser excitation
at $800~{}\mathrm{nm}$. The maximum contrast values, integrated over the whole
frequency range, were: $R_{\phi}=5.2$ and $R_{\phi}=4.6$, for the two shown
hotspots, which should be large enough to be experimentally tested. Such
hotspot switching could be useful for triggering local non-linear phenomena,
either in a adjacent sample object, such as a fluorescent molecule or
semiconductor nanocrystal or, more simply, in the substrate of the nanoantenna

Analyzing the spectral phase profiles found by NN-GA combination in Figure 2a)
and d) for the two different fitness functions $R_{\phi}$ and $A_{\phi}$, we
can comment on the mechanisms underlying hotspot control in the present
system. In case of the optimization of the relative flux $R_{\phi}$ at hotspot
’1’, the phase profile $\phi$ remains essentially flat between 320 and 390
THz, the spectral range of the plasmon resonance connected to this particular
hotspot. From 395 to 450 THz, on the other hand, the phase becomes very large
featuring huge higher order variations. This leads to an efficient suppression
of the SH flux of the plasmon mode in this spectral range, which is associated
with a different hotspot distribution including a peak at position ’2’.
Optimizing the relative flux at hotspot ’2’ is seen to result in the opposite
spectral phase characteristics with a flat curve between 450 to 395 THz and
strong phase variations below.

For the target function $A_{\phi}$, which maximizes the absolute SH flux at a
given hotspot position, the observed phase variations are much smaller. Here,
the dominant control mechanism will be that of local pulse compression, also
mentioned in the introduction [10, 13, 14]. In other words, the optimum
spectral phase profiles found by the NN+GA scheme will compensate the local
phase response of the nanoantenna at the particular positions thereby
maximizing the non-linear response.

An important advantage of NN-based control is expected to arise from the
adaptability of a trained network to similar problem sets. Applied to coherent
control of plasmonic nanoantennas, this means that we could train our NN on a
specific nanoantenna and then use the learned weights also for nanoantennas
with different sizes or aspect ratios. The feasibility of this idea is
demonstrated in Figure 3. Panel a) shows the prediction accuracy of a pre-
trained NN for a set of randomly generated spectral phases and for a
nanoantenna of size $90~{}\mathrm{nm}\times 250~{}\mathrm{nm}$, using a
training set of $800.000$ phases, a test set of $90.000$ phases and $1000$
epochs. The blue line refers to the true value for the relative intensity of a
target hotspot, marked with an arrow in the inset of panel c). The red line
shows the prediction by the NN. The overall root square error was about
$2.6\times 10^{-2}$. Panel b) shows the performance of the same NN but for a
nanoantenna with size $90~{}\mathrm{nm}\times 150~{}\mathrm{nm}$, and panel c)
after fine tuning only the last layer of the NN for $10$ epochs. The root mean
square error was $0.98$ with no training and $5.7\times 10^{-2}$ after fine
tuning. The NN was still able to quite accurately predict the relative
intensity of the SH hotspot. Panels d,e,f) report similar results, but for a
different hotspot, which was bright of the training antenna and dim on the
test one: also in this case the NN was able to make quite accurate
predictions, with a relative error of about $10$%.

These results indicate that it should be feasible to use a NN trained with
simulated data to control a real nanoantenna in an experimental setup. Due to
small changes in the dielectric environment and uncertainties in the
fabrication process, the plasmon frequencies could be different between the
real and the simulated nanoantennas. The ability of the NN to compensate for
the frequency shifts caused by changes in size and aspect ratio suggests that
the they could also adapt to shifts caused by the dielectric environment.

##  4 Conclusion

We introduced a novel scheme to control the near-field of plasmonic
nanostructures based on a neural network in conjunction with a genetic
algorithm. The neural network accelerates the optimization of the genetic
algorithm and stores information about the sample, which can readily be
generalized to other samples, with none or minimal further training. In order
to prove the efficacy of this approach, we showed how the algorithms can find
the optimal spectral phases for switching the position of the brightest
hotspot in an L-shaped gold nanoantenna. We also showed that a NN trained on a
specific nanoantenna provides quite accurate results also for a nanoantenna
with different size and aspect ratio, even without any further training. Our
results suggests that NNs are a powerful tool for optimal control of near-
fields at the nanoscale and, in perspective, for more complex nanoplasmonics
and nanophotonics devices. Coherent control of hotspot positions on the
nanoscale could be experimentally observed using photoemission electron
microscopy (PEEM) or scanning near-field optical microscopy with passive
probes [12, 22]. Possible further developments include the coupling of
nanoantennas to different emitters or waveguides situated near the hotspots
providing a means for all-optical ultrafast switching or the spatially
selective initiation of photochemical reactions at plasmonic nanoantennas
[23].

##  5 Appendix

###  5.1 Neural Network

The neural-network (NN) used to obtain the results shown in the article was a
multi-layer perceptron composed of four fully-connected layers, plus an input
and an output layer. The activation function was hyperbolic tangent (tanh) for
all layers except for the output layer, which had linear activation. The input
layer contained 6 neurons representing the spectral phase of the excitation
laser at 6 equidistant nodal points evenly distributed within 2 standard
deviations from the central frequency of the laser excitation pulse. The
phases were then interpolated on a finer mesh, using piecewise cubic Hermite
interpolation (pchip). The sizes of the inner layers were: 100, 100, 50, 50.
The output layer was used to predict the second-harmonic flux intensity for
the 12 corners of a L-shaped nanoantenna.

The number of layers and neurons was chosen as a trade off between prediction
accuracy and training speed. A small network requires a relatively small
training dataset, and it can potentially be trained using experimental data
within a reasonable measurement time. The network was trained using a custom
implementation of the back-propagation algorithm in Matlab. The language was
chosen for easier integration with existing toolboxes for optics and
plasmonics. The accuracy of the code was tested using a widely-used machine-
learning toolbox (Keras2 with TensorFlow back-end). The NN was trained using
using the population generated by a genetic algorithm (GA). Several training
sessions were run while changing optimization goals, e.g. relative or absolute
hotpot intensity for different choices of hotspots. The final train set size
was about 1 million, and the test set size about 0.1 million. A small
($10^{-6}$) $L_{2}$ regularization factor was used, however the regularization
choice was not crucial when training the NN using noiseless simulated data: It
will be important for training using real experimental data.

###  5.2 SHG Simulations

The second-harmonic response of the L-shaped nanoantennas was calculated using
the boundary-element method, as described by Garcia de Abajo et al. [24]. An
open-source Matlab toolbox (MNPBEM) was used to calculate the linear response
of the nanoantennas [20]. The surface density of charge and current was than
used to estimate the dipolar contributions to the SHG from each surface
elements [25, 26]. For simplicity, the bulk contribution to the SHG was not
considered. The simulation was performed in frequency domain using 200 points
between $200~{}nm$ and $1000~{}nm$, the results were then interpolated using a
finer mesh over the spectral range of the laser excitation. The simulated
laser pulse was a Gaussian with temporal full-width-half-maximum width of
$10~{}\mathrm{fs}$ and central frequency of $375~{}\mathrm{THz}$
($800~{}\mathrm{nm}$).

##  6 Funding Information

Deutsche Forschungsgemeinschaft (DFG) HA4405/8-1

##  7 Acknowledgments

The authors thanks Giovanni Piredda, Nicolas Coca Lopez, Veit Giegold and
Richard Ciesielski for helpful discussions.