Connection between inverse engineering and optimal control in shortcuts to adiabaticity
aa r X i v : . [ qu a n t - ph ] J a n Connection between inverse engineering and optimal control in shortcuts toadiabaticity
Qi Zhang,
1, 2
Xi Chen,
1, 3, ∗ and David Gu´ery-Odelin † International Center of Quantum Artificial Intelligence for Science and Technology(QuArtist) and Department of Physics, Shanghai University, 200444 Shanghai, China Laboratoire Collisions, Agr´egats, R´eactivit´e, IRSAMC, Universit´e de Toulouse, CNRS, UPS, France Department of Physical Chemistry, University of the Basque Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain
We consider fast high-fidelity quantum control by using a shortcut to adiabaticity (STA) techniqueand optimal control theory (OCT). Three specific examples, including expansion of cold atoms fromthe harmonic trap, atomic transport by moving harmonic trap, and spin dynamics in the presenceof dissipation, are explicitly detailed. Using OCT as a qualitative guide, we demonstrate how STAprotocols designed from inverse engineering method, can approach with very high precision optimalsolutions built about physical constraints, by a proper choice of the interpolation function and witha very reduced number of adjustable parameters.
I. INTRODUCTION
The last ten years witnessed the huge developmentof “shortcuts to adiabaticity” (STA) with wide appli-cations ranging from atomic, molecular, and opticalphysics (AMO) to quantum information transfer or pro-cessing [1, 2]. The concept of STA was originally pro-posed to speed up the adiabatic processes in quantumcontrol. Nowadays, STA become versatile toolboxes forcontrolling the dynamics and transformation in quan-tum physics [1, 2], statistical physics [3, 4], integratedoptics [5], and classical physics [6–9]. In this context,the most popular STA techniques are the fast-forwardscaling [10, 11], the counterdiabatic driving [12, 13](or transitionless quantum algorithm [14–17]), and theinvariant-based inverse engineering [18], and their vari-ants. These three techniques can be shown to be math-ematically equivalent [19, 20]. However, the diversity ofthe designs of shortcut protocols or their combinationmay be required for a realistic experimental implementa-tion [21]. Furthermore, some counterdiabatic hamiltoni-ans turn out to be unfeasible [22], or some systems cannotbe treated by means of invariant-based engineering.STA method provides a useful toolbox for fast and ro-bust quantum controls with applications in a wide va-riety of quantum platforms such as cold atoms [23, 24],NV center spin [25, 26] including for their use as a quan-tum sensor [27], trapped ion [28], and superconductingqubit [29–32] to name a few.Such controls have also a clear added value to quantumoptimal control in quantum information processing andquantum computing [33], in terms of analytical tools, nu-merical tools, and a combination of these two. Numericaloptimal control such as the gradient ascent pulse engi-neering (GRAPE) algorithm works to some extent as ablack box. The dynamics and the structure of the con-trol field are not easily predictable [34]. STA techniques ∗ [email protected] † [email protected] based on a clear physical picture deliver a more easily un-derstandable framework but are mostly addressing prob-lems of low complexity. However, these techniques haverecently been combined with deep machine learning formore involved physical problems [35–38].Interestingly, shortcut protocols can be readily engi-neered to accommodate for various physical constraints[39], or to mitigate an environmental noise. In this re-spect, the combination of inverse engineering methodsand optimal control theory (OCT) has been particularlyfruitful [40–48]. Most STA techniques provides solutionsthat are robust against a small variation of the durationof the parameter engineering. In Ref. [49], it is shownhow OCT solutions can be adapted to accommodate forextra boundary conditions to ensure a similar robustness.Alternatively, the STA technique of Ref. [50] provides anexplicit solution for linear control problems fulfilling theKalman criterium [51].In this article, we compare systematically the inverseengineering method with the result of optimal controltheory on three specific examples that can be addressedanalytically in both formalisms: expansion of cold atomsfrom the harmonic trap, atomic transport by moving har-monic trap, and spin dynamics in the presence of dissi-pation. Our aim is to provide a pedagogical introduc-tion and comparison between a simple if not the simplestShortcut To Adiabaticity technique, the direct inverse en-gineering of the equation of motion of the dynamical vari-ables, and the optimal control theory. STA techniquesare built about the boundary conditions while OCT in-volves the minimization of a cost function. To facilitatethe comparison we therefore discuss how inverse engi-neered (IE) solutions can be modified in order to mini-mize a cost function and mimic OCT solutions. Similarlyto the variational method in quantum mechanics, and asillustrated in the following, the family of functions overwhich the minimization is performed play a crucial role.In the following, we also show how a simple ansatz hav-ing just a few tunable parameters can approach very pre-cisely the optimal solution obtained for a given physicalconstraint. II. FAST COOLING IN TIME-VARYINGHARMONIC TRAPS
Fast frictionless coooling for ultracold and Bose-Einstein condensates belongs to the first experimentaldemonstrations of STA techniques [23, 24]. Such tech-niques have been subsequently adapted and applied tocold-atom mixtures [52], Tonks-Girardeau gas [17, 53],Fermi gases [3, 54], and many-body systems [55].In this section, we address the problem of fast atomiccooling in a time-dependent harmonic trap [18]. We de-rive the time-dependence of the trap frequency by aninverse engineer procedure on an Ermakov equation andusing OCT. We subsequently compare the two types ofsolutions. Interestingly, the tunability inherent to the in-verse engineering method provides the required flexibilityto shape the inverse-engineered trajectories to minimizea cost function. We show how such solutions can be sim-ply adapted to get results very close to optimal solutionsfor a time-averaged energy cost function [18, 40, 56].
A. Model, Hamiltonian, and the InverseEngineering Approach
More specifically, we consider in the following the fastdecompression of a one-dimensional (1D) harmonic po-tential from an initial angular frequency ω (0) = ω tothe final target one ω ( t f ) = ω f , ( ω f < ω ). The prob-lem amounts to finding the time-dependent solution ofthe Schr¨odinger equation that ensures the transformationfrom the ground state of the initial trap to the groundstate of the final trap in a finite amount of time t f : i ~ ∂ψ∂t = (cid:20) − ~ m ∂ ∂x + 12 mω ( t ) x (cid:21) ψ. (1)For this purpose, we look for a scaling solution of the form ψ ( x, t ) = exp[ − β ( t )] exp[ − α ( t ) x ] f ( ρ ( t ) = x/b ( t ) , t ).The first factor accounts for the normalization, the sec-ond factor for the evolution of the phase (we show belowthat it is purely imaginary) and the last one for the de-sired scaling dynamics. By plugging such an ansatz intothe Schr¨odinger equation, we find how the different pa-rameters are related: i ~ ∂ t f = (cid:18) i ~ ˙ β + ~ m α (cid:19) f + ~ m α + i ~ ˙ bb ! ρ∂ ρ f + (cid:18) i ~ ˙ α − ~ m + 12 mω (cid:19) b ρ f − ~ mb ∂ ρρ f. (2)By introducing the renormalized time ˜ t ( t ) = R t dt ′ /b ( t ′ ) and for the choice α = ( − im/ ~ )˙ b/b and β = ln b/
2, theeffective wave function Ψ( ρ, ˜ t ) = f ( ρ, t ) obeys a time-independent Schr¨odinger equation: i ~ ∂ Ψ ∂ ˜ t = (cid:20) − ~ m ∂ ∂ρ + 12 mω ρ (cid:21) Ψ , (3) provided that the scaling parameter b ( t ) satisfies the fol-lowing Ermakov equation¨ b + ω ( t ) b = ω b . (4)Interestingly, this latter equation is amenable to a set oflinear equations. Indeed, it is the equation of an effective2D oscillator in polar coordinates, the 1 /b − is nothingbut the centrifugal barrier which acts as a repulsive forcethat prohibits the access to a zero value of b . Alter-natively, the very same result can be obtained by usingLewis-Riesenfeld dynamical invariant [18]. The groundstate wave function in such a time-dependent harmonictrap reads ψ ( x, t ) = N √ b exp im ˙ b ~ b x ! exp (cid:18) − x a b (cid:19) , (5)where N accounts for the normalization and a = p ~ / ( mω ). The self-consistent boundary conditions fora smooth continuous interpolation function are [18]: b (0) = 1 , ˙ b (0) = 0 , ¨ b (0) = 0 ,b ( t f ) = γ = q ω /ω f , ˙ b ( t f ) = 0 , and ¨ b ( t f ) = 0 . (6)As a simple example, one can choose for the scaling factor b ( t ) a fifth order polynomial ansatz that fulfills the abovesix boundary conditions [18]: b ( τ ) = 1 + ( γ − τ − τ + 6 τ ) . (7)In view of the comparison with optimal protocols, wecalculate hereafter the mean energy associated to theground state wave function (5) [56]: E ≡ t f Z t f E ( t ) dt = 1 t f Z t f h ψ ( t ) | H ( t ) | ψ ( t ) i dt = ~ ω t f Z t f (cid:18) ˙ b + ω b (cid:19) dt, (8)where E ( t ) = K ( t ) + E p ( t ) is the sum of the kinetic en-ergy K ( t ) = ~ (˙ b + ω /b ) / (4 ω ) and the potential energy E p = ~ ω ( t ) b / (4 ω ). The mean energies obey the Virialtheorem: E p = K = E/
2. For any b ( t ) trajectory thatfulfills the boundary conditions, one can infer ω ( t ) fromEq. (4) and calculate explicitly the mean energies. Usingthe available freedom to shape the scaling factor b ( t ), theinverse-engineered solutions can be tuned so to minimizethe time-averaged energy as discussed in section II C. B. Optimal control theory
Shortcut To Adiabaticity protocols such as inverse en-gineering are built about the boundary conditions. Wehave provided an example using a polynomial interpola-tion. Optimal control theory (OCT) offers an alternativeto find a path between two states but shall be built abouta cost function. We propose hereafter to use OCT on theErmakov equation (4).For this purpose, we recast Eq. (4) into a set of first or-der nonlinear coupled equations, ˙ x = f ( x ( t ) , u ) by defin-ing the x components as x = b ( t ) and x = ˙ b/ω , and in-troducing the (scalar) control function, u ( t ) = ω ( t ) /ω :˙ x = x , (9)˙ x = − ux + 1 x . (10)In the following, we work out two OCT solutions as-sociated to the minimization of the final time and thenof the mean energy. As a result of the nonlinear charac-ter of the set of Hamiltonian equations, the Pontryaginmaximum principle only gives a necessary condition toget an extremum.
1. Time-optimal solution
The so-called time-optimal solution amounts to mini-mizing the cost function J = Z t f dt. (11)with the boundary conditions (6) which translates onthe x vector components as x (0) = 1, x (0) = 0 and x ( t f ) = γ and x ( t f ) = 0. We furthermore choose theconstraint | u | ≤ p with components( p , p , p ), fulfilling Hamilton’s equations [40, 41]: ˙ x = ∂H c /∂ p and ˙ p = − ∂H c /∂ x . With the cost function J ,the control Hamiltonian H c reads H c = p + p x + p (cid:18) − x u + 1 x (cid:19) , (12)where p is a non-zero normalization constant, and p and p are generalized Lagrange multipliers. The Pon-tryagin’s maximum principle states that at any instant(0 ≤ t ≤ t f ), the values of the control function u max-imize H c . As H c is linear in the control function u andsince x >
0, the sign of the factor in front of u , ( − p x )is fully determined by the sign of − p . This latter pa-rameter plays the role of a switching function for “bang-bang” type control as discussed in the literature [40–42, 57]. The fact that the Hamilton equations are non-linear enables the possibility to have multiple bang-bang ω ω ω ω - ω ω ω f ω - - t / t f ω ( t ) / ω FIG. 1. Fast cooling in time-varying harmonic traps: The3-jump “bang-bang” control function, u ( t ) = ω ( t ) /ω . solutions [40]. We consider in the following the simplestsolution with analytical expression. This “bang-bang”solution has a single intermediate time (see Fig. 1): u ( t ) = , t ≤ − ( ω /ω ) , < t < t ( ω /ω ) , t < t < t f ( ω f /ω ) , t ≥ t f (13)With such a control function, we infer the value of thescaling factor b ( t ) from the Ermakov equation and findthe following solution for “bang-bang” control that fulfillsthe boundary conditions (6): b ( t ) = q ω + ω ω sinh ( ω t ) , ≤ t ≤ t q γ + ω − γ ω γω sin [ ω ( t f − t )] , t ≤ t ≤ t f . (14)It is worth noticing that the Ermakov equation impliesthat the quantity x + ux + x − = c is constant. Thevalue of the constant c is fixed by the initial conditionsfor 0 < t < t and by the final conditions for t < t < t f .Using the continuity of b ( t ) at t and t f due to the secondderivative in the Ermakov equation, we find the explicitexpression for both times [40, 57]: t = 1 ω arcsinh s ω ( γ − γ ω − ω ) γ ( ω + ω )( ω + ω ) , (15) t f = t + 1 ω arcsin s ω ( γ − γ ω + ω )( ω + ω )( γ ω − ω ) . (16)As the time t shall remain real, we deduce from Eq. (15)that ω ≥ ω /γ > ω f . The last inequality is naturallysatisfied because of the cooling constraint ω > ω f . Thefirst inequality requires ω /γ ≤ ω ≤ ω . In Fig. 2,we plot the normalized final time s f = t f ω as a func-tion of ω /ω and ω /ω in their accessible domains. We FIG. 2. Fast cooling in time-varying harmonic traps: 2D colorplot of the final normalized time t f ω for a 3-jump “bang-bang” control as a function of the first pulse amplitude ω /ω ,and the second pulse amplitude ω /ω . conclude that the shortest normalized final time s f isobtained for the largest ω and ω . With the choice ω = ω = ω , we obtain the shortest time s minf = π (cid:18) ω ω f (cid:19) . (17)The lowest bound for ω , namely, ω = ω /γ providesthe upper bound for final time s f = π γ, (18)where the first period of time is reduced to t = 0, sothat only two jumps are needed.In this latter range of parameter, the scaling factorreads b ( τ ) = s γ + (1 − γ ) sin (cid:20) π (1 − τ )2 (cid:21) , (19)with τ = t/t f . In Fig. 3 (a), we plot such an example ofthe time evolution of b ( t ). The solution that correspondsto the upper bound for the final time also provides theminimum time-averaged energy. Using Eq. (8), we cal-culate this latter quantity: E p = ε (cid:18) γ (cid:19) = ε π s f ! , (20)where ε = ~ ω /
4. In Fig. 3 (b), we plot this time-averaged energy E p as a function of the final time s f . ( a ) τ = t / t f b ( τ ) ( b ) s f = t f ω E p ( s f ) / ε FIG. 3. Fast cooling in time-varying harmonic traps:(a)Example of time-optimal trajectory of b ( t ) from Eq. (14).Parameters: ω f = ω /
5. (b) The time-averaged energy asa function of the normalized final time s f = πγ/
2, obtainedfrom the time-optimal control solution.
It is worth noticing that ω f and t f are not independentsince s f = πγ/ π p ω /ω f /
2. Time-averaged energy minimization
In this section, we consider optimal control solutionassociated to the minimization of time-averaged energywith unbounded constraint [56]. The lower bound for thetime-averaged potential (total) energy in Eq. (8) reads[56, 58]: E pop = ε "(cid:18) Bs f (cid:19) − − s f arctanh B + B − s f s f ! + 2 s f arctanh (cid:18) Bs f (cid:19)(cid:21) , (21)with the following solution of b ( τ ) = q ( B − s f ) τ + 2 Bτ + 1 and B = − q s f + γ .In Fig. 4, we plot this lower bound for optimizedtime-averaged energy as a blue dashed line. ●● ★★★★ ★★ (cid:11)(cid:12)(cid:13)(cid:14)(cid:15)(cid:16)(cid:17)(cid:18)(cid:19)(cid:20)(cid:21)(cid:22) s f = t f ω E p ( s f ) / ε FIG. 4. Fast cooling in time-varying harmonic traps: Com-parison of time-averaged potential energy for different opti-mal protocols: (1) energy-minimization (blue dashed line),(2) time-optimal protocol (with ω f fixed, E p constant) (redpoint); and inversed-engineered protocols: (1) 0-freedompolynomial in Eq. (7) (black dotted line), (2) polynomialIE solution with two free parameters optimized for a givennormalized final time (stars) (see Table I): s f = 1 . s f = πγ/ s f = 4 (orange solidline). Parameters: ω f = ω / C. Comparison between IE and OCT
In the previous subsections, we have reviewed thestreamline of IE and OCT protocols to ensure a fast fric-tionless decompression in a harmonic trap whose strengthcan be time-engineered. As already discussed, there isa lot of freedom to design inverse-engineered protocolssince the only requirements concern the boundary con-ditions. However, the question of the mean energy costof such protocols may be relevant since a real potentialalways exhibits some anharmonicity when the potentialenergy becomes too large. In what follows, we proposeto design IE protocols having a minimal mean potentialenergy. We will show how we can readily approach theoptimal results.The IE solution exhibited in Eq. (7) relies on a fifth-order polynomial that fulfills the six boundary condi-tions. In Fig. 4, we plot the corresponding mean potentialenergy E p ( s f ) using a black dotted line which turns outto be quite far from the optimal solution (dashed blueline).To reduce E p ( s f ), we remove the constraints on ˙ b and ¨ b at initial and final time since they are not strictly speak-ing necessary neither fulfilled by the optimal solution. Wealso enlarge the parameter space for b ( τ ) using a third-order polynomial ansatz b ( τ ) = P n =0 a n τ n to keep somefree parameters. The two boundary conditions yields a = 1 and a = − − a − a + γ . For different normal-ized final time s f , we can therefore minimize the time-averaged energy with respect to the two parameters a and a . In Table I, we provide the optimal values a and a that minimizes the mean potential energy for the ( a ) τ = t / t f b ( τ ) ( b ) τ = t / t f b ( τ ) ( c ) τ = t / t f b ( τ ) FIG. 5. Fast cooling in time-varying harmonic traps: Com-parison of time-dependent normalized variable b ( τ ) obtainedfrom optimal control theory (averaged energy optimization)(blue dashed line) and from inverse engineered solutions op-timized to minimize the time-averaged energy (red solid line)for different final times (a) s f = 1 .
1, (b) s f = πγ/
2, and(c) s f = 4. The corresponding optimal values of polynomialfunctions are detailed in Table I. Parameters: ω f = ω / three cases with s f = 1 . s f = πγ/
2, and s f = 4. s f a a πγ/ a and a in the three-order polynomial ansatz for the IE protocol thatminimize the time-averaged energy. Parameter ω f = ω / The results are represented as stars in Fig. 4. Theynearly coincide with the result of the optimal control the-ory. This is confirmed by plotting the scaling functionsfor both protocols (see Fig. 5). We conclude that the IEtrajectories inspired by the OCT solutions can be read-ily designed to approach with an impressive accuracy theexact OCT solutions.
III. FAST TRANSPORT OF ATOMS INMOVING HARMONIC TRAPS
STA techniques have also been applied to high-fidelityfast quantum transport of neutral atoms [59] or chargedions [60, 61] using a moving trap. Such developmentshave a wide range of applications from quantum infor-mation processing [62, 63] to atom fountain clock, atomchip manipulation [64–66] or atomic interferometry [67].In recent closely related works, optimal trajectories thatminimize the excitation in ion shuttling in the pres-ence of stochastic noise have been designed by combininginvariant-based inverse engineering, perturbation theory,and optimal control [68, 69].In this section, we address the problem of the fasttransport of a single atom based on a moving 1D har-monic potential. The particle is supposed to be initiallyin the ground state and shall remain in the ground stateat the final time. We follow the same kind of presentationas previously: We first design inverse-engineered proto-cols, we then derive the OCT protocols for time [43] andmean-energy optimization [45], and eventually compareboth approaches.
A. Classical and quantum inverse-engineeredsolutions
The time-dependent Hamiltonian of atomic transportusing a moving harmonic trap reads H ( t ) = ˆ p m + 12 mω [ˆ x − x ( t )] , (22)where ω is the constant trap angular frequency, and x ( t ) the time-dependent position of the trap center.This problem amounts to finding the appropriate driv-ing of this harmonic oscillator. The exact mapping be-tween the classical and quantum solutions enables oneto solve the classical problem to get a solution validquantum mechanically [50]. The time-evolution of thecoordinate, x ( t ), of a classical particle under the time-dependent Hamiltonian (22) is given by¨ x + ω ( x − x ( t )) = 0 . (23)A smooth perfect transport i.e. a transport without anyresidual oscillations at final can be obtained using inverseengineering by imposing the six boundary conditions: x (0) = x (0) = 0 , ˙ x (0) = 0 , ¨ x (0) = 0 ,x ( t f ) = x ( t f ) = d, ˙ x ( t f ) = 0 , and ¨ x ( t f ) = 0 . (24)Any interpolation function x ( t ) that fulfills these bound-ary conditions provides a possible solution of our prob-lem. For instance, one can use the following fifth order polynomial interpolation function: x ( t ) = d " (cid:18) tt f (cid:19) − (cid:18) tt f (cid:19) + 6 (cid:18) tt f (cid:19) . (25)Once x ( t ) is known, the trajectory of the trap center x ( t ) can be directly inferred from Eq. (23). A similarresult can be derived quantum mechanically using theproperties of dynamical invariants [43, 62]. In view ofthe optimization that we will perform later on, it is worthworking out the instantaneous average potential energy h V ( t ) i = ~ ω E p ( t ) , (26)where the first term accounts for the zero-point energycontribution and E p ( t ) = mω ( x ( t ) − x ( t )) i.e. theinstantaneous potential energy for the effective classicalparticle. The time-averaged potential energy is definedby E p = 1 t f Z t f E p ( t ) dt. (27) B. Optimal control theory
To recast this problem as an optimal problem, we de-fine the variables x ( t ) = x ( t ) and x ( t ) = ˙ x , and thecontrol function u ( t ) = x ( t ) − x ( t ). The control func-tion corresponds to the relative position of the effectiveparticle with respect to the trap center. The equationof motion (23) for the effective particle can be encapsu-lated in the following set of linearly coupled first orderdifferential equations ˙ x = f [ x ( t ) , u ]:˙ x = x , (28)˙ x = − ω u. (29)Interestingly, for this linear system, the solution deducedfrom the Pontryagin formalism provides the unique con-trol solution u that minimizes the cost function.
1. Time minimization
In this section, we solve the time-optimal problem withan upper bound on the relative displacement | u | ≤ δ . Thecost function to minimize t f is J = Z t f dt. (30)The corresponding Pontryagin Hamiltonian reads H c = p + p x − ω p u , where the Lagrange multipliers p and p fulfill ˙ p = 0 and ˙ p = − p . We deduce p = c and p = − c t + c where c and c are constants tobe determined. The Hamiltonian H c is a linear functionof the bounded control function u ( t ). As a result, thesign of p sets the sign of u ( t ) to maximize H c . Theparameter p being a linear function of time, the signof p can only change once. By considering the initialand final boundary conditions, the appropriate controlsequence taking into account the upper bound for | u ( t ) | is a (three-jump) “bang-bang” control u ( t ) = , t ≤ − δ, < t < t δ, t < t < t f , t ≥ t f . (31)With such a control function, the time-optimal solutionof Eq. (23) compatible with the boundary conditions (24)reads x ( t ) = , t ≤ ω δt / , < t < t d − ω δ ( t − t f ) / , t < t < t f d. t ≥ t f . (32)The driving of the trap bottom is then given by x ( t ) =¨ x ( t ) /ω + x ( t ). By imposing, the continuity on x ( t ) and˙ x ( t ), one gets the explicit expression for the switchingand final times: t = t f , t f = 2 ω r dδ . (33)According to Eq. (27), the time-averaged potential en-ergy E p for this constrained protocol is E p = 8 md ω t f = 12 mω δ . (34)
2. Mean potential energy minimization
In this section, we work out the energy-optimal proto-col. We here provide a solution that minimizes the time-averaged potential energy for a given transport time t f and distance d , with unbounded constraint. According tothe definition of potential energy, E p = mω ( x − x ) ,the cost function for this problem is J = 12 mω Z t f u dt, (35)and the Pontryagin Hamiltonian H c = 12 mω p u + p x − p ω u. (36)The Hamilton equations give two costate equations sim-ilar to those derived in the previous section. For the ( a ) s = t / t f x (cid:23) ( s ) / d ( b ) s = t / t f x ( s ) / d FIG. 6. Fast transport of atoms in a moving harmonic trap:Comparison of the trajectories of (a) the center of mass x ( t/t f ) /d and (b) the trap center x ( t/t f ) /d , obtained fromthe OCT formalism by minimizing the time-averaged poten-tial energy (blue dashed line) and using the IE approach (redsolid line) based on a fifth-order polynomial ansatz. Parame-ters: ω = 2 π ×
50 Hz and t f = 22 ms. normalization, we can choose the constant parameter p = − /m , so that the optimal problem amounts tomaximizing the quantity − u / − p u .For convenience, we consider the unbounded case ( u ( t )is unbounded) which sets the lowest bound for time-averaged potential energy E p . The quantity − u / − p u is maximal for u = − p . This expression for the controlfonction combined to Eq. (23) and the boundary condi-tions (24) enables one to determine the optimal trajec-tory of the center of mass: x ( t ) = dt t f (cid:18) − tt f (cid:19) , (37)from which we infer the trap center trajectory x ( t ) us-ing Eq. (23) with initial and final boundary conditions x (0) = 0 and x ( t f ) = d : x ( t ) = , t ≤ (cid:18) − tt f (cid:19) dω t f + (cid:18) − tt f (cid:19) t dt f , < t < t f d, t ≥ t f . (38)In Fig. 6, we plot the OCT center of mass along with thebottom trap trajectories for some specific values using (cid:24)(cid:25)(cid:26)(cid:27)(cid:28)(cid:29)(cid:30)(cid:31) !" t f ω / π E p / ε FIG. 7. Fast transport of atoms in a moving harmonictrap: Time-averaged potential energy E p /ε (normalized to ε = mω d /
2) as a function of final time t f by using differentprotocols: time-optimal (orange dash-dotted line), energy-minimization with unbounded constraint (blue dashed line),and IE approaches with a fifth-order polynomial (black up-per solid line), a seventh-order polynomial (purple solid line),and nineteenth-order polynomial (red lower solid line). Sameparameters as Fig. 6. blue dashed lines. It is worth noticing that according toour optimal solution the trap center has to include twosudden jumps at initial and final time. With such anoptimization performed for an unbounded control func-tion, we get the following lowest time-averaged potentialenergy E p ( OCT ) = 6 md ω t f . (39)In Fig. 7, we also plot this minimal time-averaged poten-tial energy as a function of the final time t f as a bluedashed line. C. Comparison between IE and OCT
In this section, we use the freedom in the interpolationfunction that enters IE solutions to approach the solutionof the optimal control theory associated to a minimiza-tion of the time-averaged potential energy.
1. IE with polynomial ansatzs
For this purpose, we first enlarge the parameter spaceof the polynomial ansatz that fulfills the boundary condi-tions (24) and search for the optimal values of the coeffi-cients that minimize the time-averaged potential energy.To satisfy the six boundary conditions (24) the mini-mal order of the polynomial interpolation function is five(see Eq. (25)). In Fig. 6, we plot the center of mass, x ( t/t f ) /d , and trap center, x ( t ), trajectories as a func-tion of time using red solid lines. The corresponding ( a ) s = t / t f x c ( s ) / d ( b ) s = t / t f x ( s ) / d FIG. 8. Fast transport of atoms in a moving harmonic trap:Comparison of trajectories of mass of center (a) and trap cen-ter (b), calculated from the OCT formalism (blue dashed line)and the IE approach (red solid line) with a 19th order poly-nomial ansatz. Same parameters as Fig. 6. time-averaged potential energy is E p ( P = 1 . E p ( OCT ) which is significantly larger than the minimal potentialenergy given by Eq. (39). It is represented as a blacksolid line in Fig. 7.In order to further reduce the time-averaged potentialenergy, we enlarge the parameter space, while keeping thesix boundary conditions satisfied. We search for a solu-tion of the forme x ( t ) = d [ P n =0 a n ( t/t f ) n ]. By applyingthe boundary conditions (24), we have a = a = a = 0, a = 21 − a − a , a = −
35 + 8 a + 3 a , and a = 15 − a − a . The time-averaged potential energycan be explicitly worked out: E p ( a , a ) = (cid:20) a − + 4385 ( a + 70) + ( a − a + 70)11 (cid:21) md ω t f . (40)The minimization of this energy yields a = 21 and a = −
70 and E p ( P ≃ . E p ( OCT ) . This curve asa function of the final time is represented in Fig. 7 as apurple solid line. It provides a clear improvement withrespect to the fifth-order polynomial solution. A pri-ori, it is possible to further improve the optimizationusing a higher order polynomial ansatz. For instance,using a 19th order well-optimized polynomial, we havefound E p ( P ≃ . E p ( OCT ) . In Fig. 8, we have plot-ted the corresponding time-dependent trajectories x c ( t )and x ( t ). However, the increase of the polynomial or-der requires a minimization with an increasing numberof parameters. This is somehow cumbersome. In thefollowing section, we propose another type of interpolat-ing function inspired by the OCT solution and yieldingastonishing results.
2. IE with hyperbolic ansatz
In this subsection, we apply the IE approach using thefollowing hyperbolic-function x ( t ) = d (cid:26) a tan (cid:20) πa (cid:18) tt f − (cid:19)(cid:21)(cid:27) + d , (41)where a > a enables one to mimic a jumpat initial and final time. This class of solution with thepossibility of an initial and final offsest and with similar ( a ) s = t / t f x c ( s ) / d ( b ) s = t / t f x ( s ) / d FIG. 9. Fast transport of atoms in a moving harmonic trap:Comparison of trajectories of mass of center (a) and trap cen-ter (b), calculated from the OCT formalism (blue dashed line)and the IE approach with the optimized hyperbolic-functionprotocol in Eq. (41) (red solid line). The “magic” values are a = 1 . a = 1 .
25, and the other parameters are the sameas those in Fig. 6. ’()*+,-./123456 t f ω / π E p / ε FIG. 10. Fast transport of atoms in a moving harmonictrap: Time-averaged potential energy E p ( t f ) /ε (normal-ized to ε = mω d /
2) calculated from different protocols:time-optimal (orange dash-dotted line), energy-minimizationwith unbounded constraint (blue dashed line), and IE ap-proach based on a hyperbolic-function-ansatz by choosing the“magic” values a = 1 . a = 1 .
25 (marked red point).Same parameters as Fig. 7. symmetry as the optimal function provides a very perfor-mant class of functions for the optimization. The freedomprovided by the two parameters a and a enables one toreduce the time-averaged potential energy while satisfy-ing the two boundary conditions x (0) = 0 and x ( t f ) = d .Such an optimization gives a = 1 . a = 1 .
25. Thecorresponding trajectories x ( t ) and x ( t ) are plotted inFig. 9, and the mean potential energy is represented inFig. 10 with marked red points. This ansatz providesa solution that nearly coincides with the exact solution, E p ( hyp ) ≃ . E p ( OCT ) .For this transport problem, we have shown how thefreedom on the interpolation ansatz enables one to opti-mize extra constraints such as the mean energy whilst ful-filling the boundary conditions. The choice of the ansatzhas a strong impact. One could naively think that a veryhigh order polynomial could always provide a succesfulstrategy. However, we have shown on this example thatthe convergence may be quite slow with the degree of thepolynomial, and that the investigation of other shapeswith a few adjustable parameters can easily outperformthe polynomial interpolation for a give constraint. IV. SPIN DYNAMICS IN THE PRESENCE OFDISSIPATION
In contrast with the previous sections, we address inthe following an example dealing with the control of in-ternal degrees of freedom. Optimal control provides apowerful tool to solve time-optimal and energy-optimalproblems in quantum two-level and three-level systems[70–73]. Such result can be directly extended to two un-coupled [72] and coupled [74] spins with similar approach.0
FIG. 11. Spin dynamics in the presence of dissipation: Equiv-alent magnetic field (
B, B c ) of transverse magnetic field B ⊥ . Using numerical optimal algorithm, robust optimal con-trol can also be designed that accounts for inhomoge-neous boarding and/or dissipation [71, 75–77]. Inverseengineering techniques have also been used for the fastand robust control of single spin [78] and two-interactingspins [78, 79] in the presence of dissipation [80]. System-atic error or perturbation induced from the parameterfluctuatiosn, dephasing noise, bit flip can be further sup-pressed using IE and OCT in atomic population transfer[46–48] and spin flip [79].Strictly speaking, the presence of dissipation rules outthe possibility of an adiabatic evolution. However, theinverse engineering can still be applied. In the fol-lowing, we consider the control of a spin 1/2 ( S =( S x , S y , S z )through the appropriate design of the time-varying magnetic field components ( B = ( B x , B y , B z ))for the desired boundary conditions. More precisely, weaddress the dissipative evolution of this spin in the pres-ence of a strong transverse relaxation rate, R >
0. Asis commonly the case in NMR, the longitudinal relax-ation rate is supposed to be negligible compared to thetransverse one, and is here neglected [73]. Under thoseassumptions, the spin components obey the Bloch equa-tions: ˙ S x = − RS x − B y S z , (42)˙ S y = − RS y + B x S z , (43)˙ S z = B y S x − B x S y . (44)Following Ref. [73], we recast the Bloch equations us-ing spherical coordinates. For this purpose, we in-troduce the angles θ ( t ) and φ ( t ) such that S =( r sin θ cos ϕ, r sin θ sin ϕ, r cos θ ) where r denotes delength of the spin r = q S x + S y + S z . It is convenient todecompose the transverse magnetic field B ⊥ = ( B x , B y ) into B ⊥ = ( B, B c ), satisfying B k S ⊥ and B c ⊥ S ⊥ (see Fig. 11): B = ( B x /R ) cos φ − ( B y /R ) sin φ and B c = ( B x /R ) sin φ + ( B y /R ) cos φ . The Bloch equa-tions can be readily rewritten with the variables a = ln r ,tan θ = q S x + S y /S z , tan φ = S y /S x , and the normal-ized time t = Rt ′ :˙ a = − [sin θ ( t )] , (45)˙ θ = B − sin θ ( t ) cos θ ( t ) , (46)˙ φ = B c cot θ ( t ) . (47)To ensure a spin rotation from an initial spin-up stateto a given final target state, we shall use the boundaryconditions θ (0) = 0 , a (0) = 0 , θ ( t f ) = θ f , and a ( t f ) = a f = − Z t f [sin θ ( t )] dt. (48)It is worth emphasizing the fact that choosing the finalspin length a f and orientation θ f for a given final time t f may have no solution for finite resources. Indeed, if thedriving by the magnetic field is not sufficiently strong,the dissipation will set an upper limit on the final spinlength.The field component B c is always perpendicular to r ⊥ and therefore only affects the spin rotation about the z -axis. The angle θ is responsible for the partial or totalspin flip. To minimize the energy cost, the trajectorylength shall be minimal. This latter condition sets thevalue of B c to zero which means φ = constant . Basically,the IE technique amounts here to fixing the θ ( t ) functionin accordance with the boundary conditions (48), andinferring the external magnetic field B ( t ) from Eq. (46). A. Energy minimization by OCT
We consider here a given spin manipulation from( a (0) = 0 , θ (0) = 0) to ( a f , θ f ) with the minimum mag-netic field amplitude. For this purpose, we aim at mini-mizing the cost function E = Z t f B ( t ) dt. (49)Let’s first recast this problem as a control problem involv-ing a set of coupled first order equations. By defining thestate variables x = a , x = θ , and the control function u ( t ) = B ( t ), the system equations (45) and (46) is of theform ˙ x = f ( x ( t ) , u ):˙ x = − sin x , (50)˙ x = u − sin x cos x , (51)1and the cost function is J = Z t f u ( t ) dt. (52)The corresponding Pontryagin Hamiltonian reads H c = − u − p sin x + p ( u − sin x cos x ) , (53)where p and p are the Lagrange multipliers fulfilling˙ p = − ∂H c /∂ x i.e. ˙ p = 0, ˙ p = p sin(2 x ) + p cos(2 x ).The maximum Pontryagin principle states for an un-bounded control u that ∂H c /∂u = 0, i.e. u = p . Inthe absence of terminal cost, the optimal solution forthis optimization between fixed initial and final statesbut without fixing the final time gives the extra condi-tion H c [ p ( t ) , x ( t ) , u ( t )] = 0: p = ( p p + cos x + cos x ) sin x . (54)From Eq. (51), we deduce˙ x = sin x p cos x + 2 p . (55)By combining Eqs. (55) and Eq. (50), we find dx = − sin x / p p + cos x dx . After integration, this rela-tion gives r ( θ ) = cos θ + p p + cos θ √ p + 1 . (56)The (constant) value of p is deduced self-consistentlywith the boundary conditions. The final time providedby OCT for an arbitrary target r f is determined by t f = Z t f dt = Z θ f θ ( θ ) dθ. (57)We note that the dissipation has an influence on the finaltime. B. Case I: reaching the horizontal plane of theBloch sphere
In this subsection, we consider the transfer of the spinfrom the quantization axis to the horizontal plane. Theboundary conditions are thus θ (0) = 0, r (0) = 1 and θ ( t f ) = θ f = π/
2. This choice sets the value of theconstant p : p π/ = 2 r f / (1 − r f ) . To address a spe-cific example, we consider in the following the final value r ( t f ) = r f = e − . The final time obtained from Eq. (57)suffers from a logarithmic divergence. To cure this prob-lem, we shift the initial and final time by a small quantity ε ≪ θ (0) = ǫ and θ f = π − ǫ : t π/ f = 1 − r f r f " ln r f r f ! − ln ǫ = 8 . ★★ ( a ) - - - :;< - - - = >? a f E ★★ ★★ - - - - @ABCDEFGH ( b ) - - - - a f E FIG. 12. Spin dynamics in the presence of dissipation: En-ergy as a function of a f = ln r f for the same target state( θ f , r f , t f ) = ( π/ , . , . a : (a) for a second order polynomial ansatz and (b) for athird order polynomial ansatz with a = 0 .
1. The inset in(b) shows the proximity of the inverse engineering result withthat of the optimal control theory. for ε = 10 − . For this specific example, the cost functionassociated to this optimal solution (see Eq. (52)) is E ( OCT ) π/ = 11 − r f = 1 . . (59)For comparison with the inverse engineering method,we propose, for the very same t f , the following secondorder polynomial ansatz: θ ( t ) = a t − a t f − θ f t f t . (60)This ansatz fulfills the boundary conditions and has asingle free parameter. The corresponding cost func-tion, E ( P , is minimal for a = − . E ( P =1 . E ( OCT ) π/ .However, our simple polynomial ansatz provides an up-per bound on the reachable values of r f . This point isillustrated in Fig. 12 (a) where we plot the energy as a2 ( a ) t B ( t ) ( b ) t θ ( t ) FIG. 13. Spin dynamics in the presence of dissipation: For a π/ B ( t ) and (b) thecorresponding variable θ ( t ) obtained from an inverse engineer-ing technique based on an optimized third-order polynomial(red solid line) and from the optimal control theory formalismassociated to a mean energy minimization (blue dashed line). function of the logarithm of the final radius a f for dif-ferent values of the free parameter a . For this example,the reachable range of values for r f is [0 . . r f = 0 . θ ( t ) = a + a t + a t + a t , (61)where t f is determined as previously ( t f = 3 . r f = 0 .
6) and, the coefficients a = 0 and a = − ( a t f + a t f − θ f ) /t f are dictated by the boundaryconditions (48). The extra freedom provided by the a coefficient enables one to (1) reach the target and (2) min-imize the cost function. With the values a = 0 . a = 0 . E ( P , is quite closeto the optimal value: E ( P = 1 . E ( OCT ) π/ . In Fig. 12(b), we plot the energy as a function of the free param-eter a for a = 0 .
1. This curve defines a new intervalof reachable r f : [0 . . θ ( t ) and itscorresponding magnetic field B ( t ) obtained from the lat-ter IE method are depicted in Fig. 13, and the associatedspin trajectory on the Bloch sphere along with the spin t S x , S y , S z FIG. 14. Spin dynamics in the presence of dissipation: Timeevolution of the spin components S z ( t ) (red solid line), S x ( t )(blue dashed line), and S y ( t ) (black dotted line) under themagnetic field obtained from the inverse engineering method.Same parameters as Fig. 13. The inset depicts the corre-sponding spin trajectory on the Bloch sphere. components in Fig. 14. Our results can be a priori fur-ther improved using an optimization on an even largerorder polynomial. C. Case II: spin flip
In this section, we consider a spin flip ( θ f = π ) forwhich the constant p parameter is p π = 2 r f / (1 − r f ) .With the same notations as previously, the final timereads (we use r f = 0 . t πf = 1 − r f r f (cid:20) ln (cid:18) (1 + r f ) r f (cid:19) − ǫ (cid:21) = 3 . . (62)The cost function associated to the optimal solution (seeEq. (52)) is E ( OCT ) π = 1 + r f − r f = 4 . . (63)This optimal solution is plotted as a blue line in Fig. 15.The optimal solution exhibits a smooth variations of θ ( t )at initial and final and a symmetry about t f /
2. Thissuggest to add the following extra condition to the poly-nomial ansatz for θ ( t ) for the inverse engineered solution: θ (0) = 0 , θ ( t f /
2) = π/ , θ ( t f ) = π, ˙ θ (0) = ˙ θ ( t f ) = 0 , and ¨ θ (0) = ¨ θ ( t f ) = 0 . (64)We have used a ninth-order polynomial to accommodatefor the 7 boundary conditions listed above, an extra pa-rameter is fixed by the final target radius, r f . The re-maining two free parameters are used to minimize theenergy. Knowing θ ( t ), we infer the magnetic field to beapplied to drive the spin in accordance with our bound-ary conditions. As explicitly shown in Fig. 15, we find a3 ( a ) t B ( t ) ( b ) t θ ( t ) FIG. 15. Spin dynamics in the presence of dissipation: (a)The magnetic field B ( t ) and (b) the variable θ ( t ) as a functionof time for a minimal-energy spin flip. An optimal ninth-order polynomial has been used for θ ( t ) to apply the inverseengineering method (red solid line). The optimal solution isplotted as a blue dashed line. bell shape for the magnetic field B ( t ) associated to this θ ( t ). However, the curves remain relatively far from theoptimal result. We find E ( P = 1 . E ( OCT ) . The rip-ples in the polynomial ansatz increase the energy andare difficult to remove by increasing the polynomial or-der. The convergence towards the optimal solution istherefore once again slow with the polynomial order.Alternatively, the shape obtained from OCT suggeststhat the following ansatz could be worth trying: θ ( t ) = π (cid:26) a tan (cid:20) πa t f ( t − t f (cid:21)(cid:27) + π . (65)Minimizing the energy, we find E = 1 . E ( OCT ) with a = 1 . a = 3 . V. CONCLUSION
In summary, we have investigated different implemen-tations of the inverse engineering method and comparethem with solutions deduced from the OCT for a givencost function. We have addressed in this manner the fast ( a ) t B ( t ) ( b ) t θ ( t ) FIG. 16. Spin dynamics in the presence of dissipation: In thecase of spin flip (b) obtained with magnetic field (a), com-pared with OCT (blue dashed line), an tanh ansatz (insteadof a polynomial) used in IE approach (red solid line) is chosento reduce energy to E = 4 . a = 1 . a = 3 . atomic cooling in harmonic trap, the atomic transportwith a moving harmonic trap, and the spin control in thepresence of dissipation. We have shown how the freedomon the ansatz inherent to inverse engineering techniquesprovide enough tunability to minimize a cost functionwhile fulfilling the boundary conditions. We have sys-tematically found class of functions with few adjustableparameters approaching the optimal control result witha relative excess of energy below one percent. Inverse en-gineered solutions are usually search as continuous andanalytical functions which is a priori an asset for theirpractical use. However, we have also exhibit the pos-sibility to design inverse engineered trajectories havinginitial and final jump to mimic the optimal control solu-tion yielding solutions that are nearly undistinguishablefrom their optimal counterpart. ACKNOWLEDGMENTS
This research was funded by NSFC (12075145),STCSM (2019SHZDZX01-ZX04, 18010500400 and18ZR1415500), Program for Eastern Scholar,Spanish Government PGC2018-095113-B-I00(MCIU/AEI/FEDER, UE), Basque Government4IT986-16. X. C. acknowledges Ram´on y Cajal program (RYC-2017-22482). [1] E. Torrontegui, S. Ib´anez, S. Mart´ınez-Garaot, M. Mod-ugno, A. del Campo, D. Gu´ery-Odelin, A. Ruschhaupt,X. Chen, and J. G. Muga, in Advances in atomic, molec-ular, and optical physics (Elsevier, 2013), pp. 117.[2] D. Gu´ery-Odelin, A. Ruschhaupt, A. Kiely, E. Tor-rontegui, S. Mart´ınez-Garaot, and J. G. Muga, Reviewsof Modern Physics 91, 045001 (2019).[3] D. Gu´ery-Odelin, J. Muga, M. J. Ruiz-Montero, and E.Trizac, Physical review letters 112, 180602 (2014).[4] I. A. Mart´ınez, A. Petrosyan, D. Gu´ery-Odelin, E. Trizac,and S. Ciliberto, Nature physics 12, 843 (2016).[5] H.-C. Chung, S. Mart´ınez-Garaot, X. Chen, J. Muga,and S.-Y. Tseng, EPL (Europhysics Letters) 127, 34001(2019).[6] S. Faure, S. Ciliberto, E. Trizac, and D. Gu´ery-Odelin,American Journal of Physics 87, 125 (2019).[7] C. Jarzynski, Physical Review A 88, 040101 (2013).[8] C. Jarzynski, S. Deffner, A. Patra, and Y. Suba¸sı, Phys-ical Review E 95, 032122 (2017).[9] Y. Li, L.-A. Wu, and Z. Wang, Physical Review A 83,043804 (2011).[10] S. Masuda and K. Nakamura, Physical Review A 78,062108 (2008).[11] S. Masuda and K. Nakamura, Proceedings of the RoyalSociety A: Mathematical, Physical and Engineering Sci-ences 466, 1135 (2010).[12] M. Demirplak and S. A. Rice, The Journal of PhysicalChemistry A 107, 9937 (2003).[13] M. Demirplak and S. A. Rice, The Journal of PhysicalChemistry B 109, 6838 (2005).[14] M. V. Berry, Journal of Physics A: Mathematical andTheoretical 42, 365303 (2009).[15] X. Chen, I. Lizuain, A. Ruschhaupt, D. Gu´ery-Odelin,and J. Muga, Physical review letters 105, 123003 (2010).[16] A. del Campo, Physical review letters 111, 100502 (2013).[17] S. Deffner, C. Jarzynski, and A. del Campo, PhysicalReview X 4, 021013 (2014).[18] X. Chen, A. Ruschhaupt, S. Schmidt, A. del Campo, D.Gu´ery-Odelin, and J. G. Muga, Physical review letters104, 063002 (2010).[19] X. Chen, E. Torrontegui, and J. G. Muga, Physical Re-view A 83, 062116 (2011).[20] E. Torrontegui, S. Mart´ınez-Garaot, A. Ruschhaupt, andJ. G. Muga, Physical Review A 86, 013601 (2012).[21] S. Ib´a˜nez, X. Chen, E. Torrontegui, J. G. Muga, and A.Ruschhaupt, Physical review letters 109, 100403 (2012).[22] Y.-X. Du, Z.-T. Liang, Y.-C. Li, X.-X. Yue, Q.-X. Lv, W.Huang, X. Chen, H. Yan, and S.-L. Zhu, Nature commu-nications 7, 1 (2016).[23] J.-F. Schaff, X.-L. Song, P. Vignolo, and G. Labeyrie,Physical Review A 82, 033430 (2010).[24] J.-F. Schaff, X.-L. Song, P. Capuzzi, P. Vignolo, and G.Labeyrie, EPL (Europhysics Letters) 93, 23001 (2011).[25] J. Zhang et al., Physical review letters 110, 240501(2013).[26] B. B. Zhou et al., Nature Physics 13, 330 (2017).[27] C. Munuera-Javaloy, Y. Ban, X. Chen, and J. Casanova,Physical Review Applied 14, 054054 (2020). [28] S. An, D. Lv, A. Del Campo, and K. Kim, Nature com-munications 7, 12999 (2016).[29] T. Wang et al., New Journal of Physics 20, 065003 (2018).[30] T. Wang et al., Physical Review Applied 11, 034030(2019).[31] A. Veps¨al¨ainen, S. Danilin, and G. S. Paraoanu, Scienceadvances 5, eaau5999 (2019).[32] T. Yan et al., Physical review letters 122, 080501 (2019).[33] S. J. Glaser et al., The European Physical Journal D 69,1 (2015).[34] E. Ass´emat, L. Attar, M.-J. Penouilh, M. Picquet, A.Tabard, Y. Zhang, S. Glaser, and D. Sugny, ChemicalPhysics 405, 71 (2012).[35] B. M. Henson, D. K. Shin, K. F. Thomas, J. A. Ross,M. R. Hush, S. S. Hodgman, and A. G. Truscott, Pro-ceedings of the National Academy of Sciences 115, 13216(2018).[36] J. J. W. Sørensen et al., Nature 532, 210 (2016).[37] D. Sels, Physical Review A 97, 040302 (2018).[38] Y. Ding, Y. Ban, J. D. Mart´ın-Guerrero, E. Solano, J.Casanova, and X. Chen, arXiv preprint arXiv:2009.04297(2020).[39] M. Larocca, E. Calzetta, and D. A. Wisniacki, PhysicalReview A 101, 023410 (2020).[40] D. Stefanatos, J. Ruths, and J.-S. Li, Physical Review A82, 063422 (2010).[41] D. Stefanatos and J.-S. Li, Physical Review A 86, 063602(2012).[42] X.-J. Lu, X. Chen, J. Alonso, and J. Muga, PhysicalReview A 89, 023627 (2014).[43] X. Chen, E. Torrontegui, D. Stefanatos, J.-S. Li, and J.Muga, Physical Review A 84, 043415 (2011).[44] Q. Zhang, X. Chen, and D. Gu´ery-Odelin, Physical Re-view A 92, 043410 (2015).[45] Q. Zhang, J. G. Muga, D. Gu´ery-Odelin, and X. Chen,Journal of Physics B: Atomic, Molecular and OpticalPhysics 49, 125503 (2016).[46] A. Ruschhaupt, X. Chen, D. Alonso, and J. Muga, NewJournal of Physics 14, 093040 (2012).[47] X.-J. Lu, X. Chen, A. Ruschhaupt, D. Alonso, S. Guerin,and J. G. Muga, Physical Review A 88, 033406 (2013).[48] D. Daems, A. Ruschhaupt, D. Sugny, and S. Guerin,Physical Review Letters 111, 050404 (2013).[49] V. Martikyan, D. Gu´ery-Odelin, and D. Sugny, PhysicalReview A 101, 013423 (2020).[50] D. Gu´ery-Odelin and J. Muga, Physical Review A 90,063425 (2014).[51] V. Martikyan, A. Devra, D. Gu´ery-Odelin, S. Glaser, andD. Sugny, Physical Review A 102, 053104 (2020).[52] S. Choi, R. Onofrio, and B. Sundaram, Physical ReviewA 84, 051601 (2011).[53] A. Del Campo, Physical Review A 84, 031606 (2011).[54] S. Deng, P. Diao, Q. Yu, A. del Campo, and H. Wu,Physical Review A 97, 013628 (2018).[55] W. Rohringer, D. Fischer, F. Steiner, I. E. Mazets, J.Schmiedmayer, and M. Trupke, Scientific reports 5, 9820(2015).5