Fast Driving of a Particle in Two Dimensions without Final Excitation

Controlling the motional state of a particle in a multidimensional space is key for fundamental science and quantum technologies. Applying a recently found two-dimensional invariant combined with linear invariants, we propose protocols to drive a particle in two dimensions so that the final harmonic trap is translated and rotated with respect to the initial one. These protocols realize a one-to-one mapping between initial and final eigenstates at some predetermined time and avoid any final excitations.


Introduction
Coupled harmonic oscillators are key models in physics as they describe many different systems near equilibrium. In particular, they play a fundamental role in several quantum technologies such as simulation, communication, and information processing [1][2][3][4][5][6]. Inverse engineering the dynamics of two coupled oscillators is thus a basic and important operation for controlling quantum systems. In a previous publication [7] we demonstrated that when the coupling is proportional to the product of oscillator coordinates xy, it is possible to inverse engineer the time dependence of the quadratic potential using a combination of invariants to swap the quantum numbers of any eigenstate of the initial uncoupled oscillators. The process may be faster than the adiabatic one and it is not state-specific, in other words, the initial quantum numbers need not be known. In more detail, we identified first a "family" of driving processes based on a quadratic invariant found by Simsek and Mintert [8][9][10]. This invariant is degenerate except for the ground state, so the processes in this family do not guarantee swapping in general. To remedy the "degeneracy problem" we complemented the quadratic invariant with linear invariants and found an explicit expression for the final states, which may become the desired swapped states by adjusting the value of a single control parameter (once and for all for a predetermined process time, i.e., for arbitrary quantum numbers). The adjustment was done with very little numerical effort by running classical trajectories until a boundary condition is met.
The specific application we worked out in ref. [7] was the swapping of quantum numbers describing a single particle state in a two-dimensional harmonic trap whose final configuration is rotated by π/2 with respect to the initial configuration (The intermediate driving though is not a pure trap rotation since the eigenfrequencies are also deformed along the process). Up to a phase factor, which may be manipulated, the final state was a replica, rotated by π/2, of the initial eigenstate. On a rotating basis, the protocol in [7] facilitates a one-to-one mapping between initial and final eigenstates keeping at the final time the same quantum numbers set initially. In other words, the "swapping" results from defining the quantum numbers instead in non-rotating bases for the initial and final oscillators in x and y directions, which were the principal axes directions of the initial and final trap configurations. We shall generalize these results to allow for more general trap manipulations in 2D including arbitrary rotations and displacements. These extensions provide a theoretical background for shuttling operations in 2D networks using trapped ions in multisegmented Paul traps [11] or neutral atoms driven by optical tweezers [12], but we shall leave aside, except for some comments in the final discussion, implementation issues, which should be studied separately and vary largely with the setting and system. The first generalization considered here is to allow for rotations by an arbitrary final angle θ f (in the examples we use π/3). The particular, the π/2 angle treated in [7] is somewhat simpler, because the directions of the initial and final principal axes coincide, and ref. [7] only hinted at how to deal with other angles. Moreover, ref. [7] made use of a restricted subset of the Hamiltonians allowed in ref. [8], whereas in this paper we shall make use of the whole set to combine trap rotations and displacements of the eigenstates of the initially uncoupled oscillator map, one to one, with the final eigenstates of the final oscillators.
The work in ref. [7], allowing for a fast, controlled rotation of the final trap with respect to the initial trap, together with known shortcut-to-adiabaticity (STA) [13] protocols for 1D displacements and trap expansions or compressions, form in principle a complete toolbox of elementary motional drivings in 2D. These elementary processes could drive sequentially (i.e. by means of a succession of translations and a rotation) an arbitrary 2D trap on a plane to any other trap (with different locations, orientations, and eigenfrequencies), implementing a one-to-one mapping from initial to final eigenstates. However, sequential processes may be outperformed by combined ones in which the different operations are done simultaneously [14], for example, to turn a corner. We shall design such combined protocols here. Figure 1 depicts an example of a designed state trajectory as well as the inverse engineered trap trajectory found in this work to drive that desired state evolution. We also depict schematic snapshots (ellipses) of the evolution of the trap at equal time intervals. A clarification on notation and terminology: "trap rotation" is understood here as a rotation of the principal axes of the trap with respect to fixed directions in the laboratory frame (the rotation may occur simultaneously with deformations of the principal frequencies and translations). The rotation is characterized by an angle θ(t) that is also used to define rotating coordinates q 1 , q 2 parallel to the principal axes but having the same origin as the x, y axes, see Figure 2. We also define a polar angle for the center of the trap θ L (t) in the (x, y) lab frame, see Figure 2. These two angles are, in general, different, but in numerical examples below and in Figure 1, they are made equal at initial and final times, This coincidence at boundary times is a natural choice for possible applications but it is not a necessary condition for the applicability of the inverse method worked out in the following. . Schematic representation of a harmonic trap (red ellipse) centered at some point in the x, y plane. The principal axes of the trap are rotated with respect to the x and y directions by an angle θ. This rotation angle defines the auxiliary rotating frame q 1 , q 2 with the same origin as the (fixed, or "laboratory") x, y axes. The polar angle for the center of the trap is θ L . In general, θ = θ L although they coincide at initial and final times (θ(0) = θ L (0) = 0, θ(t f ) = θ L (t f ) = θ f ) in the process represented in Figure 1.
The paper is organized as follows: in Section 2, we describe the two-dimensional invariant adapted from Simsek et al. [8]; in Section 3, we set the model and notation; Section 4 gives the ansatzes for inverse engineering, and in Section 5 the degeneracy of the invariant is discussed; in Section 6, auxiliary linear invariants are introduced to solve the degeneracy problem. The article ends with a discussion and technical appendices.

A Quadratic Invariant for Two-Dimensional Harmonic Traps
Simsek and Mintert [8] consider Hermitian Hamiltonians with a 4 × 4 matrix form whereX T = (x, y, p x , p y ), with p x and p y being conjugate momenta for x and y, 1 represents a 2 × 2 unit matrix, and the superscript T means the transpose of the matrix. The formal treatment is dimensionless here (see a discussion of the dimensionless units in ref. [7]) so that no mass orh appear explicitly in any equation. As a reminder, the unit of (angular) frequency is one of the initial eigenfrequencies of the trap, and the unit of time its inverse. The potential has a quadratic part with real coefficients in the time-dependent 2 × 2 symmetric matrix M, and a linear part characterized by the force vector F, Here M determines the eigenfrequencies and orientation of the "trap", which may be as well a saddle, or an antitrap if one or two eigenfrequencies become imaginary (When we refer to a "trap" from now on, these possibilities-a saddle or an antitrap-should also be included), whereas F determines its displacement from the origin (The exact displacement is in Equation (19) where Γ (Γ −1 /2 is the covariance matrix of the Gaussian ground state of the invariant) is given by and the dots represent time derivatives. W and γ are given by where R is a real, Hermitian, positive semidefinite 2 × 2 matrix, and L is a real vector, satisfying with A and J determined from Here [X, Y] Z = XZY − YZX and {X, Y} = XY + YX is the anti-commutator. Equation (6) is a generalized (matrix) Ermakov equation. When inverse-engineering the Hamiltonian, the designed R determines the form of the quadratic part M. R 2 /2 is physically the spatial correlation matrix for the Gaussian ground state of the invariant. In the (slow) adiabatic limit, R → M −1/4 . Equation (7) has the form of a vector Newton equation. L gives the trajectory of the center of the ground state of the invariant [8]. Using this equation inversely, the designed (state) trajectory L together with M will determine the homogeneous two-component force F.
R and L are designed by setting their boundary conditions and interpolating between them. Analogously to one dimensional (1D) configurations [15,16], the boundary conditions are chosen so that which impliesR with t b = 0, t f being the boundary times.

Model for Rotation and Transport
The Hamiltonian (1) takes the explicit form It is useful to express M, and later the whole Hamiltonian, in terms of rotated coordinates centered at the origin of the laboratory frame, see Figure 2, where p 1,2 are the corresponding momenta and The angle θ is chosen so that q 1 and q 2 become the natural coordinates along the principal axes for the quadratic part, i.e., with (squared) eigenfrequencies given by The inverse relations are In terms of these rotated coordinates and momenta (time-dependent functions of the lab frame coordinates and momenta), the lab-frame Hamiltonian may be rewritten as where and, to simplify notation, we have dropped the explicit dependences of q 1 (x, y; t), q 2 (x, y; t), p 1 (p x , p y ; t), p 2 (p x , p y ; t), and the dependences on time of F 1 , F 2 , ω 1 , ω 2 . The trap-center trajectory, in the lab frame, is given by

Ansatzes for Inverse Engineering
We assume the following initial and final traps which are depicted in Figure 1. Assuming w 1 = ω 1 (0) = ω 1 (t f ) and w 2 = ω 2 (0) = ω 2 (t f ), the matrix M and F at initial and final time are where θ b = 0, θ f are the initial and final angles at the boundary times. At the initial time, U(0) becomes the unit matrix. L is chosen as an arch with radius r. R(t) and the polar angle θ L (t) of the point L(t) are interpolated using polynomials that satisfy the boundary conditions in Equations (10), (11) and (21), where This choice of R c is done for simplicity after some experimentation [7]. Note that R c is not positive semidefinite, so we numerically check our protocols to make sure that R(t) is. p r (t) satisfies the same boundary conditions as p(t) but it is chosen as a higher order polynomial in order to have some more freedom in the protocol design, Without trap translation this additional freedom was not necessary [7], but with trap translation, it will be needed later to solve the "degeneracy problem". In summary, with the described parameterizations the Hamiltonian is set by choosing the parameters λ and b 6 .

Degeneracy of the Invariant I 0
In principle, the Hamiltonian can be inversely designed from the invariant I 0 described in Section 2 and the boundary conditions in Equations (10) and (11). However, this invariant is degenerate except for the ground state. At the boundary times t b = 0, t f , and the invariant takes the form Using number operators n 1,2 for the oscillators along q 1 and q 2 , the invariant is at boundary times I 0 (t b ) = n 1 (t b ) + n 2 (t b ) + 1. We define for simplicity the shifted invariant I + (t) = I 0 (t) − 1, which at boundary times is a pure sum of number operators. The spectrum of I + (t) is constant and degenerate for all eigenvalues except zero (the ground state of the invariant). This poses a "degeneracy problem" since inverse engineering of the Hamiltonian based solely on I 0 (t) does not guarantee a one-to-one mapping between initial and final eigenstates except for the ground state. Just as an example, let us take the initial state |ψ(0) = |1, 0 i . The driving protocol corresponding to λ = 0, b 6 = 0, and t f = 3, gives the final linear combination |ψ(t f ) = 0.8238 e −0.0197i |0, 1 f + 0.5668 e 1.531i |1, 0 f due to the degeneracy of |0, 1 f and |1, 0 f . We solve this problem in the next section making use of additional invariants to design the driving protocol. (Our notational convention here for initial (subscript i) and final (subscript f ) eigenstates of the Hamiltonian is that the quantum numbers in |n, k i, f refer to the oscillators along q 1 and q 2 ).

Linear Invariants
To solve the degeneracy problem and get a more explicit expression for the final state, we need to introduce more invariants. The linear operators [17][18][19] G(t) = u x (t)p x −u x (t)x + u y (t)p y −u y (t)y + f (t) (29) are invariants for the Hamiltonian (12) provided In a "direct problem" (where M and F are given) the functions u x,y are classical trajectories that depend on the quadratic part M, and f may be found from them and the forces. Note that u x,y may be complex, with real and imaginary parts representing independent trajectories.
The linear invariant can be written in terms of the rotation coordinates and momenta as where Defining annihilation operators a z (t) = ω z (t)/2(q z − q z0 (t)) + ip z / 2ω z (t), z = 1, 2, with q 10 (t) = F 1 (t)/ω 2 1 (t), q 20 (t) = F 2 (t)/ω 2 2 (t), the linear invariant can be expressed in terms of creation and annihilation operators, in particular, at boundary times Different linear invariants may be constructed by choosing specific boundary conditions for u 1,2 and u 1r,2r . In particular, the initial conditions [19] define an invariant which is initially G(0) = a 1 (0). Instead, the initial conditions (the prime distinguishes them from the ones in Equation (35)) define a different invariant which at time zero is G (0) = a 2 (0). Now we may construct corresponding quadratic invariants to form initial number operators. In fact, I + (0) = G † (0)G(0) + G † (0)G (0) = a † 1 (0)a 1 (0) + a † 2 (0)a 2 (0). From Equation (28), the corresponding final invariant is I , which means the "final" linear invariants G(t f ) and G (t f ) have no independent terms, or terms that depend only on a † 1 (t f ) or a † 2 (t f ). A consequence is, see Appendix A, that the linear invariants at t f are given by combinations of the form where, using Equation (34), The coefficients satisfy the following relations, see Equation (A3), Considering the normalization condition at the final time, we also have the additional relations [7] c Finally, combining Equations (39) and (40), we have since the zeros of b, which can be found for special λ and b 6 values, see some examples in Table 1, achieve a remarkable simplification, namely, where the non-zero coefficients have unit modulus, according to Equation (40), and the phases are found from the final conditions of the trajectories, see Equation (38). These special values of λ and b 6 are found here by running classical trajectories in the inertial (xy) lab frame. The initial conditions are found from Equation (35). Fixing b 6 and sweeping on λ values, b = |c 2 | 2 is calculated for a given {b 6 , λ} pair from the final conditions using Equation (38). (Alternatively Equation (36) could also be used for the initial conditions with b computed as |c 1 | 2 from the final conditions.) In Figure 3, corresponding to Figure 1, the eigenfrequency ω 1,2 , the coordinate rotation angle θ, the polar angle for the center of the trap θ L , and the force F x,y are plotted for the parameters θ f = π/3, b 6 = 24.35, and λ = 34.84, which make b = 0. As shown in Figure 4, for a chosen b 6 there may be several λ that make b = 0 or extremely small (say 10 −9 ). In general, we pick up the smallest λ as larger values imply tighter traps [7]. b 6 may be as well optimized by looking for the smallest b values found this way, for example, see Table 1. The search could be done using different optimization subroutines but we are only interested in a proof of principle here. Note also in Figure 4 Table 1. They are chosen to minimize b. Table 1. Examples of parameters b 6 and λ that realize b = 0 for θ f = π/3, w 1 = 1, w 2 = 5. Since the ground state of the invariant is not degenerate, the initial state |0, 0 i will lead to |0, 0 f up to a phase factor. In a given basis the final state that dynamically evolves from |0, 0 i will generally be affected by a Lewis-Riesenfeld phase, i.e., it will be e Φ 00 |0, 0 f with a Φ 00 that depends on the specific transient protocol, is the instantaneous ground state of the invariant I 0 (t). For a given protocol this turns out to be a "global phase" common to all states, see Equation (45) below, so we may ignore it, as in [7], which amounts to redefining the final state basis absorbing the common phase. Such a phase would be relevant though in interferometric processes with different time-dependent protocols applied to the wavefunction branches [20,21]. An arbitrary initial eigenstate of H(0) may be written from the ground state as |ψ(0) = |n, k i = I n,k (0)|0, 0 i , where where we have used the boundary conditions (35) and (36) of the auxiliary trajectories to define G(0) and G (0). An invariant I n,k (t) is constructed with G(0) → G(t) and G (0) → G (t) using the corresponding trajectories. An invariant acting on a solution of the Schrödinger equation is also a solution. In particular I n,k (0)|0, 0 i will evolve dynamically to I n,k (t f )|0, 0 f , For the special parameters λ and b 6 such that b = 0, c 1 = e iφ 1 , and c 2 = e iφ 2 while the other coefficients vanish, c 2 = c 1 = 0. Thus, Equation (45) can be simplified as which realizes a one-to-one mapping. These special parameters λ and b 6 can be found in Table 1 and Figure 4. Finally, let us summarize the steps to design the protocol in practice: Boundary conditions at initial and final times (t b = 0, t f ) are set first for the quadratic and linear parts of the potential, M(t b ), and F(t b ). The auxiliary functions R(t) and L(t) are interpolated according to expressions in Equations (22) and (23), see their physical meaning below Equation (8). Using the form of R(t) in Equation (22), M(t) in Equation (6) can be solved with two unknown parameters λ and b 6 . For given λ and b 6 , the system Equation (30) can be solved with the initial conditions in Equations (35) and (36), and the coefficients c 1,2 , c 1,2 are found from Equation (38). The final state for some initial eigenstate of the Hamiltonian is given in Equation (45) in terms of these coefficients. To go from |ψ n,k (0) to |ψ n,k (t f ) up to a phase, avoiding the degeneracy problem, b = 0 must be satisfied, see Equation (42), so that λ and b 6 cannot be arbitrary. They can be chosen with different optimization subroutines, in particular, we choose to minimize b and λ, to avoid traps that are too tight. Once b 6 and λ are fixed so that b = 0, the time-dependent potential that defines the protocol is defined by the quadratic part M(t), and by the linear part F(t), which is found using Equation (7).

Conclusions and Discussion
In this paper, we have designed time-dependent driving protocols to map the eigenstates of an initial harmonic trap into corresponding eigenstates of a final trap, which is translated and rotated with respect to the initial one. The processing time for the one-to-one mapping can be chosen beforehand, in particular, the process may be faster than the adiabatic one, and the eigenstates whose combination form the initial state need not be known, in other words, the protocol is initial-state independent. Quadratic and linear invariants are used for inverse engineering the driving harmonic potential. The linear ones are needed to solve the "degeneracy problem" posed by the quadratic invariant proposed in [8].
Physical realizations may be performed for neutral atoms with optical traps or trapped ions with 2D Paul traps. Shorter times may imply larger and possibly imaginary potential eigenfrequencies, as well as trap trajectories that exceed a given spatial domain, see as an illustration of this point trap-center trajectories found for t f = 2, t f = 3, and t f = 5 in Figure 5. Note, in particular, that shorter times imply larger deviations of the wavepacket trajectory (green-line arch) with respect to the trap trajectory (rest of lines) and, therefore, larger potential energies. These aspects set in practice the limits for the applicability of the approach, which will depend on each specific setting and its technical constraints regarding allowed spatial domains, frequencies, energies, and time resolutions.
While our method uses invariant-based inverse engineering, this is just one among other shortcut-to-adiabaticity approaches [16] that may be worth exploring, for example, transitionless quantum driving [22], fast-forward approach [23], or time-rescaling methods [24,25]. We also note that trap rotations of charged particles have been studied with shortcut methods [26][27][28]. However, the present approach is also applicable to neutral particles. There are also many examples that show the usefulness of combining STA with optimal control theory [16]. Such a combined approach would be interesting in particular to find optimal b 6 and λ parameters. . Trajectories of the trap center for different total times (t f = 2, 3, 5) to drive the states along the arc. The values for b 6 and λ that make b = 0, assuring a one-to-one dynamical mapping between initial and final eigenstates are given in Table 1. Initial and final traps are set as in Figure 1, with the long axis in radial directions from the origin. w 1 = 1, w 2 = 5.