Semiclassical fast-forward shortcuts to adiabaticity

In fast forward quantum shortcuts to adiabaticity, a designed potential $U_{FF}(q,t)$ steers a wavefunction to evolve from the $n$'th eigenstate of an initial Hamiltonian $\hat H(0)$ to the $n$'th eigenstate of a final Hamiltonian $\hat H(\tau)$, in finite time $\tau$. Previously proposed strategies for constructing $U_{FF}$ are (in the absence of special symmetries) limited to the ground state, $n=0$. We develop a method that overcomes this limitation, thereby substantially expanding the applicability of this shortcut to adiabaticity, and we illustrate its effectiveness with numerical simulations. Semiclassical analysis provides insight and establishes a close correspondence with the analogous classical fast forward method.

For continuous systems with one degree of freedom, we formulate the problem of STA as follows (see Sec. V for a slightly different formulation). Consider a Hamiltonian and let |k(t) or φ k (q, t) = q|k(t) denote its k'th eigenstate, with eigenvalue E k (t). Here φ k (q, t) is a real function and H 0 (t) is time-dependent only in the interval 0 ≤ t ≤ τ . Initializing the system in an energy eigenstate ψ t≤0 = |n(0) , we want it to evolve to the final eigenstate ψ t≥τ = |n(τ ) . For slow driving (τ → ∞) the system naturally follows the adiabatic path ψ t = |n(t) at all times [44]. But for rapid driving, additional measures -"shortcuts" -must be taken to prevent non-adiabatic transitions and guide the system to the desired final state. Here and below we suppress overall time-dependent phases when writing the state of our system. In one STA approach, known as fast-forward (FF) driving [6,7,10], an auxiliary potential U F F (q, t) is added tô H 0 (t), where U F F is designed to guide the system to the final state ψ τ = |n(τ ) . This approach offers a promising tool for achieving the controlled acceleration of continuous quantum systems, towards desired ends such as the manipulation of Bose-Einstein condensates [45][46][47][48][49] or ions in trapping potentials [50,51], or population transfer between vibrational molecular states [2,52].
Unfortunately the procedure for constructing U F F (q, t) generically leads to divergent behavior -"infinities" in the potential -at locations where the eigenstate φ n (q, t) has nodes [3,22,42]. As a result this method has been restricted to ground-state wavefunctions, which lack nodes, or else to situations in which special scale-invariant symmetries eliminate the divergences [17]. Truncating divergences with finite cutoffs has proven useful when driving a system from its ground state to its first excited state (i.e. |0(0) → |1(τ ) ) [22], but this strategy has not been explored for the STA problem formulated above (|n(0) → |n(τ ) ), and it is likely to lead to difficult-to-engineer potentials for eigenstates with many nodes.
In this paper we extend the fast-forward approach to excited states, by developing a semiclassical strategy that is free of divergences. In Sec. II we derive our main result: a semiclassically motivated recipe for constructing a divergence-free U F F (q, t) for excited energy eigenstates. In Sec. III we use numerical simulations to illustrate the effectiveness of this approach. We find that our fast-forward potential steers the wavefunction to the desired final eigenstate with high accuracy, with only a small amount of probability spilling out into sideband states. In Sec. IV we compare these results with corresponding classical simulations, and we develop a semiclassical theory for these sidebands. We end with a brief discussion in Sec. V. Several technical steps of our analysis are relegated to the Appendix, for clarity of presentation.

II. DERIVATION OF MAIN RESULT
As in Refs. [37,42], our starting point is a velocity field v(q, t) (as yet undetermined) that vanishes for t / ∈ [0, τ ]. We then use v(q, t) to define an acceleration field a(q, t) from which the fast-forward potential is constructed, as follows: Given v(q, t), Eq. 2 defines U F F (q, t) up to an integration constant u F F (t), which is chosen so that U F F = 0 for t / ∈ [0, τ ]. Thus v, a and U F F all vanish outside the interval [0, τ ].
The task now is to design v(q, t) so that U F F (q, t) (given by Eq. 2) generates the desired shortcut. To this end we define a function S(q, t) via the relation where the integration constant s(t) is adjusted (see Appendix for details) so that S obeys the Hamilton-Jacobi equation By Eq. 3, S(q, t) = S − for t ≤ 0 and S(q, t) = S + for t ≥ τ , for some constants S ± . With these definitions in place we propose the ansatz (for a given n ≥ 0) as a solution to the time-dependent Schrödinger equation (TDSE) Note thatψ t≤0 = |n(0) andψ t≥τ = |n(τ ) . Substituting Eq. 5 into Eq. 6 yields (see Appendix) wherep = −i ∂ q andv(t) = v(q, t). Thus if we can find a function v(q, t) that satisfies Eq. 7, then by construction our ansatzψ (Eq. 5) will satisfy the TDSE, and will have the desired property of beginning and ending in the n'th energy eigenstate. In other words our ansatz converts the problem of constructing a potential U F F that satisfies Eq. 6, to that of constructing a velocity field v that satisfies Eq. 7, for a given choice of n. We now explore strategies for accomplishing this task. In what follows, we drop the subscript n, and simply write φ(q, t) and E(t) for the n'th eigenstate and eigenenergy ofĤ 0 (t). Eq. 7 implies (see Appendix) which is the continuity equation for a probability density φ 2 (q, t) evolving under the deterministic flow dq/dt = v(q, t).
Thus it seems we should construct v(q, t) to satisfy Eq. 8. Unfortunately the resulting velocity field generically diverges at the nodes of the eigenstate, where φ 2 = 0 [42]. In fact, constructing U F F via Eq. 8 is equivalent to previously proposed fast-forward approaches [6,7,10], and leads to the divergent behavior mentioned earlier. Let us therefore try a different strategy, motivated by semiclassical arguments. In the WKB approximation, an excited eigenstate ofĤ 0 (t) is given by [44] φ(q, t) = ρ 2 Here we treat t as a fixed parameter, and we consider a classical trajectory of energy E = E(t) evolving under the Hamiltonian H 0 = p 2 /2m + U 0 (q, t). Then q L (t) is the left turning point of the trajectory; ±p(q, t) are the momenta FIG. 1: An energy eigenstate φ at times t and t + δt (dashed and solid curves). By Eq. 11a, the node qν (t) is displaced by v(qν , t) dt when t → t + dt.
on the upper (+) and lower (−) branches of the energy shell (i.e. level set) H 0 (q, p, t) = E(t); and ρ(q, t) is the position probability density sampled by the trajectory over one period of oscillation. In Eq. 10 and below, we replace the eigenenergy E by its WKB approximation, determined by Ep dq = 2π [n + (1/2)]. The eigenstate probability distribution φ 2 (q, t) oscillates with q due to interference between φ + and φ − (Eq. 9), but after local averaging over these oscillations, φ 2 (q, t) is semiclassically equivalent to ρ(q, t) [44]. Since the troublesome nodes in φ arise from destructive interference between φ + and φ − , let us see what happens if we attempt to construct v(q, t) to satisfy Eq. 7 separately for each term: i ∂ t φ ± =Dφ ± . Upon substitution and separation of real and imaginary parts (see Appendix) we obtain To interpret these conditions, we observe that φ ± (q, t) reflects both a classical probability distribution ρ(q, t) and a quantum phase e ±i(Σ/ −π/4) . Eq. 11b is a continuity equation for ρ(q, t) under the flow dq/dt = v(q, t), whereas Eq. 11a implies that the phase remains constant, dΣ/dt = 0, under this flow. By Eq. 9 the nodes {q ν (t)} of φ(q, t) satisfy thus the value of Σ remains constant along a node q ν (t), as the parameter t is varied. Hence Eq. 11a implies that the field v(q, t) describes the nodal "velocities" (with respect to t), as illustrated in Fig. 1. Eqs. 11a and 11b are not generally satisfied by a single v(q, t). We have thus replaced one condition, Eq. 8, that suffers from divergences, with a pair of conditions, Eq. 11, that cannot be simultaneously satisfied. To make further progress, we choose to satisfy Eq. 11a rather than Eq. 11b. Thus we take and construct U F F using Eq. 2, in the hope that this fast-forward potential will accurately steer the wavefunction to the desired final adiabatic state. A posteriori justification of this choice will be provided later. First we describe a test of this approach, using numerical simulations of a model system evolving under the time-dependent Schrödinger equation.

III. NUMERICAL RESULTS
We simulated a quantum particle of unit mass evolving in a potential where λ varies from +16 to −16 as λ(t) = 4 cos(πt/τ )[5 − cos(2πt/τ )], with τ = 1.0. The system was initialized in the state ψ 0 = |n(0) , and then evolved under the time-dependent Schrödinger equation -first using the Hamiltonian H 0 (t), then usingĤ 0 (t) +Û F F (t), with the fast-forward potential determined by Eqs. 2 and 13 as described above. FIG. 3: For evolution underĤ0 +ÛF F from the initial state ψ(q, 0) = φ17(q, 0), the final probability distribution |ψ(q, τ )| 2 (magenta) and eigenstate distribution |φ17(q, τ )| 2 (blue) are plotted. Fig. 2 presents numerical results for n = 17, setting = 2. Each frame displays a histogram showing how the evolving wavefunction ψ t is distributed among the instantaneous eigenstates |φ k (t) , as quantified by the overlap p k (t) = | φ k |ψ t | 2 . Figs. 2(a) -2(c) and 2(d) -2(f) correspond to evolution underĤ 0 andĤ 0 +Û F F , respectively. In both cases the system evolves from ψ 0 = |17 (Figs. 2(a), 2(d)) to a superposition of energy eigenstates at intermediate times (Figs. 2(b), 2(e)). While the system evolving underĤ 0 (t) ends in a broad superposition of eigenstates, Fig. 2(c), the one evolving with the fast-forward shortcut reaches the final adiabatic state to a very good approximation: in Fig. 2(f), p 17 = 0.91 and p 16 + p 17 + p 18 = 0.98. Thus the potential U F F (q, t) constructed from Eqs. 2 and 13 achieves fast-forward driving with high accuracy. We have performed simulations using different parameter values, and have found that the potentialÛ F F consistently guides the system to a final state whose probability is concentrated in the target eigenstate |n(τ ) , along with a few sideband states such as k = 16, 18 in Fig. 2 For evolution underĤ 0 +Û F F , Fig. 3 shows the distributions |ψ(q, τ )| 2 and |φ 17 (q, τ )| 2 . While the minima and maxima align nicely, the values of |ψ| 2 and |φ 17 | 2 differ visibly. This is understandable, as v(q, t) was constructed to track the eigenstate's phase Σ (Eq. 11a) rather than its magnitude ρ (Eq. 11b). The mismatch in Fig. 3 is reflected in the sidebands in Fig. 2(f). We now develop intuition by comparing our results with classical simulations, and we derive a semiclassical prediction for the sidebands (Eq. 15).  , t), the classical counterpart of the Hamiltonian used to generate Fig. 2 (lower row). The closed black curves show the adiabatic energy shell E(t) of H0(q, p, t).

IV. CLASSICAL RESULTS AND SEMICLASSICAL ANALYSIS
We first note that the potential U F F (q, t) determined by Eqs. 2 and 13 is identical to the one designed in Ref. [37] (see Eqs. 8-12 therein) for a purely classical fast-forward shortcut. This agreement reflects the Correspondence Principle (CP) and provides a measure of justification for choosing to satisfy Eq. 11a rather than Eq. 11b. Indeed, by invoking the CP we could have argued for constructing U F F (q, t) from Eqs. 2 and 13 directly from the classical results of Ref. [37], rather than by proceeding via the WKB approximation as done above. These considerations further motivate us to compare our quantum simulations with classical counterparts. Fig. 4 depicts snapshots of 50 trajectories evolving under the classical counterpart of the HamiltonianĤ 0 +Û F F used to generate Figs. 2(d) -2(f) above, with τ = 1.0. Initial conditions were spaced uniformly, with respect to the angle θ of action-angle coordinates (I, θ) [53], over an energy shell H 0 (q, p, 0) = 53.86 = E 17 (0). The closed black curves depict an energy shell E(t) determined by a fixed value of the action, I = (2π) −1 p dq. We refer to E(t) as the adiabatic energy shell.
The classical action I is an adiabatic invariant: for slow driving (τ → ∞), trajectories evolving under H 0 (q, p, t) cling to the adiabatic shell E(t) at all times, whereas for rapid driving (τ = 1.0) the trajectories evolve away from the adiabatic shell [37]. Fig. 4 illustrates the effect of adding the fast-forward potential U F F when τ = 1.0: the trajectories at first stray "off shell", Fig. 4(b), but return faithfully to the adiabatic shell at t = τ , Fig. 4(c).  Fig. 4(c) are distributed non-uniformly over the energy shell E(τ ). These features are related. Using time-dependent WKB theory [54], we obtain the following prediction for the overlap between the final wavefunction ψ(q, τ ) and the eigenstate φ n+l (q, τ ) (see Appendix): where η is the classical distribution of final conditions on the energy shell E(τ ), expressed in the angle variable θ. Eq. 15 relates the sidebands in Fig. 2(f) to the non-uniformity of points in Fig. 4(c) -if the classical distribution were uniform, η = 1/2π, the sidebands would vanish: | φ n+l |ψ | 2 t=τ = δ l0 . Fig. 5 compares the histogram appearing in Fig. 2(f) with the prediction of Eq. 15, evaluated using the trajectories in Fig 4(c). The evident agreement supports our semiclassical analysis, which suggests that the non-zero sidebands will persist even in the limit → 0. Namely, imagine that for the sameĤ 0 (t) used to generate Fig. 2, we tuned the value of (keeping other parameters unchanged) so that E = 53.86 were the 1017'th, rather than 17'th, eigenenergy ofĤ(0); and imagine we then re-ran both the quantum and classical simulations. Since the latter are unaffected by the choice of , we would again obtain Fig. 4, but the quantum simulations would produce a final histogram p k (τ ) sharply peaked at k = 1017. Eq. 15 implies that the distribution of probability among the eigenstates k = 1014 − 1020 would be the same as that observed among k = 14 − 20 in Figs. 2(f) and 5. This implies that in the limit → 0, the variance of the energy distribution associated with the non-zero sidebands tends to zero (since the energy level spacing tends to zero), consistent with a classical final state in which all trajectories end on the adiabatic energy shell (Fig. 4(c)).

V. CONCLUSIONS
We have expanded the toolkit for accelerating quantum dynamics, by developing a method for constructing fastforward shortcuts that are free of divergences and thereby applicable to excited states of generic kinetic-plus-potential Hamiltonians in one degree of freedom. Our approach is the semiclassical analogue of the classical method of Ref. [37]. It will be interesting to extend our approach to systems with d > 1 degrees of freedom. For integrable systems the existence of WKB expressions for energy eigenstates are likely to prove useful to this end, but for systems with chaotic classical dynamics the task may prove challenging. As noted elsewhere, in a number of experimentally relevant situations a separation of variables may effectively reduce a d = 3 system to a d = 1 system [37]. Finally, while we have chosen to satisfy Eq. 11a in our construction of U F F (q, t), and have used the Correspondence Principle to justify our choice, it will be interesting to test numerically whether satisfying Eq. 11b instead would lead to a fast-forward potential capable of steering the wavefunction (approximately) to the desired final energy eigenstate.
From the start, we have assumed that the HamiltonianĤ 0 (t) =p 2 /2m + U 0 (q, t) is specified for all t ∈ [0, τ ]. However, in other formulations of shortcuts to adiabaticity only initial and final HamiltoniansĤ 0 (0) andĤ 0 (τ ) are given, and the problem is then to construct a HamiltonianĤ(t) under which an initial state ψ(q, 0) = φ n (q, 0) evolves to a final state ψ(q, τ ) = φ n (q, τ ). In this situation, we can simply choose a potential U 0 (q, t) that interpolates smoothly from U 0 (q, 0) to U 0 (q, τ ) and then apply our method to obtain U F F (q, t). The HamiltonianĤ(t) =p 2 /2m + U 0 (q, t) + U F F (q, t) then generates the desired evolution.
It is instructive to compare our results with those obtained for scale-invariant Hamiltonians [17], which are characterized by time-dependence of the form Here, the parameters f (t) and γ(t) describe translations and dilations of the potential energy function. [55] For HamiltoniansĤ 0 (t) of this form, Eqs. 8, 11a and 11b are exactly satisfied by the non-divergent velocity field v(q, t) = f − (γ/γ)(q − f ), where dots indicate derivatives with respect to time. Eq. 2 then leads to the fast-forward potential The same potential was obtained in Ref. [17] using a trick involving canonical or unitary transformations of variables. This comparison suggests that the success of that trick may be related to the fact that Eqs. 8, 11a and 11b are solved by a single, well-behaved field v(q, t), which is a particular feature of scale-invariant Hamiltonians. If this is the case then the approach of Ref. [17] might not readily generalize beyond scale-invariant Hamiltonians. This issue deserves further investigation. Here we provide technical details of the derivations of Eqs. 4, 7, 8, 11 and 15.

Derivation of Eq. 4
Suppose we have constructed a function S • (q, t) that satisfies ∂ q S • = mv (Eq. 3), without as yet adjusting the constant of integration s(t). From Eqs. 2 and 3 we obtain for some β(t). By setting we arrive at S(q, t) that satisfies both Eqs. 3 and 4.

Derivation of Eq. 7
Substituting ψ(q, t) given by Eq. 5 into the time-dependent Schrödinger equation we separately evaluate each side using Eqs. 1 and 4 to obtain (dropping the subscript n) where dots and primes denote ∂ t and ∂ q , respectively, and γ(t) = −(1/ ) t 0 E(t )dt . Setting the expressions in Eq. A6 equal to one another and usingĤ 0 φ = Eφ, we arrive at Using Eq. 3 we then get which is the desired result (Eq. 7).

Derivation of Eq. 8
Dividing both sides of Eq. 7 by i and usingp = −i ∂ q , we get Multiplying both sides by 2φ we obtain which is the continuity equation for the probability density φ 2 (q, t).

Derivation of Eq. 11
Writing where A = √ ρ and α = e −iπ/4 , and substituting this expression into i ∂ t φ + =Dφ + , the two sides evaluate as follows (dropping the extraneous factor α/ √ 2): Equating these two expressions, then dividing both sides of the equation by e iΣ/ , and then collecting real and imaginary terms, we obtain Dividing the first equation by A and multiplying the second by 2A, we arrive at Eq. 11. Evaluating i ∂ t φ − =Dφ − produces identical results.

Derivation of Eq. 15
Time-dependent WKB theory [54] provides a set of tools for constructing approximations to quantum wavefunctions evolving under the Schrödinger equation, from ensembles of classical trajectories obeying Hamiltonian dynamics. Applying these tools to the final wavefunction ψ(q, τ ) after evolution underĤ 0 +Û F F , we obtain (aside from an overall phase that we ignore) where α = e −iπ/4 , η ± (q, τ ) denote the classical probability distributions on the upper and lower branches of the energy shell E = E n (τ ), depicted in Fig. 4(c), and Σ n (q, t) is given by Eq. 10. At a given location q, the two terms appearing on the right side of Eq.A14 describe a right-moving and a left-moving wave train, with local momenta ±p = ±∂ q Σ n (since the classical trajectories all end on the adiabatic energy shell) and amplitudes √ η ± .
We wish to take the inner product between ψ(q, τ ) and the k'th energy eigenstate φ k (q, τ ), for k ≈ n. Timeindependent WKB theory gives (see also Eq. 9) where ρ k± = ρ k /2 and the subscript k indicates that we are evaluating Σ(q, τ ) and ρ(q, τ ) (Eq. 10) at the energy E = E k (τ ). The functions ρ k± represent equal contributions to the probability distribution ρ k from the two branches of the energy shell. Taking the inner product between φ k and ψ gives the contribution from the upper branch, and a similar contribution φ k |ψ − from the lower branch. (We discard cross-term contributions between the two branches, as they approximately vanish upon integrating over rapidly oscillating phases.) We now use action-angle coordinates (θ, I) to analyze these contributions. In what follows, we mostly suppress the dependence of various quantities on the constant τ .
The action I corresponding to energy E is given by Since the function Σ is evaluated on a particular energy shell E, it can be viewed as a function of the action, Σ = Σ(q, I), where I = I(E). Written in this way, it is the generating function for a canonical transformation from coordinates (q, p) to (θ, I) [53]. On the upper branch of the energy shell, the angle variable θ ∈ [0, π] is given by θ(q, I) = ∂Σ ∂I . (A18) Defining I j ≡ I(E j ) = [j + (1/2)] we obtain Σ n (q) − Σ k (q) = Σ(q, I n ) − Σ(q, I k ) ≈ ∂Σ ∂I (I n − I k ) = θ (n − k) .
The distribution ρ k+ (q) is uniform in the angle variable, ρ k+ (q)dq = dθ/2π, which allows us to perform a change of variables, from q to θ, in the evaluation of Eq. A16: where η(θ, τ ) = η + (q, τ ) ∂q ∂θ , 0 < θ < π (A21) is the probability distribution on the upper branch of the energy shell, given in terms of the angle variable θ.
Animated .gif files: • movie withoutFastForwardPotential.gif : shows p k (t) = | φ k (t)|ψ(t) | 2 plotted as a function of k (red bars), for a quantum system evolving in the potential given by Eq. 14 of the main text. Parameters are the same as for Fig. 2, except = 1 and ψ 0 = |35 . The inset shows U 0 (q, t).
• movie withFastForwardPotential.gif : same as above, but with the addition of U F F (q, t). The inset shows U 0 (q, t) and U F F (q, t).