Controlling Many-Body Quantum Chaos: Bose-Hubbard systems

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


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 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 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.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], 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,31,32,33,34,35,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,38,39,40,41,42,43,44,45,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,51,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 numberprojected 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 timedependent 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.

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.

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 min , the immediate neighborhood of trajectories that imitate this heteroclinic motion is of the order of h L exp (−h KS t min ), where h 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 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. [56] 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, V 1 , so that the number of quantum coarse grained cells for the full system is 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 (h KS t).This grows to the order of N L on a very short time scale given approximately by the relation where λ = 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 [57,58,59, 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 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 [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 reduceddimensional 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.

Bose-Hubbard systems
For many reasons the Bose-Hubbard model turns out to be an ideal bosonic manybody 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 ≡ 1), interaction hopping chem.pot.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 Fig. 2 for a simplified rendition of a 1D ring.
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 ∼ U N , 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, ℏ eff = L/N .As ℏ eff → 0, the mean field approximation gets better and better.A unique Hamiltonian (classical) dynamics arises for the mean field if the ratio γ = 1 is held fixed.It is possible to introduce Hermitian operators in terms of âj and â † j known as quadratures that obey position-momentum-like commutation relations, which re-expresses the Bose-Hubbard Hamiltonian as Here an irrelevant constant energy offset has been dropped (which is also effectively an O(ℏ 2 eff ) correction to the spectrum and, in any case, beyond semiclassical argumentation used ahead).Multiplying both sides of the equation by ℏ eff /J leads to a modified Hamiltonian, where µ ′ j = (µ 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(ℏ 2 eff ), just through the association of q → q and p → p.The classical analog Hamiltonian (mean field limit) is given by where the degree of chaos is determined by the value of γ.Typically, the strongest chaos is approximately in the range 1 < γ < 3; see [48,61].The dynamics of the quadratures is thus governed by Hamilton's equation of motion 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 µ ′ j 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.

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 they can be written explicitly as Their phase space representation is given by the Wigner function and is well known to be a simple symmetric 2L-dimensional Gaussian wave packet centered around the mean-field point ⃗ ϕ α = (⃗ q α +i⃗ p α )/ √ 2ℏ eff .The coherent state variance per degree of freedom is given by σ 2 = ℏ eff /2, which for a single degree of freedom has the 2σ-contour enclosing a volume 2πℏ eff .More generally, the volume per quantum state is (2πℏ 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 ≤ θ < 2π.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 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 ϕ α j = exp (2πiαj/L) 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., ϕ j = √ Lδ j,k .

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] ∂⃗ 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 ℏ eff yields the truncated Wigner approximation (TWA) [66,67] in which quantum zero-point fluctuations are encoded in Except for γ = 4 where the period diverges since ⃗ ϕ lies on a separatrix of the mean field dimer, τ DW provides a well defined time scale of the system.The time scale is marked for all the specific trajectories used throughout the work.
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 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.

Control Hamiltonian
Given that a (presumably heteroclinic) control trajectory (⃗ q(t), ⃗ p(t)) 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.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 ℏ 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, Ĥ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), Evaluated using the control trajectory, (⃗ q(t), ⃗ p(t)) 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 which is in essence a harmonic oscillator with chemical potentials that are timedependent functions of the control trajectory's site occupancies, The control Hamiltonian makes use of the chemical potentials µ j (t) as purely timedependent functions evaluated along the chosen control (heteroclinic) trajectory A minimum uncertainty coherent state (here for illustration chosen to be ℏ eff = 1/10), initially centered on the density wave ⃗ ϕ α = 2/ℏ eff (1, 0, 1, 0) (indicated as dot in the upper panel) is propagated in the Bose-Hubbard system (γ ≃ 1.965) using the TWA.With the scaled dynamics that is used, the mean-field corresponds to the phase space point ⃗ q α = (2, 0, 2, 0) and ⃗ p α = (0, 0, 0, 0), independent of ℏ eff .The state, whose volume is for t = 0 well contained in the 2σ-contour (white dashed circle), quickly spreads, already being far from its initial Gaussian shape after one period of the density wave.At t = 4τ DW the state is fully equilibrated.
(⃗ q(t), ⃗ 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, (⃗ q(t), ⃗ p(t)) c is a solution of Hamilton's equations for both Hamiltonians H (⃗ q, ⃗ p) and H c (⃗ q, ⃗ p, t).In addition, H c (⃗ q, ⃗ p, t) having a quadratic form leads to the TWA being exact for its quantized version, and there is no longer any spreading in the neighborhood of (⃗ q(t), ⃗ p(t)) c , or for that matter, anywhere else.Although, (⃗ q(t), ⃗ p(t)) c is an exponentially unstable solution of H (⃗ q, ⃗ p), it is a stable solution of H c (⃗ q, ⃗ p, t).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 [50,51,52] 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 with Ĥ given in Eq. ( 2), in terms of auxiliary fields valid in the limit of large mean occupations N ≫ 1.Here K ⃗ σγ ( ⃗ ϕ β , ⃗ ϕ α ; t) is the transition amplitude between initial | ⃗ ϕ α ⟩ and target | ⃗ ϕ α ⟩ coherent states defined by the control Hamiltonian, Eq. ( 18), with control field ⃗ µ(t) = ⃗ σ γ (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 , over the set of control fields ⃗ σ γ (s) obtained by σ γ,j (s) = ϕ γ,j (s)ϕ * γ,j (s) from the solutions of the mean field boundary problem ⃗ ϕ γ (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 Ĥc (t), implies that the coherent state | ⃗ ϕ(t)⟩ coh (possibly number projected) centered on the mean-field solution used to drive the system satisfies the time-dependent Schrödinger equation, i.e.
The explicit proof is provided in Appendix B. Given the simple form of Eq. ( 18), the actual realization of the protocol requires switching off the interactions [50,51,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 Sec.3.2.

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, (⃗ q(t), ⃗ p(t)) c , are not taken into account.The resulting effects as well as imperfections in making the residual interactions vanish both introduce an ℏ 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.

Targeting
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 meanfield 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 As already discussed, the protocol is initialized with a coherent state centered on a density wave ⃗ ϕ α = 2/ℏ eff (1, 0, 1, 0).By choosing the chemical potential offset accordingly, i.e. ⃗ µ ′ = (0, γ, 0, γ), both mean-fields share the same classical energy surface for every value of γ 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 τ DW (Fig. 3b), 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 | ⃗ ϕ α ⟩ 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 µ ′ j given in the text immediately after Eq. ( 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.
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.

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 {µ 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 {µ j (t)}.Nevertheless, the effects of uncorrelated (white) noise in time are considered as a simple starting point.

Proximity
As mentioned at the end of Sec. 2 and proved in Appendix B, the specific quadratic nature of Ĥ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 resulting in a non-perfect overlap between the target state | ⃗ ϕ β ⟩ coh and the actual time evolved final state 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 ℏ eff dependency gets introduced into the fidelity 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 δx initial and δx target bounds can be derived for the fidelity of the protocol: The exact overlap is given as where δ⃗ 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 δ⃗ 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 ℏ eff .

Residual interactions
The second mechanism mentioned above that could potentially degrade the coherent quantum targeting is related to switching off the onsite 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, Eq. ( 18), as follows: 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 (⃗ q(t), ⃗ p(t)) c would result in missing the desired target.Both effects are illustrated in Fig. 9 for a constant residual interaction ϵ(t) = ϵ, as well as the improvement using the appropriately modified {µ j (t)}.If the behavior of ϵ(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 {µ j (t)}, except with the inclusion of the ϵ term.Now resolve the equations as before, the result is the replacement, in Eq. (26).It suffices to subtract ϵ(t) from the {µ j (t)} to counteract the over or undershooting of the target.If ϵ is not compensated in the {µ j (t)}, the fidelity decay with ϵ is more significant as ℏ eff → 0. This indicates that the dominant error source is the perturbed control 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 ϵ increases.The spreading also increases with ϵ as illustrated for the 1σ-contour of the propagated Wigner density for ℏ eff = 1 (black) and ℏ eff = 1/4 (gray).For very small ϵ (= 0.01), both effects remain negligible.However for ϵ = 0.1, the control trajectory completely overshoots the target and there is a strong deformation of the coherent state.Accounting for ϵ in the µ 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σ-contour being only a cut through a higher-dimensional object.
trajectory not arriving at the target.It follows from the above argument that a nonzero 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 Ĥc by adjusting the driving by the replacement, Eq. ( 27), an overall improved robustness to the residual interactions is obtained.Moreover, notably, the ℏ eff = 1/4 coherent state is less affected than the ℏ eff = 1 one.The remaining error source is the deformation of the state, which is less significant for smaller volumes.

White noise
The final potential mechanism that might degrade the fidelity of the control protocol has to do with how perfectly the µ 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 µ j (t).Consider adding to each individual µ 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 Fig. 11.
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.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 ℏ eff case, due to the deformation being the remaining error.Since TWA is no longer exact with the interaction term present, the overlaps have been calculated in a full quantum simulation.The orange triangles mark the cases depicted in Fig. 9.An RMS average is calculated over 50 noise realizations for the final distance to the unperturbed trajectory.There is a linear dependency for low amplitudes due to the stability of the target protocol and cancellation of random changes.For strong noise magnitude (dominant with respect to the control µ j (t)), the final point is random and the distance is bounded by the phase space diameter.The errorbars display the run-to-run standard variance.
For perturbations comparable to the phase space size, i.e., occupations (∼ 1/ℏ), the actual control fields µ 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.

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,71,72,73], i.e., δ(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,78,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 manybody 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.

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 γ.One of these values is selected and illustrated in Fig. 12.For each cycle of population inversion, the phase advances by 2π/3, which after the third cycle completes the periodic trajectory.Note the simple appearance, similar to a Lissajous For γ = 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.
Figure 13: Quantum recurrences in the controlled system.When propagated in the uncontrolled system, no recurrences appear.For ℏ 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 ℏ eff = 1/64 case, together with more imperfect recurrences.Naturally, the smaller volume of the ℏ eff = 1/64 case also generates a sharper peak.
figure, of the quadrature variables for the trajectory.In essence, the initially empty site is π/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 ℏ eff .This is illustrated in Fig. 13 for the two cases ℏ 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 ℏ eff = 1 case falling off slowly with increased time, and the ℏ 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 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 γ = 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.periodic mean-field trajectory.

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. 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 ) is mapped via the patterns (a) or (c) to a 2D mean-field solution with a 2D Hamiltonian H 2D (J, U, ⃗ µ).

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 [84,61], 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 lowerdimensional 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 nearestneighbor 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 Fig. 15, that allow for the dynamics of a periodic 4 × 4 lattice to be mapped onto the dynamics of a four-site ring.In fact, the mappings work for any 4n × 4m 2D lattice with periodic boundary conditions (including n → ∞).Here the mean-field dynamics reduce to the four-site ring with double the hopping parameter J 1D = 2J 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, µ ⃗ R (t), corresponding to those of the reduced 1D system, ⃗ µ(t).The 2D-lattice control Hamiltonian, 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. 16 result from the mappings of Fig. 15 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. 16a or a shifted stripe structure as shown in Fig. 16b.
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 1D = 3J 3D .

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 meanfield/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, timedependent 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 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 {µ 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 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.
As is well known, the reason why the path integral, Eq. (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, Eq. (A.6), can be decoupled in favor of the (vector) real auxiliary field ⃗ σ(s) obtaining where the reduced propagator now has the key property of representing a free field evolving under the time-dependent external field σ j (s) coupled at the j-th site to the local occupations, i.e., a timedependent local chemical potential.Next, Gaussian path integrals, even with time-dependent parameters as in Eq. (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, where R mf , n mf j are now functionals of the auxiliary fields ⃗ σ, as indicated by [.], and functions of t, ⃗ ϕ β , ⃗ ϕ α .These objects are obtained by Finally, the van Vleck-Morette determinant A(t, [⃗ σ]) can be shown to be independent of the boundary conditions ⃗ ϕ β , ⃗ ϕ α for quadratic actions, and its precise form can be found in [88].
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. (A.5), the action functional contains a free term of the form and therefore the mean field solutions that fully determine the reduced propagator for a given auxiliary field ⃗ σ(s) are easily obtained from the variation in Eq. (A.12) to be 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, ) 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, [σ]) are smooth, the full propagator turns out to be dominated by the configurations, ⃗ σ = ⃗ σ γ , satisfying the corresponding vanishing of the first variation, namely Using the implicit dependence of R mf free and n mf α on the auxiliary field as given in Eqs.(A.10, A.11) and the mean-field equations, Eq. (A.12), this gives as the set of equations that determines the dominant auxiliary fields ⃗ σ γ .The extreme nonlinearity of this system resulting from the role of ⃗ σ as time-dependent contributions to the couplings in the linear equations (A.12, A.15) makes it natural to expect a countable, but large, number of solutions.
The implicit problem in Eq. (A. 19) can be made more familiar by invoking the dependence of the mean field occupations n mf α on ⃗ σ as given precisely by the meanfield equations (A.12).Substitution of Eq. (A. 19) into the equations obtained from Eq. (A.12) then results in the explicit solution for the Bose-Hubbard type of action.Note that the fields ⃗ σ γ inherit the boundary conditions, Eq. (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 Eq. (A.21).
The importance of this representation is, however, clear upon return to the semiclassical evaluation of the exact integral in Eq. (A.16), given now within the saddle point approximation by where K ⃗ σγ 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 ⃗ ϕ α 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 ⃗ ϕ α with ⃗ ϕ β or other physical considerations, a single one of the auxiliary fields.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 ∥ ⃗ ϕ∥ 2 = N = const.â † j {ϕ j+1 (t) + ϕ j−1 (t) + µ j (t)ϕ j (t)} | ⃗ ϕ(t)⟩ coh , which matches exactly the expression coming from the left hand side.

Remark 2
The control trajectory used to drive the chemical potentials {µ j (t)} 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 timedependent orthogonal matrix.This implies that for a coherent state initially centered off the control trajectory by ||δ⃗ x 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 Fig. 8 based solely on the knowledge of the initial state and the control trajectory.

Figure 1 :
Figure1: 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.

Figure 2 :
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, onsite interactions between atoms, as well as energetic offsets for each site given by individual chemical potentials.

Figure 3 :
Figure 3: Properties of a density wave.A density wave initially centered around ⃗ ϕ = 2/ℏ eff (1, 0, 1, 0) exhibits partial or full periodic population inversion with a period τ DW as illustrated in (a).How the period depends on dynamical parameter γ, Eq. (3), of the system is shown in (b).Except for γ = 4 where the period diverges since ⃗ ϕ lies on a separatrix of the mean field dimer, τ DW provides a well defined time scale of the system.The time scale is marked for all the specific trajectories used throughout the work.

Figure 4 :
Figure4: Spreading of a localized coherent state in a Bose-Hubbard system.A minimum uncertainty coherent state (here for illustration chosen to be ℏ eff = 1/10), initially centered on the density wave ⃗ ϕ α = 2/ℏ eff (1, 0, 1, 0) (indicated as dot in the upper panel) is propagated in the Bose-Hubbard system (γ ≃ 1.965) using the TWA.With the scaled dynamics that is used, the mean-field corresponds to the phase space point ⃗ q α = (2, 0, 2, 0) and ⃗ p α = (0, 0, 0, 0), independent of ℏ eff .The state, whose volume is for t = 0 well contained in the 2σ-contour (white dashed circle), quickly spreads, already being far from its initial Gaussian shape after one period of the density wave.At t = 4τ DW the state is fully equilibrated.

Figure 5 :
Figure5: Schematic of initial and target quantum state.Raising the unoccupied sites energetically by choosing ⃗ µ = (0, γ, 0, γ) 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.

Figure 6 :
Figure 6: Control Trajectory to target excited momentum eigenstate.The trajectory connecting the centroids of the initial density wave ⃗ ϕ α = 2/ℏ eff (1, 0, 1, 0) (red dot) with the final target state ⃗ ϕ β = 1/ℏ eff (1, i, −1, −i) (red cross) is shown for the different sites in the left panel (γ ≃ 1.965).The respective 1σ-contours for the initial and target coherent state's Wigner densities, Eq. (11), are depicted centered on their respective quadrature centroids with ℏ eff = 1 (black dashed circles ) and ℏ eff = 1/4 (gray dashed circles).Generally, the area inside the circles shrinks proportionally to ℏ 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 µ ′ j 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τ 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.

Figure 7 :
Figure 7: Coherent targeting along chaotic control trajectory.A coherent state (here for illustration ℏ eff = 1/10 was chosen) is initialized centered on the density wave ⃗ ϕ α = 2/ℏ eff (1, 0, 1, 0) (white dot).The inital phase space point is connected to the target ⃗ ϕ β = 1/ℏ 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σ-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τ DW the propagated state lands well inside the volume of the target state.

Figure 8 :
Figure8: 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 ℏ 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.

Figure 9 :
Figure9: 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 ϵ increases.The spreading also increases with ϵ as illustrated for the 1σ-contour of the propagated Wigner density for ℏ eff = 1 (black) and ℏ eff = 1/4 (gray).For very small ϵ (= 0.01), both effects remain negligible.However for ϵ = 0.1, the control trajectory completely overshoots the target and there is a strong deformation of the coherent state.Accounting for ϵ in the µ 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σ-contour being only a cut through a higher-dimensional object.

Figure 10 :
Figure 10: Partial compensation of constant residual interactions.For a constant ϵ(t) = ϵ, the fidelity of an uncompensated protocol quickly decays.The ℏ eff = 1/4 decays much faster than ℏ eff = 1.The over or undershooting of the target is the dominant difficulty.If the behavior of ϵ(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 ℏ eff case, due to the deformation being the remaining error.Since TWA is no longer exact with the interaction term present, the overlaps have been calculated in a full quantum simulation.The orange triangles mark the cases depicted in Fig.9.

Figure 11 :
Figure 11: Final deviation of the target due to white noise.An RMS average is calculated over 50 noise realizations for the final distance to the unperturbed trajectory.There is a linear dependency for low amplitudes due to the stability of the target protocol and cancellation of random changes.For strong noise magnitude (dominant with respect to the control µ j (t)), the final point is random and the distance is bounded by the phase space diameter.The errorbars display the run-to-run standard variance.

Figure 12 :
Figure12: Symmetric periodic trajectory.For γ = 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.

Figure 14 :
Figure14: 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 γ = 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.

Figure 15 :
Figure 15: 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 1D = 2J 2D .In greater detail, a 1D mean-field solution ⃗ ϕ 1D = (ϕ 1 , ϕ 2 , ϕ 3 , ϕ 4 )) with a 1D Hamiltonian H 1D (2J, U, ⃗ µ/2) is mapped via the patterns (a) or (c) to a 2D mean-field solution with a 2D Hamiltonian H 2D (J, U, ⃗ µ).

Figure 16 :
Figure 16: 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. 15a or Fig. 15c 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.