Invariant-based inverse engineering of time-dependent, coupled harmonic oscillators

Two-dimensional systems with time-dependent controls admit a quadratic Hamiltonian modelling near potential minima. Independent, dynamical normal modes facilitate inverse Hamiltonian engineering to control the system dynamics, but some systems are not separable into independent modes by a point transformation. For these"coupled systems"2D invariants may still guide the Hamiltonian design. The theory to perform the inversion and two application examples are provided: (i) We control the deflection of wave packets in transversally harmonic waveguides; and (ii) we design the state transfer from one coupled oscillator to another.

Introduction. Controlling the motional dynamics of quantum systems is of paramount importance for fundamental science and quantum-based technologies [1]. Often the external driving needs to be fast, but also gentle, to avoid excitations. Slow adiabatic driving is gentle in this sense, but it exposes the system for long times to control noise, heating, and perturbations. Shortcuts to adiabaticity (STA) are techniques to reach, via fast non-adiabatic routes, the results of slow adiabatic processes [2,3]. A distinction can be made between: STA methods that keep the structure of some Hamiltonian form and design the time dependence of the controls, e.g. using invariants [4]; and those techniques that add new terms, e.g. counterdiabatic driving [5]. Both may be useful depending on system-dependent practical considerations. A frequent problem with added terms is the difficulty to implement them, whereas a limitation of structure-preserving, invariant-based methods is that they need Hamiltonian-invariant pairs with specific forms, such as the Lewis-Leach family of Hamiltonianinvariant pairs [6], to go beyond brute-force parameter optimization [2,3].
The eigenvectors of Lewis-Riesenfeld "time-dependent invariants" [7], with appropriate phase factors, are independent solutions of the Schrödinger equation and span a basis to expand any solution with constant expansion coefficients. These invariants are useful to inverse engineer the Hamiltonian and drive some desired dynamics [4]. The multiplicity of solutions for the trajectories of the control parameters, allows for adjustments or optimization with respect to different objectives or cost functions [8]. The multiplicity is also very helpful when several oscillators have to be controlled simultaneously [9,10].
This work extends the domain of systems that can be controlled by invariant-based inverse engineering. We shall deal with two-dimensional (2D) systems with quadratic Hamiltonians, found in particular in smalloscillation regimes of ultracold atom physics. In fact quadratic Hamiltonians are ubiquitous as they represent the systems near potential minima [11]. For timeindependent Hamiltonians the dynamics is simple to de-scribe and, possibly, manipulate by finding normal modes for effective uncoupled oscillators. This decomposition though, may not be possible if the Hamiltonian parameters depend on time. Lizuain et al. [12] described the condition for which a point transformation of coordinates decouples the instantaneous modes leading to truly independent "dynamical normal modes" [9] for two timedependent harmonic oscillators: the principal axes of the potential should not rotate in the 2D space.
When the two dynamical-mode motions separate, inverse engineering the dynamics to perform some fast operation free from final excitations is relatively easy: each of the time-dependent effective oscillators implies a one-dimensional Hamiltonian-invariant "Lewis-Leach" pair [6] for which inverse engineering can be performed. The two oscillators have to be driven simultaneously with common controls but, among the plethora of parameter trajectories, it is possible to find the ones that satisfy simultaneously the boundary conditions imposed on both oscillators. This strategy has been successfully applied to design the driving of different operations on two trapped ions such as transport or expansions [9,10], separation of two equal ions in double wells [13], phase gates [14], or dynamical exchange cooling [15].
If the effective potential rotates, the motions do not separate, so inverse engineering the external driving cannot in principle be done using two independent 1D Hamiltonian-invariant pairs. Solutions to the ensuing control problem exist that depend on the system and/or the operation, such as taking refuge in a perturbative regime [14], adding terms to cancel the inertial effects [12], increasing the number of time-dependent controls to uncouple the modes [15], or using more complex, nonpoint transformations to find independent modes [16].
Here we explore instead the use of 2D dynamical invariants associated with the coupled Hamiltonian.
Hamiltonian model. Consider the Hamiltonian We use throughout dimensionless variables such that no mass factors or appear explicitly. Eq. (1) describes different physical systems, such as a single particle in a 2D potential, or two coupled harmonic oscillators on a line. Other systems different from (one or two) particles but driven by Hamiltonians of the form (1) are, e.g., coupled superconducting qubits [17][18][19][20][21] or optomechanical oscillators [22][23][24]. All these systems are analogous to each other but, arguably, the single particle in a 2D potential is easiest to visualize so we shall use a terminology (such as longitudinal and transversal directions for principal axes, rotations...) borrowed from that system. Indeed, our first example, see below, deals with a single particle.
The Hamiltonian (1) may be instantaneously diagonalized by "rotated" variables [12] where , and Subscripts l and t stand for "longitudinal" and "transversal". The original Hamiltonian, expressed in terms of the new variables, is The formal decoupling in Eq. (4) is a mirage. H is not the Hamiltonian that describes the dynamics in the rotated variables {p l , p t , q l , q t } [12,25]. In general the dependence of A(t) on time couples dynamically the "instantaneous normal modes", i.e., the normal modes that would separate the motion if the Hamiltonian kept for all times the values that the parameters have at a particular instant. In the moving frame the oscillators are coupled by a term proportional toθ = dθ/dt [12]. Some peculiar, but physically significant relations between ω 1 (t), ω 2 (t), and γ(t) can make θ(t) time independent. Here we consider instead the scenario where θ(t) changes with time. This is unavoidable if the process we want to implement implies boundary conditions for the parameters such that θ(0) = θ(t f ), as in the examples below.
2D Invariant. Urzúa et al. [11], generalizing previous results in 1D [26,27] and the work in [28] for classical coupled oscillators, see also [29], have recently found that the linear combination of operators (dots stand for time derivatives hereafter) which are classical equations of motion driven by a Hamiltonian (1).

For any state driven by H(t), G(t) is the sum of two Wronskians
, where all functions in their arguments evolve as Eq. (7). The geometrical meaning of W i (t) is an "oriented" phase-space area formed by phase-space We consider two phase spaces, i = 1, 2, one for each oscillator. W i (t) is plus or minus the triangle area A i (t) depending on whether going from U i to Q i needs an anticlockwise or clockwise displacement. For γ = 0, the two areas (and Wronskians) remain constant in time. When γ = 0 the individual Wronskians are not conserved. The conserved quanti- , the initial phase-space oriented areas. The added terms cancel each other, namely, is the sum of oriented areas and it is constant. This result is a particular case of the preservation of sums of oriented areas in classical Hamiltonian systems [30].
We construct from G a quadratic invariant that may become proportional to some relevant energy at boundary times by choosing specific boundary conditions for the u i andu i , Designing the u i we may manipulate the invariants and therefore the dynamics. From the u i we can as well get the Hamiltonian as demonstrated in the following two examples.
Controlled deflection. A single particle is launched along a potential "waveguide" which is harmonic in the transversal direction. Our goal is to deflect it, that is, manipulate the potential to change the waveguide direction, controlling the input/output scaling factor of the longitudinal velocity. To have waveguide potentials at the boundary times t b = 0, t f we impose As a consequence, Ω l (t b ) = 0 and Ω t ( Thus, at boundary times, the potential is a harmonic "waveguide" with longitudinal direction defined by the angle can take any value between 0 and π/2 for θ(t f ) ≥ θ(0). The condition (8) in Eq. (7) implies thatü 1,2 (t b ) = 0, which also gives i.e., the reference trajectories must start and end at q t (t b ) = 0, on the axis of the waveguide. If the frequencies at t b are fixed, either q l (t b ), or one of the u i (t b ) can still be chosen freely.
Rewriting the invariant G in terms of the rotated variables {q t , q l } and imposingu 1,2 (t b ) = 0 we find that Initial waveguide Final waveguide γ const.   Table I for values at t = t f .
i.e., I(t b ) is proportional to the longitudinal energy. With Eq. (10) we get where F = u2 (0) u2(t f ) sin θ(t f ) sin θ(0) and E l = p 2 l /2 . For some chosen deflection angle ∆θ and waveguide frequencies Ω t (t b ) we may impose any scaling factor by manipulating the ratio u 2 (0)/u 2 (t f ). This scaling factor will affect all wave packets. Deflection angle, velocity scaling, and waveguide compression/expansion factors can be chosen independently.
There are three external parameters, ω 1 (t), ω 2 (t) and γ(t), but two coupled equations in Eq. (7). Thus we may fix one of the external parameters or some combination. We consider two simple, not exhaustive, possibilities: i) γ constant, so initial and final Ω t coincide; and ii) ω 2 constant, which implies a compression (transverse focusing to Ωt(t f ) = 6.29, see Table I. useful to avoid transversal excitation) of the final waveguide with respect to the initial one, see Table I.
The initial state chosen for the numerical examples is a product of the ground state of the transversal harmonic oscillator and a minimum-uncertaintyproduct Gaussian in the longitudinal direction centered at q l0 , with initial momentum p l0 , ψ l (q l , t = 0) = [σ √ 2π] −1/2 e i p l0 q l e −(q l −q l0 ) 2 /(4 σ 2 ) . Firstly, we design a process that interchanges ω 1 (t) and ω 2 (t) with ∆θ = π/4 and constant γ, conceived to preserve the initial longitudinal velocities in the outgoing waveguide, E l (t f ) = E l (0), and use linear ramps for the same boundary waveguides as a benchmark to compare the performance of the invariant-based protocol. Figure 1a depicts the final longitudinal energy. For the linear ramps it oscillates with operation time. The envelope for the minima is at zero but the maximum tends for long times to some value that depends on the initial wave packet. Contrast this with the full stability of the invariant-lead processes. They guarantee a fixed result, the final longitudinal energy being identical to the initial one for any initial wave packet. The transversal excitation by the linear ramps in fast processes increases considerably as the initial wave packet deviates from the origin, while the transversal excitation in the invariantbased protocol is, in general, small and much more stable. It could be further suppressed by transverse focusing and/or optimizing the u i (t). Figure 1b verifies that, for some chosen deflection angle, we can scale the final longitudinal energy at will in both scenarios (γ or ω 2 constant). Since the invariant does not affect the transversal direction, the transversal energy may be excited, but it depends on the design of the u i (t) so it can be minimized or even suppressed. Figure 2 provides snapshots of the evolution of the 2D potential for a ω 2 -constant processes that slows down the particle by a factor of two with deflection ∆θ = π/4. State transfer. Up to now we have considered real u j (t), but the coupled Newton's equations admit purely real and purely imaginary solutions combined into complex solutions. Exploiting this complex structure, u i = u R i + iu I i , leads to interesting forms of the invariant. In particular the invariant may become proportional to the uncoupled Hamiltonians at boundary times. Let us first drop the waveguide condition (8) and go back to the laboratory frame variables {q 1 , q 2 }. Defining annihilation operators in the usual manner, a i (t) = ω i (t)/2 q i + ip i / 2ω i (t), i = 1, 2, G in Eq. (6) may become a 1 or a 2 by certain choices of the u j . Let us choose at initial time (12) and u 2 (0) =u 2 (0) = 0 with c 0 real. This implies G(0) = c 0 a 1 (0), and I(0) = c 2 0 a † 1 (0)a 1 (0)/2. Instead, at final time we impose The same constant c 0 appears in Eqs. (12) and (13) because the solutions of Eq. (7) Eigenstates of H 1 (0) may thus be mapped into eigenstates of H 2 (t f ) by proper inverse engineering of the u j (t). If ω 1 (0) = ω 2 (t f ), for all initial wavepackets. (Any other scale factor may be chosen.) The system (7), which now involves four real functions, u R 1 (t), u I 1 (t), u R 2 (t), u I 2 (t), has to be solved inversely for ω 1 (t), ω 2 (t) and γ(t). The inversion is done following techniques developed for trapped ions [14] or mechanical systems [31], see appendices for a detailed account. Figure (3)a displays the resulting evolution of the control parameters for a specific example in which the frequencies ω i swap their boundary values and γ(0) = γ(t f ). Figure (3 Discussion. In some multidimensional systems with time-dependent control there are no point transformations that lead to uncoupled normal modes. Our main point here is that in these "coupled systems", invariants of motion may still guide us to inversely design the time dependence of the controls for driving specific dynamics. This inversion procedure extends the domain of invariant-based engineering, which had been applied so far to one dimensional or uncoupled systems [3]. An important difference with respect to uncoupled systems is the diminished role of commutativity of Hamiltonian and invariant at boundary times. Commutativity, because of degeneracy, does not guarantee one-to-one mapping of eigenstates of the total Hamiltonian from initial to final configurations (see appendices). One should then focus on the invariant itself for applications, and, if required, rely on design freedom to keep other variables -e.g. the total energy-controlled. An alternative to be explored is to make use of a second invariant corresponding to a linearly independent set of classical solutions of Eq. (7), {u ′ 1 (t), u ′ 2 (t)}, linearly independent with respect to {u 1 (t), u 2 (t)} [28]. Imposing boundary conditions to the second set we would aim to control the second invariant as well, but the inversion problem becomes more demanding, as the number of conditions double while the number of (common) controls remains the same.
As for further open questions, invariant-based engineering is known to be related to other STA approaches such as counterdiabatic driving for single oscillators [32]. It would be of interest to connect the current work with CD driving for coupled oscillators [33,34].  The eigenvectors of I(t b ) in the waveguide deflection example are highly degenerate, since a longitudinal plane wave multiplied by an arbitrary function of q t is a valid eigenvector with the same eigenvalue. This means that even if I(t b ) commutes with H(t b ) and shares some eigenvectors with H(t b ) the vast majority of them are not eigenvectors of H(t b ). This phenomenon -i.e., the existence of eigenvectors of one operator not shared with the other one-is well known but, since it sets an important difference with previous applications of invariant-based inverse engineering, we shall review briefly a few relevant aspects.
Let us consider a generic observable A and an orthonormal set of eigenvectors of A that forms a basis in the state space, where i = 1, 2...g n is the index to distinguish the eigenstates in the degenerate subspace for eigenvalue a n , and g n is the degree of degeneracy of a n . Now let us introduce an operator B that commutes with A. Since ψ 1 |B|ψ 2 = 0 for any two eigenstates of A, |ψ 1 and |ψ 2 , with different eigenvalues, we find a block-diagonal matrix for B in the {|φ i n } basis, with blocks of dimension g n for each eigenvalue, where any element within each block can be nonzero [35].
If g n = 1 for all n, that is, all the eigenvalues are nondegenerate, the matrix is diagonal as all the blocks reduce to numbers 1 × 1, and therefore the elements of the basis {|φ i n } are eigenvectors of B. This applies in virtually all previous works on invariant-based shortcuts in 1D [36][37][38], in which the Hamiltonian and the invariant commute at initial and final time and are not degenerate, so they share the same eigenvectors at boundary times. Thus the inversely engineered protocol drives an eigenstate of the invariant from an eigenstate of the initial Hamiltonian H(0) to an eigenstate of the final Hamiltonian H(t f ).
If g n > 1 for some n, the corresponding block does not reduce to a number and it is not, in general, diagonal. Thus, the elements of {|u i n } are not, in general, eigenvectors of B. The consequence for inverse engineering applications is that, if I(t) is an invariant and H(t) the Hamiltonian, an initial eigenvector of both I(0) and H(0) is guaranteed to be dynamically mapped at final time into an eigenvector of I(t b ), but there is no guarantee that it will be simultaneously an eigenvector of H(t b ).
In particular, in the waveguide deflection example, the longitudinal energy may be conserved "asymptotically" at the time boundaries 1 , but the process does not necessarily conserve the total energy. A factorized initial state with some longitudinal state multiplied by the transversal ground state will have at final time the same longitudinal energy that it had initially, in a different direction, but it can be transversally excited. Avoiding only longitudinal excitations is of interest per se, but we may make use of the flexibility of the shortcut design to minimize transversal excitation as well. Similarly, in the statetransfer example the energy of oscillator 1 is transferred to oscillator 2, but oscillator 1 could be excited.
In summary, commutativity of H(t b ) and I(t b ) plays a lesser role in the 2D scenario, and may in fact be abandoned for different applications. In particular H(t b ) and I(t b ) do not commute in the state-transfer example. The following section (B) explores alternative boundary conditions for the u j (t b ) that imply different meanings for the invariant at the boundary times, and therefore different possible controlled processes.

B: Other boundary conditions.
As the commutation of I(t b ) with H(t b ) does not guarantee the mapping among initial and final eigenstates of H(t b ), we shall drop this condition and explore other boundary conditions and forms of the invariant I(t b ).
For example, keeping the waveguide condition γ(t b ) = ω 1 (t b )ω 2 (t b ), note the following alternative sets of boundary conditions and corresponding quadratic invariants: where i = 1, 2 and the invariant at the boundary time t b is proportional to the transversal kinetic energy. As well, , , where the invariant at the boundary is proportional to the transversal potential energy.
The boundary conditions imposed on the u i (t) and their derivatives do not need to be of the same type at t = 0 and t f . Designing u i (t) so as to satisfy at t = 0 and t f different boundary conditions opens several control possibilities such as, for example, driving the initial longitudinal energy into final transversal kinetic energy or viceversa.

C: Uncoupled limit (waveguide with no deflection)
In the main text the time-dependent guiding from an incoming to an outgoing waveguide implies formally time-dependent coupled oscillators. If no deflection is desired, i.e. for ∆θ = 0, the easiest approach is to keep the angle θ constant and the oscillators uncoupled, γ = 0. We may thus identify q l = q 1 and q t = q 2 , and the dynamics is separable into independent dynamical normal modes. There are linear invariants for each orthogonal direction, G i (t) = u i (t)p i −u i (t)q i , and corresponding quadratic invariants I i (t) = G † i (t)G i (t)/2. Focusing on the longitudinal direction, by imposing the boundary conditions proportional to the longitudinal momentum. Consequences and applications, e.g. for cooling, are worked out in [39]. The transversal direction evolves independently. The simplest possibility is a waveguide with ω 2 constant. Of course, compressions and or expansions may as well be designed for ω 2 (t f ) = ω 2 (0) via invariants, free from any transversal excitation as is well known [4]. In the current formal framework, making use of complex solutions as in the main text we would impose u 2 (0) = ic 0 2ω 2 (0) ,u 2 (0) = −c 0 ω 2 (0) 2 ,