# Controlling Many-Body Quantum Chaos: Bose-Hubbard systems

Lukas Beringer  Institut für Theoretische Physik, Universität Regensburg,
D-93040 Regensburg, Germany  Mathias Steinhuber  Institut für Theoretische
Physik, Universität Regensburg, D-93040 Regensburg, Germany  Juan Diego Urbina
Institut für Theoretische Physik, Universität Regensburg, D-93040 Regensburg,
Germany  Klaus Richter  Institut für Theoretische Physik, Universität
Regensburg, D-93040 Regensburg, Germany  Steven Tomsovic  Institut für
Theoretische Physik, Universität Regensburg, D-93040 Regensburg, Germany
Department of Physics and Astronomy, Washington State University, Pullman, WA
USA<EMAIL_ADDRESS>
###### Abstract

This work develops a quantum control application of many-body quantum chaos
for ultracold bosonic gases trapped in optical lattices. It is long known how
to harness exponential sensitivity to changes in initial conditions for
control purposes in classically chaotic systems. In the technique known as
targeting, instead of a hindrance to control, the instability becomes a
resource. Recently, this classical targeting has been generalized to quantum
systems either by periodically countering the inevitable quantum state
spreading or by introducing a control Hamiltonian, where both enable localized
states to be guided along special chaotic trajectories toward any of a broad
variety of desired target states. Only strictly unitary dynamics are involved;
i.e., it gives a coherent quantum targeting. In this paper, the introduction
of a control Hamiltonian is applied to Bose-Hubbard systems in chaotic
dynamical regimes. Properly selected unstable mean field solutions can be
followed quite rapidly to states possessing precise phase relationships and
occupancies. In essence, the method generates a quantum simulation technique
that can access rather special states. The protocol reduces to a time-
dependent control of the chemical potentials, opening up the possibility for
application in optical lattice experiments. Explicit applications to custom
state preparation and stabilization of quantum many-body scars are presented
in one- and two-dimensional lattices (three-dimensional applications are
similarly possible).

##  1 Introduction

The exponential instability inherent in chaotic dynamics typically leads a
system towards relaxation, equilibration, and/or thermalization. Naturally,
chaos poses fundamental challenges for controlling system dynamics as well.
Nevertheless, there are circumstances in which chaos provides a resource for
control, such as gravity assist conjectured by Ulam [1], nicely illustrated as
a control problem by the diversion of the International Sun-Earth Explorer
satellite to a Giacobini-Zinner comet flyby [2, 3], and more comprehensively,
by the theory of classical targeting [4, 5, 6, 7, 8], which falls under the
more general moniker, controlling chaos [9]. Roughly speaking, for fully
chaotic systems the origin of the resource is the existence of heteroclinic
motion (trajectories) [10]. In a nutshell, there exists a set of neighboring
initial conditions, which lead to trajectories that deviate with a maximal
exponential divergence from the trajectory of the given initial conditions.
Likewise, but the inverse process, there exists a set of trajectories which
approach the final conditions in a maximal exponential sense ending up in its
neighborhood. Those rare trajectories belonging to both sets capture the
essence of heteroclinic motion, and give rise to the opportunity to make a
tiny, but very specific, change in the initial conditions that accelerates
greatly how quickly the system dynamically evolves to the predetermined final
point. Hence, a great deal of the effort in controlling classical chaos
(targeting) is the identification of the precise and quite rare (heteroclinic)
trajectory which accomplishes the goal in an optimal way.

Recently, the translation of classical targeting into the quantum realm has
been introduced and applied to a longstanding paradigm of simple chaotic
systems [11, 12], the kicked rotor [13, 14]. For quantum systems with a well
defined classical analog, it turns out that the same heteroclinic motion can
be exploited both to arrive at a predetermined, possibly exotic, target or
final state, and greatly accelerate the way there. It can also be done with
exclusively unitary dynamics, and thus be fully coherent. The crucial new
element not found in the classical case essentially arises from the
uncertainty principle [15]. The closest a quantum state can get to the initial
conditions of a classical analog is in the form of a minimum uncertainty wave
packet and a coherent state, depending on the context. By virtue of a Wigner
transform analysis [16] of these states, they correspond most closely to a
Liouvillian density of neighboring initial conditions as opposed to a single
trajectory. Thus, a quantum analog of classical targeting must also control
the Liouvillian flow around the heteroclinic trajectory; see Fig. 1. It is
this extra requirement

Figure 1: Schematic illustration of an initially localized quantum state
guided along a classical transport trajectory. A localized quantum state, here
represented as an orange cut through its Wigner function, centered around an
initial phase space point (upper right dot) can be guided to a specific target
state (lower left cross) by following a solution of the classical limit. In a
generic system exhibiting quantum chaos the localized state would equilibrate
in a logarithmically short time period, and thus the spreading must be
counteracted in a control protocol. The blue grid schematically represents a
Hamiltonian landscape generating unstable dynamics.

that led to the introduction of two methods of optimal coherent quantum
targeting. The first technique follows the local stability analysis and
counteracts from time to time the inevitable spreading of the density [11].
The second technique is best conceived of as a form of quantum simulation [17]
in which a quite different Hamiltonian simulates the heteroclinic evolution in
a stable way [12]. In this work, the second technique is being applied to an
important, experimentally realizable many-body system, i.e., the Bose-Hubbard
model [18, 19].

Quantum control techniques comprise a large longstanding research field with
much of the early work prompted by the dream of controlling chemical reactions
[20, 21, 22, 23, 24, 25]. A survey [26] and an overview of optimal control
theory [27] provide more information on this larger field. On the whole, this
field is not directly aimed at the challenge posed by quantum chaos, where
exponential instability is converted from hindrance to resource. There are a
few works however [28, 29, 30, 31, 32, 33, 34], but none address general
approaches to optimal coherent targeting in quantum chaotic many-body systems.

As alluded to before, a quantum targeting in an isolated many-body dynamical
system may offer some important capabilities. A fully chaotic dynamics
eventually visits all possible states of a system independent of starting
point, hence there is the possibility of placing a system in an otherwise
quite unusual or difficult-to-access state. In a space of enormous volume as a
many-body system would possess, in most cases it would take an exceedingly
long time to arrive at a desired state through an ergodic wandering, however,
heteroclinic motion provides those rare, most rapid paths from particular
starting to final states, thus providing a tremendously accelerated route to a
desired state.

Many-body systems based on bosonic ultracold atoms provide ideal candidates
for implementation of coherent quantum targeting for a number of reasons. In
the main though, there is a broad variety of experiments with great control
over tunable parameters being performed, and they provide an excellent
platform for quantum simulations [35, 36, 37, 38, 39, 40, 41, 42, 43, 44].
However, in addition they possess well defined classical analogs in the form
of a mean-field limit with large particle number as required by the control
protocol. The classical dynamics is governed by the Gross-Pitaevskii equation
or its discretized version on a lattice (in optical lattices the system gets
effectively described by a Bose-Hubbard model [19]); see [45] for a detailed
review. This allows for a phase space formulation in which strongly chaotic
regimes can be identified [46, 47] and where the resulting heteroclinic
dynamics can be explored. Finally, the systems admit both a continuous $U(1)$
and discrete dynamical symmetries that can be taken advantage of in developing
control protocols for the purposes of introducing extremely helpful mappings
between different dimensional lattices and expanding the range of initial
states to which the protocol can be applied. The critical parameters needed to
be controlled in this work are the onsite interaction strengths between the
atoms, which can be magnetically tuned to arbitrary values using Feshbach
resonances [48, 49, 50], and the chemical potentials of the sites.

This paper is structured as follows: in the next section background material
and notation are introduced related to dynamical considerations, aspects of
heteroclinic motion, Bose-Hubbard models, coherent and Fock states resulting
from number-projected coherent states, truncated Wigner approximations, and
finally control protocol considerations. This is followed in Sec. 3 with
coherent quantum targeting applied to one-dimensional, $1D$, lattices. The
realization of the protocol reduces to switching off the interactions and
controlling the chemical potentials of the individual sites in a time-
dependent fashion. Examples of state preparation with non-trivial condensate
phases as well as periodic many-body states are given. Error sources and
fidelity of the protocol is discussed for explicit examples. Section 4 gives
the application to two-dimensional, $2D$, lattices. An example of a lattice of
discrete vortices is given in addition to periodic mean field solutions.
Finally, there is a summary and outlook for further research.

##  2 Background

In order to describe coherent quantum targeting and its application to Bose-
Hubbard models, it is helpful to start with some common background concerning
chaotic dynamics, stability analysis of dynamical systems, Bose-Hubbard
systems, their mean field limit and symmetries, coherent and number projected
coherent states, truncated Wigner approximations, and the protocol for the
optimal coherent quantum targeting itself.

###  2.1 Chaos, heteroclinic motion, and times scales

Within a chaotic volume of phase space, given any arbitrary initial and target
state that respects the restrictions due to constants of the motion,
ergodicity guarantees the existence of a transport pathway connecting the
respective states. In fact, almost all of the trajectories would provide a
connection. However, the time scale for making the connection with a typical
ergodic trajectory would be prohibitively long. Phase space volume grows
exponentially with system size, number of degrees of freedom $L$, and its
random exploration time scale would grow exponentially as well.

It is therefore necessary to seek special trajectories that connect the
initial and final phase points in as short a time as possible. The works of
Lyapunov [51] and Poincaré [10] intimate the solution to this optimization
problem. First, a quantum state provides a coarse grained phase space volume
of Planck’s constant raised to the power of the number of degrees of freedom,
$h^{L}$ [16]. Within such a localized neighborhood of phase space points are
initial conditions for trajectories which flee exponentially rapidly from the
central phase space point’s trajectory. Likewise, considering the target
state’s neighborhood, there exist trajectories that are approaching the
trajectory associated with the central phase space point exponentially
rapidly; i.e., for every positive Lyapunov exponent, there is a negative
exponent. There are invariant sets of phase points associated with the
positive and negative Lyapunov exponents called unstable and stable manifolds,
respectively. The intersection of these manifolds defines heteroclinic
trajectories. They have the combined features of fleeing the initial point’s
trajectory at a maximally exponential rate while simultaneously approaching
the final point’s trajectory at the maximal exponential rate. Heteroclinic
trajectories exist that can be as close as necessary to the initial and final
points.

One of the infinitely many heteroclinic trajectories gives the shortest
dynamical connection time between the initial and final phase space points
depending on the volume $h^{L}$, which determines the practical constraint of
closeness to initial and final phase space points, although $h^{L}$ can be
chosen as small as preferred. Denoting this time as $t_{min}$, the immediate
neighborhood of trajectories that imitate this heteroclinic motion is of the
order of $h^{L}\exp\left(-h_{KS}t_{min}\right)$, where $h_{KS}$ is Kolmogorov-
Sinai entropy [52, 53] or for our purposes, the sum of positive Lyapunov
exponents. For an extensive system, i.e., the mean positive Lyapunov exponent
is constant and $h_{KS}$ grows linearly with the number of degrees of freedom,
this neighborhood vanishes exponentially with system size. Parenthetically,
for the mean field Bose-Hubbard Hamiltonian introduced ahead, Eq. (7), holding
the filling factor constant ($N/L$) as the number of sites $L$ increases leads
to such an extensive system. Thus, locating a suitable solution leads to a
search for an extraordinarily rare trajectory; see Ref. [54] for a general
technique developed precisely for searching exponentially rare solutions.

Of great interest here is how short the time scale $t_{min}$ can be. Assuming
an extensive system, the available phase space volume per degree of freedom
has a geometric mean, ${\cal V}_{1}$, so that the number of quantum coarse
grained cells for the full system is $N_{L}=({\cal V}_{1}/h)^{L}$. There is a
time scale for moving to an adjacent phase space cell, $t_{0}$, which is also
the time for the quantum dynamics to move to the “closest” orthogonal state.
The time scale for an ordinary meandering chaotic trajectory to connect two
random cells is proportional to $N_{L}t_{0}$, which is an extraordinarily long
time. This is completely different from the optimal heteroclinic trajectory.
It is found somewhere inside a cell stretching as $\exp\left(h_{KS}t\right)$.
This grows to the order of $N_{L}$ on a very short time scale given
approximately by the relation

$t_{min}\sim\frac{1}{h_{KS}}\ln N_{L}=\frac{1}{\lambda}\ln\frac{{\cal
V}_{1}}{h}\ ,$ (1)

where $\lambda=h_{KS}/L$ is the average positive Lyapunov exponent. This is
another example of a so-called log-time typically found in chaotic systems.
The quintessential example is the Ehrenfest time at which quantum and
classical dynamics must diverge, which is a log-time [55, 56, 57, 58] for a
chaotic system. Depending on the purposes, the details of a particular logtime
may vary slightly in terms of Lyapunov exponents (maximum or average, etc..),
either classical actions or phase space volumes, but they are all similarly
logarithmic in Planck’s constant.

This rather simplistic logic not only gives a vastly shorter time scale than
$N_{L}t_{0}$ and one logarithmic in Planck’s constant, but surprisingly,
independent of the number of degrees of freedom. Perhaps a more realistic
calculation would lead to some weak system size dependence, after all
dynamical systems are not uniformly hyperbolic, different cells are more or
less difficult to approach, and larger systems occupy larger physical spaces,
but clearly even a more realistic estimate would still lead to the conclusion
that the optimal heteroclinic motion gives something close to an exponential
time scale speed-up compared with $N_{L}t_{0}$ and very weak dependence on
Planck’s constant. This turns out to be the remarkable way in which chaos
becomes a targeting resource. The cost of this speed up is the difficulty of
identifying the requisite exponentially rare trajectory.

Without further analysis of how to locate the optimal heteroclinic trajectory,
the search space dimension is the immense double-the-degrees-of-freedom of the
entire lattice, $2L$, such that any systematic search quickly becomes
numerically challenging if not impossible. Conveniently though, chaotic
systems usually allow for substantial reduction of the search space, as
explicitly demonstrated for a search for heteroclinic trajectories and complex
saddles in a Bose-Hubbard system in [59]: for any phase space point located
inside a chaotic region, the local dynamics can be decomposed into pairs of
stable and unstable manifolds. By disregarding all stable directions and those
corresponding to constants of the motion, any search in a Bose-Hubbard system
can be reduced to an $(L-2)$-dimensional problem. By focusing only on the most
unstable directions an even further reduction is usually possible. Even so, a
systematic search for general large or higher-dimensional systems is typically
futile.

There is a further significant reduction of the search space dimension
possible in cases where symmetries are present in the initial and final
states; details ahead. Indeed, this work focuses on Bose-Hubbard systems
possessing discrete symmetries, periodic boundary conditions, and symmetric
initial conditions leading to a reduced-dimensional symmetry plane of the mean
field dynamics [60]. The reduced-dimensional dynamics can be mapped onto
larger systems with each dimension in the lattice (one, two, or three) being
any integer multiple of the reduced-dimensional value, including infinity. In
this way, the search space remains the reduced-dimensional value independent
of the size or dimension of the lattice. Ahead for illustrative purposes, the
symmetry reduced-dimensional value of the lattice dimension is $L=4$, and thus
an optimized search can be performed in an $L-2=2$-dimensional space no matter
what the entire lattice size is. Thus, it turns out that a simple Monte Carlo
search in the reduced phase space around the initial state’s centroid is
sufficient to determine the suitable heteroclinic pathway, i.e., control
trajectory. It is assumed in this work that the technical problem of
identifying suitable control trajectories is solved and is not discussed
further. Note finally that $t_{min}$ is independent of system size in the case
of relying on discrete symmetries, just as Eq. (1) also gives a system size
independent result, i.e., thus the essential motivation of taking advantage of
discrete symmetries is the reduced search space. Furthermore, the control
scheme can also be applied to tasks in which a certain time for reaching the
target state is given, as long as it is greater than $t_{min}$. For such
larger times, there always exists a corresponding hetereoclinic trajectory
that can be used as optimal transport path.

###  2.2 Bose-Hubbard systems

For many reasons the Bose-Hubbard model turns out to be an ideal bosonic many-
body system with which to investigate optimal coherent quantum control. Not
only is this model experimentally relevant for ultracold atoms in an optical
lattice, but a classical limit can be defined as a mean-field limit of large
occupation numbers and the specific form of the interaction greatly
facilitates the control implementation as discussed in Secs. 2.2 and 2.3
ahead. A convenient property is that a quadrature phase space representation
leads to a Hamiltonian formalism with canonically conjugate real variables.
Furthermore, the control protocol works on a many-body version of minimum
uncertainty states in the form of Glauber coherent states, and crucially all
the results apply equally well to their total number projected states.
Finally, for the control system, the truncated Wigner approximation becomes a
critical analysis tool as its dynamics become essentially exact.

####  2.2.1 Hamiltonian and classical (mean field) limit

A gas of $N$ ultracold atoms subject to a sufficiently strong periodic optical
trapping potential generating a lattice may be well described by the Bose-
Hubbard model [19]; the lattice can be created in one, two, or three
dimensions. Further ahead, $2D$ lattices are treated, but here the Bose-
Hubbard Hamiltonian is only given for a periodic $1D$ lattice of $L$ sites
($L+1\equiv 1$),

$\displaystyle\begin{split}\hat{H}=-J\sum_{j=1}^{L}\big{(}\hat{a}^{\dagger}_{j}\hat{a}_{j+1}+\text{h.c.}\big{)}+\frac{U}{2}\sum_{j=1}^{L}\hat{n}_{j}\left(\hat{n}_{j}-1\right)+\sum_{j=1}^{L}\mu_{j}\hat{n}_{j}\
.\end{split}$ (2)

It energetically accounts for the possible hopping of bosons to neighboring
sites, governed by the kinetic hopping parameter $J$, as well as on-site
interactions between the bosons given by an interaction strength $U$, which is
related to the s-wave scattering length of the two-body collisions between
atoms. The individual sites can additionally be offset by chemical potentials
$\mu_{j}$; see Fig. 2 for a simplified rendition of a $1D$ ring.

interactionhoppingchem. pot. Figure 2: Schematic representation of a $1D$
Bose-Hubbard ring. Cold atoms trapped in a $1D$ periodic lattice, here
consisting of four sites. The model includes nearest-neighbor hopping, on-site
interactions between atoms, as well as energetic offsets for each site given
by individual chemical potentials.

As there are two preserved quantities, energy and total particle number, the
dimer ($L=2$) is always integrable, as are the limiting cases of $U=0$ or
$J=0$, where the absence of an interaction leads to a superfluid ground state
and the absence of hopping gives a Mott insulator phase. For $L>2$ and in the
parameter range $J\sim UN$, the quantum system is known to exhibit quantum
chaos [46, 47], which is the dynamical regime of current interest.

A mean field limit, which depends on an effective Planck constant, can be
defined equal to the inverse filling factor, $\hbar_{\rm eff}=L/N$. As
$\hbar_{\rm eff}\rightarrow 0$, the mean field approximation gets better and
better. A unique Hamiltonian (classical) dynamics arises for the mean field if
the ratio

$\gamma=\frac{1}{\hbar_{\rm eff}}\frac{U}{J}\,,$ (3)

is held fixed.

It is possible to introduce Hermitian operators in terms of $\hat{a}_{j}$ and
$\hat{a}^{\dagger}_{j}$ known as quadratures that obey position-momentum-like
commutation relations,

$\hat{a}_{j}=\frac{(\hat{q}_{j}+i\hat{p}_{j})}{\sqrt{2\hbar_{\rm
eff}}},\quad\hat{a}^{\dagger}_{j}=\frac{(\hat{q}_{j}-i\hat{p}_{j})}{\sqrt{2\hbar_{\rm
eff}}},\quad\left[\hat{q}_{j},\hat{p}_{j}\right]=i\hbar_{\text{eff}}\ ,$ (4)

which re-expresses the Bose-Hubbard Hamiltonian as

$\displaystyle\begin{split}\hat{H}=-\frac{J}{\hbar_{\text{eff}}}\sum_{j=1}^{L}\left(\hat{q}_{j}\hat{q}_{j+1}+\hat{p}_{j}\hat{p}_{j+1}\right)+\frac{U}{2\hbar_{\text{eff}}^{2}}\sum_{j=1}^{L}\left(\frac{\hat{q}_{j}^{2}+\hat{p}_{j}^{2}}{2}\right)^{2}+\sum_{j=1}^{L}\frac{\mu_{j}-U}{\hbar_{\text{eff}}}\left(\frac{\hat{q}_{j}^{2}+\hat{p}_{j}^{2}}{2}\right)\
.\end{split}$ (5)

Here an irrelevant constant energy offset has been dropped (which is also
effectively an $O(\hbar_{\text{eff}}^{2})$ correction to the spectrum and, in
any case, beyond semiclassical argumentation used ahead). Multiplying both
sides of the equation by $\hbar_{\text{eff}}/J$ leads to a modified
Hamiltonian,

$\displaystyle\begin{split}\hat{H}^{\prime}=-\sum_{j=1}^{L}\left(\hat{q}_{j}\hat{q}_{j+1}+\hat{p}_{j}\hat{p}_{j+1}\right)+\frac{\gamma}{2}\sum_{j=1}^{L}\left(\frac{\hat{q}_{j}^{2}+\hat{p}_{j}^{2}}{2}\right)^{2}+\sum_{j=1}^{L}\mu_{j}^{\prime}\left(\frac{\hat{q}_{j}^{2}+\hat{p}_{j}^{2}}{2}\right)\
,\end{split}$ (6)

where $\mu_{j}^{\prime}=(\mu_{j}-U)/J$. This Hamiltonian is in a sufficiently
symmetric operator form for a full time-dependent WKB approximation [61],
i.e., good to $O(\hbar_{\text{eff}}^{2})$, just through the association of
$\hat{q}\rightarrow q$ and $\hat{p}\rightarrow p$. The classical analog
Hamiltonian (mean field limit) is given by

$\displaystyle\mathcal{H}\left(\vec{q},\vec{p}\right)=-\sum_{j=1}^{L}\left(q_{j}q_{j+1}+p_{j+1}p_{j}\right)+\frac{\gamma}{2}\sum_{j=1}^{L}\left(\frac{q_{j}^{2}+p_{j}^{2}}{2}\right)^{2}+\sum_{j=1}^{L}\mu_{j}^{\prime}\left(\frac{q_{j}^{2}+p_{j}^{2}}{2}\right)\
,$ (7)

where the degree of chaos is determined by the value of $\gamma$. Typically,
the strongest chaos is approximately in the range $1<\gamma<3$; see [46, 59].
The dynamics of the quadratures is thus governed by Hamilton’s equation of
motion

$\displaystyle\dot{\vec{x}}=\mathbb{\Sigma}\,\vec{\nabla}_{\vec{x}}\mathcal{H}\
,\quad\text{ with }\vec{x}=(\vec{q},\vec{p})\text{ and
}\mathbb{\Sigma}=\begin{pmatrix}0&\mathbb{1}\\\ -\mathbb{1}&0\end{pmatrix}\ ,$
(8)

in a 2$L$-dimensional phase space. The classical dynamics conserves energy and
total number of particles along each trajectory, i.e., each mean-field
solution. For any targeting application using the Bose-Hubbard model this
implies that the space of possible targeting states is restricted to those
associated to phase space points sharing the same energy and particle number
surface as the corresponding initial state. Choosing the initial chemical
potentials $\mu_{j}^{\prime}$ of the individual sites however allows for quite
a bit of freedom in putting initial and target states on the same energy
surface, so the constraint posed by energy conservation can be largely evaded.

####  2.2.2 Localized states

Similar to [11, 12], the protocol for Bose-Hubbard systems is designed to
control states that are initially localized in the quadrature phase space.
Glauber coherent states provide a natural and optimal choice with which to
embark as they are minimum uncertainty states and thus the most classical
possible. Defined as eigenstates of the annihilation operator

$\hat{\vec{a}}|{\vec{\phi}}\rangle_{\rm
coh}={\vec{\phi}}|{\vec{\phi}}\rangle_{\rm coh}\ ,$ (9)

they can be written explicitly as

$\displaystyle|{\vec{\phi}}\rangle_{\rm
coh}=\exp{-\frac{\lVert{{\vec{\phi}}}\rVert^{2}}{2}+\vec{\phi}\cdot\hat{\vec{a}}^{\dagger}}\ket{0}\
.$ (10)

Their phase space representation is given by the Wigner function and is well
known to be a simple symmetric $2L$-dimensional Gaussian wave packet

$\mathcal{W}\left(\vec{q},\vec{p}\right)=\left(\pi\hbar_{\text{eff}}\right)^{-L}\exp{-\frac{\left(\vec{q}-\vec{q}_{\alpha}\right)^{2}}{\hbar_{\text{eff}}}-\frac{\left(\vec{p}-\vec{p}_{\alpha}\right)^{2}}{\hbar_{\text{eff}}}}\
,$ (11)

centered around the mean-field point
$\vec{\phi}_{\alpha}=(\vec{q}_{\alpha}+i\vec{p}_{\alpha})/\sqrt{2\hbar_{\text{eff}}}$.
The coherent state variance per degree of freedom is given by
$\sigma^{2}=\hbar_{\text{eff}}/2$, which for a single degree of freedom has
the $2\sigma$-contour enclosing a volume $2\pi\hbar_{\text{eff}}$. More
generally, the volume per quantum state is $(2\pi\hbar_{\text{eff}})^{L}$
[16], as mentioned above.

The Bose-Hubbard system possesses a $U(1)$ dynamical symmetry, which in mean
field quadratures equates to an $SO(2)$ symmetry. Each trajectory maps
perfectly onto another by simultaneously rotating each pair $(q_{j},p_{j})$
through any given angle $0\leq\theta<2\pi$. This symmetry has the consequence
that it implies a generalization of the control protocol in a straightforward
manner.

The same scheme as for a Glauber coherent state also works for a number
projected coherent state

$\displaystyle|{\vec{\phi}}\rangle_{\rm
proj}^{N}=\frac{1}{\sqrt{\lVert{{\vec{\phi}}}\rVert^{2N}N!}}\big{(}\vec{\phi}\cdot\hat{\vec{a}}^{\dagger}\big{)}^{N}\ket{0}\
,$ (12)

which contains a fixed number of particles, and is a Fock state in some basis;
states of this kind can be constructed [62]. This additional class of states
is of great physical significance just by virture of fixing the particle
number, but it also contains interesting cases such as
$\phi^{\alpha}_{j}=\exp\left(2\pi i\alpha j/L\right)$ with $\alpha=0$
corresponding to the non-interacting ground state of the superfluid phase and
$\alpha=1$ being the first excited state of the momentum operator. Another
class of states that can be represented as number projected coherent states
are all Fock states concentrated on a single site $k$, e.g.,
$\phi_{j}=\sqrt{L}\delta_{j,k}$.

####  2.2.3 Truncated Wigner Approximation

Due to their simpler representation, the discussion follows coherent states,
however all results are equally applicable to the number projected
counterparts. In quadratures, the propagation of the state is given by the
Moyal bracket [63]

$\displaystyle\begin{split}\frac{d}{dt}\mathcal{W}\left(\vec{q},\vec{p}\right)&=\left\\{\mathcal{H}\left(\vec{q},\vec{p}\right),\mathcal{W}\left(\vec{q},\vec{p}\right)\right\\}_{\text{MB}}=\frac{2}{\hbar_{\text{eff}}}\mathcal{H}\left(\vec{q},\vec{p}\right)\sin\left[\frac{\hbar_{\text{eff}}\overleftrightarrow{\Lambda}}{2}\right]\mathcal{W}\left(\vec{q},\vec{p}\right)\
,\end{split}$ (13)

with
$\overleftrightarrow{\Lambda}=\frac{\overleftarrow{\partial}}{\partial\vec{q}}\frac{\overrightarrow{\partial}}{\partial\vec{p}}-\frac{\overleftarrow{\partial}}{\partial\vec{p}}\frac{\overrightarrow{\partial}}{\partial\vec{q}}$
denoting the symplectic operator, which acts like the Poisson bracket. An
important consequence of Eq. (13) is that for any quadratic Hamiltonian all
higher order terms vanish, which dovetails perfectly with the control protocol
introduced ahead.

Dropping third and higher order terms in $\hbar_{\text{eff}}$ yields the
truncated Wigner approximation (TWA) [64, 65] in which quantum zero-point
fluctuations are encoded in the Wigner transforms, but otherwise quantum
interference between mean field paths is dropped. The distribution of initial
conditions then evolves under a classical Liouville equation

$\frac{d}{dt}\mathcal{W}\left(\vec{q},\vec{p}\right)=\left\\{\mathcal{H}\left(\vec{q},\vec{p}\right),\mathcal{W}\left(\vec{q},\vec{p}\right)\right\\}\
.$ (14)

The Gaussian form for coherent states of Eq. (11) is straightforward to
implement within the TWA. This renders the TWA an excellent analysis tool for
this work.

####  2.2.4 Density Wave

Density waves corresponding to coherent states centered on mean-fields of the
form $\vec{\phi}=\sqrt{2/\hbar_{\text{eff}}}\left(1,0,1,0...\right)$,
alternating sites for higher-dimensional lattices, are special solutions to
Hamilton’s equations given by initially populating every second site equally,
while leaving the rest unoccupied. The mean field solutions, using initial
conditions based on their centroid values, can be mapped to and from the dimer
[60], which is an integrable system. The resulting trajectories exhibit
periodic population inversion for $\gamma<4$ with a period
$\tau_{DW}(\gamma)$, see Fig. 3(a). For $\gamma>4$, the system is in a self-
trapping regime where the solution is still periodic, but only a partial
population imbalance is achieved. If $\vec{\phi}$ lies on the symmetry plane
separatrix of the dimer, $\gamma=4$, the period diverges as shown in Fig. 3.

(a)

(b)

Figure 3: Properties of a density wave. A density wave initially centered
around $\vec{\phi}=\sqrt{2/\hbar_{\text{eff}}}\left(1,0,1,0\right)$ exhibits
partial or full periodic population inversion with a period $\tau_{DW}$ as
illustrated in (a). How the period depends on dynamical parameter $\gamma$,
Eq. (3), of the system is shown in (b). Except for $\gamma=4$ where the period
diverges since $\vec{\phi}$ lies on a separatrix of the mean field dimer,
$\tau_{DW}$ provides a well defined time scale of the system. The time scale
is marked for all the specific trajectories used throughout the work.

Due to these special properties, it serves as a convenient initial state for
the applications discussed throughout this work, even though theoretically any
other state with a centroid located in a chaotic phase space region could be
used. This makes $\tau_{DW}$ a natural choice for a dynamical time scale, as
long as the systems considered are sufficiently removed from the singularity
at $\gamma=4$, as is done ahead. The corresponding time scales $\tau_{DW}$ of
this work’s specific cases are depicted in Fig. 3(b).

###  2.3 Control Hamiltonian

Given that a (presumably heteroclinic) control trajectory
$\left(\vec{q}(t),\vec{p}(t)\right)_{c}$ connecting the mean-field centers of
the initial and target state at time $t_{min}$, respectively, is identified,
the initial goal is to guide the coherent state along this solution.

$t=0$

$t=0.5\tau_{\text{DW}}$

$t=\tau_{\text{DW}}$

$t=4\tau_{\text{DW}}$

Figure 4: Spreading of a localized coherent state in a Bose-Hubbard system. A
minimum uncertainty coherent state (here for illustration chosen to be
$\hbar_{\text{eff}}=1/10$), initially centered on the density wave
$\vec{\phi}_{\alpha}=\sqrt{2/\hbar_{\text{eff}}}\left(1,0,1,0\right)$
(indicated as dot in the upper panel) is propagated in the Bose-Hubbard system
($\gamma\simeq 1.965$) using the TWA. With the scaled dynamics that is used,
the mean-field corresponds to the phase space point
$\vec{q}_{\alpha}=(2,0,2,0)$ and $\vec{p}_{\alpha}=(0,0,0,0)$, independent of
$\hbar_{\text{eff}}$. The state, whose volume is for $t=0$ well contained in
the $2\sigma$-contour (white dashed circle), quickly spreads, already being
far from its initial Gaussian shape after one period of the density wave. At
$t=4\tau_{DW}$ the state is fully equilibrated.

Due to the quartic polynomial in quadratures present in the interaction term
of Eq. (7), the initially localized state spreads exponentially in time as the
dynamics are strongly chaotic in the dynamical regimes of interest. This
spreading is so rapid that without some control mechanism, a quantum targeting
is impossible beyond an exceedingly short logarithmic time scale in
$\hbar_{\text{eff}}$. Figure 4 demonstrates this fast equilibration for the
example of a coherent state initially centered on a density wave. The full
time evolution is provided in Animation 1 in the supplementary material.

It is theoretically possible to proceed as in [11] and apply an additional
unitary transformation to the system sufficiently often that is designed to
counteract the spreading. However, this generally introduces terms in the time
evolution that are not particle number conserving and thus challenging, if not
impossible, with respect to experimental realization. There exists a much
simpler approach developed in [12] based on introducing a new time-dependent
simulation control Hamiltonian, $\hat{H}_{c}$, that is uniquely designed for
some particular targeting problem with fixed initial and final quadrature
conditions; note that a minimal time related to heteroclinic motion could be
the desired goal or any specific time longer than that.

Consider Hamilton’s equations for each component using Eq. (7),

$\displaystyle\begin{split}&\dot{q}_{j}=\frac{\partial{\mathcal{H}}\left(\vec{q},\vec{p}\right)}{\partial{p_{j}}}\biggr{\rvert}_{\vec{q}(t),\vec{p}(t)}=\left[\mu_{j}^{\prime}+\gamma\left(\frac{q_{j}^{2}+p_{j}^{2}}{2}\right)\right]p_{j}-p_{j+1}-p_{j-1}\,,\\\
&\dot{p}_{j}=-\frac{\partial\mathcal{H}\left(\vec{q},\vec{p}\right)}{\partial{q_{j}}}\biggr{\rvert}_{\vec{q}(t),\vec{p}(t)}=-\left[\mu_{j}^{\prime}+\gamma\left(\frac{q_{j}^{2}+p_{j}^{2}}{2}\right)\right]q_{j}+q_{j+1}+q_{j-1}\
.\end{split}$ (15)

Evaluated using the control trajectory,
$\left(\vec{q}(t),\vec{p}(t)\right)_{c}$, the terms in large square
parentheses can be replaced by time-dependent values that are no longer
functions of $\\{q_{j},p_{j}\\}$. With that replacement, it is straightforward
to see that a control Hamiltonian can be created with the interaction term
switched off given by

$\mathcal{H}_{c}\left(\vec{q},\vec{p},t\right)=-\sum_{j=1}^{L}\left(q_{j}q_{j+1}+p_{j+1}p_{j}\right)+\sum_{j=1}^{L}\mu_{j}(t)\left(\frac{q_{j}^{2}+p_{j}^{2}}{2}\right)\
,$ (16)

which is in essence a harmonic oscillator with chemical potentials that are
time-dependent functions of the control trajectory’s site occupancies,

$\mu_{j}(t)=\mu_{j}^{\prime}+\gamma\left(\frac{q_{j}^{2}(t)+p_{j}^{2}(t)}{2}\right)=\mu_{j}^{\prime}+\gamma
n_{j}(t)\ .$ (17)

The control Hamiltonian makes use of the chemical potentials $\mu_{j}(t)$ as
purely time-dependent functions evaluated along the chosen control
(heteroclinic) trajectory $(\vec{q}(t),\vec{p}(t))_{c}$, which effectively
replaces the originally cubic terms in quadratures with time-dependent linear
terms. In this way, the control trajectory and only the control trajectory,
$(\vec{q}(t),\vec{p}(t))_{c}$ is a solution of Hamilton’s equations for both
Hamiltonians $\mathcal{H}\left(\vec{q},\vec{p}\right)$ and
$\mathcal{H}_{c}\left(\vec{q},\vec{p},t\right)$. In addition,
$\mathcal{H}_{c}\left(\vec{q},\vec{p},t\right)$ having a quadratic form leads
to the TWA being exact for its quantized version,

$\displaystyle\begin{split}\hat{H}_{\text{c}}(t)=-\sum_{j=1}^{L}&\left(\hat{a}^{\dagger}_{j}\hat{a}_{j+1}+\hat{a}^{\dagger}_{j+1}\hat{a}_{j}-\mu_{j}(t)\hat{n}_{j}\right)\
,\end{split}$ (18)

and there is no longer any spreading in the neighborhood of
$(\vec{q}(t),\vec{p}(t))_{c}$, or for that matter, anywhere else. Although,
$(\vec{q}(t),\vec{p}(t))_{c}$ is an exponentially unstable solution of
$\mathcal{H}\left(\vec{q},\vec{p}\right)$, it is a stable solution of
$\mathcal{H}_{c}\left(\vec{q},\vec{p},t\right)$.

Up to this point there is a broad freedom in the choice of a particular
chaotic Hamiltonian that may serve as a black box to produce a desired control
solution. The apparently natural choice of a Bose-Hubbard type of system,
however, points to a deeper connection between the dispersionless dynamics
generated by Eq. (18) and the precise form of the interaction (quartic) term
in Eq. (7). In addition, this choice is quite favorable from the point of view
of a practical implementation, since among the infinity of nonlinear systems
that one might consider for producing a control trajectory, a key selection
criteria is the desire to have a minimal alteration of a realistic
experimental set up. The Bose-Hubbard Hamiltionian, describing a broad range
of self-interacting bosonic systems, satisfies this condition thanks to the
technical possibilities to tune the interaction [48, 49, 50] and chemical
potential parameters, and to implement Eq (18).

There is however another deep conceptual reason for the special connection
between the control and Bose-Hubbard Hamiltonians. It relies on the coherent
representation of the time evolution operator for Bose-Hubbard systems

$K(\vec{\phi}_{\beta},\vec{\phi}_{\alpha};t)=\langle\vec{\phi}_{\beta}\left|{\rm
e}^{-\frac{i}{\hbar}\hat{H}t}\right|\vec{\phi}_{\alpha}\rangle,$ (19)

with $\hat{H}$ given in Eq. (2), in terms of auxiliary fields

$K^{\rm sc}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha};t)=\sum_{\gamma}\frac{{\cal
D}_{\gamma}(\vec{\phi}_{\alpha},\vec{\phi}_{\beta},t)}{Z({\bf v})}{\rm
e}^{\frac{i}{2}\int_{0}^{t}ds\vec{\sigma}_{\gamma}(s){\bf
v}^{-1}\vec{\sigma}_{\gamma}(s)}\times
K_{\vec{\sigma}_{\gamma}}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha};t)$ (20)

valid in the limit of large mean occupations $N\gg 1$. Here
$K_{\vec{\sigma}_{\gamma}}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha};t)$ is the
transition amplitude between initial $|\vec{\phi}_{\alpha}\rangle$ and target
$|\vec{\phi}_{\alpha}\rangle$ coherent states defined by the control
Hamiltonian, Eq. (18), with control field
$\vec{\mu}(t)=\vec{\sigma}_{\gamma}(t)$. As shown in A, the total transition
amplitude for the full interacting problem is given as a coherent sum,
weighted by the factors $\frac{{\cal
D}_{\gamma}(\vec{\phi}_{\alpha},\vec{\phi}_{\beta},t)}{Z({\bf v})}{\rm
e}^{\frac{i}{2}\int_{0}^{t}ds\vec{\sigma}_{\gamma}(s){\bf
v}^{-1}\vec{\sigma}_{\gamma}(s)}$, over the set of control fields
$\vec{\sigma}_{\gamma}(s)$ obtained by
$\sigma_{\gamma,j}(s)=\phi_{\gamma,j}(s)\phi^{*}_{\gamma,j}(s)$ from the
solutions of the mean field boundary problem $\vec{\phi}_{\gamma}(s)$. This
representation leads to a natural interpretation of the control Hamiltonian as
a contribution to the full interacting problem, that is selected based on
further physical considerations. One such physical consideration, relevant for
the approach discussed here, is the real character of the control fields.
Since, generically, the solutions of the mean field problem linking initial
and target coherent states admit only solutions with complexified quadratures
and the corresponding complexified auxiliarly fields, one instead looks for
nearby states admiting one solution with real quadratures, making the control
Hamiltonian Hermitian.

Furthermore, by its very construction, the transition amplitude for a given
auxiliary field is dispersionless due to the quadratic form of the control
Hamiltonian. Specifically, this linearity of the control evolution
$\hat{H}_{\text{c}}(t)$, implies that the coherent state
$|{\vec{\phi}(t)}\rangle_{\rm coh}$ (possibly number projected) centered on
the mean-field solution used to drive the system satisfies the time-dependent
Schrödinger equation, i.e.

$\displaystyle i\frac{\partial}{\partial t}|{\vec{\phi}(t)}\rangle_{\rm
coh}=\hat{H}_{\text{c}}(t)|{\vec{\phi}(t)}\rangle_{\rm coh}\ .$ (21)

The explicit proof is provided in B. Given the simple form of Eq. (18), the
actual realization of the protocol requires controlling the chemical
potentials of each individual site in a time-dependent fashion, and using
Feshbach resonance to tune the interaction strength to a sufficiently
vanishing value [48]. Some aspects of the protocol robustness to imperfections
are addressed in Sec. 3.2.

##  3 Coherent quantum targeting in 1D lattices

Here, a few representative example applications of the protocol on a 1D
periodic lattice are presented. Preparing states with well-defined phase
relations between the condensates on different sites is covered in Sec. 3.1
and the creation of non-trivial periodic evolution in Sec. 3.3. As it turns
out, even though the idea of placing quantum states on classical trajectories
is semiclassical in nature, the control system is harmonic and thus behaves
purely classical. This implies that the method should work perfectly well for
both the large particle number limit, as well as in dilute quantum regime.
However, in this treatment the small shifts between initial and target
quadrature centers to the initial and final points of the control trajectory,
$(\vec{q}(t),\vec{p}(t))_{c}$, are not taken into account. The resulting
effects as well as imperfections in making the residual interactions vanish
both introduce an $\hbar_{\text{eff}}$ dependency to the fidelity of the
protocol and are addressed in Sec. 3.2. Both the TWA and full quantum
simulations of the time evolution are given.

###  3.1 Targeting

(a)

${\pi/2}$${\pi/2}$${\pi/2}$${\pi/2}$ (b)

Figure 5: Schematic of initial and target quantum state. Raising the
unoccupied sites energetically by choosing $\vec{\mu}=(0,\gamma,0,\gamma)$
puts the density wave, (a), and excited momentum eigenstate (as target), (b),
on the same classical energy surface. This opens up the possibility of finding
a control trajectory, connecting both mean fields in quadrature phase space.
The corresponding initial coherent quantum state can be guided along the
control trajectory into the target quantum state by driving the control system
with the associated time-dependent occupation numbers.

As discussed in Sec. 2.1, the great advantage of using heteroclinic pathways
of the chaotic mean-field limit as guiding control trajectories is their
exponential speed-up of connecting any two points in phase space provided by
ergodicity. In principle, any state of the form of Eqs. (10) or (12) can be
prepared by targeting the corresponding mean-field starting from an initial
field sharing the same energy and particle number surfaces. This includes
states where the condensates on each site carry well defined, non-trivial
phases [50]. As a proof of principle, the targeted many-body state is chosen
to have the bosons occupying the first excited momentum eigenstate on a four-
site ring, which corresponds to the homogeneous mean-field
$\vec{\phi}_{\beta}=\sqrt{1/\hbar_{\text{eff}}}(1,i,-1,-i)$ with
$\pi/2$-phases between neighboring sites.

As already discussed, the protocol is initialized with a coherent state
centered on a density wave
$\vec{\phi}_{\alpha}=\sqrt{2/\hbar_{\text{eff}}}(1,0,1,0)$. By choosing the
chemical potential offset accordingly, i.e.
$\vec{\mu}^{\prime}=(0,\gamma,0,\gamma)$, both mean-fields share the same
classical energy surface for every value of $\gamma$ and can be connected with
a transport trajectory, see Fig. 5. Figure 6 shows a heteroclinic trajectory
with a control time of less then five oscillations of the density wave
$\tau_{DW}$ (Fig. 3(b)), which is quite short. If greater precision is
required, then ergodicity and the nature of unstable and stable manifolds
guarantees the existence of heteroclinic trajectories closer to the initial
and target mean-fields at the expense of a longer control time.

Figure 6: Control Trajectory to target excited momentum eigenstate. The
trajectory connecting the centroids of the initial density wave
$\vec{\phi}_{\alpha}=\sqrt{2/\hbar_{\text{eff}}}(1,0,1,0)$ (red dot) with the
final target state $\vec{\phi}_{\beta}=\sqrt{1/\hbar_{\text{eff}}}(1,i,-1,-i)$
(red cross) is shown for the different sites in the left panel ($\gamma\simeq
1.965$). The respective $1\sigma$-contours for the initial and target coherent
state’s Wigner densities, Eq. (11), are depicted centered on their respective
quadrature centroids with $\hbar_{\text{eff}}=1$ (black dashed circles ) and
$\hbar_{\text{eff}}=1/4$ (gray dashed circles). Generally, the area inside the
circles shrinks proportionally to $\hbar_{\text{eff}}$. The time-dependent
control chemical potentials shown in the right panel are determined by the
respective site occupations by Eq. (17) and the offsets of $\mu_{j}^{\prime}$
given in the text in Sec. 3.1. Initially, somewhat similar to the integrable
behavior of the density wave, the chaotic nature of the control trajectory
quickly dominates and it reaches the target at $t=4.9\tau_{DW}$. Additionally,
the expectation values of the occupations are shown for the quantum
propagation of the coherent state evolving under the original Bose-Hubbard
Hamiltonian (uncontrolled, grey dotted), which deviate rapidly from the
classical solution. Under the control Hamiltonian (controlled, black dotted)
the expectation values shadow the control trajectory, as they must.

The initial state $|{\vec{\phi}_{\alpha}}\rangle_{\rm coh}$ then evolves under
the control Hamiltonian, Eq. (18), by driving the system with the time
dependent chemical potentials that are determined by the respective site
occupations shown in the right panel of Fig. 6. The time-dependent site
chemical potentials follow using Eq. (17) and the definition of
$\mu_{j}^{\prime}$ given in the text immediately after Eq. (6). A TWA
simulation is used to calculate the propagation of the state, being exact due
to the quadratic nature of the control Hamiltonian.

Figure 7 shows the initial and final state of the protocol. The full time
evolution is provided in Animation 2 of the supplementary material. In
contrast to the spreading observed in Fig. 4, the coherent state keeps its
minimum uncertainty form throughout the entire propagation, successfully
arriving close the target mean-field. For completeness, a full quantum
simulation of the time evolution is performed. It gives the same results.
Figure 6 shows that the expectation value of the occupations follow the
classical mean-field solution perfectly in the control system, whereas they
almost immediately deviate when evolving in the original interacting Bose-
Hubbard system; see Fig. 4.

$t=0$

$t=4.9\tau_{\text{DW}}$

Figure 7: Coherent targeting along chaotic control trajectory. A coherent
state (here for illustration $\hbar_{\text{eff}}=1/10$ was chosen) is
initialized centered on the density wave
$\vec{\phi}_{\alpha}=\sqrt{2/\hbar_{\text{eff}}}(1,0,1,0)$ (white dot). The
inital phase space point is connected to the target
$\vec{\phi}_{\beta}=\sqrt{1/\hbar_{\text{eff}}}(1,i,-1,-i)$ (red cross)
through a control trajectory found in the Bose-Hubbard mean-field limit
(indicated as blue line). The $2\sigma$-curves of the initial and target state
are represented with white dashed circles. A TWA simulation shows that driving
the control system with the occupations of the control trajectory results in
the coherent state evolving along the predesigned path, keeping its minimum
uncertainty form. After the control time $t_{c}=4.9\tau_{\text{DW}}$ the
propagated state lands well inside the volume of the target state.

###  3.2 Protocol robustness

Three mechanisms that might potentially degrade the quality of a particular
application of coherent quantum targeting are: i) the proximity of initial and
target state centroids from initial and final conditions of the control
trajectory, respectively, ii) the precision with which the interactions can be
made to vanish, and iii) the faithfulness with which the time-dependencies of
the set of $\\{\mu_{j}(t)\\}$ match the occupancies of the control trajectory.
An analysis of the latter item is presumably model dependent, and amongst
other issues, would depend on the type of noise inherent in creating the
$\\{\mu_{j}(t)\\}$, and is not discussed further here. The first two items are
discussed below.

####  3.2.1 Proximity

As mentioned at the end of Sec. 2 and proved in B, the specific quadratic
nature of $\hat{H}_{c}$ guarantees that a coherent state initially (exactly)
centered on the control trajectory used to drive the system follows this
trajectory perfectly without any spreading, and satisfies the associated time-
dependent Schrödinger equation. In addition, any coherent state centered
anywhere remains a coherent state for all times, only they do not follow the
control trajectory. In general however, as mentioned at the beginning of this
section, the initial and final conditions of any control trajectory, Sec. 3.1,
are slightly shifted from the initial and desired target states centroids,
respectively. Letting $t_{c}$ denote the time at which the control
(heteroclinic) trajectory arrives closest to the target point, the shifts can
be expressed as the norms

$\delta x_{\text{initial}}=||\vec{x}_{\alpha}-\vec{x}(0)||,\qquad\delta
x_{\text{target}}=||\vec{x}_{\beta}-\vec{x}(t_{c})||\ ,$ (22)

resulting in a non-perfect overlap between the target state
$|{\vec{\phi}_{\beta}}\rangle_{\rm coh}$ and the actual time evolved final
state

$|{\vec{\phi}_{\alpha}(t_{c})}\rangle_{\rm
coh}\coloneqq\exp{-\frac{i}{\hbar}\int\limits_{0}^{t_{c}}dt\hat{H}_{\text{c}}\left(t\right)}|{\vec{\phi}_{\alpha}}\rangle_{\rm
coh}\ .$ (23)

The shifts are minimized in the search for a control trajectory while
simultaneously keeping the protocol time $t_{c}$ near its minimum. Since this
distance only has a meaning relative to the volume of the coherent state, an
$\hbar_{\text{eff}}$ dependency gets introduced into the fidelity
$\mathcal{F}$ of the control process. States with larger filling factors and
thus smaller volume are more affected by the shifts as shown in Fig. 8. Based
on the knowledge of $\delta x_{\text{initial}}$ and $\delta x_{\text{target}}$
bounds can be derived for the fidelity of the protocol:

$\mathcal{F}_{\text{min}}=\exp{-\frac{\left(\delta x_{\text{initial}}+\delta
x_{\text{target}}\right)^{2}}{2\hbar_{\text{eff}}}},\quad\mathcal{F}_{\text{max}}=\exp{-\frac{\left(\delta
x_{\text{initial}}-\delta x_{\text{target}}\right)^{2}}{2\hbar_{\text{eff}}}}\
.$ (24)

The exact overlap is given as

$\mathcal{F}=\left|{}_{\text{coh}}\bra{\vec{\phi}_{\beta}}e^{-\frac{i}{\hbar}\int\limits_{0}^{t_{c}}dt\hat{H}_{\text{c}}\left(t\right)}\ket{\vec{\phi}_{\alpha}}_{\text{coh}}\right|^{2}=\exp{-\frac{\delta\vec{x}^{2}}{2\hbar_{\text{eff}}}}\
,$ (25)

where $\delta\vec{x}$ is the distance between the centroids of the target and
final state. It could be derived using an adapted version of Heller’s
linearized wave packet dynamics [66], which would be rather involved due to
the dependencies on the stability matrix. Here, it is simpler to extract
$\delta\vec{x}$ numerically. Interestingly, even though the control method is
semiclassical in nature, the harmonic nature of the Hamiltonian causes the
protocol to have almost perfect fidelity in the dilute quantum limit and then
fall off for states with larger filling factor, i.e., inverse
$\hbar_{\text{eff}}$.

Figure 8: Fidelity of the targeting protocol. The protocol depends on the
specific control trajectory used for targeting. The initial and final distance
from the particular trajectory to the centroids of the respective states
introduce a implicit $\hbar_{\text{eff}}$ dependence into the fidelity. Upper
and lower bounds depending on those distances are depicted. The exact fidelity
(black curve) depends on the distance of the final states centroid to the
actual target state and can be obtained numerically. The results obtained from
both TWA and quantum simulations match the fidelity curve. Note that the
fidelity would remain unity for all $N$ if shift operators were additionally
used to translate the initial and final coherent states.

####  3.2.2 Residual interactions

The second mechanism mentioned above that could potentially degrade the
coherent quantum targeting is related to switching off the on-site
interactions. Despite the fact that there can be great experimental control
over interactions [48], imagine nevertheless that the strength of the
interaction cannot be tuned to vanish precisely, and there is a small
perturbative residual interaction, which might also have some slow or weak
time dependence. Also, imagine that the residual strength can be measured or
known by some means. This would add a perturbative term to the control
Hamiltonian, Eq. (18), as follows:

$\displaystyle\hat{H}_{\rm
c}^{\rm\epsilon}\left(t\right)=-\sum_{j=1}^{L}\left(\hat{a}^{\dagger}_{j}\hat{a}_{j+1}+\hat{a}^{\dagger}_{j+1}\hat{a}_{j}-\mu_{j}(t)\hat{n}_{j}\right)+\frac{\epsilon(t)}{2}\sum_{j=1}^{L}\hat{n}_{j}\left(\hat{n}_{j}-1\right)\
.$ (26)

The effect of such residual interactions would be twofold: the interactions
would again lead to a deformation of the initially minimum uncertainty
coherent states, and constructing the control Hamiltonian with
$(\vec{q}(t),\vec{p}(t))_{c}$ would result in missing the desired target.

Figure 9: Effect of constant weak residual interactions. The trajectory no
longer ends near the target and is no longer perfectly stable, which distorts
the contours of the Wigner transform. Shown for only the first site, the final
point of the control trajectory ends at the red cross, which is increasingly
far from the target (orange cross) as $\epsilon$ increases. The spreading also
increases with $\epsilon$ as illustrated for the $1\sigma$-contour of the
propagated Wigner density for $\hbar_{\text{eff}}=1$ (black) and
$\hbar_{\text{eff}}=4$ (gray). For very small $\epsilon$ $(=0.01)$, both
effects remain negligible. However for $\epsilon=0.1$, the control trajectory
completely overshoots the target and there is a strong deformation of the
coherent state. Accounting for $\epsilon$ in the $\mu_{j}(t)$, the overshoot
is corrected, however the spreading remains, as shown in the most right panel.
The different state contour areas are a result of the depicted
$1\sigma$-contour being only a cut through a higher-dimensional object.

Both effects are illustrated in Fig. 9 for a constant residual interaction
$\epsilon(t)=\epsilon$, as well as the improvement using the appropriately
modified $\\{\mu_{j}(t)\\}$. If the behavior of $\epsilon(t)$ is known, it is
quite simple to take into account. Consider Hamiltonian equations, just as
before in Eq. (15), which were used to determine the $\\{\mu_{j}(t)\\}$,
except with the inclusion of the $\epsilon$ term. Now resolve the equations as
before, the result is the replacement,

$\mu_{j}(t)\Longrightarrow\mu_{j}(t)-\epsilon(t)$ (27)

in Eq. (26). It suffices to subtract $\epsilon(t)$ from the $\\{\mu_{j}(t)\\}$
to counteract the over or undershooting of the target.

If $\epsilon$ is not compensated in the $\\{\mu_{j}(t)\\}$, the fidelity decay
with $\epsilon$ is more significant as $\hbar_{\text{eff}}\rightarrow 0$. This
indicates that the dominant error source is the perturbed control trajectory
not arriving at the target. It follows from the above argument that a non-zero
distance of the final state to the target state affects more significantly the
overlap of states with smaller volume (larger mean particle number). If
correcting the $\hat{H}_{c}$ by adjusting the driving by the replacement, Eq.
(27), an overall improved robustness to the residual interactions is obtained.
Moreover, notably, the $\hbar_{\text{eff}}=1/4$ coherent state is less
affected than the $\hbar_{\text{eff}}=1$ one. The remaining error source is
the deformation of the state, which is less significant for smaller volumes.

Figure 10: Partial compensation of constant residual interactions. For a
constant $\epsilon(t)=\epsilon$, the fidelity of an uncompensated protocol
quickly decays. The $\hbar_{\text{eff}}=1/4$ decays much faster than
$\hbar_{\text{eff}}=1$. The over or undershooting of the target is the
dominant difficulty. If the behavior of $\epsilon(t)$ is known, it can be
accounted for following the replacement in Eq. (27). This leads to a
significant overall improvement of the fidelity, even favoring the smaller
$\hbar_{\text{eff}}$ case, due to the deformation being the remaining error.

###  3.3 Stabilizing periodic trajectories

Nearly forty years ago, the concept of “quantum scarring of eigenstates” was
introduced by Heller [67]. He argued for and observed excess intensity in the
neighborhood of short unstable periodic trajectories in the quantum
eigenstates of a strongly chaotic dynamical system. The critical point was
that there were expectations that the eigenstates would behave similarly to
random wave functions subject only to energy surface constraints [68, 69, 70,
71], i.e., $\delta(E-H)$ must somehow be respected. Furthermore, there were
quantum ergodic theorems regarding equidistribution and expectation values of
smooth operators giving classical values for nearly all eigenstates [72]. The
excess intensity observed can be considered a weak form of dynamical
localization, and yet somehow must be consistent with the mathematical ergodic
theorems, which is a key element of the surprise of scarring’s existence.

In the many-body context, the Eigenstate Thermalization Hypothesis is the
closest analog of the quantum ergodicity discussed above [73, 74]. Unexpected
deviations from this hypothesis in certain cases might therefore provide the
notion of a many-body quantum scarring analogous to Heller’s quantum scarring.
Recently , experimental effects seen in a 51-atom Rydberg chain were
introduced as many-body scarring [75, 76, 77], although it is not necessarily
clear that it is conceptually the same as in the original context. In fact, an
emergent local regularity of the dynamics, at least on short to medium time
scales, may underlie this many-body scarring phenomenon [78]. Thus, in some
cases many-body quantum scarring may be a result of at least some partial
local integrable or near-integrable dynamics. Nevertheless, there exists at
least one example in the many-body context in which the explicit connection to
unstable periodic mean-field solutions has been made [79]. With the keen
current interest in many-body periodic or partially periodic dynamics, it is
worth looking at how optimal quantum coherent control can be utilized to
stabilize periodic dynamics in a many-body system. The method is applied to
simpler symmetric and much more complicated asymmetric unstable trajectories
next.

####  3.3.1 Symmetric perodic mean-field control trajectory

In Sec. 2.2.4, it is mentioned that the density wave exhibits periodic
population inversion. This does not constitute a periodic trajectory unless
the phase is also periodic at the same time. It turns out that there are truly
periodic trajectories very close to the density wave initial condition for
certain values of $\gamma$. One of these values is selected and illustrated in
Fig. 11. For each cycle of population inversion,

Figure 11: Symmetric periodic trajectory. For $\gamma=2.0607$ there is a
periodic mean-field solution very close to the density wave that is symmetric
for all even and odd sites. After three periods of the occupations (or three
driving periods) the solution returns to its original phase.

the phase advances by $2\pi/3$, which after the third cycle completes the
periodic trajectory. Note the simple appearance, similar to a Lissajous
figure, of the quadrature variables for the trajectory. In essence, the
initially empty site is $\pi/2$ out of phase with the occupied site. There is
an Animation 3 in the supplementary material illustrating the continuous time
evolution.

For the same reasons that the fidelity of the control protocol is not perfect
in Fig. 8, i.e., the trajectory initial conditions do not perfectly match the
initial state centroids, the achievement of perfect periodicity degrades with
increasing time and decreasing $\hbar_{\text{eff}}$. This is illustrated in
Fig. 12

Figure 12: Quantum recurrences in the controlled system. When propagated in
the uncontrolled system, no recurrences appear. For $\hbar_{\text{eff}}=1$,
the recurrences are almost perfect due to the large size of the state in phase
space. In addition, smaller side peaks are visible because the volume is so
large. Those disappear for the smaller inherent volume for the
$\hbar_{\text{eff}}=1/64$ case, together with more imperfect recurrences.
Naturally, the smaller volume of the $\hbar_{\text{eff}}=1/64$ case also
generates a sharper peak.

for the two cases $\hbar_{\text{eff}}=1,1/64$ up to three periods of the
periodic motion. The tall peaks indicate the initial state is recurring at
integer multiples of the periodic control trajectory with the
$\hbar_{\text{eff}}=1$ case falling off slowly with increased time, and the
$\hbar_{\text{eff}}=1/64$ case falling off much quicker. Perfect recurrences
would be restored by an initial shift of the initial state centroids to the
initial conditions of the periodic mean-field trajectory.

####  3.3.2 Asymmetric periodic mean-field control trajectory

Although more challenging, it is possible to locate periodic mean-field
trajectories with far less symmetry that are in the close neighborhood of a
density wave. An example is shown in Fig. 13. In either the quadrature
variables, which do not even faintly resemble a Lissajous figure, or even just
the occupancies, the control trajectory follows a highly nontrivial evolution
before it returns. The complete time evolution can again be found for this
case as Animation 4 in the supplementary material. The quality of the
recurrence for this longer, more complicated trajectory, is not quite as good
as for the symmetric trajectory in the previous subsection.

Figure 13: Asymmetric periodic trajectory. There exist also a dense set of
highly nontrivial periodic mean-field trajectories, which could serve as
control trajectories. Here, a trajectory is shown for $\gamma=1.3$. Its
evolution is extremely close to a periodic mean-field trajectory. Due to
technical complications, the exact periodic trajectory’s initial conditions
are not sought to greater precision, and in any case, the density wave is not
perfectly centered on the control trajectory either.

##  4 Coherent quantum targeting in $2D$ lattices

A considerable utility of cold atom optical lattices lies in their capacity to
be created in any number of dimensions. A priori, the control protocol is not
limited to $1D$ systems, however higher-dimensional systems lead to
exponentially expanding search spaces for heteroclinic trajectories. Thus, to
target arbitrary states requires sophisticated search algorithms [80, 59], and
even so it may become effectively impossible. Nevertheless, there do exist
certain classes of symmetric, periodic lattice configurations for which it is
possible to map the mean-field dynamics to that of $1D$ periodic rings [60].
In this way, the protocol can be reduced to identifying the corresponding
trajectory in the lower-dimensional space, which renders the search feasible.
Here, $2D$ lattice examples are treated in which the mean-field dynamics are
explicitly mapped to that of a $1D$ ring, reducing the problem again to
finding control trajectories in a numerically accessible search space.

There is more than one way to generate an appropriate mapping for these
purposes. For the one relied on below, the sites are mapped into higher-
dimensional lattices such that the nearest-neighbor associations are
preserved. There is a change in nearest-neighbor multiplicity with increasing
dimension that requires a renormalization of the hopping parameter, which once
taken into account generates a dynamics in the larger lattice that appears as
multiple copies of the $1D$ ring used in the mapping [60]. As long as the
initial and target states respect the implied discrete symmetries, the optimal
coherent quantum targeting in higher-dimensional lattices reduces to again
identifying optimal control trajectories in the $1D$ ring.

Following this particular prescription, two possible lattice arrangements are
shown in Fig. 14, that allow for the dynamics

(a)

(b)

(c)

Figure 14: Two possible schematic mappings of $2D$ mean-field dynamics to a
$1D$ ring. By initializing the condensates on a $2D$ periodic lattice in the
arrangement shown in panels (a) or (c), the mean-field dynamics of each
lattice site can be mapped to the dynamics of the corresponding site on a $1D$
ring, (b), with double the hopping parameter, i.e.
$J_{\text{1D}}=2J_{\text{2D}}$. This necessitates the relation between the
$1D$ ring and $2D$-lattice chemical potentials to be
$\vec{\mu}(t)/2\Rightarrow\mu_{\vec{R}}(t)$.

of a periodic $4\times 4$ lattice to be mapped onto the dynamics of a four-
site ring. In fact, the mappings work for any $4n\times 4m$ $2D$ lattice with
periodic boundary conditions (including $n\rightarrow\infty$). Here the mean-
field dynamics reduce to the four-site ring with double the hopping parameter
$J_{\text{1D}}=2J_{\text{2D}}$ [60]. The two examples given are a lattice
consisting of four-site ring building blocks or of shifted stripes,
respectively. By initializing the lattice in one of those symmetric ways, it
is possible to guide a many-body state along a control trajectory by driving
the $2D$-lattice system with its mean-field occupations, $\mu_{\vec{R}}(t)$,
corresponding to those of the reduced $1D$ system, $\vec{\mu}(t)$. The
$2D$-lattice control Hamiltonian,

$\displaystyle\hat{H}_{\text{c}}^{\text{2D}}\left(t\right)=-\sum_{\left<\vec{R},\vec{R}^{\prime}\right>}\big{(}\hat{a}^{\dagger}_{\vec{R}}\hat{a}_{\vec{R}^{\prime}}+\text{h.c.}\big{)}+\sum_{\vec{R}}\mu_{\vec{R}}(t)\hat{n}_{\vec{R}},$
(28)

is analogous to the $1D$ case after accounting for the doubled number of
neighboring sites in the chemical potentials.

The two target configurations shown in Fig. 15 result from the mappings of
Fig. 14 applied to the $1D$-ring control trajectory of Sec. 3.1. The initial
condition consists of alternating occupied and unoccupied sites satisfying
either mapping. The full evolutions are shown in the supplementary material,
Animation 5 and Animation 6, respectively, where the sites evolve according to
the symmetry of the mappings. This way it is possible, for example, to target
the homogeneous lattice state where the respective condensate phases form a
vortex structure as shown in Fig. 15a or a shifted stripe structure as shown
in Fig. 15b.

(a)

(b)

Figure 15: Targeting in $2D$ lattices. Mapping the mean-field dynamics of a
symmetrically arranged $2D$ lattice to a four-site ring allows for targeting
certain classes of states by identifying their respective lower-dimensional
control trajectory. By initializing the lattice using the arrangements
presented in Fig. 14(a) or Fig. 14(c) with the initial conditions of the
control trajectories previously discussed and propagating the mean-field
dynamics, it is possible to prepare the homogeneous lattice states shown in
(a) and (b), respectively. The occupation of each site is encoded in the arrow
length and the phases in the angle of the corresponding arrows read clockwise
from the right horizontal.

Applying the same logic it is possible to also create periodic many-body
states on a $2D$ lattice by for example using the periodic trajectories
discussed in Sec. 3.3. The entire time evolution is again presented in the
supplementary material, Animation 7 and Animation 8, respectively. Naturally
there also exist multiple mappings that reduce the dynamics of a $3D$ lattice
to that of a four-site ring with adjusted hopping parameter
$J_{\text{1D}}=3J_{\text{3D}}$.

##  5 Conclusion

The presence of chaos poses fundamental challenges for controlling the
evolution of a classical dynamical system, nevertheless it is long known that
it is possible to convert the inherent difficulties caused by instability and
ergodicity into a resource [4]. Recently, it was shown that with respect to
quantum systems, an additional requirement is needed, i.e. the suppression of
quantum state spreading [11]. Done originally in a Schrödinger or single
particle context, a feasible extension into the many-body domain, here for
ultracold atoms in an optical lattice, requires the identification of a time-
dependent control Hamiltonian, and in this way build a quantum simulator
controlling many-body quantum chaos.

The extension of chaotic classical targeting to coherent quantum targeting is
facilitated through the use of semiclassical methods in which trajectories in
the mean-field/classical limit and their stability analysis provide the basis
for determining the necessary perturbations or time-dependent modifications of
the system. The quadrature phase space formulation of a many-body system of
ultracold bosons results in a mean field limit, i.e. large filling factor, of
ordinary Hamiltonian dynamics. Thus, the adaptation to a many-body system
parallels closely the prior work in a Schrödinger dynamics context. In
particular, the same class of trajectories, heteroclinic ones, provides the
time-dependent chemical potentials that determine the control Hamiltonian.

As in [11, 12], the protocol is designed for guiding localized quantum states.
All of the results are shown for minimum uncertainty coherent states, but they
apply immediately to the wider and important class that includes number
projected coherent states. To control the spreading of the many-body state a
control system is introduced that is essentially a time-dependent harmonic
oscillator driven with the on-site occupations along the specific control
trajectory. In practice this means that interactions between the bosons is
turned off in an experiment, which is something that in optical lattice
systems can be achieved to high precision using Feshbach resonances. Together
with the fact that the protocol reduces to controlling the chemical potentials
of the individual sites in a smooth, time-dependent fashion we believe this
protocol to be suitable for practical experimental implementation.

Due to the harmonic nature of the system and the $U(1)$ symmetry of the
dynamics any (number projected) coherent state placed on the control
trajectory used to drive the Hamiltonian satisfies the time evolution of the
system exactly. This means that the protocol provides a tool to propagate a
broad range of states from arbitrary initial and target points in quadratures
given that the underlying dynamics is sufficiently ergodic, opening up new
possibilities in state preparation for optical lattices. An example is shown
of a relatively short trajectory that can be used to target a homegeneous
lattice state with well defined, non-trivial phases between the condensates.

There are a few aspects of quantum targeting that may impact any assessment of
how robust it can be against imperfections. First, the ideal would be to make
the tiny shifts of the initial state centroids to the initial conditions of
the optimal control (heteroclinic) trajectories. Likewise, the final state
centroids on those trajectory endpoints would be slightly shifted to the
target state’s values. In fact, the initial shift is the single most important
element in classical targeting, whereas suppressing the state spreading is the
critical element in the quantum case. By virtue of the uncertainty principle,
the optimal trajectory is already well represented inside the quantum state’s
Wigner transform. By not performing the shifts, a decrease of the fidelity is
introduced into the protocol, which gets worse as the filling factor
increases. This is observed numerically by calculating the overlap between the
propagated and desired target state.

A second imperfection would be due to the inability to turn off the
interactions exactly, which left uncompensated quickly breaks down the control
protocol. Although, it is anticipated that the on-site interaction strength
can, in many circumstances, be controlled experimentally to high precision, if
any nonvanishing residual interaction can be characterized, then a simple
rectification expressed in Eq. (27) solves the over or undershooting effect,
which would be the main contributor to a fidelity decay. The local deformation
of the initially circular Wigner density contours at the target state remains,
but improves with increasing filling factor. A third issue, that of imperfect
control over the $\\{\mu_{j}(t)\\}$, requires appropriate modeling of noise or
other imperfections and is left for future consideration.

For the examples given, it is necessary to introduce chemical potential
offsets in order to place an initial state and a target state on the same
energy surface such that chaotic trajectories can provide transport paths
between them. If it were desired not to have these offsets at the end,
although not discussed here, it would be possible to consider say
adiabatically turning off the offsets and seeking some modified chaotic
transport path between the initial and target states. For example, the
momentum target state of Sec. 3 would become a fixed state of the system
without interactions and chemical potential offsets, which might be a desired
goal as well. Additional time-dependence in the chemical potential offsets
does not add much, if any, additional complications with respect to solving
Hamilton’s equations, but it adds an extra dimension to the complexity of the
search for an optimal control trajectory. This issue of possible unwanted
final chemical potential offsets does not arise if the initial and targets
states are identical, as for the periodic trajectory cases, as they are
automatically on the same energy surface without offsets.

The presented $1D$ examples are of greater utility and versatility than they
might appear at first sight. By mapping the sites into higher-dimensional
lattices such that the nearest-neighbor associations are preserved, it is
possible to control much larger and higher-dimensional lattices without
suffering from the exponential growth of search volumes in phase space. For
initial and target states that respect certain symmetric lattice
configurations, optimal coherent quantum targeting in higher-dimensional
lattices is no more complicated than finding control trajectories in the $1D$
ring to which the mean-field dynamics can be mapped. On the other hand, for
arbitrary lattice configurations in high dimensions, a systematic search for
control trajectories becomes exponentially challenging with system size.

There are a number of interesting directions that future optimal coherent
quantum targeting research might follow from the technical/pratical side to
more general extensions. This includes extending the methods to spin chains
and fermionic systems, which present several new challenges. It might be
advantageous to replace the current methods of identifying control
trajectories with machine learning methods. Additionally, more work is needed
to understand how the exponential speed up of heteroclinic pathways relates to
quantum speed limits. Finally, a direction worth mentioning is the
investigation of many-body interference using the control protocol.

### Acknowledgments

We thank Monika Aidelsburger for useful conversations. We acknowledge support
by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation),
project Ri681/15-1, within the Reinhart-Koselleck Programme and Vielberth
Foundation for financial support. MS is funded through a fellowship by the
Studienstiftung des Deutschen Volkes.

##  Appendix A A formulation of quantum control from the perspective of
auxiliary fields

The purpose of this appendix is to provide a picture of the control protocol
discussed in the main text from the point of view of the semiclassical
approximation implemented within the auxiliary field formulation of
interacting field theories. The auxiliary field method, see for example [81],
is an exact representation of the path integral form of the propagator where
interactions are accounted for by the response of the system to all possible
time-dependent external potentials. It is based on the functional identity

$\int\frac{{\cal D}[x]}{Z(M)}{\rm
e}^{\frac{i}{2}\int_{0}^{t}ds\int_{0}^{t}ds^{\prime}x(s)M(s,s^{\prime})x(s^{\prime})}{\rm
e}^{i\int_{0}^{t}dsx(s)y(s)}={\rm
e}^{\frac{i}{2}\int_{0}^{t}ds\int_{0}^{t}ds^{\prime}y(s)(M^{-1})(s,s^{\prime})y(s^{\prime})}$
(29)

also known in the field theory literature as the Hubbard-Stratanovich (HS)
transformation. Here $M(s,s^{\prime})$ is a kernel with a functional inverse
defined by

$\int_{0}^{t}ds^{\prime}M(s,s^{\prime})(M^{-1})(s^{\prime},s^{\prime\prime})=\delta(s-s^{\prime\prime}),$
(30)

and $Z(M)$ is a normalization constant defined by comparing both sides of Eq.
(29) with $y(s)=0$. The generalization for multi-component fields,
($\vec{x}(s),\vec{y}(s))$, is straightforward by also introducing a matrix
structure $M(s,s^{\prime})\to{\bf M}(s,s^{\prime})$.

To apply the HS transformation to the case of interest here, start with the
coherent state form of the propagator for interacting bosonic fields,

$K(\vec{\phi}_{\beta},\vec{\phi}_{\alpha};t):=\langle\vec{\phi}_{\beta}|{\rm
e}^{-\frac{i}{\hbar}\hat{H}t}|\vec{\phi}_{\alpha}\rangle$ (31)

with $\hat{H}$ given by Eq. (2). This propagator admits a path integral
representation, found using standard techniques generically given by

$K(\vec{\phi}_{\beta},\vec{\phi}_{\alpha};t)=\int{\cal
D}[\vec{\phi},\vec{\phi}^{*}]{\rm e}^{i(R_{\rm
free}[\vec{\phi},\vec{\phi}^{*}]+R_{\rm int}[\vec{\phi},\vec{\phi}^{*}])},$
(32)

and defined by the non-interacting (free) $R_{\rm free}$ and interacting
$R_{\rm int}$ terms in the action. As a rule, the free term contains
functionals that are at most quadratic in the complex fields $\phi_{j}(s)$.
Interaction terms, for the case of particular interest here, are in turn given
as linear combinations of terms bilinear in the local occupations

$\hat{V}_{\rm
int}=\frac{1}{2}\sum_{i,j}v_{ij}\hat{n}_{i}\left(\hat{n}_{j}-\delta_{ij}\right){\rm\
\ \ \ with\ \ \ }\hat{n}_{i}=\hat{a}^{\dagger}_{i}\hat{a}_{i},{\rm\ \ and\
}v_{ij}=v_{ji}.$ (33)

Therefore, neglecting issues of ordering in defining the classical action
functional (see [82] and [83] for details)

$R_{\rm
int}[\vec{\phi},\vec{\phi}^{*}]=\frac{1}{2}\int_{0}^{t}ds\sum_{i,j}v_{ij}|\phi_{i}(s)|^{2}\left(|\phi_{j}(s)|^{2}-\delta_{ij}\right).$
(34)

As is well known, the reason why the path integral, Eq. (32), cannot be
evaluated in closed form is that the interaction term is not of Gaussian type.
However, by means of the HS transformation, the fourth-order terms, Eq. (34),
can be decoupled in favor of the (vector) real auxiliary field
$\vec{\sigma}(s)$ obtaining

$K(\vec{\phi}_{\beta},\vec{\phi}_{\alpha};t)=\int\frac{{\cal
D}[\vec{\sigma}]}{Z({\bf v})}{\rm
e}^{\frac{i}{2}\int_{0}^{t}ds\vec{\sigma}(s){\bf
v}^{-1}\vec{\sigma}(s)}K_{\vec{\sigma}}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha};t)$
(35)

where the reduced propagator

$K_{\vec{\sigma}}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha};t)=\int{\cal
D}[\vec{\phi},\vec{\phi}^{*}]{\rm e}^{i(R_{\rm
free}[\vec{\phi},\vec{\phi}^{*}]+\int_{0}^{t}ds\sum_{j}\sigma_{j}(s)n_{j}(s))},{\rm\
\ with\ }n_{j}(s)=|\phi_{j}(s)|^{2},$ (36)

now has the key property of representing a free field evolving under the time-
dependent external field $\sigma_{j}(s)$ coupled at the $j$-th site to the
local occupations, i.e., a time-dependent local chemical potential.

Next, Gaussian path integrals, even with time-dependent parameters as in Eq.
(36) preventing us to find their solution in closed form, can be exactly
expressed by means of the solution of the corresponding classical (mean-field)
problem, which in this case is unique. In other words,

$K_{\vec{\sigma}}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha};t)=A(t,[\vec{\sigma}]){\rm
e}^{i(R_{\rm free}^{\rm
mf}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha},[\vec{\sigma}])+\int_{0}^{t}ds\sum_{j}\sigma_{j}(s)n_{j}^{\rm
mf}(s,\vec{\phi}_{\beta},\vec{\phi}_{\alpha},[\vec{\sigma}])))}$ (37)

where $R^{\rm mf},n_{j}^{\rm mf}$ are now functionals of the auxiliary fields
$\vec{\sigma}$, as indicated by $[.]$, and functions of
$t,\vec{\phi}_{\beta},\vec{\phi}_{\alpha}$. These objects are obtained by

$\displaystyle R_{\rm free}^{\rm
mf}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t,[\vec{\sigma}])$ $\displaystyle=$
$\displaystyle R_{\rm free}[\vec{\phi}(s)=\vec{\phi}^{\rm
mf}(s,\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t,[\vec{\sigma}]),\vec{\phi}^{*}(s)=\vec{\phi}^{\rm
mf}(s,\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t,[\vec{\sigma}])^{*}]$

and

$n_{j}^{\rm
mf}(s,\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t,[\vec{\sigma}])=|\phi_{j}^{\rm
mf}(s,\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t,[\vec{\sigma}])|^{2}$ (39)

from the unique solution $\vec{\phi}^{\rm
mf}(s,\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t,[\vec{\sigma}])$ of the Euler-
Lagrange (mean field) equations

$\frac{\delta}{\delta\vec{\phi}^{*}}\left(R_{\rm
free}[\vec{\phi},\vec{\phi}^{*}]+\int_{0}^{t}ds\sum_{\alpha}\sigma_{\alpha}(s)n_{\alpha}(s)\right)_{\vec{\phi}=\vec{\phi}^{\rm
mf}}=0$ (40)

satisfying

$\vec{\phi}(s=0)=\vec{\phi}_{\alpha}{\rm\ \ and\
}\vec{\phi}(s=t)=\vec{\phi}_{\beta}.$ (41)

Finally, the van Vleck-Morette determinant $A(t,[\vec{\sigma}])$ can be shown
to be independent of the boundary conditions
$\vec{\phi}_{\beta},\vec{\phi}_{\alpha}$ for quadratic actions, and its
precise form can be found in [84].

It is instructive to consider the specific situation for a system described by
the Bose-Hubbard Hamiltonian. In this case, besides the interaction term, Eq.
(33), the action functional contains a free term of the form

$R_{\rm
free}[\vec{\phi},\vec{\phi}^{*}]=\int_{0}^{t}ds\sum_{j}\imaginary\left(\phi_{j}(s)\frac{d\phi_{j}^{*}(s)}{ds}\right)+\int_{0}^{t}ds\sum_{i,j}h_{ij}\phi_{i}^{*}(s)\phi_{j}(s)$
(42)

and therefore the mean field solutions that fully determine the reduced
propagator for a given auxiliary field $\vec{\sigma}(s)$ are easily obtained
from the variation in Eq. (40) to be

$i\frac{d}{ds}z_{i}(s)=\sum_{j}h_{ij}z_{j}(s)+\sigma_{i}(s)z_{i}(s),$ (43)

namely the mean field equations for a bosonic free field under a time-
dependent chemical potential, as expected.

Up to here, the original interacting problem is transformed into a coherent
weighted sum of reduced propagators over all possible time-dependent auxiliary
fields, each of them representing linear evolution,

$\displaystyle K(\vec{\phi}_{\beta},\vec{\phi}_{\alpha};t)=\int\frac{{\cal
D}[\vec{\sigma}]}{Z({\bf v})}$ $\displaystyle{\rm
e}^{\frac{i}{2}\int_{0}^{t}ds\vec{\sigma}(s){\bf
v}^{-1}\vec{\sigma}(s)}A(t,[\vec{\sigma}])$ $\displaystyle\times{\rm
e}^{i(R_{\rm free}^{\rm
mf}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t,[\vec{\sigma}])+\int_{0}^{t}ds\sum_{j}\sigma_{j}(s)n_{j}^{\rm
mf}(s,\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t,[\vec{\sigma}])))},$

which is an exact expression. It is at this stage that the semiclassical
approximation, essential for the construction of the control protocol, enters.
Evaluating the integral over auxiliary fields in a saddle point approximation,
that in turn requires the assumption that the prefactors $A(t,[\sigma])$ are
smooth, the full propagator turns out to be dominated by the configurations,
$\vec{\sigma}=\vec{\sigma}_{\gamma}$, satisfying the corresponding vanishing
of the first variation, namely

$\displaystyle\frac{\delta}{\delta\vec{\sigma}}\left(\frac{1}{2}\int_{0}^{t}ds\vec{\sigma}(s){\bf
v}^{-1}\vec{\sigma}(s)\right.$ $\displaystyle+$ $\displaystyle R_{\rm
free}^{\rm mf}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t,[\vec{\sigma}])$
$\displaystyle+$
$\displaystyle\left.\left.\int_{0}^{t}ds\sum_{j}\sigma_{j}(s)n_{j}^{\rm
mf}(s,\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t,[\vec{\sigma}])\right)\right|_{\vec{\sigma}_{\gamma}}=0.$

Using the implicit dependence of $R_{\rm free}^{\rm mf}$ and $n_{\alpha}^{\rm
mf}$ on the auxiliary field as given in Eqs. (A, 39) and the mean-field
equations, Eq. (40), this gives

${\bf v}^{-1}\vec{\sigma}_{\gamma}(s)=\vec{n}^{\rm
mf}(s,\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t,[\vec{\sigma}_{\gamma}])$ (46)

or

$\sigma_{\gamma,i}(s)=\sum_{j}v_{ij}n_{j}^{\rm
mf}(s,\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t,[\vec{\sigma}_{\gamma}])$ (47)

as the set of equations that determines the dominant auxiliary fields
$\vec{\sigma}_{\gamma}$. The extreme nonlinearity of this system resulting
from the role of $\vec{\sigma}$ as time-dependent contributions to the
couplings in the linear equations (40, 43) makes it natural to expect a
countable, but large, number of solutions.

The implicit problem in Eq. (47) can be made more familiar by invoking the
dependence of the mean field occupations $n_{\alpha}^{\rm mf}$ on
$\vec{\sigma}$ as given precisely by the mean-field equations (40).
Substitution of Eq. (47) into the equations obtained from Eq. (40) then
results in the explicit solution

$\sigma_{\gamma,i}(s)=\sum_{j}v_{ij}|\phi_{\gamma,j}^{\rm clas}|^{2}$ (48)

where $\vec{\phi}^{\rm~{}class}_{\gamma}$ solves the nonlinear equation

$i\frac{d}{ds}\phi_{i}(s)=\sum_{j}h_{ij}\phi_{j}(s)+\sum_{j}v_{ij}n_{i}(s)\phi_{j}(s)$
(49)

for the Bose-Hubbard type of action. Note that the fields
$\vec{\sigma}_{\gamma}$ inherit the boundary conditions, Eq. (41), finally
expressing the equivalence of the semiclassical approximation at the mean
field level with a nonlinear classical field equation. This is reassuring as
it is expected that the introduction of auxiliary fields at an intermediate
step should not change the classical limit of the original interacting
problem, precisely given by Eq. (49).

The importance of this representation is, however, clear upon return to the
semiclassical evaluation of the exact integral in Eq. (A), given now within
the saddle point approximation by

$\displaystyle K^{\rm sc}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha};t)$
$\displaystyle=$ $\displaystyle\sum_{\gamma}\frac{{\cal
D}_{\gamma}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t)}{Z({\bf v})}\left.{\rm
e}^{\frac{i}{2}\int_{0}^{t}ds\vec{\sigma}(s){\bf
v}^{-1}\vec{\sigma}(s)}\right|_{\vec{\sigma}=\vec{\sigma}_{\gamma}}$
$\displaystyle\times$ $\displaystyle\left.A(t,[\vec{\sigma}]){\rm e}^{i(R_{\rm
free}^{\rm
mf}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t,[\vec{\sigma}])+\int_{0}^{t}ds\sum_{j}\sigma_{j}(s)n_{j}^{\rm
mf}(s,\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t,[\vec{\sigma}])))}\right|_{\vec{\sigma}=\vec{\sigma}_{\gamma}},$

where the prefactors ${\cal D}_{\gamma}$ result from the usual integration
over Gaussian fluctuations around the classical configurations
$\vec{\sigma}_{\gamma}$. The second line of Eq. (A) contains the reduced
propagator for the quadratic action, Eqs. (36, 37) in its semiclassical (but
exact) form. Thus,

$K^{\rm sc}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha};t)=\sum_{\gamma}\frac{{\cal
D}_{\gamma}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha},t)}{Z({\bf v})}{\rm
e}^{\frac{i}{2}\int_{0}^{t}ds\vec{\sigma}_{\gamma}(s){\bf
v}^{-1}\vec{\sigma}_{\gamma}(s)}\times
K_{\vec{\sigma}_{\gamma}}(\vec{\phi}_{\beta},\vec{\phi}_{\alpha};t)$ (51)

where $K_{\vec{\sigma}_{\gamma}}$ is the exact coherent state propagator for
the $\gamma$-th classical auxiliary field that, given the quadratic nature of
the corresponding action, acts invariantly on the initial coherent state,
i.e., it transports $\vec{\phi}_{\alpha}$ without dispersion along the
classical phase space trajectory.

Equation (51) is the final result of this analysis. It expresses the
semiclassical approximation to the coherent state propagator of the
interacting theory as a coherent sum over a countable set of contributions.
Each term in the sum represents the exact and dispersionless quantum
propagation of the initial coherent state under auxiliary fields playing the
role of time-dependent chemical potentials that are in turn determined by the
solution of the interacting problem. In the language of the main text, the
full interacting propagator is given as a sum over all possible classical
(mean field) control protocols. The control protocol proposed in the main text
amounts then to choosing, based on minimal time $t$ to dynamically connect
$\vec{\phi}_{\alpha}$ with $\vec{\phi}_{\beta}$ or other physical
considerations, a single one of the auxiliary fields.

##  Appendix B Coherent Time-Evolution in Control Systems

A proof is outlined showing that a coherent state
$|{\vec{\phi}(t)}\rangle_{\rm coh}$ centered on any mean-field solution
$\vec{\phi}(t)$ of the control Hamiltonian Eq. 18 solves the controlled
quantum time evolution

$\displaystyle|{\vec{\phi}(t)}\rangle_{\rm coh}=\hat{U}_{\rm
c}(t)|{\vec{\phi}_{0}}\rangle_{\rm
coh}=\hat{\mathcal{T}}\exp\\{i\int_{0}^{t}\hat{H}_{\rm
c}(s)\operatorname{d}\\!s\Big{\\}}|{\vec{\phi}_{0}}\rangle_{\rm coh},$

where $\vec{\phi}_{0}=\vec{\phi}(0)$. More explicitly, the coherent state is
shown to be a solution to the corresponding Schrödinger equation, i.e.,

$\displaystyle
i\frac{\operatorname{d}\\!}{\operatorname{d}\\!t}|{\vec{\phi}(t)}\rangle_{\rm
coh}=\hat{H}_{\rm c}(t)|{\vec{\phi}(t)}\rangle_{\rm coh}\ ,$

which is equivalent to the previous equation. To this end, recall that the
classical control Hamiltonian, Eq. 16, can also be written in terms of complex
mean-fields

$\mathcal{H}_{\text{c}}(\vec{\phi},\vec{\phi}^{*},t)=-\sum\limits_{j=1}^{L}\left\\{\phi_{j+1}^{*}+\phi_{j-1}^{*}-\mu_{j}(t)\phi_{j}^{*}\right\\}\phi_{j},$

that evolve according to the equations motion

$\displaystyle i\dot{\phi}_{j}$
$\displaystyle=\frac{\partial\mathcal{H}_{\text{c}}}{\partial\phi_{j}^{*}}=-\left\\{\phi_{j+1}+\phi_{j-1}-\mu_{j}(t)\phi_{j}\right\\}$
$\displaystyle i\dot{\phi}_{j}^{*}$
$\displaystyle=-\frac{\partial\mathcal{H}_{\text{c}}}{\partial\phi_{j}}=\left\\{\phi_{j+1}^{*}+\phi_{j-1}^{*}-\mu_{j}(t)\phi_{j}^{*}\right\\}.$

Start by assuming that there is a coherent state solution and look for a
contradiction or consistency. Taking the left hand side time derivative of the
coherent state reduces to taking the derivative of the mean-field since the
mean particle number conservation along any classical solution implies
$\lVert{\vec{\phi}}\rVert^{2}=N=\text{const.}$

$\displaystyle
i\frac{\operatorname{d}\\!}{\operatorname{d}\\!t}|{\vec{\phi}(t)}\rangle_{\rm
coh}=$ $\displaystyle
ie^{-\lVert{\vec{\phi}}\rVert^{2}/2}e^{\vec{\phi}(t)\cdot\hat{\vec{a}}^{\dagger}}\left(\sum\limits_{j=1}^{L}\dot{\phi}_{j}(t)\hat{a}_{j}^{\dagger}\right)\ket{0}.$

Next use the classical mean-field equations of motion from above

$\displaystyle
i\frac{\operatorname{d}\\!}{\operatorname{d}\\!t}|{\vec{\phi}(t)}\rangle_{\rm
coh}$
$\displaystyle=-ie^{-\lVert{\vec{\phi}}\rVert^{2}/2}e^{\vec{\phi}(t)\cdot\hat{\vec{a}}^{\dagger}}\sum_{j=1}^{L}\left(\left\\{\phi_{j+1}(t)+\phi_{j-1}(t)-\mu_{j}(t)\phi_{j}(t)\right\\}\hat{a}_{j}^{\dagger}\right)\ket{0}$
$\displaystyle=-\sum_{j=1}^{L}\left(\left\\{\phi_{j+1}(t)+\phi_{j-1}(t)-\mu_{j}(t)\phi_{j}(t)\right\\}\hat{a}_{j}^{\dagger}\right)|{\vec{\phi}(t)}\rangle_{\rm
coh}.$

With respect to the right hand side, use the definition given in Eq. 9 of a
coherent state being an eigenstate of the annihilation operator to obtain

$\displaystyle\hat{H}_{\rm c}(t)|{\vec{\phi}(t)}\rangle_{\rm coh}$
$\displaystyle=-\sum_{j=1}^{L}\left(\left\\{\hat{a}_{j+1}^{\dagger}+\hat{a}_{j-1}^{\dagger}+\mu_{j}(t)\hat{a}_{j}^{\dagger}\right\\}\hat{a}_{j}\right)|{\vec{\phi}(t)}\rangle_{\rm
coh}$
$\displaystyle=-\sum_{j=1}^{L}\left(\hat{a}_{j}^{\dagger}\left\\{\phi_{j+1}(t)+\phi_{j-1}(t)+\mu_{j}(t)\phi_{j}(t)\right\\}\right)|{\vec{\phi}(t)}\rangle_{\rm
coh},$

which matches exactly the expression coming from the left hand side.

##### Remark 1

The proof for number projected states is analogous using the property

$\hat{a}_{j}|{\vec{\phi}(t)}\rangle_{\rm
proj}^{N}=\phi_{j}|{\vec{\phi}(t)}\rangle_{\rm proj}^{N-1}.$

##### Remark 2

The control trajectory used to drive the chemical potentials
$\left\\{\mu_{j}(t)\right\\}$ is a special case since it is a solution to both
the classical equations of motion of the control system, as well as the Bose-
Hubbard mean-field limit. With that information, the above proof implies that
a coherent state initially centered on the control trajectory arrives at the
end of the trajectory with perfect fidelity.

##### Remark 3

For any specified control Hamiltonian of the form of Eq. (18), recalling the
equation for calculating the stability matrix, every initial mean-field
condition or trajectory leads to the exact same time-dependent stability
matrix. Thus, all trajectories have to rotate around the control trajectory
and the stability matrix solution is a time-dependent orthogonal matrix. This
implies that for a coherent state initially centered off the control
trajectory by $||\delta\vec{x}_{\text{init}}||$, this distance stays constant.
By using the properties of the stability matrix, and Heller’s linearized wave
packet dynamics [66], it would be possible to provide an analytical formula
for the fidelity, for example shown in Fig. 8 based solely on the knowledge of
the initial state and the control trajectory.

## References

  * [1] Stanislaw M. Ulam.  On the possibility of extracting energy from gravitational systems by navigating space vehicles.  Technical Report LAMS-2219, Los Alamos National Laboratory, 1958.
  * [2] R. Farquhar, D. Muhonen, and S. A. Davis.  Trajectories and orbit maneuvers for the ISEE-3/ICE comet mission (1985).  J. Astronaut. Sci., 33:235–254, 1985.
  * [3] Robert W. Farquhar.  The flight of ISEE-3/ICE: Origins, mission history, and a legacy.  J. Astronaut. Sci., 49:23–73, 2001.
  * [4] Edward Ott, Celso Grebogi, and James A. Yorke.  Controlling chaos.  Phys. Rev. Lett., 64:1196–1199, 1990.
  * [5] Troy Shinbrot, Edward Ott, Celso Grebogi, and James A. Yorke.  Using chaos to direct trajectories to targets.  Phys. Rev. Lett., 65:3215–3218, 1990.
  * [6] Eric J. Kostelich, Celso Grebogi, Edward Ott, and James A. Yorke.  Higher-dimensional targeting.  Phys. Rev. E, 47:305–310, 1993.
  * [7] Erik M. Bollt.  Controlling chaos, targeting, and transport.  PhD thesis, University of Colorado, 1995.
  * [8] Christian G. Schroer and Edward Ott.  Targeting in Hamiltonian systems that have mixed regular/chaotic phase spaces.  Chaos, 7:512–519, 1997.
  * [9] Edward Ott.  Controlling chaos.  Scholarpedia, 1(8):1699, 2006.  revision #91167.
  * [10] Henri Poincaré.  Les méthodes nouvelles de la mécanique céleste, volume 3.  Gauthier-Villars et fils, Paris, 1899.
  * [11] Steven Tomsovic, Juan Diego Urbina, and Klaus Richter.  Controlling quantum chaos: Optimal coherent targeting.  Phys. Rev. Lett., 130:020201, 2023.  arXiv:2211.07408 [nlin.CD].
  * [12] Steven Tomsovic, Juan Diego Urbina, and Klaus Richter.  Controlling quantum chaos: time-dependent kicked rotor.  Phys. Rev. E, 108:044202, 2023.  arXiv:2305.14187 [quant-ph].
  * [13] B. V. Chirikov.  A universal instability of many-dimensional oscillator systems.  Phys. Rep., 52:263–379, 1979.
  * [14] Felix M. Izrailev.  Simple models of quantum chaos: Spectrum and eigenfunctions.  Phys. Rep., 196:299–392, 1990.
  * [15] Werner Heisenberg.  Über den anschaulichen inhalt der quantentheoretischen kinematik und mechanik.  Zeit. der Phys., 43:172–198, 1927.
  * [16] Eugene P. Wigner.  On the quantum correction for thermodynamic equilibrium.  Phys. Rev., 40:749–759, 1932.
  * [17] Tomi H. Johnson, Stephen R. Clark, and Dieter Jaksch.  What is a quantum simulator.  EPJ Quant. Tech., 1:10, 2014.
  * [18] H. A. Gersch and G. C. Knollman.  Quantum cell model for bosons.  Phys. Rev., 129:959–967, 1963.
  * [19] Dieter Jaksch, Cristoph Bruder, J. Ignatio Cirac, Crispin W. Gardiner, and Peter Zoller.  Cold bosonic atoms in optical lattices.  Phys. Rev. Lett., 81:3108–3111, 1998.
  * [20] David J. Tannor and Stuart A. Rice.  Control of selectivity of chemical reaction via control of wave packet evolution.  J. Chem. Phys., 83:5013–5018, 1985.
  * [21] Paul Brumer and Moshe Shapiro.  Control of unimolecular reactions using coherent light.  Chem. Phys. Lett., 126:541–546, 1986.
  * [22] Richard S. Judson and Herschel Rabitz.  Teaching lasers to control molecules.  Phys. Rev. Lett., 68:1500–1503, 1992.
  * [23] Warren S. Warren, Herschel Rabitz, and Mohammed Dahleh.  Coherent control of quantum dynamics: The dream is alive.  Science, 259:1581–1589, 1993.
  * [24] Anthony P. Peirce, Mohammed A. Dahleh, and Herschel Rabitz.  Optimal control of quantum-mechanical systems: Existence, numerical approximation, and applications.  Phys. Rev. A, 37:4950–4964, 1988.
  * [25] Ronnie Kosloff, Stuart A. Rice, Pierre Gaspard, S. Tersigni, and David J. Tannor.  Wavepacket dancing: Achieving chemical selectivity by shaping light pulses.  Chem. Phys., 139:201–220, 1989.
  * [26] Daoyi Dong and Ian R. Petersen.  Quantum control theory and applications: A survey.  IET Control Theory Appl., 4:2651–2671, 2010.
  * [27] M. R. James.  Optimal quantum control theory.  Annu. rev. control robot. auton. syst., 4:343–367, 2021.
  * [28] Steven Tomsovic and Julie H. Lefebvre.  Can wave packet revivals occur in chaotic quantum systems?  Phys. Rev. Lett., 79:3629–3632, 1997.
  * [29] Jiangbin Gong and Paul Brumer.  Quantum chaos meets coherent control.  Annu. Rev. Phys. Chem., 56:1–23, 2005.
  * [30] Martin Gruebele and Peter G. Wolynes.  Quantizing ulam’s control conjecture.  Phys. Rev. Lett., 99:060201, 2007.
  * [31] Martin Bitter and Valery Milner.  Experimental demonstration of coherent control in quantum chaotic systems.  Phys. Rev. Lett., 118:034101, 2017.
  * [32] Guillaume Vanhaele and Peter Schlagheck.  NOON states with ultracold bosonic atoms via resonance- and chaos-assisted tunneling.  Phys. Rev. A, 103:013315, 2021.
  * [33] Guillaume Vanhaele, Arnd Bäcker, Roland Ketzmerick, and Peter Schlagheck.  Creating triple-NOON states with ultracold atoms via chaos-assisted tunneling.  Phys. Rev. A, 106:L011301, 2022.
  * [34] Javier Madronero, Alexey Ponomarev, Andre R. R. Carvalho, Sandro Wimberger, Carlos Viviescas, Andrey Kolovsky, Klaus Hornberger, Peter Schlagheck, Andreas Krug, and Andreas Buchleitner.  Quantum chaos, transport, and control - in quantum optics.  Adv. in Atom. Mol. Opt. Phys., 53:33–73, 2006.
  * [35] Immanuel Bloch, Jean Dalibard, and Wilhelm Zwerger.  Many-body physics with ultracold gases.  Rev. Mod. Phys., 80:885–964, 2008.
  * [36] Jean Dalibard, Fabrice Gerbier, Gediminas Juzeliūnas, and Patrik Öhberg.  Artificial gauge potentials for neutral atoms.  Rev. Mod. Phys., 83:1523–1543, 2011.
  * [37] Immanuel Bloch, Jean Dalibard, and Sylvain Nascimbène.  Quantum simulations with ultracold gases.  Nature Physics, 8:267–276, 2012.
  * [38] Chih-Chun Chien, Sebastiano Peotta, and Massimiliano Di Ventra.  Quantum transport in ultracold atoms.  Nature Physics, 11:998–1004, 2015.
  * [39] Tim Langen, Remi Geiger, and Jörg Schmiedmayer.  Ultracold atoms out of equilibrium.  Annu. Rev. Condens. Matter Phys., 6:201–217, 2015.
  * [40] Florian Schäfer, Takeshi Fukuhara, Seiji Sugawa, Yosuke Takasu, and Yoshiro Takahashi.  Tools for quantum simulation with ultracold atoms in optical lattices.  Nature Reviews Physics, 2:411–425, 2020.
  * [41] Bing Yang, Hui Sun, Robert Ott, Han-Yi Wang, Torsten V. Zache, Jad C. Halimeh, Zhen-Sheng Yuan, Philipp Hauke, and Jian-Wei Pan.  Observation of gauge invariance in a 71-site Bose-Hubbard quantum simulator.  Nature, 587:392–396, 2020.
  * [42] Ehud Altman, Kenneth R. Brown, Giuseppe Carleo, Lincoln D. Carr, Eugene Demler, Cheng Chin, Brian DeMarco, Sophia E. Economou, Mark A. Eriksson, Kai-Mei C. Fu, Markus Greiner, Kaden R.A. Hazzard, Randall G. Hulet, Alicia J. Kollár, Benjamin L. Lev, Mikhail D. Lukin, Ruichao Ma, Xiao Mi, Shashank Misra, Christopher Monroe, Kater Murch, Zaira Nazario, Kang-Kuen Ni, Andrew C. Potter, Pedram Roushan, Mark Saffman, Monika Schleier-Smith, Irfan Siddiqi, Raymond Simmonds, Meenakshi Singh, I.B. Spielman, Kristan Temme, David S. Weiss, Jelena Vučković, Vladan Vuletić, Jun Ye, and Martin Zwierlein.  Quantum simulators: Architectures and opportunities.  PRX Quantum, 2:017003, Feb 2021.
  * [43] Christoph Braun, Raphaël Saint-Jalm, Alexander Hesse, Johannes Arceri, Immanuel Bloch, and Monika Aidelsburger.  Real-space detection and manipulation of topological edge modes with ultracold atoms.  arXiv:2304.01980 [cond-mat.quant-gas], 2023.
  * [44] Dolev Bluvstein, Simon J. Evered, Alexandra A. Geim, Sophie H. Li, Hengyun Zhou, Tom Manovitz, Sepehr Ebadi, Madelyn Cain, Marcin Kalinowski, Dominik Hangleiter, J. Pablo Bonilla Ataides, Nishad Maskara, Iris Cong, Xun Gao, Pedro Sales Rodriguez, Thomas Karolyshyn, Giulia Semeghini, Michael J. Gullans, Markus Greiner, Vladan Vuletić, and Mikhail D. Lukin.  Logical quantum processor based on reconfigurable atom arrays.  Nature, 2023.
  * [45] Lev P. Pitaevskii and Sandro Stringari.  Bose-Einstein Condensation.  Oxford University Press, Oxford, 2003.
  * [46] Andrey R. Kolovsky and Andreas Buchleitner.  Quantum chaos in the Bose-Hubbard model.  Europhys. Lett., 68:632, 2004.
  * [47] Lukas Pausch, Edoardo G. Carnio, Alberto Rodríguez, and Andreas Buchleitner.  Chaos and ergodicity across the energy spectrum of interacting bosons.  Phys. Rev. Lett., 126:150601, Apr 2021.
  * [48] Cheng Chin, Rudolf Grimm, Paul Julienne, and Eite Tiesinga.  Feshbach resonances in ultracold gases.  Rev. Mod. Phys., 82:1225–1286, 2010.
  * [49] Guo-Xian Su, Hui Sun, Ana Hudomal, Jean-Yves Desaules, Zhao-Yu Zhou, Bing Yang, Jad C. Halimeh, Zhen-Sheng Yuan, Zlatko Papić, and Jian-Wei Pan.  Observation of many-body scarring in a bose-hubbard quantum simulator.  Phys. Rev. Research, 5:023010, 2023.
  * [50] Alexander Impertro, Simon Karch, Julian F. Wienand, SeungJung Huh, Christian Schweizer, Immanuel Bloch, and Monika Aidelsburger.  Local readout and control of current and kinetic energy operators in optical lattices.  arXiv:2312.13268 [cond-mat.quant-gas], 2023.
  * [51] Alexandr MIklailovich Lyapunov.  The general problem of the stability of motion.  PhD thesis, University of Kharkov, 1892.  in Russian.
  * [52] A. N. Kolmogorov.  New metric invariant of transitive dynamical systems and endomorphisms of lebesgue spaces.  Doklady of Russian Academy of Sciences, 119(5):861–864, 1958.
  * [53] Yakov Grigor’evich Sinai.  On the notion of entropy of a dynamical system.  Doklady of Russian Academy of Sciences, 124:768–771, 1959.
  * [54] Jorge C. Leitão, João M. Viana Parente Lopes, and Eduardo G. Altmann.  Importance sampling of rare events in chaotic systems.  Eur. Phys. J. B, 90:181, 2017.
  * [55] Paul Ehrenfest.  Bemerkung über die angenäherte gültigkeit der klassischen mechanik innerhalb der quantenmechanik.  Zeit. der Phys., 45(7-8):455–457, 1927.
  * [56] Gennady P. Berman and George M. Zaslavsky.  Condition of stochasticity of quantum nonlinear systems.  Physica A, 91:450–460, 1978.
  * [57] George M. Zaslavsky.  Stochasticity in quantum systems.  Phys. Rep., 80:157–250, 1981.
  * [58] Klaus Richter, Juan Diego Urbina, and Steven Tomsovic.  Semiclassical roots of universality in many-body quantum chaos.  J. Phys. A: Math. Theor., 55:453001, 2022.  arXiv:2205.02867 [quant-phys].
  * [59] Steven Tomsovic.  Complex saddle trajectories for multidimensional quantum wave packet/coherent state propagation: application to a many-body system.  Phys. Rev. E, 98:023301, 2018.  arXiv:1804.10511 [cond-mat.stat-mech].
  * [60] Mathias Steinhuber.  Semiclassical approach to the high-dimensional self-trapping of interacting cold atoms in optical lattices.  Master’s thesis, Universität Regensburg, 2020.
  * [61] V. P. Maslov and M. V. Fedoriuk.  Semiclassical approximation in quantum mechanics.  Reidel Publishing Company, Dordrecht, 1981.
  * [62] Markus Greiner, Olaf Mandel, Theodor W. Hänsch, and Immanuel Bloch.  Collapse and revival of the matter wave field of a bose-einstein condensate.  Nature, 419:51–54, 2002.
  * [63] J. E. Moyal.  The monte carlo method.  Proc. Camb. Philos. Soc., 45:99–124, 1949.
  * [64] M. J. Steel, M. K. Olsen, L. I. Plimak, P. D. Drummond, S. M. Tan, M. J. Collett, D. F. Walls, and R. Graham.  Dynamical quantum noise in trapped bose-einstein condensates.  Phys. Rev. A, 58:4824–4835, 1998.
  * [65] Alice Sinatra, Carlos Lobo, and Yvan Castin.  The truncated wigner method for bose-condensed gases: limits of validity and applications.  J. Phys. B: Atom. Molec. Opt. Phys., 35:3599–3631, 2002.
  * [66] Eric J. Heller.  Time-dependent approach to semiclassical dynamics.  J. Chem. Phys., 62:1544–1555, 1975.
  * [67] Eric J. Heller.  Bound-state eigenfunctions of classically chaotic hamiltonian systems: Scars of periodic orbits.  Phys. Rev. Lett., 53(16):1515–1518, Oct 1984.
  * [68] Michel V. Berry.  Regular and irregular semiclassical wavefunctions.  J. Phys. A, 10:2083–2091, 1977.
  * [69] Andre Voros.  Semiclassical ergodicity of quantum eigenstates in the wigner representation.  In Giulio Casati and Joseph Ford, editors, Stochastic Behaviour in Classical and Quantum Hamiltonian Systems, pages 326–333. Springer-Verlag, Berlin, 1979.
  * [70] Steven W. McDonald and Allan N. Kaufman.  Spectrum and eigenfunctions for a hamiltonian with stochastic trajectories.  Phys. Rev. Lett., 42(18):1189–1191, Apr 1979.
  * [71] Steven William McDonald.  Wave dynamics of regular and chaotic rays.  PhD thesis, University of California, Lawrence Berkeley Laboratory, 1983.  [Report No. LBL-14837].
  * [72] Alexander I. Shnirelman.  Ergodic properties of eigenfunctions.  Usp. Math. Nauk., 29:181–182, 1974.
  * [73] Josh M. Deutsch.  Quantum statistical mechanics in a closed system.  Phys. Rev. A, 43:2046, 1991.
  * [74] Mark Srednicki.  Chaos and quantum thermalization.  Phys. Rev. E, 50:888–901, Aug 1994.
  * [75] Hannes Bernien, Sylvain Schwartz, Alexander Keesling, Harry Levine, Ahmed Omran, Hannes Pichler, Soonwon Choi, Alexander S. Zibrov, Manuel Endres, Markus Greiner, Vladan Vuletić, and Mikhail D. Lukin.  Probing many-body dynamics on a 51-atom quantum simulator.  Nature, 551:579–584, 2017.
  * [76] C. J. Turner, A. A. Michailidis, Dmitry A. Abanin, M. Serbyn, and Z. Papić.  Weak ergodicity breaking from quantum many-body scars.  Nature Physics, 14:745–749, 2018.
  * [77] Maksym Serbyn, Dmitry A. Abanin, and Zlatko Papić.  Quantum many-body scars and weak breaking of ergodicity.  Nature Physics, 17:675–685, 2021.
  * [78] Vedika Khemani, Chris R. Laumann, and Anushya Chandran.  Signatures of integrability in the dynamics of rydberg-blockaded chains.  Phys. Rev. B, 99:161101, 2019.
  * [79] Quirin Hummel, Klaus Richter, and Peter Schlagheck.  Genuine many-body quantum scars along unstable modes in bose-hubbard systems.  Phys. Rev. Lett., 130:250402, 2023.
  * [80] Ehud Altman.  Many-body localization and quantum thermalization.  Nature Physics, 14:979–983, 2018.
  * [81] John W. Negele and Henri Orland.  Quantum Many-particle Systems.  Advanced Books Classics. Taylor & Francis, 1998.
  * [82] Justin H. Wilson and Victor Galitski.  Breakdown of the coherent state path integral: Two simple examples.  Phys. Rev. Lett., 106:110401, Mar 2011.
  * [83] Falk Bruckmann and Juan Diego Urbina.  Rigorous construction of coherent state path integrals through dualization.  arXiv:1807.10462v1 [quant-ph], 2018.
  * [84] Lawrence S. Schulman.  Techniques and Applications of Path Integration.  John Wiley and Sons, Inc., New York, 1981.