Time-dependent harmonic potentials for momentum or position scaling

Cooling methods and particle slowers as well as accelerators are basic tools for fundamental research and applications in different fields and systems. We put forward a generic mechanism to scale the momentum of a particle, regardless of its initial position and momentum, by means of a transient harmonic potential. The design of the time-dependent frequency makes use of a linear invariant and inverse techniques drawn from ''shortcuts to adiabaticity''. The timing of the process may be decided beforehand and its influence on the system evolution and final features is analyzed. We address quantum systems but the protocols found are also valid for classical particles. Similar processes are possible as well for position scaling.

Introduction. Particle slowers and accelerators are basic tools for fundamental research and applications in different fields covering a huge range of systems, from high energy physics to atomic and molecular physics. Zeeman [1] or Stark slowers [2], optical slowers [3], magnetic inverse coil-guns [4,5], and delta-kick cooling (DKC) [6], for example, have played a central role to develop cold and ultracold physics, while accelerators are needed to launch beams for controlled collisions, deposition [7], or implantation [8] at chosen speeds. For such a vast domain of systems and conditions many different techniques have been developed. A broad family of methods applies electromagnetic fields adapted to the particle type and the operation, taking into account if the particle is charged, its magnetic moment, its dipole moment, its polarizability, or if it allows for cyclic transitions. The results often depend heavily on the initial states, initial location, velocity, or spreads, and methods that could suppress or mitigate these dependences are of general interest.
In this paper we find a simple, generic mechanism, and work out protocols, to scale the momentum of a classical particle, or of a quantum wave packet. The scaling can speed up or slow down the particle by a predetermined factor; this factor could even be negative, to produce a "momentum mirror". The main features of this mechanism are system independent, the only formal requirement is that the particle is subjected to a transient harmonic potential with time-dependent frequency during a prearranged duration. The specific system will determine the practical details on how the harmonic potential is implemented, using optical, magnetic, electrical or mechanical means. An astonishing property of the protocols described below is that the scale factor is the same for all initial conditions, i.e. for arbitrary quantum wave packets or for all initial positions and momenta of the classical particles. While, in principle, information on the exact initial condition is not needed to perform the scaling, practical considerations will of course set limits. These limits are not fundamental, but depend on the spatial, energetic, and temporal domain in which the needed harmonic potential can be effectively implemented in some specific setting.
The theory behind the time-dependent protocols for the harmonic potential makes use of an invariant of motion linear in position and momentum. Basically, we deal with an inverse problem, where the Hamiltonian is found from the desired dynamics encoded in the invariant, along the lines of the set of inverse techniques known as "shortcuts to adiabaticity" [9,10]. The theory is worked out here for a quantum particle represented by a wave packet but the resulting protocols apply equally well to classical particles since, as it is well known, harmonic potentials lead to classical equations of motion for the expectation values of position and momentum. In fact the dynamics of an arbitrary wave packet can be exactly reproduced by swarms of classical particles using the Wigner representation to fix the, possibly negative, "weighting factors" [11].
We shall present first the theory and deduce the protocols. Then we provide expressions for the time dependence of expectation values of position and momentum for a chosen scale factor, as well as expressions for second order moments for positions and momenta in terms of the initial values. This is valuable information to set both practical limits and design details depending on the intended target and resources available. We end the paper by considering related processes, in particular the scaling of positions, i.e., focusing or antifocusing.
Lewis-Riesenfeld invariants. Lewis-Riesenfeld "timedependent invariants" are operators whose expectation values remain constant for states driven by the associated time-dependent Hamiltonian [12]. The time-dependent eigenvectors of the invariant conform a convenient basis, since their probabilities remain constant along the evolution. The phases can be chosen to make each eigenvector a solution of the time-dependent Schrödinger equation. This structure has been used systematically to inverse engineer Hamiltonians from desired faster-than-adiabatic dynamics since [13], for operations to control internal or motional states. Specifically in harmonic systems, most applications have made use of quadratic invariants in positions and momenta. The existence of linear invariants was known [14][15][16] but has not been exploited for inverse engineering. The bias towards quadratic invariants in most inverse engineering applications is in part explained by the fact that an "Ermakov" quadratic invariant may be set to commute with the harmonic oscillator Hamiltonian at initial and final process times [13]. Thus fast expansions, transport, rotations, or splittings between initial and final traps can be designed so that the final energy is the same as if the process had been very slow, i.e., adiabatic [10]. Instead, the linear invariants offer the possibility to control (scale) other observables, such as the momentum, the position, and therefore kinetic or potential energies. The linear invariant eigenvectors provide continuum representations well adapted to processes where the initial and final harmonic frequencies vanish, a challenging limit for the discrete representations associated with the conventional Ermakov invariant.
Linear invariant. The Hamiltonian of a particle subjected to a harmonic potential with its center fixed at the origin and time-varying frequency is given by (1) Here we consider q and p noncommuting operators, but the same symbols may represent c-numbers in wave function representations, or as conjugate variables of a classical particle. The context should avoid any confusion. The linear combination of operators (dots stand for time derivatives hereafter) [14,15,17] satisfies the invariant equation i ∂G/∂t − [H, G] = 0, provided the reference trajectory u satisfies which is a classical equation of motion for a particle driven by a "classical" Hamiltonian (1). For any quantum state evolving with H(t), the expectation value of times m, where both functions in the argument evolve classically, i.e. following a harmonic oscillator equation (3), according to Ehrenfest's theorem. G is indeed invariant asẆ (t) = 0 using Eq. (3). This result does not depend on the particular state so the expectation values can be substituted by operators in Eq. (2). Here we shall consider only real solutions u. A corresponding quadratic invariant takes the form (To get the Ermakov quadratic invariant [10,12] u has to be made complex, see e.g. [17].) By imposing the boundary conditions at initial and final times t b = 0, t f , which also implyü(t b ) = 0, see Eq.
(3), we find G(t b ) = u(t b )p, proportional to the momentum. Thus, the final and initial momenta are proportional to each other for any wave packet, p f = (u 0 /u f ) p 0 , with a corresponding relation for kinetic energies due to the associated quadratic invariant, , and generally subscripts f and 0 for final and initial times. The scaling does not only affect expectation values but also each momentum component as we shall see. To design a harmonic slower or accelerator we first choose the scaling factor u 0 /u f , and a u(t) that satisfies the boundary conditions (5) and the scaling factor. ω(t) is found from Eq. (3) as With the chosen boundary conditions the eigenvectors of G(t b ) or I(t b ) are plane waves, i.e., not square integrable, but they form a valid and useful basis. The (constant-intime) eigenvalues of G(t) can be conveniently computed at time 0 as λ = u 0 p 0 . The initial plane wave momentum p 0 will play the role of integration variable to expand the wave functions. At an arbitrary time the eigenvectors of G(t), G(t)|φ p0 (t) = u 0 p 0 |φ p0 (t) , may be calculated as The phase ϕ p0 (t) is chosen so that Eq. (7) represents a solution of the time-dependent Schrödinger equation, and it is found by inserting Eq. (7) into the Schrödinger equation, where I t = t 0 dt ′ u 2 0 /u 2 t ′ , and in general we use the subscript t as a shorthand for the argument (t). The factor h −1/2 in the eigenvector (7) is chosen to have delta-normalized momentum plane waves at time t = 0, q|φ p0 (0) = q|p 0 , i. e., p 0 |p ′ 0 = δ(p 0 − p ′ 0 ). Instead, at final time, x|φ p0 (t f ) = e iϕp 0 (t f ) q|p 0 u 0 /u f . The invariant eigenstate that starts as a plane wave with momentum p 0 ends being proportional to a plane wave with momentum p f = p 0 u 0 /u f . An arbitrary wave function may be expanded in the basis of functions (7) as Integrating first over q in the (implicit) triple integral dq|ψ(q, t)| 2 gives a delta function in momentum so dq|ψ(q, t)| 2 = dp 0 | p 0 |ψ(0) | 2 , i.e., the norm is conserved at all times.
Here we choose polynomial trajectories for simplicity, with the coefficients fixed so thatu(t b ) =ü(t b ) = 0, where s = t/t f . See Fig. 1 for examples of this function and the corresponding ω(t) 2 . u(t) in Eq. (10) goes from u 0 to u f monotonously and posseses the symmetry The following first order moments are calculated from Eq. (9) by using triple integrals and delta-function derivatives. Since u(t) appears only in the form of the ratio U t = u t /u 0 we can work out all expressions in terms of U t , Similarly, the second order moments are pq The above first and second order moments are consistent with the invariants G and I, in the sense that the expectation values of G and I are indeed constant with them. The variances for position and momentum take the form Considering that A f = 1, we get at t f that (∆p) 2 f = (∆p) 2 0 /U 2 f for any state. Moreover I f = t f term is neglected performs the momentum scaling preserving the uncertainty product ∆p f ∆q f ≈ ∆p 0 ∆q 0 . This comes at a price, as the maximal transient value of |ω 2 | (and thus of the absolute value of the potential energy) scales as ∼ t −2 f for short times. In other words, demanding shorter and shorter process times requires the ability to implement the harmonic oscillator potential for energies growing as t −2 f . The practical limitations of the opposite, large time limit are due to the the first term in q 2 t , which grows as t 2 f . Thus, large process times need a potential implemented over a large spatial range. Similar limitations concern the first moments, in particular q t should not exceed the region where the potential may be implemented. In a realistic setting the harmonic potential will be realized within a temporal, spatial and energetic domain, which will determine the range of val-ues allowed for the initial (first or second) moments so that the final and/or transient moments do not exceed the set limits.
"Cooling", conserving phase-space volume, is an obvious application of the above by setting a large factor U f . Notice that some of the constraints of delta-kick cooling do not apply here, specifically, in DKC [6,[18][19][20] the initial state must be centered at the origin in phase space, so that a free expansion elongates the state along a given well defined angle (phase line) and a transient harmonic trap rotates the state to the horizontal (position) axis. The present method, instead, does not require any condition for the initial state, other than those imposed by the geometry of the actual setting and technical limits to implement the harmonic potential. Figure 2 (left) shows the evolution of a state in phase space, initially a minimum uncertainty product state which is initially "off center". In the simulations and figures we use dimensionless variables for coordinates, times, or momenta, defined from dimensional ones as Q = q/l, s = t/t f , P = pl/ , where l = ( t f /m) 1/2 . The Schrödinger equation becomes i∂Ψ(Q, s)/∂s = [P 2 /2 + Ω(s) 2 Q 2 /2]Ψ(Q, s) where Ω(s) = t f ω(t), Ψ(Q, s) = l −1/2 ψ(q, t) and P = −i∂/∂Q.
Position focusing, momentum mirrors, and more. We may consider as well negative scaling factors with a u(t) designed to avoid singularities in ω(t). The simplest case is U f = −1 which inverts all momenta regardless of their initial sign and the initial state. The u(t) function in Eq. (10) is valid for this purpose as the zero of u at t f /2 is canceled by a zero ofü(t f /2), see Fig. 1.
Zeros of u(t) at some intermediate time t 0 > 0 might seem to imply singularities in the wave function ψ(q, t) even if ω 2 (t) remains finite. A detailed analysis though shows that cancellations occur, e.g. due to the asymptotic property lim t→t0 u(t)I t = −u 2 0 /u t0 , so that the singularities are in fact avoided. A simple example is a Gaussian state for which the momentum integral in Eq. (9) can be done formally.
A second extension of the current methodology is "position focusing" or antifocusing, namely to scale positions rather than momenta. Formally the procedure is very similar, with a different design for u(t) so that u(t b ) = 0. Thus, the linear invariant (2) is at initial and final times proportional to q. The process scaling is of the form q f = q 0u0 /u f . In parallel with Eqs. (7,9) we work out the eigenvectors of G(t) in momentum representation with eigenvalues −mu 0 q 0 , where J t = t 0 dt ′ ω 2 t ′u 2 0 /u 2 t ′ , and a corresponding representation for arbitrary wave functions, ψ(p, t) = dq 0 φ q0 (p, t) q 0 |ψ(0) . The invariant eigenvector and solution of the Schrödinger equation φ q0 (t) evolves from an eigenvector of po-sition, φ q0 (p, 0) = p|q 0 , to a scaled version φ q0 (p, t f ) = (u 0 /u t ) 1/2 exp[−imq 2 0 J f /(2 )] p|q 0u0 /u t . For completness, the first moments are p t = (−mJ t q 0 + p 0 )(u t /u 0 ), and q t = p 0 u t /(mu 0 ) + q 0 [u 0 /u t − (u t /u 0 )J t ]. These processes may lead to position focusing or to position expansions that can be combined with side inversions if the scaling factoru 0 /u f is made negative, see Fig. 1. Again, the initial state is arbitrary. A process for focusing with side inversion is depicted in Fig. 2 for an initially off-center state.
So far we have considered, in all examples and boundary conditions, processes from free motion to free motion, i.e., ω(t b ) = 0. However the frequencies at the boundaries may have any desired value by choosing u(t b ) and its derivatives consistently. Specifically for momentum scaling, G(t b ) = u(t b )p is valid as long asu(t b ) = 0, so ω(t b ) = 0 is not necessary. Thus the approach can be adapted to scale the momenta from a trap with ω 0 = ω(0) to a trap with ω f = ω(t t ). Also the kinetic energy is scaled but not necessarily the total energy. A possible application could be to control the temperature if its final desired value does not correspond to the one of an adiabatic process. The momentum does not commute with H(t b ) for nonzero ω(t b ) so the final momenta will not be conserved for t > t f unless the trap is switched off abruptly at t f . As for position scaling, its combination with a nonzero ω(t b ) = 0 provides a way to scale the potential energy at will, since the quadratic invariant I becomes proportional to the potential energy at boundary times in these protocols.
Discussion. Spreads of momentum or velocity of initial particles often lead to particle loss and inefficiencies in focusing, slowing or acceleration processes. Shortcuts to adiabaticity techniques can be made very robust with respect to initial conditions or protocol imperfections. This feature and the possibility to choose and shorten the process time make them powerful tools to design cooling [13,21,22], even for open systems [23][24][25][26], launching [27], or compression and expansion protocols [9,10]. This work, in particular, demonstrates that, making use of linear invariants, momentum or position scaling, irrespective of initial conditions of the particle, can be achieved. The proposed methodology can be adapted to sequential interactions for beam control, or for trapped particles, for example providing a robust alternative to DKC to reach picokelvin temperatures. This work was supported by the Basque Country Government (Grant No. IT986-16), and by Grants PGC2018-101355-B-I00 (MCIU/AEI/FEDER,UE) and FIS2016-80681P.