Controlling many-body quantum chaos: Bose–Hubbard systems

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–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 figure 1. It is this extra requirement 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].

Zoom In Zoom Out Reset image size

Quantum control techniques comprise a large longstanding research field with much of the early work prompted by the dream of controlling chemical reactions [20–25], which more recently have been applied to many-body systems [26, 27]. A survey [28] and an overview of optimal control theory [29] 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. Nevertheless, there are a few works [30–36], 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 [37–46]. 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 [47] for a detailed review. This allows for a phase space formulation in which strongly chaotic regimes can be identified [48, 49] 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. For this work it is critical to be able to switch off the on-site interaction between the atoms [50–52], and control 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 (TWAs), and finally control protocol considerations. This is followed in section 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.

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, TWAs, 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 [53] 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_\textrm{min}$, the immediate neighborhood of trajectories that imitate this heteroclinic motion is of the order of $h^L\exp\left(-h_\textrm{KS}t_\textrm{min}\right)$, where $h_\textrm{KS}$ is Kolmogorov–Sinai entropy [54, 55] 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_\textrm{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, equation (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 [56] for a general technique developed precisely for searching exponentially rare solutions.

Of great interest here is how short the time scale $t_\textrm{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₀, 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_Lt_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_\textrm{KS}t\right)$. This grows to the order of N_L on a very short time scale given approximately by the relation

$\begin{equation} t_\textrm{min} \sim \frac{1}{h_\textrm{KS}}\ln N_L = \frac{1}{\lambda}\ln \frac{{\cal V}_1}{h} \ ,\end{equation} \tag{ 1 }$

where $\lambda = h_\textrm{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 [57–60] 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_Lt_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_Lt_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 [61]: 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 [62]. 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_\textrm{min}$ is independent of system size in the case of relying on discrete symmetries, just as equation (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_\textrm{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 sections 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 TWA 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$),

$\begin{align} \begin{split} \hat{H} = -J \sum_{j = 1}^{L} \left( {\hat{a}^{\dagger}_{j}}{\hat{a}_{j+1}} + \text{h.c.} \right) + \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}\end{align} \tag{ 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 µ_j ; see figure 2 for a simplified rendition of a 1D ring.

Zoom In Zoom Out Reset image size

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 [48, 49], 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_\textrm{eff}} = L/N$. As ${\hbar_\textrm{eff}} \rightarrow 0$, the mean field approximation gets better and better. A unique Hamiltonian (classical) dynamics arises for the mean field if the ratio

$\begin{equation} \gamma = \frac{1}{{\hbar_\textrm{eff}}}\frac{U}{J} \, ,\end{equation} \tag{ 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,

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

which re-expresses the Bose–Hubbard Hamiltonian as

$\begin{align} \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}\end{align} \tag{ 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,

$\begin{align} \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}\end{align} \tag{ 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 [63], 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

$\begin{align} \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)\ ,\end{align} \tag{ 7 }$

where the degree of chaos is determined by the value of γ. Typically, the strongest chaos is approximately in the range $1 \lt \gamma \lt 3$; see [48, 61]. The dynamics of the quadratures is thus governed by Hamilton's equation of motion

$\begin{align} \dot{\vec{x}} = \mathbb{\Sigma} \,\vec{\nabla}_{\vec{x}} \mathcal{H}\ , \quad \text{with } \vec{x} = \left(\vec{q},\vec{p}\right) \text{and } \mathbb{\Sigma} = \begin{pmatrix} 0 &\unicode{x1D7D9} \\ -\unicode{x1D7D9} & 0 \end{pmatrix}\ ,\end{align} \tag{ 8 }$

in a 2L-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

$\begin{equation} \hat{\vec{a}}{\vert{{\vec{\phi}}}\rangle_\textrm{coh}} = {{\vec{\phi}}}{\vert{{\vec{\phi}}}\rangle_\textrm{coh}}\ ,\end{equation} \tag{ 9 }$

they can be written explicitly as

$\begin{align} {\vert{{\vec{\phi}}}\rangle_\textrm{coh}} = \exp\left\{{-\frac{\lVert{{{\vec{\phi}}}}\rVert^{2}}{2} + {\vec{\phi}} \cdot {\hat{\vec{a}}^{\dagger}}}\right\} \vert 0\rangle\ .\end{align} \tag{ 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

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

centered around the mean-field point $\vec{\phi}_{\alpha} = (\vec{q}_{\alpha} + \textrm{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σ-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\unicode{x2A7D} \theta \lt 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

$\begin{align} {|{{\vec{\phi}}}\rangle_\textrm{proj}^{N}} = \frac{1}{\sqrt{\lVert{{{\vec{\phi}}}}\rVert^{2N} N!}} \left({\vec{\phi}} \cdot {\hat{\vec{a}}^{\dagger}}\right)^{N} \vert 0\rangle\ ,\end{align} \tag{ 12 }$

which contains a fixed number of particles, and is a Fock state in some basis; states of this kind can be constructed [64]. 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 α = 0 corresponding to the non-interacting ground state of the superfluid phase and α = 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 [65]

$\begin{align} \begin{split} \frac{\textrm{d}}{\textrm{d}t} \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}\end{align} \tag{ 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 equation (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 TWA [66, 67] 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

$\begin{equation} \frac{\textrm{d}}{\textrm{d}t} \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\}\ .\end{equation} \tag{ 14 }$

The Gaussian form for coherent states of equation (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{\ldots}\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 [62], which is an integrable system. The resulting trajectories exhibit periodic population inversion for γ < 4 with a period $\tau_\textrm{DW}(\gamma)$, see figure 3(a). For γ > 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, γ = 4, the period diverges as shown in figure 3. 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 τ_DW a natural choice for a dynamical time scale, as long as the systems considered are sufficiently removed from the singularity at γ = 4, as is done ahead. The corresponding time scales τ_DW of this work's specific cases are depicted in figure 3(b).

Zoom In Zoom Out Reset image size

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_\textrm{min}$, respectively, is identified, the initial goal is to guide the coherent state along this solution. Due to the quartic polynomial in quadratures present in the interaction term of equation (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.

Zoom In Zoom Out Reset image size

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 equation (7),

$\begin{align} \dot{q}_j & = \frac{\partial{\mathcal{H}}\left( \vec{q}, \vec{p} \right) }{\partial{p_j}}\biggr\rvert_{\vec{q}\left(t\right), \vec{p}\left(t\right)} = \left[\mu_j^{\prime} + \gamma \left(\frac{q_j^2+p_j^2}{2}\right) \right] p_j - p_{j+1} - p_{j-1} \, , \nonumber\\ \dot{p}_j & = -\frac{\partial\mathcal{H}\left(\vec{q}, \vec{p}\right)}{\partial{q_j}}\biggr\rvert_{\vec{q}\left(t\right), \vec{p}\left(t\right)} = -\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{align} \tag{ 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

$\begin{equation} \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\left(t\right) \left(\frac{q_{j}^{2} +p_{j}^{2}}{2}\right) \ ,\end{equation} \tag{ 16 }$

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

$\begin{equation} \mu_j\left(t\right) = \mu_j^{\prime} + \gamma \left(\frac{q_{j}^{2}\left(t\right) +p_{j}^{2}\left(t\right)}{2}\right) = \mu_j^{\prime} + \gamma n_j\left(t\right) \ .\end{equation} \tag{ 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,

$\begin{align} \begin{split} \hat{H}_{\text{c}}\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}\left(t\right)\hat{n}_j \right) \ , \end{split}\end{align} \tag{ 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 equation (18) and the precise form of the interaction (quartic) term in equation (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 [50–52] and chemical potential parameters, and to implement equation (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

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

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

$\begin{equation} K^\textrm{sc}\left({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha};t\right) = \sum_{\gamma} \frac{{\cal D}_{\gamma}\left({\vec{\phi}}_{\alpha},{\vec{\phi}}_{\beta},t\right)}{Z\left(\textbf{v}\right)} \textrm{e}^{\frac{i}{2}\int_{0}^{t}ds\vec{\sigma}_{\gamma}\left(s\right)\textbf{v}^{-1}\vec{\sigma}_{\gamma}\left(s\right)} \times K_{\vec{\sigma}_{\gamma}} \left({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha};t\right)\end{equation} \tag{ 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, equation (18), with control field $\vec{\mu}(t) = \vec{\sigma}_{\gamma}(t)$. As shown in appendix 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(\textbf{v})} \textrm{e}^{\frac{\textrm{i}}{2}\int_{0}^{t}\textrm{d}s\vec{\sigma}_{\gamma}(s)\textbf{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 ${\vert{{\vec{\phi}}(t)}\rangle_\textrm{coh}}$ (possibly number projected) centered on the mean-field solution used to drive the system satisfies the time-dependent Schrödinger equation, i.e.

$\begin{align} \textrm{i}\frac{\partial}{\partial t} {\vert{{\vec{\phi}}\left(t\right)}\rangle_\textrm{coh}} = \hat{H}_{\text{c}}\left(t\right) {\vert{{\vec{\phi}}\left(t\right)}\rangle_\textrm{coh}}\ .\end{align} \tag{ 21 }$

The explicit proof is provided in appendix B. Given the simple form of equation (18), the actual realization of the protocol requires switching off the interactions [50–52] and controlling the chemical potentials of each individual site in a time-dependent fashion. Some aspects of the protocol robustness to imperfections are addressed in section 3.2.

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 section 3.1 and the creation of non-trivial periodic evolution in section 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 section 3.2. Both the TWA and full quantum simulations of the time evolution are given.

3.1. Targeting

As discussed in section 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 equation (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 [52]. 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 γ and can be connected with a transport trajectory, see figure 5. Figure 6 shows a heteroclinic trajectory with a control time of less then five oscillations of the density wave τ_DW (figure 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. The initial state ${\vert{\vec{\phi}_{\alpha}}\rangle_\textrm{coh}}$ then evolves under the control Hamiltonian, equation (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 figure 6. The time-dependent site chemical potentials follow using equation (17) and the definition of $\mu_j^{\prime}$ given in the text immediately after equation (6). A TWA simulation (quantum calculation) is used to calculate the propagation of the state. Note that the TWA calculation exactly matches the quantum calculation due to the quadratic nature of the control Hamiltonian, and they cannot be distinguished in this case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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. These mechanisms are discussed below individually. Note that an analysis of the latter item is presumably model dependent, and amongst other issues, might depend on the type of noise inherent in creating the $\{\mu_j(t)\}$. Nevertheless, the effects of uncorrelated (white) noise in time are considered as a simple starting point.

3.2.1. Proximity

As mentioned at the end of section 2 and proved in appendix 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, section 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

$\begin{equation} \delta x_{\text{initial}} = ||\vec{x}_{\alpha} - \vec{x}\left(0\right)||, \qquad \delta x_{\text{target}} = ||\vec{x}_{\beta} - \vec{x}\left(t_c\right)||\ ,\end{equation} \tag{ 22 }$

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

$\begin{equation} {\vert{\vec{\phi}_\alpha\left(t_c\right)}\rangle_\textrm{coh}} : = \exp\left\{-\frac{\textrm{i}}{\hbar}\int\limits_0^{t_c} \textrm{d}t \hat{H}_{\text{c}}\left(t\right)\right\}{\vert{\vec{\phi}_{\alpha}}\rangle_\textrm{coh}} \ .\end{equation} \tag{ 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 figure 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:

$\begin{equation} \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}}}} \ .\end{equation} \tag{ 24 }$

The exact overlap is given as

$\begin{equation} \mathcal{F} = \left|_{\text{coh}}\langle \vec{\phi}_{\beta}\vert \textrm{e}^{-\frac{\textrm{i}}{\hbar}\int\limits_0^{t_c} \textrm{d}t \hat{H}_{\text{c}}\left(t\right)}\vert \vec{\phi}_{\alpha} \rangle_{\text{coh}}\right|^2 = \exp{-\frac{\delta\vec{x}^2}{2\hbar_{\text{eff}}}}\ ,\end{equation} \tag{ 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 [68], 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}}$.

Zoom In Zoom Out Reset image size

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 [50], imagine nevertheless that the strength of the interaction cannot be turned off 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, equation (18), as follows:

$\begin{align} \hat{H}_\textrm{c}^\mathrm{\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}\left(t\right)\hat{n}_j \right) + \frac{\epsilon\left(t\right)}{2} \sum_{j = 1}^{L}\hat{n}_j\left(\hat{n}_j-1\right) \ .\end{align} \tag{ 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. Both effects are illustrated in figure 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 equation (15), which were used to determine the $\{\mu_j(t)\}$, except with the inclusion of the ε term. Now resolve the equations as before, the result is the replacement,

$\begin{equation} \mu_j\left(t\right) \Longrightarrow \mu_j\left(t\right) - \epsilon\left(t\right)\end{equation} \tag{ 27 }$

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

Zoom In Zoom Out Reset image size

If ε is not compensated in the $\{\mu_j(t)\}$, the fidelity decay with ε shown in figure 10 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, equation (27), an overall improved robustness to the residual interactions is obtained (see again figure 10). 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.

Zoom In Zoom Out Reset image size

3.2.3. White noise

The final potential mechanism that might degrade the fidelity of the control protocol has to do with how perfectly the $\mu_j(t)$ can be produced experimentally. As a simple modeling starting point, imagine that there is a small amount of white (uncorrelated) noise generated in the $\mu_j(t)$. Consider adding to each individual $\mu_j(t)$ white noise (uncorrelated over time and site index) of some given root mean square (RMS) magnitude. This changes the mean field trajectory and its endpoint differs from the ideal. As a function of the white noise strength (i.e. the RMS of a white noise realization), the RMS distance δ_RMS between noisy and ideal trajectory endpoints are calculated for 50 realizations. The results are shown in figure 11.

Zoom In Zoom Out Reset image size

It turns out that the RMS deviation depends linearly on the strength of the white noise up to a level beyond which it saturates. The control trajectory is a stable solution and the perturbations are random with a vanishing mean. This gives a perturbed trajectory that stays close to the ideal. As the time evolution is integrated, there are significant cancellations of the uncorrelated deviations. The trajectory with noise only slowly diffuses away from the ideal.

For perturbations comparable to the phase space size, i.e. occupations ($\sim 1/\hbar$), the actual control fields $\mu_j(t)$ are more or less negligible and the motion is random, but nevertheless the deviation to the target point is bounded by the phase space size. Therefore, the curve saturates for strong enough white noise.

3.3. Stabilizing periodic trajectories

Nearly forty years ago, the concept of 'quantum scarring of eigenstates' was introduced by Heller [69]. 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 [70–73], 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 [74]. 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 [75, 76]. 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. Experimental effects seen in a 51-atom Rydberg chain, and recently even in tilted Bose–Hubbard systems, were introduced as many-body scarring [51, 77–79], 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 [80]. 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 [81]. With the keen current interest in many-body periodic or partially periodic dynamics [82, 83], 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 section 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 γ. One of these values is selected and illustrated in figure 12. For each cycle of population inversion, 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.

Zoom In Zoom Out Reset image size

For the same reasons that the fidelity of the control protocol is not perfect in figure 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 figure 13 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.

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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 [61, 84], 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 [62]. 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 [62]. 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 figure 15, that allow for the dynamics 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}}$ [62]. 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,

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

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

Zoom In Zoom Out Reset image size

The two target configurations shown in figure 16 result from the mappings of figure 15 applied to the 1D-ring control trajectory of section 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 figure 16(a) or a shifted stripe structure as shown in figure 16(b).

Zoom In Zoom Out Reset image size

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 section 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}}$.

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 usually be achieved to high precision, e.g. using Feshbach resonances as a basis for one possible method. 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 non-vanishing residual interaction can be characterized, then a simple rectification expressed in equation (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)\}$, may require appropriate modeling of noise or other imperfections in a particular physical setup, which is left for future consideration. However, in the case of white noise (uncorrelated fluctuations in the individual chemical potentials over time), the strength of this type of noise has a linear dependence on the distance deviation in the targeting protocol and the control protocol is not particularly sensitive.

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

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. M S is funded through a fellowship by the Studienstiftung des Deutschen Volkes.

The data that support the findings of this study are openly available at the following URL/DOI: https://epub.uni-regensburg.de/58463/ [85].

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 [86], 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

$\begin{equation} \int \frac{{\cal D}\left[x\right]}{Z\left(M\right)}\textrm{e}^{\frac{\textrm{i}}{2}\int_{0}^{t}\textrm{d}s \int_{0}^{t}\textrm{d}s^{^{\prime}} x\left(s\right)M\left(s,s^{^{\prime}}\right) x\left(s^{^{\prime}}\right)}\textrm{e}^{\textrm{i}\int_{0}^{t}\textrm{d}sx\left(s\right)y\left(s\right)} = \textrm{e}^{\frac{\textrm{i}}{2}\int_{0}^{t}\textrm{d}s \int_{0}^{t}\textrm{d}s^{^{\prime}} y\left(s\right)\left(M^{-1}\right)\left(s,s^{^{\prime}}\right) y\left(s^{^{\prime}}\right)}\end{equation} \tag{ A.1 }$

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

$\begin{equation} \int_{0}^{t} \textrm{d}s^{^{\prime}} M\left(s,s^{^{\prime}}\right)\left(M^{-1}\right)\left(s^{^{\prime}},s^{^{\prime\prime}}\right) = \delta\left(s-s^{^{\prime\prime}}\right),\end{equation} \tag{ A.2 }$

and Z(M) is a normalization constant defined by comparing both sides of equation (A.1) 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 \textbf{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,

$\begin{equation} K\left({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha};t\right): = \langle {\vec{\phi}}_{\beta}|\textrm{e}^{-\frac{\textrm{i}}{\hbar}\hat{H}t}|{\vec{\phi}}_{\alpha}\rangle\end{equation} \tag{ A.3 }$

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

$\begin{equation} K\left({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha};t\right) = \int{\cal D}\left[{\vec{\phi}},{\vec{\phi}}^{*}\right] \textrm{e}^{\textrm{i} \left(R_\textrm{free}\left[{\vec{\phi}},{\vec{\phi}}^{*}\right]+R_\textrm{int}\left[{\vec{\phi}},{\vec{\phi}}^{*}\right]\right)},\end{equation} \tag{ A.4 }$

and defined by the non-interacting (free) $R_\textrm{free}$ and interacting $R_\textrm{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

$\begin{align} \hat{V}_\textrm{int} = \frac{1}{2}\sum_{i,j}v_{ij}\hat{n}_{i}\left(\hat{n}_{j}-\delta_{ij}\right)\,\,\,\,\, \textrm{with}\,\,\,\,\, \hat{n}_{i} = \hat{a}^{\dagger}_{i}\hat{a}_{i},\,\,\, \textrm{and}\,\,\,\,\,v_{ij} = v_{ji}.\end{align} \tag{ A.5 }$

Therefore, neglecting issues of ordering in defining the classical action functional (see [87, 88] for details)

$\begin{equation} R_\textrm{int}\left[{\vec{\phi}},{\vec{\phi}}^{*}\right] = \frac{1}{2}\int_{0}^{t}\textrm{d}s \sum_{i,j}v_{ij} |\phi_{i}\left(s\right)|^{2} \left(|\phi_{j}\left(s\right)|^{2}-\delta_{ij}\right).\end{equation} \tag{ A.6 }$

As is well known, the reason why the path integral, equation (A.4), 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, equation (A.6), can be decoupled in favor of the (vector) real auxiliary field $\vec{\sigma}(s)$ obtaining

$\begin{equation} K\left({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha};t\right) = \int \frac{{\cal D}\left[\vec{\sigma}\right]}{Z\left(\textbf{v}\right)}\textrm{e}^{\frac{\textrm{i}}{2}\int_{0}^{t}\textrm{d}s\vec{\sigma}\left(s\right)\textbf{v}^{-1}\vec{\sigma}\left(s\right)}K_{\vec{\sigma}}\left({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha};t\right)\end{equation} \tag{ A.7 }$

where the reduced propagator

$\begin{equation} K_{\vec{\sigma}}\left({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha};t\right) = \int{\cal D}\left[{\vec{\phi}},{\vec{\phi}}^{*}\right] \textrm{e}^{\textrm{i} \left(R_\textrm{free}\left[{\vec{\phi}},{\vec{\phi}}^{*}\right]+\int_{0}^{t}\textrm{d}s\sum_{j}\sigma_{j}\left(s\right) n_{j}\left(s\right)\right)},\,\,\,\,\textrm{with} \,\,\, n_{j}\left(s\right) = |\phi_{j}\left(s\right)|^{2},\end{equation} \tag{ A.8 }$

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

Next, Gaussian path integrals, even with time-dependent parameters as in equation (A.8) 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,

$\begin{equation} K_{\vec{\sigma}}({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha};t) = A(t,\left[\vec{\sigma}\right])\textrm{e}^{\textrm{i} (R_\textrm{free}^\textrm{mf}({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},\left[\vec{\sigma}\right])+\int_{0}^{t}\textrm{d}s\sum_{j}\sigma_{j}(s) n_{j}^\textrm{mf}(s,{\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},\left[\vec{\sigma}\right])))}\end{equation} \tag{ A.9 }$

where $R^\textrm{mf},n_{j}^\textrm{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

$\begin{align} R_\textrm{free}^\textrm{mf}\left({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},t,\left[\vec{\sigma}\right]\right) = R_\textrm{free}\left[{\vec{\phi}}\left(s\right) = {\vec{\phi}}^\textrm{mf}\left(s,{\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},t,\left[\vec{\sigma}\right]\right),{\vec{\phi}}^{*}\left(s\right) = {\vec{\phi}}^\textrm{mf}\left(s,{\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},t,\left[\vec{\sigma}\right]\right)^{*}\right]\end{align} \tag{ A.10 }$

and

$\begin{equation} n_{j}^\textrm{mf}\left(s,{\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},t,\left[\vec{\sigma}\right]\right) = \left|\phi_{j}^\textrm{mf}\left(s,{\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},t,\left[\vec{\sigma}\right]\right)\right|^{2}\end{equation} \tag{ A.11 }$

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

$\begin{equation} \frac{\delta}{\delta {\vec{\phi}}^{*}}\left(R_\textrm{free}\left[{\vec{\phi}},{\vec{\phi}}^{*}\right]+\int_{0}^{t}\textrm{d}s\sum_{\alpha}\sigma_{\alpha}\left(s\right) n_{\alpha}\left(s\right)\right)_{{\vec{\phi}} = {\vec{\phi}}^\textrm{mf}} = 0\end{equation} \tag{ A.12 }$

satisfying

$\begin{equation} {\vec{\phi}}\left(s = 0\right) = {\vec{\phi}}_{\alpha}\,\,\,\,\,\textrm{and}\,\,\,\,\, {\vec{\phi}}\left(s = t\right) = {\vec{\phi}}_{\beta}.\end{equation} \tag{ A.13 }$

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 [89].

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

$\begin{equation} R_\textrm{free}\left[{\vec{\phi}},{\vec{\phi}}^{*}\right] = \int_{0}^{t}\textrm{d}s\sum_{j}\operatorname{Im} \left(\phi_{j}\left(s\right)\frac{\textrm{d}\phi_{j}^{*}\left(s\right)}{\textrm{d}s}\right)+ \int_{0}^{t} \textrm{d}s \sum_{i,j}h_{ij}\phi_{i}^{*}\left(s\right)\phi_{j}\left(s\right)\end{equation} \tag{ A.14 }$

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 equation (A.12) to be

$\begin{equation} i\frac{\textrm{d}}{\textrm{d}s}z_{i}\left(s\right) = \sum_{j}h_{ij}z_{j}\left(s\right)+\sigma_{i}\left(s\right)z_{i}\left(s\right),\end{equation} \tag{ A.15 }$

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,

$\begin{align} K({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha};t) = \int \frac{{\cal D}\left[\vec{\sigma}\right]}{Z(\textbf{v})} \textrm{e}^{\frac{\textrm{i}}{2}\int_{0}^{t}\textrm{d}s\vec{\sigma}(s)\textbf{v}^{-1}\vec{\sigma}(s)}A(t,\left[\vec{\sigma}\right]) \times \textrm{e}^{\textrm{i} (R_\textrm{free}^\textrm{mf}({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},t,\left[\vec{\sigma}\right])+\int_{0}^{t}\textrm{d}s\sum_{j}\sigma_{j}(s) n_{j}^\textrm{mf}(s,{\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},t,\left[\vec{\sigma}\right])))},\end{align} \tag{ A.16 }$

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

$\begin{align} &\frac{\delta}{\delta \vec{\sigma}}\left(\frac{1}{2}\int_{0}^{t}\textrm{d}s\vec{\sigma}\left(s\right)\textbf{v}^{-1}\vec{\sigma}\left(s\right) +R_\textrm{free}^\textrm{mf}\left({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},t,\left[\vec{\sigma}\right]\right)+ \left.\int_{0}^{t}\textrm{d}s\sum_{j}\sigma_{j}\left(s\right) n_{j}^\textrm{mf}\left(s,{\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},t,\left[\vec{\sigma}\right]\right)\right)\right|_{\vec{\sigma}_{\gamma}} = 0. \end{align} \tag{ A.17 }$

Using the implicit dependence of $R_\textrm{free}^\textrm{mf}$ and $n_{\alpha}^\textrm{mf}$ on the auxiliary field as given in equations (A.10) and (A.11) and the mean-field equations, equation (A.12), this gives

$\begin{equation} \textbf{v}^{-1}\vec{\sigma}_{\gamma}\left(s\right) = \vec{n}^\textrm{mf}\left(s,{\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},t,\left[\vec{\sigma}_{\gamma}\right]\right)\end{equation} \tag{ A.18 }$

or

$\begin{equation} \sigma_{\gamma,i}\left(s\right) = \sum_{j}v_{ij}n_{j}^\textrm{mf}\left(s,{\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},t,\left[\vec{\sigma}_{\gamma}\right]\right)\end{equation} \tag{ A.19 }$

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 (A.12) and (A.15) makes it natural to expect a countable, but large, number of solutions.

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

$\begin{equation} \sigma_{\gamma,i}\left(s\right) = \sum_{j}v_{ij}|\phi_{\gamma,j}^\textrm{clas}|^{2}\end{equation} \tag{ A.20 }$

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

$\begin{equation} i\frac{\textrm{d}}{\textrm{d}s}\phi_{i}\left(s\right) = \sum_{j}h_{ij}\phi_{j}\left(s\right)+\sum_{j}v_{ij}n_{i}\left(s\right) \phi_{j}\left(s\right)\end{equation} \tag{ A.21 }$

for the Bose–Hubbard type of action. Note that the fields $\vec{\sigma}_{\gamma}$ inherit the boundary conditions, equation (A.13), 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 equation (A.21).

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

$\begin{align} K^\textrm{sc}({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha};t)& = \sum_{\gamma} \frac{{\cal D}_{\gamma}({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},t)}{Z(\textbf{v})} \left.\textrm{e}^{\frac{\textrm{i}}{2}\int_{0}^{t}\textrm{d}s\vec{\sigma}(s)\textbf{v}^{-1}\vec{\sigma}(s)}\right|_{\vec{\sigma} = \vec{\sigma}_{\gamma}}\nonumber\\\ &\quad\times \left.A(t,\left[\vec{\sigma}\right])\textrm{e}^{\textrm{i} (R_\textrm{free}^\textrm{mf}({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},t,\left[\vec{\sigma}\right])+\int_{0}^{t}\textrm{d}s\sum_{j}\sigma_{j}(s) n_{j}^\textrm{mf}(s,{\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},t,\left[\vec{\sigma}\right])))}\right|_{\vec{\sigma} = \vec{\sigma}_{\gamma}}, \nonumber\end{align} \tag{ A.22 }$

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 equation (A.22) contains the reduced propagator for the quadratic action, equations (A.8) and (A.9) in its semiclassical (but exact) form. Thus,

$\begin{equation} K^\textrm{sc}\left({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha};t\right) = \sum_{\gamma} \frac{{\cal D}_{\gamma}\left({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha},t\right)}{Z\left(\textbf{v}\right)} \textrm{e}^{\frac{\textrm{i}}{2}\int_{0}^{t}\textrm{d}s\vec{\sigma}_{\gamma}\left(s\right)\textbf{v}^{-1}\vec{\sigma}_{\gamma}\left(s\right)} \times K_{\vec{\sigma}_{\gamma}} \left({\vec{\phi}}_{\beta},{\vec{\phi}}_{\alpha};t\right)\end{equation} \tag{ A.23 }$

where $K_{\vec{\sigma}_{\gamma}}$ is the exact coherent state propagator for the γ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 (A.23) 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.

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

$\begin{align*} {\vert{\vec{\phi}\left(t\right)}\rangle_\textrm{coh}} = \hat{U}_\textrm{c}\left(t\right) {\vert{\vec{\phi}_{0}}\rangle_\textrm{coh}} = \hat{\mathcal{T}}\exp\left\{-\textrm{i}\int_{0}^{t} \hat{H}_\textrm{c} \left(s\right){\operatorname{d}\!} s\right\} {\vert{\vec{\phi}_{0}}\rangle_\textrm{coh}},\end{align*}$

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.

$\begin{align*} \textrm{i}\frac{{\operatorname{d}\!} }{{\operatorname{d}\!} t} {\vert{\vec{\phi}\left(t\right)}\rangle_\textrm{coh}} = \hat{H}_\textrm{c}\left(t\right) {\vert{\vec{\phi}\left(t\right)}\rangle_\textrm{coh}}\ ,\end{align*}$

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

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

that evolve according to the equations motion

$\begin{align*} \textrm{i}\dot{\phi}_j & = \frac{\partial \mathcal{H}_{\text{c}}}{\partial \phi_j^*} = - \left\{\phi_{j+1} + \phi_{j-1} - \mu_j\left(t\right)\phi_{j}\right\}\\ \textrm{i}\dot{\phi}_j^* & = -\frac{\partial \mathcal{H}_{\text{c}}}{\partial \phi_j} = \left\{\phi_{j+1}^* + \phi_{j-1}^* - \mu_j\left(t\right)\phi_{j}^*\right\}.\end{align*}$

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

$\begin{align*} \textrm{i}\frac{{\operatorname{d}\!} }{{\operatorname{d}\!} t} {\vert{\vec{\phi}\left(t\right)}\rangle_\textrm{coh}} = \, \textrm{i} \textrm{e}^{-\lVert{{\vec{\phi}}}\rVert^{2}/2} \textrm{e}^{{\vec{\phi}}\left(t\right) \cdot {\hat{\vec{a}}^{\dagger}}} \left(\sum\limits_{j = 1}^L \dot{\phi}_j\left(t\right) \hat{a}_{j}^{\dagger} \right)\vert 0\rangle.\end{align*}$

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

$\begin{align*} \textrm{i}\frac{{\operatorname{d}\!} }{{\operatorname{d}\!} t} {\vert{\vec{\phi}\left(t\right)}\rangle_\textrm{coh}} & = -\textrm{i} \textrm{e}^{-\lVert{{\vec{\phi}}}\rVert^{2}/2} \textrm{e}^{{\vec{\phi}}\left(t\right) \cdot {\hat{\vec{a}}^{\dagger}}} \sum_{j = 1}^{L}\left(\left\{\phi_{j+1}\left(t\right)+\phi_{j-1}\left(t\right) - \mu_{j}\left(t\right) \phi_{j}\left(t\right)\right\}\hat{a}_{j}^{\dagger}\right)\vert 0\rangle\\ & = -\sum_{j = 1}^{L}\left(\left\{\phi_{j+1}\left(t\right)+\phi_{j-1}\left(t\right) - \mu_{j}\left(t\right) \phi_{j}\left(t\right) \right\}\hat{a}_{j}^{\dagger} \right){\vert{\vec{\phi}\left(t\right)}\rangle_\textrm{coh}}.\end{align*}$

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

$\begin{align*} \hat{H}_\textrm{c}\left(t\right) {\vert{\vec{\phi}\left(t\right)}\rangle_\textrm{coh}} & = -\sum_{j = 1}^{L} \left(\left\{\hat{a}_{j+1}^{\dagger} + \hat{a}_{j-1}^{\dagger} + \mu_j\left(t\right)\hat{a}_{j}^{\dagger}\right\} \hat{a}_{j} \right) {\vert{\vec{\phi}\left(t\right)}\rangle_\textrm{coh}}\\ & = -\sum_{j = 1}^{L} \left(\hat{a}_{j}^{\dagger}\left\{\phi_{j+1}\left(t\right)+\phi_{j-1}\left(t\right) + \mu_{j}\left(t\right) \phi_j\left(t\right)\right\}\right) {\vert{\vec{\phi}\left(t\right)}\rangle_\textrm{coh}},\end{align*}$

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

Remark 1. The proof for number projected states is analogous using the property

$\begin{equation*} \hat{a}_{j}{|{\vec{\phi}\left(t\right)}\rangle_\textrm{proj}^{N}} = \phi_{j}|{\vec{\phi}\left(t\right)}\rangle_\textrm{proj}^{N-1}.\end{equation*}$

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 equation (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 [68], it would be possible to provide an analytical formula for the fidelity, for example shown in figure 8 based solely on the knowledge of the initial state and the control trajectory.