Fast transitionless expansions of cold atoms in optical Gaussian beam traps
E. Torrontegui, Xi Chen, M. Modugno, A. Ruschhaupt, D. Guéry-Odelin, J. G. Muga
aa r X i v : . [ qu a n t - ph ] O c t Fast transitionless expansions of cold atoms in optical Gaussian beam traps
E. Torrontegui, Xi Chen,
1, 2
M. Modugno,
3, 4
A. Ruschhaupt, D. Gu´ery-Odelin, and J. G. Muga
1, 7 Departamento de Qu´ımica F´ısica, Universidad del Pa´ıs Vasco - Euskal Herriko Unibertsitatea, Apdo. 644, Bilbao, Spain Department of Physics, Shanghai University, 200444 Shanghai, P. R. China Departamento de F´ısica Te´orica e Historia de la Ciencia,Universidad del Pa´ıs Vasco - Euskal Herriko Unibertsitatea, Apdo. 644, Bilbao, Spain IKERBASQUE, Basque Foundation for Science Institut f¨ur Theoretische Physik, Leibniz Universit¨at Hannover, Appelstra β e 2, 30167 Hannover, Germany Laboratoire Collisions Agr´egats R´eactivit´e, CNRS UMR 5589, IRSAMC,Universit´e Paul Sabatier, 118 Route de Narbonne, 31062 Toulouse CEDEX 4, France Max Planck Institute for the Physics of Complex Systems, N¨othnitzer Str. 38, 01187 Dresden, Germany
We study fast expansions of cold atoms in a three-dimensional Gaussian-beam optical trap. Threedifferent methods to avoid final motional excitation are compared: inverse engineering using Lewis-Riesenfeld invariants, which provides the best overall performance, a bang-bang approach, and afast adiabatic approach. We analyze the excitation effect of anharmonic terms, radial-longitudinalcoupling, and radial-frequency mismatch. In the inverse engineering approach these perturbationscan be suppressed or mitigated by increasing the laser beam waist.
PACS numbers: 03.75.-b,03.65.-w,03.65.Nk
I. INTRODUCTION: DRIVEN EXPANSIONS OFCOLD ATOMS
Driving the expansion of a Bose-Einstein condensateor more generally of a cold atom cloud in a controlledway, by implementing a designed time-dependence of thetrap frequency, is a common and basic operation in acold-atom laboratory. The expansion may be aimed atdifferent goals, such as decreasing the temperature [1, 2],adjusting the density to avoid three body losses [3], fa-cilitating temperature and density measurements [3], orchanging the size of the cloud for further manipulations[4]. Expansions, isolated [5] or as part of refrigerationcycles [6], are also important from a fundamental pointof view to quantify the third principle of thermodynam-ics. In general, changes of the confining trap will excitethe state of motion of the atoms unless they are donevery slowly or, in the usual quantum-mechanical sense ofthe word, “adiabatically”, but slow-change processes arealso prone to perturbations and decoherence, or imprac-tical for performing many cycles. Engineering fast ex-pansions without final excitation is thus receiving muchattention recently both theoretical and experimental [5–22]. Most theoretical treatments so far are for idealizedone-dimensional (1D) systems, but the implementationrequires a three-dimensional (3D) trap [11, 14, 15], inprinciple with anharmonicities and couplings among dif-ferent directions.In this paper we address these important aspects to im-plement the expansion in practice. Specifically we shallmodel a simple physical realization based on an elongated Compressions may also play a role in many experiments but theirtreatment is similar to expansions so we shall not deal with themhere explicitly. zx O y FIG. 1: (Color online) Schematic view of the optical trap. cigar-shaped optical dipole trap with cylindrical symme-try, see Fig. 1. This trap is formed by a single laserwhich is red detuned with respect to an atomic transi-tion to make the potential attractive (Section II), and ischaracterized in the harmonic approximation by longitu-dinal and radial frequencies. While magnetic traps allowfor an independent control of longitudinal and radial fre-quencies [11, 14, 15], this is not the case for a simple lasertrap that therefore requires a special study. We assumethat the time-dependence of the longitudinal frequencyis engineered to avoid final excitations with a simple 1Dharmonic theory (Section III) and analyze the final fi-delity in the actual trap. Even though for full 3D-resultswe resort to a purely numerical calculation in Section V,an understanding of the effects involved is achieved firstby analyzing separately longitudinal and radial motionsin Section IV. Conclusions and open questions are drawnin the final Section VI.
II. THE MODEL
The intensity profile of a Gaussian laser beam in theparaxial approximation is given by I ( r, z, t ) = I ( t ) e − r /w ( z )
11 + z /z R , (1)where r and z are the radial and longitudinal coordinatesrespectively, and the variation of the spot size w with z is given by w ( z ) = w s (cid:18) zz R (cid:19) , where z R = πw /λ is the Rayleigh range, w the waist,and λ the laser wavelength. The paraxial approximationis valid for waists larger than about 2 λ/π .If an atom is placed in this field with a detuning δ which is large with respect to the Rabi frequency Ω andthe inverse of the lifetime of the excited state Γ = τ − ,but small with respect to the transition frequency, inter-nal excitations and counter-rotating terms are negligible,and the potential that the ground state atoms feel dueto the dipole force is proportional to the laser intensity, V ( r, z, t ) ≃ ~ Ω δ = ~ Γ8 Γ δ II sat , | δ | ≫ Ω , Γ , (2)and inversely proportional to δ . I sat = πhc/ (3 λ τ ) is thesaturation intensity. Combining Eqs. (1) and (2), andadding for convenience in later expressions the physicallyirrelevant term V ( t ), the potential takes finally the form V ( r, z, t ) = − V ( t ) e − r /w ( z )
11 + z /z R + V ( t ) , (3)where V ( t ) = I ( t ) ~ Γ / (8 δI sat ). The analysis carriedout here is not restricted to a two-level atom. Indeed,the dipole force is always proportional to the intensityfor a multi-level atom but with a coefficient that dependson the level structure and the light polarization [23]. Asan example of the magnitudes involved in practice, forrubidium-87 atoms, a Gaussian beam linearly polarizedof waist 8 µ m, power P = 30 W and wavelength 1060 nmprovides a longitudinal trap frequency ω z / π ≃ Indeed other regimes could be used as well [23]. For larger detun-ings counter-rotating terms become important, but the dipole po-tential is still proportional to the laser intensity. In the extremelimit of quasi-electrostatic traps the frequency of the trappinglight is much smaller than the resonance frequency precludingthe possibility of a potential sign change from attractive to re-pulsive; this change is in principle allowed by Eq. (2) by changingthe sign of the detuning across resonance, but see the discussionbelow. This requires that the vertical shift of the trap center in thegravity field is small with respect to the waist, or g ≪ w ω R ,where ω R is the radial angular frequency, and g the accelerationdue of gravity. To solve the time dependent Schr¨odinger equation as-sociated with the potential in Eq. (3), i ~ ∂ Ψ ∂t = − ~ m ∇ Ψ + V Ψ , (4)we use cylindrical coordinates, ∇ ≡ r ∂∂r (cid:18) r ∂∂r (cid:19) + 1 r ∂ ∂φ + ∂ ∂z . With the standard separation of variables Ψ = F ( r, z, t )Φ( φ ), we find − ~ m (cid:20) r (cid:18) F r F + r F rr F (cid:19) + 1 r Φ φφ Φ + F zz F (cid:21) + V = i ~ F t F , (the simple or double subscripts denote first or secondderivatives) whereas the equation for Φ readsΦ φφ Φ = − ν ν ∈ Z , (5)where ν is a magnetic quantum number determining theconserved angular momentum component along the z di-rection. This is solved by e iνφ so we shall concentrate on F ( r, z, t ). Defining now ˜Ψ( r, z, t ) = √ r F ( r, z, t ) we get i ~ ∂ ˜Ψ ∂t = − ~ m ˜ ∇ ˜Ψ + ˜ V ˜Ψ , (6)where ˜ ∇ ≡ ∂ ∂r + ∂ ∂z , (7)˜ V ≡ ~ m (cid:18) ν − / r (cid:19) + V ( r, z, t ) , (8)and V ( r, z, t ) is given by Eq. (3). The task is now todesign V ( t ) so as to achieve a fast expansion withoutfinal excitation. III. EXPANSION PROTOCOLS
In this section we briefly present three expansion pro-tocols for a one-dimensional harmonic trap. They arechosen mostly because they have been realized experi-mentally [1, 11, 14, 24]. Their simplicity will help usto identify relevant physical phenomena that will affectalso other protocols. No claim of completeness or globaloptimization is made, and there is certainly room for in-vestigating other protocols.The objective is to expand the trap from an initial fre-quency ω / π to a final frequency ω f / π in a time t f ,without inducing any final excitation in the adiabatic orinstantaneous basis for which the Hamiltonian is diago-nal. Note that transient excitations are generally allowed. A. Inverse engineering with Lewis-Riesenfeldinvariants
Lewis-Riesenfeld invariants may be used to engineerefficient expansion protocols as explained in [5, 7, 8, 11].We refer the reader to these works for the details andgive here only the elements required for a practical ap-plication. For any harmonic oscillator expansion, and infact for any potential with the structure V ( q, t ) = m ω ( t ) q + 1 b ( t ) U (cid:18) qb ( t ) (cid:19) , (9)where the function U is arbitrary, there is a quadratic-in-momentum invariant operator of the form I ( q, p, t ) = 12 m [ b p − m ˙ bb ( qp + pq ) + m ˙ b q ]+ 12 mω (cid:16) qb (cid:17) + U (cid:16) qb (cid:17) . (10)where q and p are generic position and momentum op-erators (particularized later for longitudinal or radial di-rections), and b = b ( t ) is a scaling function that satisfiesthe Ermakov equation¨ b + ω ( t ) b = ω /b . (11)To inverse engineer the trap frequency, b ( t ) is designed tosatisfy the boundary conditions that guarantee no finalexcitations and continuity of the trap frequency, b (0) = 1 , ˙ b (0) = 0 , ¨ b (0) = 0 ,b ( t f ) = γ, ˙ b ( t f ) = 0 , ¨ b ( t f ) = 0 , (12)where γ = p ω /ω f . Finally ω ( t ) is deduced from Eq.(11). When these conditions are satisfied the n -th eigen-state of H (0) = p / (2 m ) + V ( q,
0) evolves as an “ex-panding mode”, which is the time-dependent n -th eigen-vector of the invariant times a “Lewis-Riesenfeld” phasefactor e iα n , where α n ( t ) = − ( n + 1 / ω R t dt ′ /b , untilit becomes an n -th eigenvector of the final Hamiltonian H ( t f ) = p / (2 m ) + V ( q, t f ).Discontinuities in ¨ b may in principle be allowed, butthey amount to perform finite jumps of the trap fre-quency. This is an idealization but it can be approachedup to technical limits. Discontinuities in ˙ b have more se-rious consequences, as they imply delta functions for ¨ b and infinite trap frequencies at the jumps. Thus theirapproximate physical realization is even more difficult,but they are useful to find bounds for physical variablesof interest as in [5] or in Appendix A. There are many b ( t ) that satisfy Eq. (12). A simple smooth function is aquintic polynomial with six coefficients, b ( t ) = 6 ( γ − s −
15 ( γ − s +10 ( γ − s +1 , (13)where s := t/t f . More sophisticated choices are possibleto minimize the process time, the average energy, or other variables with or without imposed constraints (“boundedcontrol”) on the allowed trap frequencies [12, 25].In principle this method allows for arbitrarily small val-ues of t f but for very small process times the potentialbecomes transitorily repulsive [8]. With a laser beam, arepulsive potential could be achieved by means of a pos-itive (blue) detuning. However the validity of Eq. (3)requires a large detuning, so care should be exercisedwhen crossing the resonance. In this work the numericalexamples are restricted to positive trap frequencies sincethe transient passage from attractive to a repulsive inter-action could involve unwanted effects such as radiationpressure. We shall in any case point out in the figuresthe minimum value of t f for which the potential staysattractive for all times with the quintic b ( t ). B. Bang bang
Bang-bang methods are those with a stepwise constantbehavior of some variable. The simplest one to avoid finalexcitations consists of using a constant intermediate trapfrequency, the geometric average of the initial and finalfrequencies, ω bang = ω , t ≤ ω ω f ) / , < t < t bangf ω f , t ≥ t bangf = π/ [2( ω ω f ) / ] . (14)during a fourth of the corresponding period. Optimalbang-bang trajectories with two intermediate frequenciesand arbitrary final times are also possible [6, 8] but theyare not considered here. To arrive at Eq. (14), one mayuse a classical argument [24] or proceed as follows: thesolution of the Ermakov equation for t > ω with initialboundary conditions b (0) = 1, ˙ b (0) = 0 is b ( t ) = q [( ω − ω ) /ω ] sin ( ω t ) + 1 . (15)Then, if we solve for t f and ω the two equations impliedby the final boundary conditions b ( t f ) = γ, ˙ b ( t f ) = 0 , (16)the values in Eq. (14) are found. C. Fast adiabatic protocols “Fast” and “adiabatic” may appear as contradictoryconcepts. A fast adiabatic protocol seeks to perform aprocess as quickly as possible keeping it adiabatic at alltimes. For harmonic oscillator expansions the adiabatic-ity condition reads ˙ ω/ω ≪
1, so making ˙ ω/ω constantfrom t = 0 to t f [1], and solving the resulting differentialequation for ω ( t ), we get [8] ω adi ( t ) = ω − ( ω f − ω ) tt f ω f . (17) IV. LONGITUDINAL AND RADIAL MOTIONS
By expanding the potential in Eq. (3) in a doubleTaylor series around ( z = 0 , r = 0), V ( r, z, t ) ≃ − V ( t ) (cid:18) − r w − z z R + 2 r w + z z R + 4 r z w z R + · · · (cid:19) , (18)we see that the first coupling term between radial andlongitudinal motions is of fourth order, proportional to r z . If we could neglect this and higher order couplingterms, longitudinal and radial motions would be indepen-dent. The approximation that considers longitudinal andradial motions completely uncoupled is useful to analyzeseveral effects separately and gain insight. Not only that,as we shall confirm later with full 3D calculations, thisis also a good approximation numerically for low energylevels. A. Longitudinal motion
To study the motion in the longitudinal direction weconsider first the full longitudinal Hamiltonian (putting r = 0 in Eq. (3)) H ( z, t ) = − ~ m ∂ ∂z − V ( t ) (cid:18)
11 + z /z R − (cid:19) . (19)To characterize this potential in terms of a longitudinalfrequency consider the harmonic approximation H har ( z, t ) = − ~ m ∂ ∂z + V ( t ) z z R = p z m + mω z ( t ) z , (20)where p z = − i ~ ∂/∂z and ω z ( t ) = 2 V ( t ) mz R . (21)We may now apply the methods described in Sec. III.We shall add a subscript z , for “longitudinal”, to thetrap frequencies and take q → z . For an imposed ω z ( t )and fixed waist and laser frequency, Eq. (21) determinesthe time dependence of the laser intensity.Figure 2 shows for the ground state, n = 0, the “lon-gitudinal fidelity” F L = |h Z n ( t f ) | U z ( t f , | Z n (0) i| , where Z n ( z,
0) and Z n ( z, t f ) are the initial and final n -th eigen-states of the full longitudinal Hamiltonian, Eq. (19),for frequencies ω z and ω fz respectively, and U z ( t f , ω z / π = 2500 Hz, and the wavelength is taken as λ = 1060 nm, characteristic of Neodymium-doped lasers. çç óó (a) t f H ms L F L çç óó (b) t f H ms L F L FIG. 2: (Color online) Longitudinal fidelity for Eq. (19) ver-sus final time t f for two different final frequencies: in (a) ω fz / π = 250 Hz; in (b) ω fz / π = 25 Hz. In both figures ω z / π = 2500 Hz, λ = 1060 nm, and the initial state is theground state. All curves have been calculated for two differentwaists, w = 10 µm and w = 3 µm , but the results are indis-tinguishable. V ( t ) is chosen according to the three protocols:inverse engineering, using Eq. (13) in Eq. (11) (solid red line);bang-bang (circles for w = 10 µm and triangles for w = 3 µm ); and fast adiabatic method (long-dashed blue line). Thegreen vertical lines denote the minimum t f for which V ( t )remains positive for all t in the inverse engineering approachwith a quintic b ( t ). The vertical green lines in this and the following figurescorrespond to the minimum t f value for which ω z ( t ) ispositive for all t using the quintic b ( t ). Of course thislimit only applies to the inverse engineering approach.The three methods are compared for two different fi-nal frequencies, in (a) ω fz / π = 250 Hz and in (b) ω fz / π = 25 Hz, and two different waists, see the captionfor details. Actually the effect of the waist change is neg-ligible in the scale of these figures. Globally the inverseengineering method outperforms the others. The bang-bang approach provides a good fidelity, but only for aspecific final time, and the fast adiabatic method fails infact to be adiabatic at short times. This is very evidentfor the smaller final frequency in Fig. 2 (b): the openingof the trap is faster and the level spacings smaller thanfor the larger final frequency, so the state cannot followan adiabatic behavior. ó ó ó ó ó ó ç ç ç ç ç ç ò ò ò ò ò òè è è è è è x x x x x x + + + + + + F L FIG. 3: (Color online) Longitudinal fidelity versus level num-ber n . For w = 3 µm (blue symbols): fidelity F L = |h Z n ( t f ) | U z ( t f , | Z n (0) i| with inverse engineering computedfor Eq. (19) and a quintic b ( t ) (empty triangles); first-orderbound in Eq. (A9) (x); second-order approximate fidelitycomputed with Eqs. (A11) and (A12) for the quintic b ( t )(filled triangles). For a larger waist, w = 10 µm (red sym-bols), the corresponding values (empty circles, crosses, andfilled circles), are very nearly one and indistinguishable. Inall cases ω fz / π = 25 Hz, ω z / π = 2500 Hz, λ = 1060 nm,and t f = 2 . For small t f and/or large γ the transient excitationenergy during the inverse engineering protocol may behigh, giving the atoms access to the anharmonic part ofthe potential. This could lead to a decay of fidelity as t f decreases, but the effect is not seen in the scale of Fig.2. The small effect of anharmonicity is enhanced by in-creasing the vibrational number, see Fig. 3. It can beavoided by increasing the waist, which reduces the anhar-monic terms, as demonstrated in Fig. 3. The AppendixA provides a perturbation theory analysis of longitudinalanharmonicity. A first order approach gives analyticallower bounds for the fidelity, and the possibility to maxi-mize them by other choices of b ( t ). More accurate resultsare found in second order, see Fig. 3. B. Radial motion
To study radial motion we define the radial Hamilto-nian by setting z = 0 in Eq. (3), and add the “centrifugalterm” (an attractive one for ν = 0), H ( r, t ) = − ~ m ∂ ∂r − V ( t ) (cid:16) e − r /w − (cid:17) + ~ m (cid:18) ν − / r (cid:19) , (22)In the harmonic approximation we take − V ( t ) (cid:16) e − r /w − (cid:17) ∼ mω R ( t ) r , (23) where ω R ( t ) ≡ V ( t ) mw , (24)defines the radial frequency of the trap ω R ( t ) / (2 π ).We assume that V ( t ) is set to satisfy the designed lon-gitudinal expansion according to Eq. (21). Substitutingthis into Eq. (24) gives the relation between radial andlongitudinal frequencies, ω R ( t ) = √ πw λ ω z ( t ) , (25)which is key to understand the behavior of the radialwavefunction. The waist value w = λ/ ( √ π ) wouldmake both frequencies equal but for such a small waistthe paraxial approximation fails and the present theorycould not be applied [26].Now we define the ground state “radial fidelity” as F R = |h φ ( t f ) | U R ( t f , | φ (0) i| , where U R ( t f ,
0) is theradial evolution operator for H ( r, t ), and φ ( r,
0) and φ ( r, t f ) are ground states of the initial and final radialtraps for Eq. (22), with V ( t ) given by Eq. (21). InFig. 4 the radial fidelity is depicted for the three pro-tocols explained above and ν = 0. The bang-bang tra-jectory causes much excitation and gives a rather poorfidelity in the radial direction (see the big empty sym-bols). This is because, even though the “correct” tran-sient radial frequency ω bangR ( t ) = ( √ πw /λ ) ω bangz ( t ) isimplemented (as the geometric average of initial and fi-nal radial frequencies), the time t f is adjusted for thelongitudinal, not for the radial frequency, compare t bangfz to t bangfR := λt bangfz / ( √ πw ).The behavior of the other two methods is much morerobust as shown by their fidelities in Fig. 4, where theinset amplifies the small- t f region. The fast adiabaticapproach, blue long-dashed lines of Fig. 4, gives an ex-cellent radial fidelity. A detailed calculation shows thatthe adiabaticity condition for the radial Hamiltonian,see the Appendix B, takes, in spite of the centrifugalterm, the same form as for the ordinary harmonic os-cillator, namely ˙ ω R /ω R ≪
1. If ˙ ω z /ω z is constant, theratio ˙ ω R /ω R will be constant too, but smaller by a factor λ/ ( √ πw ), so that the radial dynamics will be in gen-eral more adiabatic. Indeed, comparing blue long-dashedlines of Figs. 2 and 4, we see that the radial fidelity isbetter than the longitudinal one, as the confinement istighter radially than longitudinally, and the energy lev-els are more separated.For the inverse engineering protocol with a quintic b ,red solid and black short-dashed lines in Fig. 4, the radialexcitation is small, the fidelity being nearly one for t f > ω Rinv ( t ) = 2 π w λ ω z b − ¨ bb , (26) FIG. 4: (Color online) Radial fidelity for different final times.Parameters: ν = 0, ω z / π = 2500 Hz, ω fz / π = 250 Hz,and λ = 1060 nm. Protocols: inverse engineering, quintic b , w = 10 µ m (solid red line); inverse engineering, quintic b , w = 3 µ m (short-dashed black line); bang bang, w = 10 µ m(empty circle); bang bang, w = 3 µ m (empty triangle); fastadiabatic, w = 10 µ m (long-dashed blue line); fast adiabatic, w = 3 µ m (dotted orange line). In the inset, corresponding(red) triangles, (black) squares, (blue) circles, and (orange)diamonds are the approximate fidelities from adiabatic per-turbation theory using Eq. (28). The green vertical line near0 . t f for which the quintic- b in-verse engineering protocol implies positive frequencies at alltransient times. where we have used the Ermakov equation and γ =( ω z /ω fz ) / = ( ω R /ω fR ) / , so the function b is thesame for the longitudinal and radial directions. Instead,the actual radial frequency varies according to Eq. (25), ω R ( t ) = 2 π w λ ω z b − ¨ bb ! . (27)The difference between Eqs. (26) and (27) may be quitesignificant, as the example in Fig. 5 shows, so the ra-dial state does not really follow the expanding modescorresponding to Eq. (26). As a consequence a pertur-bative approach based on Eq. (26) as the zeroth orderfails completely, even for qualitative guidance. (Furtherevidence is provided in Fig. 6.) Instead, a valid pertur-bation theory analysis may be based on the zeroth orderof the adiabatic modes, which are instantaneous eigen-states of the radial Hamiltonian. The physical reason forthe relatively good radial behavior of the invariant-basedinverse engineering protocol is thus the adiabaticity inthat direction. When t f approaches the critical lowervalue, ω R ( t ) becomes too steep at small values of ω R , seeFig. 5, and adiabaticity breaks down, which explains thefidelity decrease there. To mitigate this problem, a largerwaist may be used to increase ω R ( t ) for a given ω z ( t ), seeEq. (25), making the radial process more adiabatic andimproving the fidelity, as shown in Fig. 4. Keep in mindthat to implement a given ω z ( t ) an increase of waist hasto be compensated by an increase of laser intensity. H ms L Ω R (cid:144) H Π L H z FIG. 5: (Color online) (Square of) radial frequency versustime: Eq. (27) based on a quintic b (red solid line); fastadiabatic protocol, Eq. (25) with Eq. (17) for ω z (blue longdashed line); Eq. (26) (black short dashed line). Parameters: λ = 1060 nm, ω z / π = 2500 Hz, ω fz / π = 250 Hz, w = 3 µ m and t f = 0 .
36 ms.
Inserting a wave function expansion in an adi-abatic basis (instantaneous eigenstates) ψ ( r, t ) = P k a k ( t ) h r | φ k ( t ) i e − i R t E k ( t ′ ) dt ′ , into the radialSchr¨odinger equation provides a set of coupled dif-ferential equations. Integrating formally for a k ( t ) andassuming that in zeroth order a (0) k ( t ) = δ ,k , we get infirst order [27, 28] a (1)1 ( t ) = − Z t dt ′ h φ ( t ′ ) | ˙ φ ( t ′ ) i e − i ~ R t ′ dt ′′ [ E ( t ′′ ) − E ( t ′′ )] = − Z t dt ′ ˙ ω R ( t ′ )2 ω R ( t ′ ) e i R t ′ dt ′′ ω R ( t ′′ ) , (28)where the second line follows by applying, in addition,the harmonic approximation. All other first order am-plitudes (for k ≥
2) are zero so the radial fidelity at t f can be calculated in this approximation as F R ( t f ) =[1 − | a (1)1 ( t f ) | ] / . This is in fact correct up to second or-der and it provides good agreement with the exact resultsas shown by the symbols in the inset of Fig. 4. V. FULL 3D ANALYSIS
Finally we study numerically the actual, coupled dy-namics driven by the exact full potential, Eq. (3),combining all the effects considered so far and thelongitudinal-radial coupling. The 3D fidelities are alsocompared with simpler 1D fidelities.Before discussing the results a remark on the technical-ities of the numerical calculations of the time-dependentwave functions is in order. To solve the longitudinal timedependent Schr¨odinger equation we approximate the evo-lution operator U L using the Split Operator Method(SOM) [29] and Fast Fourier Transform (FFT) [30, 31].For the radial time dependent Schr¨odinger equation, wehave used instead the finite differences Crank-Nicholsonscheme [31]. In this way the singular point r = 0 can H ms L ov e r l a p FIG. 6: (Color online) (Modulus of) the overlap integral be-tween the state evolving from the ground state with the radialHamiltonian, using the actual radial frequency in Eq. (27),and (a) the instantaneous eigenstate in harmonic approxima-tion (red solid line), or (b) the expanding mode in harmonicapproximation for the radial frequency in Eq. (26) (dashedblue line). Parameters: ν = 0, ω z / π = 2500 Hz, λ = 1060nm, ω fz / π = 250 Hz, w = 3 µ m, t f = 0 . ν = 0, ω z / π = 2500Hz, λ = 1060 nm, ω fz / π = 250 Hz. Protocols: inverseengineering, quintic b , w = 10 µm (red solid line), and radialfidelity (red filled triangles); inverse engineering, quintic b , w = 3 µm (black short-dashed line), and radial fidelity (filledsquares); bang bang, w = 10 µm (empty circle in inset), andradial fidelity (filled circle); bang bang, w = 3 µm (emptytriangle in inset), and radial fidelity (filled triangle in inset);fast adiabatic, w = 3 µm and w = 10 µm (blue long-dashedline), and longitudinal fidelity (blue diamonds). The greenvertical line marks the threshold for positive trap frequenciesin the inverse engineering, quintic- b protocol. be excluded imposing a zero of the wavefunction therefor all t . For the 3D calculation we use similarly theCrank-Nicholson scheme with operator splitting for mul-tidimensions (the radial and longitudinal directions) [31].We calculate, starting from the ground state of theinitial trap for ν = 0, the fidelity of the final state with respect to the ground state of the final trap using thethree expansion protocols. Figure 7 clearly shows thatthe inverse engineering protocol provides the best over-all performance, with a fidelity decay at smaller timesthat can be avoided by increasing the waist. The radial-longitudinal quartic coupling does not play any signifi-cant role for the parameters of the example, and in anycase it may also be suppressed by a waist increase.As for the simple bang-bang approach, it fails in 3D,as expected, due to its poor radial behavior. (The radialfidelity alone is close to the 3D fidelity.) The results forfast adiabatic and inverse engineering methods are nottoo different from each other for the chosen expansion,but the fast adiabatic method will fail sooner for moredemanding expansions with smaller final frequencies, asclearly illustrated in Fig. 2 (b), due to its inability to re-main really adiabatic for the longitudinal direction. Themain limitation of the inverse engineering method is in-stead the possible failure of adiabaticity in the tighterradial direction, so it is intrinsically more robust thanthe fast adiabatic one. Note that the 3D fidelity is mim-icked accurately by the 1D radial fidelity for inverse en-gineering and by the 1D longitudinal fidelity for the fastadiabatic approach. VI. DISCUSSION AND OUTLOOK
In this work we have analyzed the practical implemen-tation of a transitionless expansion in a simple Gaussianoptical dipole trap, taking into account anharmonicities,radial-longitudinal couplings, and the radial-longitudinalfrequency mismatch. The main conclusion of the studyis that the transitionless expansions in optical traps arefeasible under realistic conditions. Despite the relationbetween the longitudinal and transversal trapping fre-quencies through the intensity, the different timescalesenable us to design fast expansions with high fidelitieswith respect to the ideal results using the invariant-basedinverse engineering method, which is particularly suitablecompared to the two other approaches examined. Ourdetailed analysis of radial and longitudinal motions re-veals the weakest points of each approach: for the inverseengineering, the main perturbation is due to the possi-ble adiabaticity failure in the radial direction, which canbe suppressed or mitigated by increasing the laser waist.This waist increase would also reduce smaller perturb-ing effects due to longitudinal anharmonicity or radial-longitudinal coupling. The simple bang-bang approachfails because the time for the radial expansion is badlymismatched with respect to the ideal time, and the fastadiabatic method fails at short times because of the adi-abaticity failure in the longitudinal direction.Complications such as perturbations due to differentnoise types, and consideration of condensates, gravity ef-fects, or the transient realization of imaginary trap fre-quencies are left for separate works. Other extensionsof the present analysis could involve the addition of asecond laser for further control of the potential shape,or alternative traps. Optical traps based on Bessel laserbeams, for example, may be interesting to decouple longi-tudinal and radial motions. Transport protocols are alsoamenable to a similar 3D analysis. Finally, the combina-tion of invariant-based inverse engineering with optimalcontrol theory to optimize the frequency design is also apromising venue of research [12, 25].
Acknowledgments
We thank G. C. Hegerfeldt for useful discussions. Weacknowledge funding by the Basque Government (GrantNo. IT472-10) and Ministerio de Ciencia e Innovaci´on(FIS2009-12773-C02-01). X. C. thanks the Juan de laCierva Programme, the National Natural Science Foun-dation of China (Grant Nos. 60806041 and 61176118)and the Shanghai Leading Academic Discipline Program(Grant No. S30105); E. T. the Basque Government(Grant No. BFI08.151); and D. G. O. the Agence Na-tional de la Recherche, the R´egion Midi-Pyr´en´ees, theUniversity Paul Sabatier (OMASYC project) and the In-stitut Universitaire de France.
Appendix A: Perturbative treatment of longitudinalanharmonicity
To quantify the effect of anharmonicity in the inverseengineering method let us first expand the longitudinalHamiltonian retaining the quartic term, V ( z ) = − V ( t )1 + z /z R + V ( t ) ≈ V z z R − V z z R = 12 mω z ( t ) z | {z } := V u − mω z ( t ) z z R | {z } := V . (A1)The harmonic term is the unperturbed potential and thequartic term the perturbation. A further simplificationto evaluate the fidelities F L = |h Z n ( t f ) | U z ( t f , | Z n (0) i| is to substitute the initial and final exact eigenstates | Z n (0) i and | Z n ( t f ) i by corresponding eigenstates of theharmonic oscillator | ψ n (0) i and | ψ n ( t f ) i , see the expres-sions below. This is an excellent approximation for thelower eigenstates. In Fig. 3, for example, the over-lap probability between exact and approximate states is0.997 in the worst possible case ( n = 5, w = 3 µ m, and ω fz = 25 Hz). In other words, the perturbation does notimply any significant deformation in the initial and finalstates. Its potential impact is instead at transient timeswhere the energy may be high.We thus expand the overlap h ψ n ( t f ) | U L ( t f , | ψ n (0) i as 1 + f (1) n,n + f (2) n,n + ... in powers of V to write F L = (cid:12)(cid:12)(cid:12) f (1) n,n + f (2) n,n + ... (cid:12)(cid:12)(cid:12) . (A2) The first order correction is given by f (1) n,n = − i ~ Z t f dt ′ h ψ n ( t ′ ) | V ( t ′ ) | ψ n ( t ′ ) i = im ~ z R Z t f dt ′ ω L ( t ′ ) h ψ n ( t ′ ) | z | ψ n ( t ′ ) i , (A3)and the unperturbed wavefunction corresponds to theknown harmonic evolution of an expanding mode underthe time-dependent frequency ω z ( t ) from an eigenstateof the initial harmonic oscillator [8], h z | ψ n ( t ′ ) i = e im ~ ˙ bz b ( b n n !) / e − i ( n +1 / ω z R t ′ dt ′′ b t ′′ ) × (cid:16) mω z π ~ (cid:17) / e − mω zz ~ b H n (cid:18)r mω z ~ zb (cid:19) . (A4)This is nothing but the eigenvector of the quadratic in-variant times the Lewis-Riesenfeld phase factor.Using the triangular inequality | x + y | ≥ (cid:12)(cid:12) | x | − | y | (cid:12)(cid:12) , F L ≥ − | f (1) n,n + f (2) n,n + ... | , and assuming that the per-turbative corrections satisfy | f (1) n,n | >> | f (2) n,n | then F L ≥ − | f (1) n,n | . (A5)The relevant matrix element in Eq. (A3) can be calcu-lated explicitly, h ψ n ( t ′ ) | z | ψ n ( t ′ ) i = (cid:18) ~ mω z (cid:19) b n + 1) + n ] . (A6)Using this result and the Ermakov equation, | f (1) n,n | = 3 ~ mz R [( n + 1) + n ] (cid:18) t f − ω z Z t f ¨ bb dt (cid:19) . (A7)To bound − R t f ¨ bb dt subjected to the imposed con-straints on b at t = 0 and t f we follow closely themethod applied in [5] to bound the average energy. Weintegrate first by parts applying ˙ b (0) = ˙ b ( t f ) = 0 andminimize the resulting positive integral R t f b ˙ b dt , usingthe Euler-Lagrange equation ¨ bb + ˙ b = 0 for b (0) = 1, b ( t f ) = γ . These two conditions are fulfilled for a set offunctions larger than the subset of functions that satisfyall boundary conditions in Eq. (12), so the minimizationsets a lower bound for the later subset. The solution is b = [ t ( γ − /t f + 1] / and with it the integral becomes Z t f b ˙ b dt = ( γ − t f . (A8)This b is not as such a physically valid solution, as thederivatives would jump at the edges, but it maximizesthe fidelity lower bound, F L ≥ − (cid:26) ~ λ mπ w [( n + 1) + n ] × (cid:20) t f + 3( γ − t f ω z (cid:21)(cid:27) . (A9)Take note that there is an optimal t f value, and that thefidelity deteriorates when increasing the quantum num-ber and improves by increasing the waist.Using for b the quintic polynomial in Eq. (13) we get,instead of Eq. (A8), Z t f b ˙ b dt = 10( γ − t f (1101 + 1351 γ + 1101 γ ) . (A10)This tends to 0 . γ /t f for γ >>
1, which, up to theconstant factor, is the same dependence found for theoptimal function.The first-order analysis provides analytical results andsome qualitative guidance but for more accurate resultswe should resort to a second order calculation. Insteadof Eq. (A2) we may also write the fidelity for the n − th state in terms of the corresponding probability as F L =(1 − P n = n ′ P n ′ ) / , where P n ′ is the probability to findthe system in the n ′ -th state at t f . This is approximatedto second order using first order non-diagonal terms, F L = vuut − X n = n ′ (cid:12)(cid:12)(cid:12)(cid:12) f (1) n,n ′ (cid:12)(cid:12)(cid:12)(cid:12) , (A11)where f (1) n,n ′ = − i ~ Z t f dt ′ h ψ n ( t ′ ) | V ( t ′ ) | ψ n ′ ( t ′ ) i = i ~ λ π mw ω z α n,n ′ β n,n ′ ( t ) √ π n + n ′ n ! n ′ ! , (A12) α n,n ′ = Z ∞−∞ dy e − y H n ( y ) H n ′ ( y ) y , (A13)and β n,n ′ ( t ) = Z t dt b ( t ) ω z ( t ) e − i ( n ′ − n ) ω z R t dt b t . (A14) Appendix B: The condition of adiabaticity
Consider a time-dependent Hamiltonian ˆ H ( t ), withinstantaneous eigenstates and energies given byˆ H ( t ) | i ( t ) i = E i ( t ) | i ( t ) i . (B1)The adiabaticity condition is given by [27] |h i ( t ) | ∂ t j ( t ) i| ≪ ~ | E i ( t ) − E j ( t ) | , i = j. (B2) Now we apply this general expression to get the adi-abaticity condition on the longitudinal and radial fre-quencies of the Gaussian beam trap. In the harmonicapproximation the longitudinal Hamiltonian is given byEq. (20), H har ( z, t ) = − ~ m ∂ ∂z + mω z ( t ) z , (B3)whereas the radial Hamiltonian is given, according toEqs. (22), (23), and (24), by H har ( r, t ) = − ~ m ∂ ∂r + mω R ( t ) r ~ m (cid:18) ν − / r (cid:19) . (B4)The instantaneous eigenstates and energies of the Hamil-tonian, Eq. (B3), are h z | n ( t ) i = (cid:16) mω z π ~ (cid:17) / r n n ! e − − mωzz ~ H n (cid:18)r mω z ~ z (cid:19) ,E n ( t ) = (cid:18) n + 12 (cid:19) ~ ω z , n = 0 , , ... where ω z = ω z ( t ) and H n is the Hermite polynomial. Toget the adiabaticity condition for the ground state i = 0of the longitudinal Hamiltonian we replace these expres-sions for the eigenstates and energies into Eq. (B2). Theoverlap of the ground state with the first excited state j = 1 is 0 so the first non-vanishing overlap of the groundstate is with the j = 2 state. The adiabatic condition issatisfied if √ ω z ω z ≪ . (B5)For the radial Hamiltonian, Eq. (B4), the instantaneouseigenstates and energies for ν = 0 are [32] h r | k ( t ) i = r mω R r ~ e − − mωRr ~ L k (cid:16) mω R ~ r (cid:17) ,E k ( t ) = (2 k + 1) ~ ω R , k = 0 , , ...k being the radial quantum number and L k the general-ized Laguerre polynomial. Again ω R = ω R ( t ). To get theadiabaticity condition for the radial direction we replacethese expressions for i = 0 and j = 1 into Eq. (B2),˙ ω R ω R ≪ . (B6) [1] A. Kastberg, W. D. Phillips, S. L. Rolston, and R. J. C.Spreeuw, Phys. Rev. Lett. , 1542 (1995). [2] A. E. Leanhardt et al., Science
01, 1513 (2003).[3] W. H¨ansel, P. Hommelhoff, T. W. H¨ansch, and J. Reichel, Nature , 498 (2001).[4] T. Kinoshita, T. Wenger, and D. S. Weiss, Phys. Rev. A , 011602 (2005).[5] X. Chen and J. G. Muga, Phys. Rev. A , 053403(2010).[6] P. Salamon, K. H. Hoffmann, Y. Rezek, and R. Kosloff,Phys. Chem. Chem. Phys. , 1027 (2009).[7] J. G. Muga, X. Chen, A. Ruschhaupt, and D. Gu´ery-Odelin, J. Phys. B , 241001 (2009).[8] X. Chen, A. Ruschhaupt, S. Schmidt, A. del Campo, D.Gu´ery-Odelin, and J. G. Muga, Phys. Rev. Lett. ,063002 (2010).[9] J. G. Muga, X. Chen, S. Ib´a˜nez, I. Lizuain, and A.Ruschhaupt, J. Phys. B , 085509 (2010).[10] S. Masuda and K. Nakamura, Proc. R. Soc. A , 1135(2010).[11] J. F. Schaff, X. L. Song, P. Vignolo, and G. Labeyrie,Phys. Rev. A , 033430 (2010); , 059911(E) (2011).[12] D. Stefanatos, J. Ruths, and Jr-Shin Li, Phys. Rev. A , 063422 (2010).[13] S. Masuda and K. Nakamura, Physical Review A ,043434 (2011).[14] J. F. Schaff, X. L. Song, P. Capuzzi, P. Vignolo, and G.Labeyrie, Europhys. Lett. , 23001 (2011).[15] J. F. Schaff, P. Capuzzi, G. Labeyrie and P. Vignolo,arXiv:1105.2119v1.[16] B. Andresen, K. H. Hoffmann, J. Nulton, A. Tsirlin, andP. Salamon, Eur. J. Phys. , 827 (2011).[17] Y. Li, L.-A. Wu, and Z.-D. Wang, Phys. Rev. A ,043804 (2011). [18] A. del Campo, Phys. Rev. A , 031606(R) (2011).[19] A. del Campo, Eur. Phys. Lett. (accepted),arXiv:1010.2854.[20] X. Chen, E. Torrontegui, and J. G. Muga, Phys. Rev. A , 062116 (2011).[21] S. Ib´a˜nez, S. Mart´ınez-Garaot, X. Chen, E. Torrontegui,and J. G. Muga, Phys. Rev. A , 023415 (2011).[22] S. Choi, R. Onofrio, and B. Sundaram, Phys. Rev. A(accepted), arXiv:1109.4908.[23] R. Grimm, M. Weidem¨uller, and T. B. Ovchinnikov, Adv.At. Mol. Opt. Phys. , 95 (2000).[24] J. M. Vogels, personal communication.[25] X. Chen, E. Torrontegui, D. Stefanatos, Jr-Shin Li, andJ. G. Muga, Phys. Rev. A , 043415 (2011).[26] S. Nemoto, Appl. Opt. , 1940 (1990).[27] L. I. Schiff, Quantum Mechanics (McGraw Hill, NewYork, 1949).[28] M. Demirplak and S. A. Rice, J. Phys. Chem. , 9937(2003).[29] E. Torrontegui, J. Echanobe, A. Ruschhaupt, D. Gu´ery-Odelin, and J. G. Muga, Phys. Rev. A. , 043420 (2010).[30] R. Kosloff, J. Phys. Chem , 2087 (1988).[31] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P.Flanery, Numerical Recipies The Art of Scientific Com-puting rd ed. (Cambridge University Press, Cambridge,2007).[32] P. Camiz A. Gerardi, C. Marchioro, E. Presutti, and E.Scacciatelli, J. Math. Phys.12