Shortcuts to adiabaticity for trapped ultracold gases

We study, experimentally and theoretically, the controlled transfer of harmonically trapped ultracold gases between different quantum states. In particular we experimentally demonstrate a fast decompression and displacement of both a non-interacting gas and an interacting Bose-Einstein condensate which are initially at equilibrium. The decompression parameters are engineered such that the final state is identical to that obtained after a perfectly adiabatic transformation despite the fact that the fast decompression is performed in the strongly non-adiabatic regime. During the transfer the atomic sample goes through strongly out-of-equilibrium states while the external confinement is modified until the system reaches the desired stationary state. The scheme is theoretically based on the invariants of motion and scaling equations techniques and can be generalized to decompression trajectories including an arbitrary deformation of the trap. It is also directly applicable to arbitrary initial non-equilibrium states.


Introduction
In quantum mechanics, the evolution of a system described by a time-dependent Hamiltonian H (t) is adiabatic when the transition probabilities between the instantaneous eigenstates of H are negligible. This happens either when H is time independent or when its rate of change is slow compared with the typical time scales involved [1][2][3]. Nevertheless, thinking in terms of instantaneous eigenstates is often much easier than looking for the solutions of time-dependent problems. In the field of atomic physics, going from the semi-classical approach of atom-field interaction to the celebrated dressed state picture [4] illustrates the convenience of such adiabatic representations.
For this reason, many adiabatic schemes to prepare interesting quantum states have been proposed. For instance, non-classical states [5,6] or new strongly correlated states [7] can be prepared by adiabatic passage. Quantum adiabatic computation has recently been demonstrated [8]. Yet adiabatic techniques are typically slow [3], while experimentalists are often constrained by finite lifetimes or coherence times of their samples. This has motivated the search for fast schemes reproducing or approaching adiabatic transformations. Some methods use minimization techniques to optimize the transition to a target state [9][10][11][12], whereas others yield the exact same state that would have been reached after an adiabatic transformation [13,14]. The latter are referred to as shortcuts to adiabaticity. In this paper, we detail how such 3 methods can be used on the motional degrees of freedom of ultracold gases confined in timedependent harmonic traps and experimentally demonstrate the validity of the approach. Two direct applications of the procedure are the fast cooling of atomic samples, and the suppression (or reduction) of any parasitic excitations that occur in experiments on ultracold gases when the trap geometry or the interactions are modified. Since the method is not restricted to equilibrium states it could be used in a variety of situations as discussed at the end of the paper.
The first part is theoretical and recalls how harmonically confined gases react to the variation of the trap. Both the one-dimensional (1D) non-interacting gas and the 3D Bose-Einstein condensate (BEC) with repulsive contact interaction between particles are treated. In the second part, the methods to realize shortcuts to adiabaticity are detailed for these two systems, and examples are given. The third part focuses on the experimental realization of these methods. Rapid decompressions have been performed on both a non-interacting gas and a BEC. The practical limitations that degrade the results are discussed. In the last part of this paper, we attempt to generalize the problem to an arbitrary variation of the 3D harmonic potential and give other examples of shortcuts which may be of experimental relevance.

Scaling properties of harmonically confined ultracold gases: two examples
In this section, we recall how the density and velocity distributions of a 1D non-interacting gas are affected by a change of the harmonic confinement. In 1D, the harmonic trap is fully described by its time-dependent angular frequency ω(t) and the position of its minimum q 0 (t). We show that the dynamics is fully described by two scaling functions, one associated with the cloud's size and the other with its center-of-mass position and exhibit the exact solutions of the Schrödinger equation. This will be used in the rest of the paper to realize shortcuts to adiabaticity (cf section 3). Similar scaling properties are also recalled for BECs with strong interactions in the Thomas-Fermi regime. The analogy between the invariant method used for the non-interacting gas [15] and the scaling often used for BECs [16][17][18] is underlined.

The non-interacting gas
We consider a 1D non-interacting gas confined in the most general time-dependent harmonic potential, described by the one-particle Hamiltonian where q and p are conjugate variables and m is the mass of a particle. We first recall how dynamical invariants can be used to find the general solutions of the Schrödinger equation.

Definition and properties of dynamical invariants.
In 1969 Lewis and Riesenfeld [15] generalized the concept of the invariant of motion to the case of explicitly time-dependent Hamiltonians H (q, p, t). Such Lewis invariants (also called dynamical invariants or first integrals) can be used to solve the Schrödinger equation In this case, the following properties hold [15]: (i) if |t is a solution of (2), then I |t is also a solution of (2), (ii) the eigenvalues λ(t) and associated eigenvectors |λ; t of I are a priori time dependent. We assume they form a complete set. It turns out that the eigenvalues are actually constant: λ(t) = λ. They are real because I is Hermitian.
(iv) If we assume that I does not contain the operator ∂/∂t, then the set of vectors {e iα λ (t) |λ; t , α λ (t) ∈ R(t)} is also a complete set of eigenvectors of I . If these functions are chosen to solve the equations then equation (4) also holds for λ = λ. Using the fact that the set is complete, this gives the general solutions of the time-dependent Schrödinger equation as where the c λ 's are constant complex numbers.
The solutions of the Schrödinger equation are thus given by the knowledge of an invariant I (q, p, t), any set of its time-dependent eigenvectors, and the phases α λ (t) that must solve equations (5).

Derivation of a dynamical invariant.
In this section, we give a simple derivation of the invariants of a 1D time-dependent harmonic oscillator (HO) described by (1). We use the classical formalism to derive the invariant, which is also an invariant of the corresponding quantum system. The canonical equations of motion associated with the Hamiltonian (1) are where are the Poisson brackets of two observables A and B.
When the angular frequency ω(t) and trap center q 0 (t) vary, one expects the cloud to be displaced and to change its size; thus one can introduce a canonical change of variables leading to a new Hamiltonian H . One has to derive conditions on the real dimensionless function b and the displacement function q cm such that the transformation is canonical. For this, we look for a new Hamiltonian of the form where ω 0 is a constant angular frequency. The Hamiltonian explicitly depends on time only through the function f (τ ) (which does not contain the variables Q and P). The transformation (8) From equation (10a), one deduces that and that where˙denotes the derivation with respect to time t. From equation (10b), one finds that the functions b and q cm must obey the two differential equations When these two equations are satisfied, the quantity corresponds to a gauge transformation which changes the phase of the wavefunction in the following manner: where F is a primitive of f .  (1). This is detailed in appendix A. The result is that the wavefunction associated with the eigenvalue λ n of the invariant I is expressed in terms of the nth Hermite polynomial H n as Here where and q 0 , b, q cm and their derivatives are functions of t (t when they are under an integral symbol) and are linked by equations (13) and (14). a ho = √h /mω 0 is the HO length of I .
From this expression, we see the physical meaning of the two scaling functions: q cm (t) is the center of the wavefunction (center of mass of a cloud that was initially at equilibrium), and a ho b(t) is the width of the wavefunction.

The case of an interacting Bose-Einstein condensate (BEC)
For the corresponding 3D interacting system of N particles, the Hamiltonian is The potential U is supposed to be a time-dependent 3D HO and the rotation of this harmonic confinement is excluded for the moment (the trap keeps the same eigenaxes): V is the interaction potential between two particles, which is well approximated by a delta function for ultracold gases [19]. The procedure described in section 2.1 cannot be easily adapted, because it would require the knowledge of an invariant of this many-body system. But when dealing with a BEC, the dynamics is well described by a single-particle wavefunction, whose evolution obeys a nonlinear Schrödinger equation, the Gross-Pitaevskii equation (GPE) [19].

Scaling approach.
Let us consider a quantum system described by the wavefunction ψ(r, t), whose time evolution is given by the GPE with m the mass of particles, N the number of particles andṼ = 4πh 2 a s /m the interaction coupling constant generated by s-wave scattering between particles, characterized by the scattering length a s . Analogously to the non-interacting case, we wish to write the solution of the time-dependent GPE as a function of the solution of a time-independent one expressed in a suitable frame of reference. Following this line, a strategy to solve equation (24) is to find a change of variables ρ(r, {b i (t)}, {r cm i (t)}) where the b i 's and the r cm i 's are scaling and translation functions such that equation (24) can be written as a time-independent equation (i.e. a GPE with a time-independent potential) on the wavefunction χ(ρ, τ ), defined by the relation A(t) being a time-dependent normalization factor and φ(r, t) a space-and time-dependent phase. All the dynamics induced by the time-dependent potential is transferred to the functions {b i (t)} and {r cm i (t)}, and the differential equations they have to satisfy (similar to equations (13) and (14)). If one can solve the new time-independent equation on χ, one solves equation (24) and knows the wavefunction ψ(r, t).
Equation (24) is invariant under the transformation in any of the following cases: (i) in the non-interacting limit [16,20]: in this case the system is equivalent to three independent HO of the kind treated in section 2.1; (ii) for a suitable driving of the interaction termṼ [20], that is, assuming one can controlṼ (t) at will (for cold gases, this can be done using Feshbach resonances); (iii) in the TF limit [17].
This third case, which is detailed in the following section, is used in the rest of the paper.

Condensate wavefunction in the Thomas-Fermi approximation.
Given a timedependent GPE, the TF approximation consists in neglecting the kinetic-energy-like term in the ρ-frame of reference, i.e. −h 2 /(2m) i b −2 i ∂ 2 χ/∂ρ 2 i , supposed to be small compared to the interaction term [17,18]. In this regime, provided that A(t) = ( i b i ) −1/2 and that where the scaling and translation functions are solutions of 8 one obtains the following equation on χ: where we defined a rescaled time If at t = 0 the condensate is at equilibrium, the solution of equation (31) is µ being the chemical potential. This gives the typical inverted parabola density profile whose being the initial TF radii.

Shortcuts to adiabaticity
In this section, the definition of a shortcut to adiabaticity is given, and the results of section 2 are used to derive angular frequency trajectories realizing such shortcuts, for both non-interacting gases and interacting BECs confined in time-dependent harmonic traps.

Shortcut to adiabaticity based on an invariant of motion
For a system described by a Hamiltonian H (t), a shortcut to adiabaticity is realized when another Hamiltonian H (t) can be found, such that the state obtained after a finite time of evolution with H (t) is identical (up to a global phase factor) to the final state of the adiabatic evolution with H (t). Shortcuts to adiabaticity are not adiabatic; only the final state is identical to that obtained after an adiabatic evolution.
The possibility of such shortcuts has been known for a long time. For instance, in the case of a HO with a time-dependent frequency ω(t) treated in [15], when discussing the transition probability P sm between two instantaneous eigenstates |s; t and |m; t , the authors noticed that some trajectories ω(t) could lead to the same result as the adiabatic case, namely P sm = δ sm . Such shortcuts to adiabaticity can thus be realized simply by engineering the time-dependent parameters of the Hamiltonian.
A practical method to find a class of appropriate ω(t) was detailed by Chen et al [14]. In this case, the Hamiltonian is chosen to be time independent (but with different frequencies) outside the time interval t ∈ [0, t f ]. An invariant is engineered to commute with the Hamiltonian outside of this interval. This yields a specific ω(t) for which all the eigenstates of H (t 0) are exactly mapped to the corresponding ones of H (t t f ) after the evolution for t ∈ [0, t f ]. Up to a global phase and a rescaling of energies and lengths, the final state (at time t = t f ) is identical to the initial one (t = 0), i.e. if the initial state was where {|n; t , n ∈ N} is a basis of instantaneous eigenstates of H (t) and {hω n (t)} the corresponding eigenvalues, and n |c n | 2 = 1, then the final state is where This is true even if the initial state is not an equilibrium state.

Frequency trajectory for a non-interacting gas. The Hamiltonian is assumed to have the form
which is identical to (1), with the additional constraint q 0 (t) = −g/ω 2 (t) (and a gauge transformation consisting in adding −m 2 (t) q 2 0 (t)/2 to H ). It describes a single particle in a harmonic trap subject to a constant force, which, in the experiments presented in section 4, comes from gravity. The angular frequency ω(t) is assumed to be constant outside the interval t ∈ [0, t f ]. During this interval, the problem is to find the appropriate frequency trajectory ω(t), connecting the initial trap of initial frequency ω(0) to a final trap of frequency ω(t f ), for the decompression (or the compression if ω(0) < ω(t f )) to implement a shortcut to adiabaticity. Figure 1 shows the initial and final situations assuming a decompression [ω(t f ) < ω(0)].
We use the strategy introduced by Chen et al [14]. If the invariant commutes with the Hamiltonian for t 0 and t t f , and provided that the functions b and q cm are sufficiently continuous, the stationary states of H (t 0) will be transferred to the corresponding ones of H (t t f ). It is convenient to use the dimensionless function instead of q cm , and to rewrite equation (14) using the rescaled time τ (equation (11)). Equation (14) becomes If one chooses to set ω 0 = ω(0), and the conditions then where h is a function of time only. These boundary conditions (BCs) thus fulfill the condition (37). Since the functions b and c must be solutions of equations (13) and (39), four additional BCs must be satisfied: In order to construct the functions b and c satisfying these BCs and the two differential equations (13) and (39), it is convenient to write all the BCs on the function c and its derivatives with respect to the rescaled time τ . Using equations (11) and (13) and differentiating equation (39) twice with respect to τ , one finds the ten conditions and, for all k ∈ {1, 2, 3, 4}, which are sufficient for the 12 BCs (40). τ f is the rescaled time corresponding to t f : Under this form, the BCs are well suited to use a polynomial ansatz for c(τ ), deduce b(τ ) with equation (39), compute τ (t) by numerically integrating equation (11) and obtain b(t). The final step consists in using equation (13) to obtain the time-dependent trap frequency as ω 2 (t) = ω 2 0 /b 4 −b/b. An example of this procedure is given in figure 2 for particular values of the initial and final frequencies. The final rescaled time τ f can be chosen at will, it can be arbitrarily small, but one important constraint on the function c is that it must not lead to vanishing values of b which give infinite ω 2 (t). Additional constraints on c arise from experimental requirements, such as positive ω 2 (t) (attractive potentials), maximal and minimal frequencies attainable with a given setup, speed at which ω(t) can be varied, etc. Since all this depends on a particular experimental setup, no mathematical analysis of the best ansatz to use was done.
For the experiments presented in section 4 and in [21,22], a polynomial of order 15 was used: The first coefficient is fixed at 1 by equation (41a) and c 1 , . . . , c 4 are fixed at 0 by equations (41c). We arbitrarily impose c 5 = c 6 = · · · = c 10 = 0, which leaves five coefficients which are uniquely determined by the remaining BCs (41b) and (41d). The calculation of these remaining coefficients is done by inverting the linear system corresponding to these five equations. In principle, since there are ten BCs, a 9th order polynomial can be used, which yields a unique solution for the ten coefficients of c. Nevertheless, the obtained trajectory was not well behaved enough to be realized experimentally (the frequency was decreasing too fast in the beginning compared to what could be achieved by the apparatus). This is the reason why a higher-order polynomial was used and six coefficients were fixed at 0.
Since the polynomial can be of any order greater than 9 and the BCs only impose a linear relation between ten of its coefficients, there is obviously an infinite number of different solutions connecting the two given initial and final states. Moreover, other functions than polynomials could be used for c, as long as they provide enough free parameters.
The obtained non-zero coefficients of (42) are given in table 1.

Example.
In this section, we determine the trajectory used in section 4.2 and in [21]. The parameters are given in table 2. Figure 2 shows the functions c(τ ), b(τ ), t (τ ) and ω(t)/2π corresponding to this decompression. Since the exact wavefunctions are known, all the properties of the atomic cloud can be calculated during decompression. For instance, figure 3 displays the size and center-of-mass Table 1. Non-zero coefficients of the polynomial ansatz for c(τ ) calculated from the BCs (41).
position of a cloud initially at equilibrium in the compressed trap. These are compared to the same values if the decompression was done very slowly as in the adiabatic theorem. The clear difference between the plain and dashed curves illustrates the fact that the decompression is not adiabatic.

Shortcut to adiabaticity for an interacting BEC in the Thomas-Fermi limit
Let us suppose that ψ(r, t 0) is a stationary state of equation (24). We can engineer the parameters of the potential U (r, t) such that ψ(r, t f ) is also a stationary state for t t f . This implies that χ (ρ, τ τ f ), with τ f = τ (t f ), must be a stationary state of equation (31) and that ∇ r φ(r, t f ) = 0. If these two conditions hold, ψ(r, t) can evolve during the time interval [0, t f ] between two stationary states even being strongly different from the adiabatic stationary state during the evolution for 0 < t < t f . In our experiment, the time-dependent trapping potential has a cylindrical symmetry of the form with initial and final angular frequencies ω ⊥, (0) and ω ⊥, (t f ) = ω ⊥, (0)/γ 2 ⊥, , respectively. This case corresponds to fix ∀t, r 0 x,y (t) = 0 in equation (23) and r 0 z (t) = −g/ω 2 ⊥ (t). By introducing the dimensionless function the differential equations (29) and (30) take the form The final state is an equilibrium state if the final TF radii verify that R ⊥, (t f )/R ⊥, (0) = γ 2 ⊥, , if the vertical center-of-mass position fulfills the condition r cm . These latter imply thatb ⊥, (0) =b ⊥, (t f ) = 0 must hold as well, giving 16 independent BCs.
Our procedure to engineer ω ⊥, (t) is to reduce the dimensionality of the problem by only considering the trajectories that lead to a constant axial size. This corresponds to keeping b (t) = b (0) for any t, fixing a trap decompression with γ ⊥ = γ 2 . In this case, equations (45)-(47) reduce tob Equation (48) is identical to equation (13) and equation (50) is nothing but equation (39) expressed with the real time (the rescaled time being given by equation (32) instead of equation (11)). Thus we can exploit for b ⊥ (t) and c(t) the solutions obtained for the non-interacting gas, provided that the axial frequency is varied according to equation (49).

Example.
As an example of the procedure described above, we determine the trajectories used in section 4.3 and in [22]. The decompression parameters are given in table 3. The radial frequency is reduced by a factor of 9 and the axial frequency by a factor of 3. The obtained trajectories are represented in figure 4.

Validity of the Thomas-Fermi approximation.
To check the validity of the Thomas-Fermi approximation that led to the trajectories of figure 4, 3D Gross-Pitaevskii simulations have been performed and compared to the analytical results of section 3.2. In the numerical solution, we use a split step operator in time combined with a fast Fourier transformation in space. The results are presented in figure 5 and show that this approximation is well justified for our experimental parameters (decompression of figure 4, the number of atoms N ∼ 10 5 , scattering length of 87 Rb of a s ∼ 100 a 0 , a 0 being the Bohr radius).

Experimental realization of shortcuts to adiabaticity
The procedure described above was tested experimentally by quickly decompressing a trapped ultracold gas of 87 Rb atoms. In this section, we describe the experimental steps involved in the preparation of the cold sample (cold thermal gas or BEC) and then explain how the decompression is controlled, monitored and compared to simpler (non-optimal) schemes.

The apparatus
4.1.1. Production of ultracold clouds. Our Bose-Einstein condensation apparatus is based on a double magneto-optical trap (MOT) design. The upper chamber, operated at a relatively high pressure (∼10 −9 mbar), contains a large MOT with 10 11 87 Rb atoms, which constitutes the primary source of cold atoms. These are then transferred with a pushing laser beam to a second MOT located in the lower, ultra-high vacuum chamber (∼10 −11 mbar). After the standard steps of polarization gradient cooling and optical pumping, atoms in |5 2 S 1/2 , F = 2, m F = 2 are transferred to a quadrupole magnetic trap. After a phase of adiabatic compression, this trap is converted into a Ioffe-Pritchard trap by ramping up the current in a third coil, following the QUIC trap design of [23]. Once the cloud is in the Ioffe-Pritchard trap, radio-frequency (rf) evaporative cooling is performed to reach BEC within 10 s. We are able to produce almost pure BECs (no discernible thermal fraction) containing up to 5 × 10 5 atoms. The fast decompression of such a BEC with sizeable interactions is studied in section 4.3. In contrast, to produce a thermal cloud with negligible interactions such as that used in section 4.2, we reduce the loading time of the secondary MOT to start the evaporation with a lower initial collision rate.

Control of the trapping frequencies.
Implementing shortcuts to adiabaticity requires a precise control of the trapping frequencies, in a dynamical fashion. In our QUIC magnetic trap, this can be achieved by varying the current i Q running through the three coils, and the current i B 0 running through an additional pair of Helmoltz coils disposed along the long (axial) dimension of the trap (compensation coils). The resulting potential is where µ/ h = 1.4 MHz G −1 for atoms in |5 2 S 1/2 , F = 2, m F = 2 . B is the radial magnetic field gradient, while B corresponds to its curvature along y. The harmonic approximation of equation (51) describes accurately the potential seen by cold enough atoms, i.e. k B T µB 0 [24]. Then, the radial and axial angular frequencies are These expressions show that the radial and axial frequencies can be controlled independently to some extent. The dynamical control of the currents is achieved using homemade, computercontrolled electronic circuits. The experimental realization of shortcut trajectories requires a careful preliminary calibration of the frequencies versus currents, which is achieved by monitoring the cloud's center-of-mass oscillations after a small excitation. Due to the finite-time response of the controlling circuit, it is also necessary to check the behavior of the frequency during an actual trajectory. This is illustrated in figure 6, where we compare the theoretical decompression trajectory of figure 2 (line) and measured experimental values (circles). In this example, the deviation is below 5%.

Shortcut to adiabaticity for a non-interacting gas
For this first experiment, we use the procedure described in section 4.1 to produce a thermal gas (N 10 5 , T 0 = 1. ω y (0)/2π = 22.2 Hz and ω z (0)/2π = 235.8 Hz. The initial cloud-averaged collision rate per particle is γ el 8 Hz, which corresponds to a collision time of ∼125 ms. This is 30 times the oscillation period and more than 3 times the decompression time, which justifies the non-interacting approximation.
We use here the decompression trajectory discussed in section 3.1.2, adapted to the vertical axis (Oz), with the parameters of table 2. To maximize the decompression factor γ 2 = ω z (0)/ω z (t f ), the compensation coil current i B 0 is increased from i B 0 (t = 0) 0 A to i B 0 (t f ) = 3.0 A, while the QUIC current is decreased from i Q (t = 0) = 26.7 A to i Q (t f ) = 3.6 A (see the resulting trajectory in figure 6). The decompression duration is t f = 35 ms.
In theory, starting from a gas at equilibrium and temperature T 0 in the compressed trap, a shortcut to adiabaticity should lead to an equilibrium state in the final trap, with a temperature T f = T 0 ω(t f )/ω(0). This corresponds to a situation where entropy has not increased. In contrast, for a non-optimal decompression, one expects to observe oscillations of the cloud's size and center of mass in the decompressed trap, once the decompression is completed. To evaluate the efficiency of our shortcut, we thus perform the fast decompression and hold the cloud in the decompressed trap for a variable amount of time. The trap is then abruptly switched off, and an absorption image is taken after a constant time of free expansion (6 ms). The amplitude of the dipole (oscillation of the center of mass) and breathing modes (oscillation of the size) give access to the excess energy provided to the cloud, as compared to an adiabatic modification of the potential. If the cloud is reasonably at equilibrium after decompression, one can also directly measure the final temperature by measuring the evolution of the size during a free expansion.
In the following, we compare four decompression trajectories: 1. the shortcut, given in figures 2(d) and 6; 2. a linear decompression of the same duration (35 ms); 3. an abrupt decompression, which, somehow, corresponds to a worst case scenario (in practice, the decompression time is 0.1 ms and ω(t) is not controlled; it is imposed by the response of the magnetic trap control electronics); 4. a 6 s long linear decompression, which can be considered nearly adiabatic. Compared to the linear decompression in 35 ms, the shortcut reduces the amplitude of the dipole mode by a factor of 7.2 (obtained from the sine fits) and the amplitude of the breathing mode by a factor of 3 (comparison of the standard deviations of the two sets of data). The excess energy, which is dominated by the center-of-mass energy, is thus reduced by a factor of ∼52. In the case of the 6 s long ramp, we measured a final temperature of the cloud of 130 nK, a factor 12.5 below the initial one. This is consistent with the expected value of 15. The small difference may arise from a small heating rate due to the fluctuations of the magnetic trap.
The fact that the shortcut decompression still produces sizeable excitations is due to experimental imperfections. Several possible causes can be invoked. Firstly, as seen in figure 6, there are still small deviations from the ideal trajectory. These may have an impact, especially in the last phase of the trajectory where the cloud is subject to a large acceleration (see figure 3). Secondly, as can again be seen in figure 3, during the trajectory the cloud wanders quite far (several hundreds of µm) from the trap center and feels the non-harmonic part of the potential. This effect is difficult to quantify since our knowledge of the potential shape is not sufficiently accurate (however, the anharmonicity could be inferred from variations of the oscillation frequency with amplitude). In principle, it could be avoided by designing other shortcut trajectories keeping the cloud closer to the trap center at all times.  Since the shortcut trajectory was designed only for the radial dimensions, the resulting axial breathing mode is of the same magnitude as for the linear decompression.
We compare in figure 9 the results of the shortcut decompression to linear ramps of various durations. The vertical axis in this figure represents amplitudes of oscillations after trap decompression, either of the center-of-mass position (filled symbols) or of the cloud radius (open symbols), scaled by their values for an abrupt decompression (t f ∼ 0.1 ms). The horizontal axis is the duration of the decompression t f (notice the logarithmic scale). The circles correspond to linear decompressions, while the stars are the shortcut results. As can be seen, fulfilling the adiabaticity criterion is easier for the breathing mode (size oscillation) than for the dipole mode (center-of-mass oscillation): the oscillation amplitude is reduced by a factor of 2 for t f = 20 ms for the earlier and for t f 150 ms for the latter. Using the amplitude of the dipole mode as a criterion to compare the linear and shortcut schemes, one sees that the decompression time is reduced by a factor of 37.

Shortcut to adiabaticity for an interacting condensate
As opposed to the previous case of non-interacting atoms, the decompression of a BEC is an intrinsically 3D problem because of the interactions. As a result, both the radial and axial frequencies have to be varied following equations (48) and (49) in order to realize a shortcut to adiabaticity. In the present section, we describe a decompression experiment based on the trajectories discussed in section 3.2.1 and represented in figure 4. In this scheme, the radial frequency is decreased by a factor of 9, while the axial frequency is adjusted to maintain the axial size of the BEC fixed during the whole trajectory. Accordingly, the axial frequency is decreased by a factor of 3.
We start from an initial BEC containing 1.3 × 10 5 atoms in the condensed fraction and 7 × 10 4 non-condensed atoms at a temperature of 130 nK. The experimental scheme is similar to that employed for the thermal cloud. Here, we use a longer time of flight of 28 ms to characterize Contrary to the previous case of a thermal cloud, the BEC cannot be held for more than 150 ms in the compressed magnetic trap because of a relatively high heating rate. Thus, here we cannot compare our scheme to the adiabatic limit corresponding to a slow linear decompression. Figure 10 shows the temporal behavior of the cloud following the linear and shortcut decompressions. These absorption images are taken in the (y, z) plane, after a certain holding time in the decompressed trap (indicated in the figure) plus a 28 ms long time of flight. The field of view is 545 µm × 545 µm. The center-of-mass motion has been subtracted from these data for better clarity. In the linear case, the BEC (yellow central part) experiences large deformations and oscillations of its aspect ratio, whereas in the shortcut case it remains nearly perfectly stationary. Surprisingly, in the case of the linear decompression the BEC also oscillates angularly. We attribute this to an uncontrolled tilt of the trap axes during the decompression. This will be discussed in more detail later. The nearly isotropic aspect of the BEC after the shortcut decompression is due to the value of the time of flight, which is close to the critical time of inversion of the aspect ratio. The thermal component surrounding the BEC (red halo) is also visible. Its temporal evolution is discussed at the end of this section.
To provide a more quantitative analysis, the column densities obtained from the absorption images were fitted with a 2D bimodal distribution consisting of a Gaussian component, accounting for the thermal fraction, plus a 3D inverted parabola integrated along one dimension, accounting for the condensed atoms. The fitting parameters were the cloud center, two angles, one for each couple of eigenaxes of each component, and the two widths of each component. In figure 11(a) is reported the center-of-mass oscillations (dipole mode) for the abrupt (squares), linear (diamonds) and shortcut (circles). Figure 11(b) shows the oscillations of the BEC aspect ratio (breathing mode). All measurements are carried out after a 28 ms time of flight. As in the case of the non-interacting cloud, the shortcut scheme reduces the amplitude of the dipole mode compared to a standard linear decompression, here by a factor of 4.3. For our relatively long time of flight, the measured positions reflect the atomic velocities. Thus, using the shortcut scheme reduces the kinetic energy associated with the dipole mode by a factor of 18.5 compared to the linear one (and 36 compared to the abrupt). The residual energy after the shortcut decompression is 580 nK. As can be seen in figure 11(b), both non-optimal schemes induce very large oscillations of the BEC aspect ratio, with a rather complicated dynamics. A Fourier analysis reveals a main oscillation frequency of 47 Hz, consistent with a radial breathing mode at 2ω ⊥ [25][26][27]. A smaller contribution at 12.5 Hz corresponds to the axial breathing mode at √ 5/2ω [27]. The shortcut scheme suppresses strikingly these breathing oscillations, yielding a BEC close to the targeted equilibrium state.
As emphasized in section 3.2, the shortcut trajectory employed in this experiment is also valid for the thermal fraction, in the radial dimensions only. This is demonstrated in figure 12, where we compare the oscillations of the radial (open symbols) and axial (filled symbols) sizes of (a) the BEC and (b) the thermal fraction, after the shortcut decompression. The BEC TF radius is stationary with an average value of 46.8 µm close to the theoretical value (43 µm). As can be observed in figure 12(b), the radial size of the thermal fraction is also quite stationary as expected from a shortcut trajectory. Thus, this experiment demonstrates that both a non-interacting thermal gas and an interacting BEC can be decompressed simultaneously using an appropriate shortcut trajectory. The observed behavior is also qualitatively consistent with our initial assumption that the BEC and thermal fraction are independent. However, we expect that ultimately the validity of this approach will be limited by the interaction between the condensed and non-condensed fractions. The temperature inferred from the radial size of the thermal component is 22 nK, a factor of 6 below the initial one. This factor is smaller than the expected one [ω ⊥ (0)/ω ⊥ (t f ) = 9], and even if we improve the experimental setup to realize the ideal frequency trajectory we would probably be limited by the transfer of energy from the axial breathing mode via the interaction with the condensate. Indeed, the axial size of the thermal fraction presents clear breathing oscillations, reflecting the fact that the shortcut trajectory ω (t) is not valid in this case, as expected.
A striking feature of figure 10(a) was the large angular oscillation of the BEC after the linear decompression. This unexpected effect is due to a slight tilt of the QUIC trap eigenaxes (3 • ) in the (y, z) plane as the trap center moves downwards due to gravity. Because of this, an angular momentum is imparted to the atoms during the decompression, exciting a 'scissors mode' [28,29]. Our nearly critical time of flight then results in a magnification and a deformation of the scissors oscillations [30,31]. Figure 13 shows an example of these oscillations, together with a GPE simulation (red line).

Other possible applications
In this section, we attempt to generalize the shortcut decompression of BECs to other situations which may find applications in experiments where a fast and large modification of the width of the velocity distribution or of the chemical potential is required.

Arbitrary variation of a harmonic potential
Let us consider the time evolution of a condensate in the time-dependent harmonic potential of the form where the symmetric matrix W (t) = R −1 (t)W (t)R(t) represents the harmonic potential of stiffnessW rotated by a rotation matrix R(t). The column vectors r and u, respectively, represent the position and a spatially homogeneous force which may depend on time. The superscript t indicates the transpose of vectors or matrices.

23
To solve equation (24) we look for a linear change of variables ρ(r, {b i j (t)}, {r cm i (t)}) where the b i j 's are scaling and rotation functions for the r i 's. Let B be a 3 × 3 matrix whose elements are the functions b i j . The transformation is + a(t). (55) In the TF limit, and if the matrixḂ B −1 is symmetric, equation (24) is invariant under this transformation. The full derivation is given in appendix B, but we give here the key elements. The TF approximation consists in neglecting the kinetic-energy-like term χ(ρ, τ ) being defined as in equation (25). In this regime, the condensate wavefunction χ(ρ, τ ) verifies the equation of motion (31), under the action of the time-independent potential if the generic scaling functions satisfÿ It is worthwhile recalling that, as shown by the above equations, the evolution of B is decoupled from the center-of-mass motion which evolves with the net external force. The phase of the wavefunction is chosen as The wavefunction normalization is and the time τ is defined by The derivation of the scaling equations (appendix B) relies on the particular choice of the above phase φ, which verifies v(r) being the velocity field of the condensate, and on the assumption that the matrixḂ B −1 is symmetric. The first condition consists in imposing that there are no terms linear in momentum in the GPE in the ρ-coordinate frame; if the first condition is fulfilled the second imposes that the velocity field is irrotational, namely that the condensate is a superfluid everywhere as well. This implies that our scaling ansatz does not take into account the presence of quantized vortices and thus can describe the dynamics of a rotated condensate only below the critical angular velocityα c 0.7ω x for a slightly anisotropic confinement [32], or in general, for a metastable configuration [33]. Nevertheless, a slightly modified ansatz may be found to incorporate the possibility of quantized vortices. It is also possible to relax the first condition and allow for terms in the GPE that contain for instance the angular momentum components. These extensions are deferred for future studies.
Equations (58) and (59) can be used to determine the dipolar, compressional and scissors modes for a harmonically trapped superfluid condensate (see appendix C). Replacing det B with (det B) β in equation (58), the same equation describes the compression and the scissors dynamics of a superfluid characterized by an equation of state µ(n) ∝ n β , as has already been shown for the quadrupolar modes [34] and as can be easily deduced from equation (B.8) of appendix B. In the following we present three possible shortcut trajectories based on these scaling equations and adapted to compress or decompress and rotate a BEC in the absence and in the presence of gravity.

Uniform decompression or compression of a condensate
We now consider the particular case of u = 0 and W diagonal. If one wants to compress or decompress the condensate without modifying the condensate aspect ratio, the condition ω i (t f ) = ω i (0)/γ 2 must hold for any i. The BCs for the shortcut solution are: with c 0 = 1, c 1 = c 2 = 0, c 3 = 10(γ 4/5 − 1), c 4 = −15(γ 4/5 − 1), c 5 = 6(γ 4/5 − 1). The time evolution of the trap frequencies ω i (t) will be given by the equation If the kinetic energy term (56) is negligible during the whole decompression, the final state is a BEC at equilibrium with a chemical potential that has been divided by a factor of γ 16/5 (because µ ∝ ( i ω i ) 2/5 ).

General compression or decompression in the presence of gravity
We now consider the case when , ω y (t) = ω (t) and u z = mg. A general compression or decompression of a condensate confined in this axially symmetric trap (43) can be realized in two steps: (i) in the first step (t ∈ [0,t ]), b is fixed as outlined in section 3.2, while the desired final value of t f ]) b ⊥ is fixed and b evolves according to the set of equations: where c(t) = −ω 2 ⊥ (t)r cm z (t)/[gb ⊥ (t)] as in equation (44). Also in this case one can write the function c(t) as a polynomial of order 9 (see equation (42)) with the first coefficient fixed at one and the following four coefficients fixed at zero. The other coefficients are fixed by the BCs at the time t f of the function c(t) and of the function b (t), which from equation (69) can be written as and by the BCs of their derivatives at the same time t f .

Rotation of the BEC in the presence of gravity
Now we propose a shortcut trajectory to rotate an axially symmetric BEC of an angleᾱ, in the presence of gravity. In this case, and Let us suppose, for instance, ω ⊥ (0) < ω (0), with ω (0) = λω ⊥ (0). The tilted ground state for the potential W (t f ) can be obtained in two steps: (i) during a timet, fixing b , decompressing the BEC in the radial direction up to the value b ⊥ (t) = λ −1 . At t =t the trap is spherical with frequencyω = λω (0) and the BEC is spherical with a TF radius equal to R (0). (ii) Fixing b along the direction y , compressing in the direction x and z , where the axis r is defined by r = Rᾱr. Using the new coordinate reference frame and setting c z (t) = −ω 2 r cm z (t)/[gb ⊥ (t) cosᾱ] and c y (t) = −ω 2 r cm y (t)/(g sinᾱ), we obtain the set of equationsb the latter describing the center-of-mass motion in the y direction. The BCs for such a problem are: and that all the first and the second derivatives with respect to time are null at t =t and t f . In this case a finite-order polynomial ansatz in τ for c i was found to be inadequate as a solution of the scaling equations due to the coupling of c y and c z . A full numerical solution of the dynamical equation using, e.g., a shooting method [35] or following a strategy such as that implemented in optimal control [9] may be needed for finding a shortcut trajectory in this case.

Conclusion
We have experimentally demonstrated the controlled transfer of trapped ultracold atoms between two stationary states using a faster-than-adiabatic process which reduces the transfer time down to a few tens of milliseconds. The transfer is achieved by engineering specific trajectories of the external trapping frequencies that take explicitly into account the spatial shift introduced by gravity. This scheme was successfully applied both to a thermal gas of atoms and to an almost pure BEC. The scheme used is flexible enough to be adapted to both situations even though in the thermal gas interactions do not play a significant role, while the BEC is strongly affected by the s-wave scattering of atoms. The residual excitations observed after the shortcut decompressions in the present demonstration experiments are due to our imperfect control over the time-varying magnetic trapping potential and could be substantially reduced in future realizations. Theoretically, the design of the transfer process was based on the invariants of motion and scaling equation techniques which turned out to be possible thanks to the harmonic shape of the external potential. In our scheme, the invariant of motion technique (for non-interacting particles) and the scaling equation technique (valid for both the non-interacting and the interacting gas) are tightly connected. The invariant of motion we used is a time-independent HO Hamiltonian that can be obtained by a time-dependent canonical transformation of position and momentum. In the scaling equation technique, we looked for a scaling plus shift transformation of the coordinate that allowed the equation of motion for the system to be time independent (except for terms that are not coordinate dependent). In both cases the whole dynamics is included in the new set of (canonical) coordinates, which depend on the trap frequencies. We also showed that these techniques can be generalized to include the rotation of the eigenaxes without much effort.
Very often, in cold-atom experiments, samples are prepared by transferring atoms from some confinement to another, e.g. from a MOT to a magnetic quadrupolar trap, from a quadrupolar trap to a Ioffe-Pritchard trap, from a harmonic confinement to an optical lattice, etc, the main limitation being that, for short transfer times, parasitic excitations may show up. The main application of our scheme is to guide this transfer in order to prepare a very cold sample in a very short time with the desired geometry and without exciting unwanted modes. The shortcutto-adiabaticity scheme proposed here could be applied to non-interacting particles such as cold gases or ultracold spin-polarized fermions, to normal or superfluid (bosonic or fermionic as well) gases in the hydrodynamic regimes, and to strongly correlated systems such as the Tonks gas. In this paper we focused on explicit solutions to transfer atoms between two stable states, but the same strategy could be applied to control the generation of metastable states, vortex states or some exotic out-of-equilibrium states. We plan to explore these possibilities in future studies.
We look for the conditions that A, B and a have to verify aiming to simplify equation By imposing the quadratic term in ρ to be equal to m 2 ρ t W 0 ρ, we obtain the condition (iv), i.e. equation (58); the fifth condition is that the linear term in ρ vanishes and thus leads to (59); finally, by requiring that the ρ-independent term be null, we obtain (61) for φ 0 .