Superfluidity and relaxation dynamics of a laser-stirred 2D Bose gas
Vijay Pal Singh, Christof Weitenberg, Jean Dalibard, Ludwig Mathey
SSuperfluidity and relaxation dynamics of a laser-stirred 2D Bose gas
Vijay Pal Singh,
1, 2, 3
Christof Weitenberg, Jean Dalibard,
4, 3 and Ludwig Mathey
1, 2, 3 Zentrum f¨ur Optische Quantentechnologien, Universit¨at Hamburg, 22761 Hamburg, Germany Institut f¨ur Laserphysik, Universit¨at Hamburg, 22761 Hamburg, Germany The Hamburg Centre for Ultrafast Imaging, Luruper Chaussee 149, Hamburg 22761, Germany Laboratoire Kastler Brossel, Coll`ege de France, ENS-PSL Research University,CNRS, UPMC-Sorbonne Universit´es, 11 place Marcelin Berthelot, 75005 Paris, France (Dated: November 5, 2018)We investigate the superfluid behavior of a two-dimensional (2D) Bose gas of Rb atoms usingclassical field dynamics. In the experiment by R. Desbuquois et al. , Nat. Phys. , 645 (2012), a 2Dquasicondensate in a trap is stirred by a blue-detuned laser beam along a circular path around thetrap center. Here, we study this experiment from a theoretical perspective. The heating induced bystirring increases rapidly above a velocity v c , which we define as the critical velocity. We identify thesuperfluid, the crossover, and the thermal regime by a finite, a sharply decreasing, and a vanishingcritical velocity, respectively. We demonstrate that the onset of heating occurs due to the creationof vortex-antivortex pairs. A direct comparison of our numerical results to the experimental onesshows good agreement, if a systematic shift of the critical phase-space density is included. We relatethis shift to the absence of thermal equilibrium between the condensate and the thermal wings,which were used in the experiment to extract the temperature. We expand on this observation bystudying the full relaxation dynamics between the condensate and the thermal cloud. I. INTRODUCTION
Frictionless flow is one of the defining features of su-perfluidity [1]. For a moving obstacle with velocity v in asuperfluid, the frictionless nature of the superfluid nearthe obstacle breaks down when v exceeds a certain criticalvelocity v c . According to Landau’s criterion this criticalvelocity is estimated as v c = min k [ (cid:15) ( k ) / (cid:126) k ], where (cid:15) ( k )is the excitation spectrum, (cid:126) is the Planck constant, and k is the wave vector, with k = | k | , see Refs. [1–3]. Anobject moving with a velocity above v c dissipates energyvia the creation of elementary excitations, for example,vortices or phonons. Superfluidity was first observed inliquid helium 4 and helium 3. Since then, superfluidityhas been studied in quantum gas systems of bosons [4–8],fermions [9–11], as well as of Bose-Fermi mixtures [12].The phenomenon of superfluidity is closely relatedto the Bose-Einstein condensation (BEC) of interactinggases. Interestingly, a uniform two-dimensional (2D) sys-tem cannot undergo the BEC transition because the for-mation of long-range order is precluded by thermal fluc-tuations [13, 14]. However, it forms a superfluid withquasi-long range order via the Berenzinskii-Kosterlitz-Thouless (BKT) mechanism [15]. The quasi-long rangeorder of this state refers to the algebraic decay of thesingle-particle correlation function. The algebraic expo-nent of this correlation function increases smoothly withtemperature. At the critical temperature, the superfluiddensity of the system undergoes a universal jump of 4 /λ ,where λ is the de Broglie wavelength. Experiments on2D bosonic systems, such as a liquid helium film [16],and trapped Bose gases [17–21] have shown indicationsof the BKT transition. Furthermore, a trapped 2D sys-tem can form a BEC due to the modified density of states[22, 23] and leads to an interesting interplay of the twophase transitions [24]. Quasi-long range order in 2D bosonic systems can bedetected via interference and time-of-flight techniques[17–21, 25–27]. However, as a direct method, superflu-idity of ultracold atomic gases was probed using a localperturbation, in particular via laser stirring. For exam-ple, superfluidity of 3D BECs was tested via laser stirringin Refs. [4, 11]. In the experiment [26], thermal relax-ation of a perturbed 2D quasicondensate is studied.Ref. [8] reported on stirring a trapped 2D Bose gasof Rb atoms with a blue-detuned laser, moving on acircular path around the trap center. The circular mo-tion ensures that the harmonically trapped 2D gas isprobed at a fixed phase-space density. By choosing differ-ent radii of the circular motion, the superfluid transitionwas explored. In this paper, we provide a quantitativeunderstanding of the experiment using a c-field simula-tion method. We demonstrate that a blue-detuned laserof intensity comparable to the mean-field energy causesdissipation due to the creation of vortex-antivortex pairs.This is in contrast to laser stirring with a red-detunedlaser [11], where dissipation occurs via phonons [28]. Fur-thermore, we study the relaxation dynamics of the stirredgas following the stirring process, which shows a slow en-ergy transport between the condensate and the thermalcloud. We identify the origin of this slow relaxation to bevortex recombination and diffusion. We show that thiseffect can explain quantitatively the shift of the measuredcritical phase-space density in the experiment.This paper is organized as follows. In Sec. II we de-scribe the simulation method that we use. In Sec. III wedetermine the critical velocity v c of the stirred gas, basedon which we identify the superfluid to thermal transition.In Sec. IV we discuss the dissipation via vortex pairs. InSec. V we compare the simulation results with the ex-periment. In Sec. VI we analyze the relaxation of thestirred gas, and in Sec. VII we conclude. a r X i v : . [ c ond - m a t . qu a n t - g a s ] M a r II. SIMULATION METHOD
We simulate the stirring dynamics of a weakly inter-acting 2D bosonic system using the c-field simulationmethod that we used for a 3D system in Ref. [28]. Wedescribe this method in the following. We start out withthe Hamiltonian of the unperturbed system,ˆ H = (cid:90) d r (cid:104) − ˆ ψ † ( r ) (cid:126) m ∇ ˆ ψ ( r ) + V ( r ) ˆ ψ † ( r ) ˆ ψ ( r )+ g ψ † ( r ) ˆ ψ † ( r ) ˆ ψ ( r ) ˆ ψ ( r ) (cid:105) . (1)ˆ ψ and ˆ ψ † are the bosonic annihilation and creation oper-ator, respectively. The 2D coupling parameter g is givenby g = ˜ g (cid:126) /m , where ˜ g = √ πa s /l z is the dimensionlessinteraction, m is the atomic mass, a s is the 3D s -wavescattering length, and l z = (cid:112) (cid:126) / ( mω z ) is the harmonicoscillator length of the confining potential mω z z / z direction. ω z is the trap frequency along the z di-rection. V ( r ) describes the external potential, which isa harmonic trap, V h ( r ) = mω r r / ω r is the trap fre-quency in the radial direction and r = ( x + y ) / is theradial coordinate. We introduce a time-dependent termto describe laser stirring,ˆ H s ( t ) = (cid:90) d r V ( r , t )ˆ n ( r ) , (2)where V ( r , t ) is the time-dependent stirring potential andˆ n ( r ) is the density operator at the location r = ( x, y ).The stirring potential is a Gaussian with a width σ anda strength V , V ( r , t ) = V exp (cid:16) − (cid:0) r − r s ( t ) (cid:1) σ (cid:17) , (3)which is centered at r s ( t ) = (cid:0) x s ( t ) , y s ( t ) (cid:1) . We move( x s , y s ) along a circular path as a function of time t .We perform numerical simulations by mapping thissystem on a lattice system, which also introduces a short-range cutoff; see Appendix A. This short-range cutoff isof the order of the healing length ξ = (cid:126) / √ mgn , with n being the density. We describe both the equations ofmotion and the initial state within a c-number repre-sentation, which corresponds to formally replacing theoperators ˆ ψ by complex numbers ψ . Furthermore, weapproximate the initial ensemble by a classical ensemble,within a grand-canonical ensemble of temperature T andchemical potential µ . We sample the initial states via aclassical Metropolis algorithm.The simulation setup consists of a disc-shaped 2D cir-cular condensate of Rb atoms. This choice of the2D circular condensate is inspired by the experimen-tal setup of Ref. [8]. In the simulations we consider N = 38 , − , Rb atoms confined by the har-monic potential in both the radial and transverse direc-tions. The trap frequencies are ω r = 2 π ×
25 Hz and ω z = 2 π × . a s = 5 . g = √ πa s /l z = 0 . T = 63 −
85 nK.The simulation parameters that we use, are in the typi-cal range of the experimental parameters of Ref. [8]. Forsimulations of a quasi- and a pure-2D trap geometry weuse a lattice of 180 × × ×
200 sites, withthe lattice discretization length l = 0 . µ m, respectively.We choose l such that it is smaller than, or comparableto, the healing length ξ and the de Broglie wavelength λ = (cid:112) π (cid:126) /mk B T , see Ref. [29]. The trapped gas is inthe pure-2D regime if k B T, µ < (cid:126) ω z . When k B T and µ are comparable to (cid:126) ω z , it is in the quasi-2D regime.After initializing the trapped system at temperature T ,we switch on the stirring potential described by Eq. 3. Inthe experiment [8] the trapped gas is stirred with a blue-detuned laser beam moving on a circular path around thetrap center. For the circular motion of stirring we choose( x s , y s ) = R (cid:0) cos( ω m t ) , sin( ω m t ) (cid:1) , where R and ω m arethe stirring radius and frequency, respectively. For thestirring potential we use the strength V = k B ×
80 nK andthe width σ = 1 µ m, in accordance with the experiment.The stirring sequence is the following: We linearly switchon the stirring potential over 5 ms, let it stir the systemfor 200 ms, and then switch it off over 5 ms. This is againinspired by the experimental choices. We repeat this forvarious stirring velocities v = Rω m by changing both R and ω m . By choosing different R we stir the differentregimes of the trapped gas, the superfluid, the thermal,and the crossover regime.After stirring we calculate the total energy E = (cid:104) H (cid:105) using the unperturbed Hamiltonian in Eq. 1, where weuse ψ instead of ˆ ψ . From this energy we determine theequilibrium temperature T eq of the stirred gas. We inferthis temperature by numerically inverting the tempera-ture dependence of the equilibrium state, E eq = E eq ( T ).We elaborate on this in Appendix A. From the tem-perature difference between the stirred and initial sys-tem, the heating ∆ T eq = T eq − T is determined. Wealso calculate the local energy, as well as the vortexand anti-vortex distribution. We define the local en-ergy as E i = − J (cid:80) j ( ψ ∗ i ψ j + ψ ∗ j ψ i ) + U n i + V i n i , where j refers to the nearest neighbor sites. ψ i , n i = | ψ i | ,and V i are the complex-valued field, the density, and thetrap potential at site i , respectively. J and U are theBose-Hubbard parameters, see Appendix A. For the vor-tex distribution, we calculate the phase winding aroundthe lattice plaquette of size l × l , using (cid:80) (cid:3) δφ ( x, y ) = δ x φ ( x, y ) + δ y φ ( x + l, y ) + δ x φ ( x + l, y + l ) + δ y φ ( x, y + l ),where the phase differences between sites is taken to be δ x/y φ ( x, y ) ∈ ( − π, π ]. φ is the phase field of ψ . We iden-tify a vortex and an antivortex by a phase winding of2 π and − π , respectively. By counting all vortices andantivortices we determine the total number of vortices.We restrict this counting to the the superfluid region ofthe gas as we describe below. R r
BKT n ( a t o m s / µ m ) r ( µ m) n BKT (a)050100
R r
BKT n ( a t o m s / µ m ) r ( µ m) n BKT (a) 048 v B ∆ T e q ( n K ) v (mm / s)(b) v c v B ∆ T e q ( n K ) v (mm / s)(b) v c r BKT
10 15 20 v c ( mm / s ) R ( µ m) Superfluid Thermal (c)
FIG. 1.
Determining the superfluid response and critical velocity.
In panel (a) we show the simulated density profile(dots) of the trapped 2D gas. We stir the gas with a repulsive Gaussian potential (green Gaussian beam of width σ = 1 µ m) ona circular path around the trap center at a stirring radius R . The phase-space density D of the gas decreases with increasing R .The prediction for the BKT transition in a uniform gas [31] yields the critical phase-space density D c = 8 .
3, which correspondsto the cloud density n BKT = 19 . /µ m at T = 85 nK. Using local density approximation, the border of the superfluid region isexpected to be at r BKT = 16 . µ m. In panel (b) we show the simulated heating ∆ T eq (circles) as a function of stirring velocity v = Rω m , at R = 14 . µ m. ∆ T eq is determined from the equilibrium temperature T eq of the stirred gas. We determine a criticalvelocity v c of 0 . ± . / s using the fitting function in Eq. 4. We show v c and the fitted curve by the vertical dashed andthe green continuous line, respectively. The Bogoliubov estimate of the phonon velocity at R is v B = 1 . / s. In panel (c)we show v c (circles) determined at various R . The y errorbars show the standard deviation. The x errorbars denote the size ofthe stirring potential (1 / √ e of diameter 2 σ ). III. SUPERFLUID RESPONSE
To study the superfluid behavior we stir a 2D quasicon-densate with a repulsive Gaussian potential. We preparea trapped 2D quasicondensate of N = 93 , Rb atomsat temperature T = 85 nK. We show the simulated den-sity profile of the trapped gas in Fig. 1(a). We stir thegas with a circularly moving, repulsive stirring potentialat stirring radius R = 14 . µ m. As mentioned in Sec.II, we use the strength V = k B ×
80 nK and the width σ = 1 µ m for the stirring potential. This strength V iswell above the local mean-field energy µ loc ≈ k B ×
21 nKat the stirring location. After stirring we determine theinduced heating ∆ T eq = T eq − T from the equilibriumtemperature T eq of the stirred gas, see Sec. II for details.By varying the stirring frequency ω m we determine ∆ T eq as a function of stirring velocity v = Rω m . We show∆ T eq determined for various v in Fig. 1(b). The inducedheating is almost negligible at low v , its onset occurs ata velocity v c , and for v > v c it increases rapidly. Wequantify the onset of heating using a fitting function,(∆ T ) fit = A · max[0 , ( v − v c )] + B, (4)which is discussed in Ref. [30], with the free param-eters A , B , and v c . For the simulated heating shownin Fig. 1(b), this function gives a critical velocity of v c = 0 . ± . / s. We compare this critical veloc-ity to the Bogoliubov estimate of the phonon velocity v B = 1 . / s at the stirrer location. The Bogoliubovvelocity is determined by v B = (cid:126) √ ˜ gn/m . The observedcritical velocity is v c ≈ . v B . This is notably differentfrom the case of an attractive stirring potential, where v c ≈ v B [28]. We explain this reduction of v c for a repul-sive stirring potential in Sec. IV. By choosing different radii R we explore the variousregimes of the trapped gas. We use the same strength V and the same width σ as above. For each R , we first de-termine the induced heating ∆ T eq as a function of v , andthen by using the fitting function given in Eq. 4 we de-termine v c . We show v c determined at various R in Fig.1(c). The stirring radii are in the range R = 10 − µ m.For R = 10 − µ m, there is no significant change of v c .As R reaches the crossover regime, v c is reduced sharplyand for R above the crossover regime, v c is zero. Accord-ing to the BKT prediction in a uniform system [31] with˜ g = 0 .
093 combined with local density approximation,the crossover regime should occur at r BKT = 16 . µ m.This prediction is in good agreement with the crossoverregime identified by the simulated v c . Thus, we clearlyidentify the superfluid, the crossover, and the thermalregimes by the finite, the sharply decreasing, and thezero critical velocities v c , respectively. We note that inthe crossover region the decrease of v c is as sharp as thesize of the stirrer allows. Furthermore, we note that theobserved almost constant v c for R < r
BKT can be due tothe accelerated circular motion and the large strength ofthe stirring potential [28].
IV. DISSIPATION MECHANISM
The observed critical velocities are in the range v c =0 . − . v B . To understand what leads to this reductionof the critical velocity with regard to the phonon velocity,we investigate the time evolution of the phase field φ ( r , t )of a single realization of the thermal ensemble. We ob-tain this phase field from the complex field ψ ( r ) via thephase-density representation ψ = √ n exp( iφ ). In Fig. 2we show the phase evolution of the trapped 2D quasi- FIG. 2.
Dissipation due to vortex-antivortex pairs.
The phase evolution of a single realization of the trapped gas stirredwith a repulsive Gaussian potential along a circular path at the times t = 5 , , , ,
150 ms. The stirring potential (greendisc of diameter 2 σ ) moves counterclockwise in a circle at R = 12 µ m with a velocity v > v c . This creates strong phase gradients,which result in the creation of vortex-antivortex pairs. We identify a vortex (red circle) and an antivortex (blue triangle) bya phase winding of 2 π and − π around the lattice plaquette, respectively. The white line indicates the superfluid-thermalboundary circle, based on r BKT . In the thermal region (white shaded region outside of r BKT ) the phase fluctuates strongly.The field of view in each figure is of size 35 µ m × µ m. condensate stirred at R = 12 µ m. We use the velocity v ≈ . v B , which is above the steep onset of dissipationrelated to the breakdown of superfluidity. The phase evo-lution of the unperturbed gas shows rather weak phasegradients. As stirring is switched on, the phase fieldaround the stirring potential starts to fluctuate. Thesefluctuations develop into strong phase gradients, whichresult in the creation of vortex-antivortex pairs. This canbe confirmed by calculating the phase winding aroundeach plaquette of our numerical grid, as described in Sec.II. We show the calculated phase winding in Fig. 2,where vortices and antivortices are shown as circles andtriangles, respectively. This indeed shows the creationof vortex-antivortex pairs during stirring. We recall thatthe stirring strength V = k B ×
80 nK is much larger thanthe mean-field energy µ loc ≈ k B ×
30 nK at the stirring lo-cation, which results in a strong reduction of the densityat the stirrer location. This density reduction serves as anucleation site for the creation of vortex pairs. We notethat this mechanism of vortex-pair-induced dissipation issuppressed for an attractive stirring potential, as shownin Ref. [28]. This scenario of dissipation induced by vor- tex pair creation is consistent with a recent experiment[32].
V. COMPARISON TO EXPERIMENT
We now compare the results of our simulation withthe experiment [8]. We first show the comparison be-tween the experiment and simulation for the heating as afunction of v . In the superfluid regime, we stir the qua-sicondensate at the radius R = 14 . µ m. The simulateddensity profile is shown in the inset of Fig. 3(a). Afterstirring we let the stirred gas relax for 100 ms of relax-ation time and then determine the induced heating fromthe temperature of the wings of the cloud. We fit thesewings to the Hartree-Fock prediction, n ( r ) = − mk B T π (cid:126) ln (cid:104) − exp (cid:16) µ − V h ( r ) − gn ( r ) k B T (cid:17)(cid:105) , (5)with the fitting parameters T and µ . This method isadopted according to experiment, in which the temper-ature of the stirred gas is determined in the same way, R n r ( µ m)050 R n r ( µ m)040 R n r ( µ m)040 R n r ( µ m)048 ∆ T w / e q ( n K ) (a)048 ∆ T w / e q ( n K ) (a)0612 0 0.5 1 1.5 ∆ T w / e q ( n K ) v (mm / s) (b) -0.1 0 0.1 0.2 0.3 0.4 0.500.51 v c ( mm / s ) µ loc /k B T -0.1 0 0.1 0.2 0.3 0.4 0.500.51 v c ( mm / s ) µ loc /k B TT w , exp T w , sim T eq , sim (c)(c) v c, w , exp v c, w , sim v c, eq , sim FIG. 3.
Comparison of simulation and experimental results.
In panel (a) we stir the quasicondensate region at R = 14 . µ m (inset), whereas in panel (b) we stir the thermal region at R = 16 . µ m (inset). We compare the measured heating∆ T w (solid circles) with the simulated ∆ T w (diamonds) and ∆ T eq (open circles) at varying stirring velocity v . ∆ T w and ∆ T eq are determined from the wing temperature T w and equilibrium temperature T eq , respectively. The y errorbars are the standarddeviation. The x errorbars indicate the spread of velocities along the size of the stirring potential. In panel (c) we comparethe measured v c, w (solid circles) with the simulated v c, w (diamonds) and v c, eq (open circles) across the superfluid-thermaltransition. According to the BKT prediction in a uniform gas [31] this transition occurs at ( µ loc /k B T ) c ≈ .
13. The thermalstate by this prediction is indicated by the gray shaded area. The x errorbars denote the region of µ loc /k B T , which is probedby the stirring potential due to its size. The y errorbars show the standard deviation. The experimental data is from Ref. [8]. following a relaxation of 100 ms as well. We denote thisheating determined from the wing temperature T w by∆ T w = T w − T , with T being the initial temperature.We show the simulated ∆ T w and their comparison withthe experimental ones for various v in Fig. 3(a). Themeasured and simulated heating are found to be in goodagreement if we base the comparison on ∆ T w . We alsocompare the measured ∆ T w with the simulated ∆ T eq de-termined from the equilibrium temperature T eq of thestirred gas. We show ∆ T eq as open circles in Fig. 3(a).They show agreement at low and intermediate velocities v , whereas they differ at large v . This noticeable differ-ence at large v is due to the absence of global thermalequilibrium of the stirred gas. As explained in Sec. VI,the stirred gas relaxes by transporting the excess energybetween the superfluid in the central part and the ther-mal cloud in the periphery, which is a slow process. Theabsence of global thermal equilibrium leads to a smallerwing temperature than the equilibrium temperature.The results shown in Fig. 3(a) indicate that the onsetof heating occurs at a velocity v c , and for v > v c heat-ing increases rapidly. Both in experiment and simulation v c is determined using the fitting function in Eq. 4. InFig. 3(b) we show the comparison between the experi-ment and simulation for stirring the thermal region of thetrapped gas of N = 38 ,
162 atoms. The simulated den-sity profile of the gas is given in the inset of Fig. 3(b).Both the measured and simulated heating ∆ T w are in good agreement. The simulated ∆ T eq determined fromthe equilibrium temperature of the system are below themeasured ∆ T w at large v . As we will explain in Sec.VI, this is again due to the absence of global thermalequilibrium. As the stirred thermal cloud has more ex-cess energy than the condensate, the wing temperatureis larger than the equilibrium temperature. The resultsshown in Fig. 3(b) indicate that heating occurs at all v ,which results in a zero v c .Next, we show in Fig. 3(c) the comparison between theexperiment and simulation for v c that are determined bystirring the superfluid, the crossover, and the thermalregime. In the experiment [8] v c is measured for differentconfigurations of the total number of atoms N , the tem-peratures T , and the stirring radii R . We compare themeasured v c with the simulated v c determined by stirringthe 2D gas in Sec. III. We show both the measured andsimulated v c as a function of the dimensionless parameter µ loc ( r ) /k B T . The parameter µ loc /k B T characterizes thedegree of degeneracy of the cloud and is the relavant pa-rameter in the sense that the thermodynamic propertiesof the gas depend only on the ratio µ/k B T [31, 33, 34].We refer to v c determined from the wing temperature T w and from the equilibrium temperature T eq as v c, w and v c, eq , respectively. Both the measured and simulated v c, w show good agreement. The measured v c, w and the sim-ulated v c, eq agree in the superfluid and thermal regime,while they differ in the crossover regime. For the mea- ∆ ˜ E i ( n K ) xy(a) pure-2Dxy(a) pure-2Dxy(a) pure-2D 0 . . . . ∆ ˜ E i ( n K ) xy(b) quasi-2Dxy(b) quasi-2Dxy(b) quasi-2D 0 . . . . FIG. 4.
Relaxation dynamics.
Energy flow of the stirred trapped gas for the relaxation times t relax = 0 . , . , . , . R = 12 µ m (circle indicated by the green dashed line). We show theevolution of the excess energy ∆ ˜ E i = (cid:0) E i ( t ) − E eq i (cid:1) /n max for the pure-2D gas in the upper panels (a) and for the quasi-2D gasin the lower panels (b). E i is the local energy of the stirred gas at time t and E eq i is its final equilibrium local energy. n max isthe maximum density of the system. The top panels of (a) and (b) are the one-dimensional cut through the trap center. Thesuperfluid-thermal boundary circle is indicated by the white line. ∆ ˜ E i in the far wings of the cloud is negative as E eq i is largerthan E i . The field of view is 90 µ m × µ m. sured and simulated v c, w , the crossover regime occurs at µ loc /k B T ≈ .
24 and 0 .
22, respectively. However, forthe simulated v c, eq , it occurs at µ loc /k B T ≈ .
11. Thetheoretical prediction for the BKT transition in a uni-form gas [31] with ˜ g = 0 .
093 occurs at ( µ/k B T ) c ≈ . v c, eq , whereas its comparison withthe crossover regimes identified by the measured and sim-ulated v c, w shows a shift. This shift was observed in Ref.[8], but could not be explained. We conclude that theexperiments of Ref. [8] can be reproduced quantitativelyif the wing temperature is used, rather than T eq . Thissuggests that the system has not relaxed to thermal equi-librium after the waiting time of 100 ms. We confirm andelaborate on this point and the underlying mechanism inthe following sections. VI. RELAXATION DYNAMICS
We now investigate the relaxation of the system, fol-lowing the stirring process in the superfluid regime. Thisincludes a discussion of the influence of the confinementof the system in the z direction. For strong confinement,the system approaches a purely 2D limit, while it is quasi-2D for intermediate confinement. A. Energy-flow dynamics
We first analyze the energy-flow dynamics of a stirredtrapped gas in the purely 2D limit, and then compare thisdynamics with a quasi-2D gas. For a pure-2D trappedgas, we consider a gas of N = 64 , Rb atoms, whichis strongly confined in the transverse direction by the har-monic potential. The temperature T = 85 nK is smallerthan the transverse trap energy (cid:126) ω z /k B = 144 nK, sothat the gas is in the ground state in this direction. As . . . . . . . . . n i /n max . n i /n max . n i /n max . . . . . . . . . . n i /n max . n i /n max . n i /n max . FIG. 5.
Relaxation of the density and vortices.
We show the evolution of the density and vortices of a single realization ofthe stirred gas for t relax = 0 . , . , . , . z = 0) xy plane.The maximum density n max is 71 /µ m for pure-2D gas and the maximum column density n max is 95 /µ m for quasi-2D gas.We show vortices (red circles) and antivortices (blue triangles) in the superfluid region at the detection radius R det = 14 . µ m(circle indicated by the white line), which is below r BKT [35]. The field of view is 40 µ m × µ m. the width of the condensate in the z direction is smallerthan the lattice discretization length l , we simulate thissystem using a single xy -layer of lattice only, see Sec.II. We stir the gas at R = 12 µ m for 200 ms at a velocity v > v c . After that we switch off the stirring potential andlet the gas relax. We calculate the local energy E i of thestirred gas and its final equilibrated local energy E eq i , asdescribed in Sec. II. We show the evolution of the excessenergy ∆ ˜ E i = (cid:0) E i ( t ) − E eq i (cid:1) /n max for various relaxationtimes t relax in Fig. 4(a). n max is the maximum densityof the system. The evolution of ∆ ˜ E i after t relax = 0 . . xy -layers of lattice in the z direction.We stir the gas using the same stirring parameters as forthe pure-2D case. We show the evolution of the excessenergy ∆ ˜ E i of the stirred gas for various t relax in Fig.4(b). In this case, E i , E eq i and n max are the column (i.e.integrated along the z axis) quantities. The evolutionof ∆ ˜ E i after t relax = 0 . . . B. Vortex dynamics
To understand what causes this slow relaxation for thesystem, we now examine the evolution of the density andvortices of the system. We calculate the local density as n i = | ψ i | and vortices as described in Sec. II. We showthe density and vortex evolution of a single realizationof the stirred pure- and quasi-2D gas after various t relax in Fig. 5. For both systems, the density relaxation ishard to recognize, whereas the vortex evolution clearlyexhibits decay of vortices. Thus, the system relaxes viadecay of the induced vortices. Vortices can decay viaboth annihilation of a vortex with an antivortex, anddrifting out to the thermal region of the cloud. For the t on t off N v time (s)pure-2Dquasi-2D bi-exponential fit FIG. 6.
Vortex relaxation.
We show the vortex number N v as a function of time t for the pure- and quasi-2D gas. Wefit the vortex decay with the bi-exponential function in Eq. 6and determine the decay times τ d , τ d (see text). The fittedcurves are shown as the dashed lines. pure-2D gas the number of vortices after t relax = 0 . R det = 14 . µ m, and average it over200 realizations. We show the averaged vortex number N v as a function of time t for both systems in Fig. 6. Asstirring is switched on at time t on , N v starts to increaseapproximately linearly. It reaches its maximum N max at t off . After the stirring is switched off, it decays ap-proximately exponentially. For both systems the natureof vortex growth and decay are the same, but the rateswith which they grow and decay are different. For thepure-2D gas the growth and decay rates are larger andsmaller than those for the quasi-2D gas, respectively. Theenhanced growth and the suppressed decay rate for thepure-2D gas can be due to the suppression of vortex an-nihilation, as mentioned above, and a slow vortex drift.We quantify the vortex decay rate using the function, f ( t ) = N e − t/τ d + N e − t/τ d + N , (6)with the free parameters N , N , N , τ d , τ d . From thefit, we determine τ d , τ d ≈ . , .
87 s and 0 . , .
87 sfor pure- and quasi-2D gas, respectively. These decaytimes are similar to those determined from the mean ex-cess energy in Appendix B 1. The fast decay τ d andthe slow decay τ d are essentially connected to the vor-tex annihilation and drift lifetime, respectively. For thepure-2D gas τ d and τ d are larger than and equal tothose for the quasi-2D gas, respectively.We show in Table I N max and the extracted τ d , τ d at varying R det , for both systems. τ d and τ d increase TABLE I. N max and the extracted τ d , τ d for different R det .Pure-2D Quasi-2D R det ( µ m) N max τ d (s) τ d (s) N max τ d (s) τ d (s)14.0 172 0.129 0.809 127 0.050 0.80614.5 196 0.135 0.866 151 0.055 0.86915.0 224 0.142 0.920 182 0.059 0.93315.5 251 0.147 0.959 214 0.066 0.989 weakly as R det is increased. However, the following con-clusions are essentially independent of the choice of R det .Overall, N max and τ d are larger for pure-2D gas thanthose for quasi-2D gas, respectively, while τ d are similarfor both systems. We compare the simulated τ d , τ d ofquasi-2D gas to the waiting time of 0 . τ d , whereas it is smaller than the slow decay τ d . Thissuggests that most vortex recombination processes haveoccurred at the time of the measurement. However, thevortex drift to the thermal cloud has not occurred, andthe system is in a metastable state, not in the equili-brated state. This is the mechanism that is responsiblefor the difference between the wing temperature and theequilibrium temperature.We note that in Ref. [39] an estimate for the timeof a vortex line drifting to the thermal cloud was given.While this estimate was for a three dimensional system,we find that the analytical estimate of Ref. [39] givesa timescale that is consistent with our simulation. Wealso note that the vortex lifetime is suppressed at hightempearures [40–42]. VII. CONCLUSIONS
We have studied the superfluid to thermal transition ofa trapped 2D Bose gas of Rb atoms by stirring it witha repulsive stirring potential on a circular path aroundthe trap center. The superfluid transition was probedby choosing different radii of the circular motion. Wehave identified the superfluid, the crossover, and the ther-mal regime by the finite, the sharply decreasing, and thezero critical velocity, respectively. The superfluid regionof the gas yields critical velocities that are in the range v c = 0 . − . v B , where v B is the phonon velocity. Wehave demonstrated that the onset of dissipation is dueto the creation of vortex-antivortex pairs. The compar-ison of the simulation with the experiment shows goodagreement if the temperature measurement of the exper-iment is imitated in the simulation, i.e. by extracting thewing temperature. However, we confirm the systematicshift that was observed in experiment, if thermal equi-librium is assumed. We have demonstrated that the ab-sence of thermal equilibrium after the waiting time thatwas used in experiment is due to a remarkably slow re-laxation mechanism: The energy transport across the su-perfluid to thermal interface occurs only on timescales ofseconds. This slow transport mechanism is due to theslow drift of vortices out of the superfluid into the ther-mal wings of the system. We emphasize that this mech-anism is relevant for many on-going experiments in thefield of ultracold atoms, and their temperature measure-ments. Furthermore, this effect of suppressed transportacross critical interfaces is in itself intriguing, and couldbe studied in a future cold atom experiment with clarity. ACKNOWLEDGEMENTS
We acknowledge support from the DeutscheForschungsgemeinschaft through Grants No. MA5900/1-1 and No. SFB 925, the Hamburg Centre for Ul-trafast Imaging, and from the LandesexzellenzinitiativeHamburg, supported by the Joachim Herz Stiftung. JDthanks J´erˆome Beugnon for many fruitful discussionsand acknowledges support from ERC (Synergy grantUQUAM).
Appendix A: Simulated heating
In this section we show how we determine the equi-librium temperature T eq of a stirred gas using the c-fieldmethod described in Sec. II. We discretize the continuumHamiltonian in Eq. 1 by the Bose-Hubbard Hamiltonian[43] on a 2D square lattice, H = − J (cid:88) (cid:104) ij (cid:105) ( ψ ∗ i ψ j + ψ ∗ j ψ i ) + U (cid:88) i n i + (cid:88) i V i n i . (A1) ψ i and n i = | ψ i | are the complex-valued field and thedensity at site i , respectively. (cid:104) ij (cid:105) indicates nearest-neighbor bonds. For a lattice discretization length l ,the Bose-Hubbard parameters are related to the contin-uum parameters via J = (cid:126) / (2 ml ) and U = gl − . The2D coupling parameter g is given by g = ˜ g (cid:126) /m , where˜ g = √ πa s /l z is the dimensionless interaction, m is theatomic mass, a s is the 3D s-wave scattering length, and l z = (cid:112) (cid:126) / ( mω z ) is the harmonic oscillator length of theconfining potential mω z z / z direction. ω z is thetrap frequency in the z direction. The harmonic trappingpotential is V i = mω r r / ω r is the trap frequency inthe radial direction and r = ( x + y ) / is the radialcoordinate.We first initialize the system in a thermal state at tem-perature T via classical Monte Carlo, and then calculateits energy E = (cid:104) H (cid:105) using the Hamiltonian in Eq. A1.By varying the temperature of the system T while keep-ing the total number of atoms N fixed, we calculate theenergy E as a function of T . In Fig. 7 we show the tem-perature dependence of the energy E for the pure- andquasi-2D gas that are described in Sec. VI.To determine the heating, we first stir the gas withthe repulsive stirring potential as described in Sec. II E ( n K ) T (nK)quasi-2Dpure-2D FIG. 7. Temperature dependence of the energy per atom forthe pure- and quasi-2D gas. and then after stirring calculate its energy E using Eq.A1. We numerically invert this energy E to the equilib-rium temperature T eq using the temperature dependenceshown in Fig. 7. Finally, from the temperature differ-ence between the stirred and initial system, the heating∆ T eq = T eq − T is determined. Appendix B: Relaxation dynamics
In this section we elaborate on the relaxation dynamicsof the stirred trapped gas that we discuss in Sec. VI. Wefirst discuss the energy flow dynamics in Sec. B 1 andthen the vortex dynamics in Sec. B 2.
1. Energy-flow dynamics
Here we elaborate on the energy flow dynamics for thepure- and quasi-2D trapped system. We stir both sys-tems with the stirring potential at velocity v ≈ . v B .After stirring we calculate the excess energy ∆ ˜ E i = (cid:0) E i ( t ) − E eq i (cid:1) /n max as described in Sec. VI and by av-eraging this energy over the superfluid region of the gas,we calculate the mean energy ∆ E mean . We show the evo-lution of ∆ E mean for both systems in Fig. 8(a). ∆ E mean decays approximately exponentially as the excess energy∆ ˜ E i within the superfluid region outflows to the ther-mal cloud. We quantify the energy decay time using thefitting function in Eq. 6. From the fit, we determinethe decay times τ d , τ d ≈ . , .
98 s and 0 . , .
02 sfor pure- and quasi-2D gas, respectively. These decaytimes are similar to the ones that we determine from thevortex decay in Sec. VI B. In addition to the stirring at v ≈ . v B , we show ∆ E mean corresponding to stirring at v ≈ . v B for the pure- and quasi-2D gas in Fig. 8(b)and (c), respectively.0 ∆ E m e a n ( n K ) (a) v ≈ . v B ∆ E m e a n / E m a x (b) pure-2D00.51 t off ∆ E m e a n / E m a x time (s)(c) quasi-2D pure-2Dquasi-2Dbi-exponential fit v ≈ . v B , E max = 12 . v ≈ . v B , E max = 5 . v ≈ . v B , E max = 11 .
03 nK v ≈ . v B , E max = 5 .
03 nK
FIG. 8. In panel (a) we show the evolution of the mean ex-cess energy ∆ E mean for the pure- and quasi-2D gas stirredat velocity v ≈ . v B . We fit the energy decay with the bi-exponential function in Eq. 6 to determine the decay times τ d , τ d (see text). The fitted curves are shown by the dashedlines. We show ∆ E mean normalized by its maximum meanenergy E max corresponding to v ≈ . v B and 0 . v B for thepure- and quasi-2D gas in panels (b) and (c), respectively. N v / N m a x (a) pure-2D00.51 t on t off N v / N m a x time (s) (b) quasi-2D R = 12 µ m , V = k B ×
80 nK R = 12 µ m , V = k B ×
56 nK R = 10 µ m , V = k B ×
80 nK R = 10 µ m , V = k B ×
56 nK R = 12 µ m , V = k B ×
80 nK R = 12 µ m , V = k B ×
56 nK R = 10 µ m , V = k B ×
80 nK R = 10 µ m , V = k B ×
56 nK
FIG. 9. We show the averaged vortex number N v normalizedby its maximum vortex number N max for the different stirringradii R and strengths V , for the pure- and quasi-2D gas inpanels (a) and (b), respectively.
2. Vortex relaxation
Next, we turn to the vortex relaxation of the stirredgas. We stir the pure- and quasi-2D gas using differentstirring radii R and stirrer strengths V . We calculate thetotal number of vortices within the superfluid region ofthe cloud, as described in Sec. VI B, and average it over128 realizations. We show the averaged vortex number N v normalized by its maximum vortex number N max asa function of time t for the pure- and quasi-2D gas inFigs. 9(a) and 9(b), respectively. For the pure-2D gasthe relaxation of vortices is slower than that for the quasi-2D gas, as shown in Sec. VI B. [1] A. J. Leggett, Quantum Liquids: Bose Condensation andCooper Pairing in Condensed-Matter Systems (OxfordUniv. Press, 2006).[2] L. Pitaevskii and S. Stringari,
Bose-Einstein Condensa-tion (Oxford Univ. Press, 2003).[3] C. J. Pethick and H. Smith,
Bose-Einstein Condensationin Dilute Gases (2nd ed. Cambridge University Press,2008).[4] C. Raman, M. K¨ohl, R. Onofrio, D. S. Durfee, C. E. Kuk-lewicz, Z. Hadzibabic, and W. Ketterle, Phys. Rev. Lett. , 2502 (1999); R. Onofrio, C. Raman, J. M. Vogels, J.R. Abo-Shaeer, A. P. Chikkatur, and W. Ketterle. Phys.Rev. Lett. , 2228 (2000).[5] K. W. Madison, F. Chevy, W. Wohlleben, and J. Dal-ibard, Phys. Rev. Lett. , 160405 (2007). [6] P. Engels and C. Atherton, Phys. Rev. Lett. , 160405(2007).[7] T. W. Neely, E. C. Samson, A. S. Bradley, M. J. Davis,and B. P. Anderson, Phys. Rev. Lett. , 160401 (2010).[8] R. Desbuquois, L. Chomaz, T. Yefsah, J. L´eonard, J.Beugnon, C. Weitenberg, and J. Dalibard, Nat. Phys. ,645 (2012).[9] M. W. Zwierlein, J. R. Abo-Shaeer, A. Schirotzek, C. H.Schunck, and W. Ketterle, Nature , 1047 (2005).[10] D. E. Miller, J. K. Chin, C. A. Stan, Y. Liu, W. Setiawan,C. Sanner, and W. Ketterle, Phys. Rev. Lett. , 070402(2007).[11] W. Weimer, K. Morgener, V. P. Singh, J. Siegl, K. Hueck,N. Luick, L. Mathey, and H. Moritz, Phys. Rev. Lett. , 095301 (2015). [12] I. Ferrier-Barbut, M. Delehaye, S. Laurent, A.T. Grier,M. Pierce, B. S. Rem, F. Chevy, and C. Salomon, Science , 1035 (2014); M. Delehaye, S. Laurent, I. Ferrier-Barbut, S. Jin, F. Chevy, and C. Salomon, Phys. Rev.Lett. , 265303 (2015).[13] N. D. Mermin and H. Wagner, Phys. Rev. Lett. , 1133(1966).[14] P. C. Hohenberg, Phys. Rev. , 383 (1967).[15] V. L. Berezinskii, Sov. Phys. JETP , 610 (1972); J. M.Kosterlitz and D. J. Thouless, J. Phys. C , 1181 (1973);P. Minnhagen, Rev. Mod. Phys. , 1001 (1987).[16] D. J. Bishop and J. D. Reppy, Phys. Rev. Lett. , 1727(1978).[17] Z. Hadzibabic, P. Kru¨ger, M. Cheneau, B. Battelier, andJ. Dalibard, Nature (London) , 1118 (2006).[18] P. Clad´e, C. Ryu, A. Ramanathan, K. Helmerson, andW. D. Phillips, Phys. Rev. Lett. , 170401 (2009).[19] S. Tung, G. Lamporesi, D. Lobser, L. Xia, and E. A.Cornell, Phys. Rev. Lett. , 230408 (2010).[20] T. Plisson, B. Allard, M. Holzmann, G. Salomon, A. As-pect, P. Bouyer, and T. Bourdel, Phys. Rev. A , 061606(2011).[21] J. Choi, S. W. Seo, and Y. Shin, Phys. Rev. Lett. ,175302 (2013).[22] V. Bagnato and D. Kleppner, Phys. Rev. A, , 5 (2004); H. Holzmann, M.Chevallier, and W. Krauth, Europhys. Lett. , 30001(2008).[23] Z. Hadzibabic and J. Dalibard, Rivista del Nuovo Ci-mento , 389 (2011).[24] Richard J. Fletcher, Martin Robert-de-Saint-Vincent,Jay Man, Nir Navon, Robert P. Smith, Konrad G. H.Viebahn, and Zoran Hadzibabic, Phys. Rev. Lett. ,255302 (2015).[25] A. Polkovnikov, E. Altman, and E. Demler, Proc. Natl.Acad. Sci. USA , 6125 (2006). [26] J. Choi, S. W. Seo, and Y. Shin, Phys. Rev. Lett. ,125301 (2012).[27] V. P. Singh and L. Mathey, Phys. Rev. A , 053612(2014).[28] V. P. Singh, W. Weimer, K. Morgener, J. Siegl, K. Hueck,N. Luick, H. Moritz, and L. Mathey, Phys. Rev. A ,023634 (2016).[29] C. Mora and Y. Castin, Phys. Rev. A , 053615 (2003).[30] G. E. Astrakharchik, and L. P. Pitaevskii, Phys. Rev. A , 013608 (2004).[31] N. Prokof’ev, O. Ruebenacker, and B. Svistunov, Phys.Rev. Lett. , 270402 (2001); N. Prokof’ev and B. Svis-tunov, Phys. Rev. A , 043608 (2002).[32] Woo Jin Kwon, Sang Won Seo, and Yong-il Shin, Phys.Rev. A , 033613 (2015).[33] T. Yefsah, R. Desbuquois, L. Chomaz, K. J. G¨unter, andJ. Dalibard, Phys. Rev. Lett. , 130401 (2011).[34] C. L. Hung, X. Zhang, N. Gemelke, and C. Chin, Nature , 236 (2011).[35] For the pure- and quasi-2D gas the border of the su-perfluid region with the final equilibrium temperature isexpected to be at r BKT = 15 . µ m and 16 . µ m, respec-tively.[36] S. J. Rooney, P. B. Blakie, B. P. Anderson, and A. S.Bradley, Phys. Rev. A , 023637 (2011).[37] P. C. Haljan, B. P. Anderson, I. Coddington, and E. A.Cornell, Phys. Rev. Lett. , 2922 (2001).[38] V. Bretin, P. Rosenbusch, F. Chevy, G.V. Shlyapnikov,and J. Dalibard, Phys. Rev. Lett. , 100403 (2003).[39] P. O. Fedichev and G. V. Shlyapnikov, Phys. Rev. A ,R1779(R) (1999).[40] J. R. Abo-Shaeer, C. Raman, and W. Ketterle, Phys.Rev. Lett. , 070409 (2002).[41] B. Jackson, N. P. Proukakis, C. F. Barenghi, and E.Zaremba, Phys. Rev. A , 053615 (2009).[42] G. Moon, W. J. Kwon, H. Lee, and Y. Shin, Phys. Rev.A , 051601(R) (2015).[43] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, andP. Zoller, Phys. Rev. Lett.81