Superfluidity and relaxation dynamics of a laser-stirred 2D Bose gas

TL;DR

This paper theoretically investigates the superfluid behavior and vortex dynamics of a 2D Bose gas under stirring, aligning well with experimental results and highlighting the role of vortex-antivortex pairs in the onset of heating.

Contribution

It provides a detailed theoretical analysis of superfluidity, critical velocity, and vortex creation in a 2D Bose gas, explaining experimental observations and the effects of thermal non-equilibrium.

Findings

Critical velocity marks the transition from superfluid to thermal regimes.

Vortex-antivortex pairs are responsible for the onset of heating.

Good agreement with experiment when accounting for thermal non-equilibrium.

Abstract

We investigate the superfluid behavior of a two-dimensional (2D) Bose gas of $^{87}$Rb atoms using classical field dynamics. In the experiment by R. Desbuquois \textit{et al.}, Nat. Phys. \textbf{8}, 645 (2012), a 2D quasicondensate in a trap is stirred by a blue-detuned laser beam along a circular path around the trap center. Here, we study this experiment from a theoretical perspective. The heating induced by stirring increases rapidly above a velocity $v_c$, which we define as the critical velocity. We identify the superfluid, the crossover, and the thermal regime by a finite, a sharply decreasing, and a vanishing critical velocity, respectively. We demonstrate that the onset of heating occurs due to the creation of vortex-antivortex pairs. A direct comparison of our numerical results to the experimental ones shows good agreement, if a systematic shift of the critical phase-space density is included. We relate this 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 by studying the full relaxation dynamics between the condensate and the thermal cloud.

Figures

Tables

Table 1. Table 1: N max subscript 𝑁 max N_{\mathrm{max}} and the extracted τ d 1 , τ d 2 subscript 𝜏 subscript 𝑑 1 subscript 𝜏 subscript 𝑑 2 \tau_{d_{1}},\,\tau_{d_{2}} for different R det subscript 𝑅 det R_{\mathrm{det}} .

Pure-2D			Quasi-2D
$R_{\det }\left(\mu m\right)$	$N_{\max }$	${\tau }_{d_1}\left(s\right)$	${\tau }_{d_2}\left(s\right)$	$N_{\max }$	${\tau }_{d_1}\left(s\right)$	${\tau }_{d_2}\left(s\right)$
14.0	172	0.129	0.809	127	0.050	0.806
14.5	196	0.135	0.866	151	0.055	0.869
15.0	224	0.142	0.920	182	0.059	0.933
15.5	251	0.147	0.959	214	0.066	0.989

Equations

$${\hat{H}}_0$$

$$\displaystyle\quad+\frac{g}{2}\hat{\psi}^{\dagger}({\bf r})\hat{\psi}^{\dagger}({\bf r})\hat{\psi}({\bf r})\hat{\psi}({\bf r})\Big{]}.$$

$${\hat{H}}_s(t)=\int \mathrm{d}\mathbf{r}V(\mathbf{r},t)\hat{n}(\mathbf{r}),$$

$$V({\bf r},t)=V_{0}\exp\Bigl{(}-\frac{\bigl{(}{\bf r}-{\bf r}_{s}(t)\bigr{)}^{2}}{2\sigma^{2}}\Bigr{)},$$

$$(\Delta T{)}_{\mathrm{fit}}=A\cdot \max [0,(v^2-v_c^2)]+B,$$

$$n({\bf r})=-\frac{mk_{\mathrm{B}}T}{2\pi\hbar^{2}}\ln\Bigl{[}1-\exp\Bigl{(}\frac{\mu_{0}-V_{h}({\bf r})-2gn({\bf r})}{k_{\mathrm{B}}T}\Bigr{)}\Bigr{]},$$

$$f(t)=N_1{e}^{-t/{\tau }_{d_1}}+N_2{e}^{-t/{\tau }_{d_2}}+N_0,$$

$$H_0=-J\sum _{⟨ij⟩}({\psi }_i^{\ast }{\psi }_j+{\psi }_j^{\ast }{\psi }_i)+\frac{U}{2}\sum _i{n}_i^2+\sum _i{V}_i{n}_i.$$

Full text

I Introduction

Frictionless flow is one of the defining features of superfluidity Leggett . For a moving obstacle with velocity $v$ in a superfluid, the frictionless nature of the superfluid near the obstacle breaks down when $v$ exceeds a certain critical velocity $v_{c}$. According to Landau’s criterion this critical velocity is estimated as $v_{c}=\mathrm{min}_{\bf k}[\epsilon({\bf k})/\hbar k]$, where $\epsilon({\bf k})$ is the excitation spectrum, $\hbar$ is the Planck constant, and ${\bf k}$ is the wave vector, with $k=|{\bf k}|$, see Refs. Leggett ; PitaevskiiStringari ; Pethick2008 . An object moving with a velocity above $v_{c}$ dissipates energy via the creation of elementary excitations, for example, vortices or phonons. Superfluidity was first observed in liquid helium $4$ and helium $3$. Since then, superfluidity has been studied in quantum gas systems of bosons Ketterle1999 ; Dalibard2000 ; Atherton2007 ; Anderson2010 ; Dalibard2012 , fermions Zwierlein2011 ; Miller2007 ; Weimer2015 , as well as of Bose-Fermi mixtures Salomon2014 .

The phenomenon of superfluidity is closely related to the Bose-Einstein condensation (BEC) of interacting gases. Interestingly, a uniform two-dimensional (2D) system cannot undergo the BEC transition because the formation of long-range order is precluded by thermal fluctuations Mermin1966 ; Hohenberg1967 . However, it forms a superfluid with quasi-long range order via the Berenzinskii-Kosterlitz-Thouless (BKT) mechanism Minnhagen1987 . The quasi-long range order of this state refers to the algebraic decay of the single-particle correlation function. The algebraic exponent of this correlation function increases smoothly with temperature. At the critical temperature, the superfluid density of the system undergoes a universal jump of $4/\lambda^{2}$, where $\lambda$ is the de Broglie wavelength. Experiments on 2D bosonic systems, such as a liquid helium film Bishop1978 , and trapped Bose gases Dalibard2006 ; Clade2009 ; Tung2010 ; Plisson2011 ; Shin2013 have shown indications of the BKT transition. Furthermore, a trapped 2D system can form a BEC due to the modified density of states Petrov2004 ; Hadzibabic2011 and leads to an interesting interplay of the two phase transitions Hadzibabic2015 .

Quasi-long range order in 2D bosonic systems can be detected via interference and time-of-flight techniques Polkovnikov2006 ; Dalibard2006 ; Clade2009 ; Tung2010 ; Plisson2011 ; Shin2013 ; Shin2012 ; Singh2014 . However, as a direct method, superfluidity of ultracold atomic gases was probed using a local perturbation, in particular via laser stirring. For example, superfluidity of 3D BECs was tested via laser stirring in Refs. Ketterle1999 ; Weimer2015 . In the experiment Shin2012 , thermal relaxation of a perturbed 2D quasicondensate is studied.

Ref. Dalibard2012 reported on stirring a trapped 2D Bose gas of 87Rb atoms with a blue-detuned laser, moving on a circular path around the trap center. The circular motion ensures that the harmonically trapped 2D gas is probed at a fixed phase-space density. By choosing different radii of the circular motion, the superfluid transition was explored. In this paper, we provide a quantitative understanding of the experiment using a c-field simulation method. We demonstrate that a blue-detuned laser of intensity comparable to the mean-field energy causes dissipation due to the creation of vortex-antivortex pairs. This is in contrast to laser stirring with a red-detuned laser Weimer2015 , where dissipation occurs via phonons Singh2016 . Furthermore, we study the relaxation dynamics of the stirred gas following the stirring process, which shows a slow energy transport between the condensate and the thermal cloud. We identify the origin of this slow relaxation to be vortex recombination and diffusion. We show that this effect can explain quantitatively the shift of the measured critical phase-space density in the experiment.

This paper is organized as follows. In Sec. II we describe the simulation method that we use. In Sec. III we determine the critical velocity $v_{c}$ of the stirred gas, based on which we identify the superfluid to thermal transition. In Sec. IV we discuss the dissipation via vortex pairs. In Sec. V we compare the simulation results with the experiment. In Sec. VI we analyze the relaxation of the stirred gas, and in Sec. VII we conclude.

II Simulation method

We simulate the stirring dynamics of a weakly interacting 2D bosonic system using the c-field simulation method that we used for a 3D system in Ref. Singh2016 . We describe this method in the following. We start out with the Hamiltonian of the unperturbed system,

[TABLE]

$\hat{\psi}$ and $\hat{\psi}^{\dagger}$ are the bosonic annihilation and creation operator, respectively. The 2D coupling parameter $g$ is given by $g=\tilde{g}\hbar^{2}/m$, where $\tilde{g}=\sqrt{8\pi}a_{s}/l_{z}$ is the dimensionless interaction, $m$ is the atomic mass, $a_{s}$ is the 3D $s$-wave scattering length, and $l_{z}=\sqrt{\hbar/(m\omega_{z})}$ is the harmonic oscillator length of the confining potential $m\omega_{z}^{2}z^{2}/2$ in the $z$ direction. $\omega_{z}$ is the trap frequency along the $z$ direction. $V({\bf r})$ describes the external potential, which is a harmonic trap, $V_{h}({\bf r})=m\omega_{r}^{2}r^{2}/2$. $\omega_{r}$ is the trap frequency in the radial direction and $r=(x^{2}+y^{2})^{1/2}$ is the radial coordinate. We introduce a time-dependent term to describe laser stirring,

[TABLE]

where $V({\bf r},t)$ is the time-dependent stirring potential and $\hat{n}({\bf r})$ is the density operator at the location ${\bf r}=(x,y)$. The stirring potential is a Gaussian with a width $\sigma$ and a strength $V_{0}$,

[TABLE]

which is centered at {\bf r}_{s}(t)=\bigl{(}x_{s}(t),y_{s}(t)\bigr{)}. We move $(x_{s},y_{s})$ along a circular path as a function of time $t$.

We perform numerical simulations by mapping this system on a lattice system, which also introduces a short-range cutoff; see Appendix A. This short-range cutoff is of the order of the healing length $\xi=\hbar/\sqrt{2mgn}$, with $n$ being the density. We describe both the equations of motion and the initial state within a c-number representation, which corresponds to formally replacing the operators $\hat{\psi}$ by complex numbers $\psi$. Furthermore, we approximate the initial ensemble by a classical ensemble, within a grand-canonical ensemble of temperature $T$ and chemical potential $\mu$. We sample the initial states via a classical Metropolis algorithm.

The simulation setup consists of a disc-shaped 2D circular condensate of 87Rb atoms. This choice of the 2D circular condensate is inspired by the experimental setup of Ref. Dalibard2012 . In the simulations we consider $N=38,162-93,267$ 87Rb atoms confined by the harmonic potential in both the radial and transverse directions. The trap frequencies are $\omega_{r}=2\pi\times 25\,\mathrm{Hz}$ and $\omega_{z}=2\pi\times 1.5\,\mathrm{kHz}$. Here the scattering length is $a_{s}=5.3\,\mathrm{nm}$, which yields $\tilde{g}=\sqrt{8\pi}a_{s}/l_{z}=0.093$. The temperature of the trapped gas is in the range $T=63-85\,\mathrm{nK}$. The simulation parameters that we use, are in the typical range of the experimental parameters of Ref. Dalibard2012 . For simulations of a quasi- and a pure-2D trap geometry we use a lattice of $180\times 180\times 5$ and $200\times 200$ sites, with the lattice discretization length $l=0.5\,\mu\mathrm{m}$, respectively. We choose $l$ such that it is smaller than, or comparable to, the healing length $\xi$ and the de Broglie wavelength $\lambda=\sqrt{2\pi\hbar^{2}/mk_{\mathrm{B}}T}$, see Ref. Mora2003 . The trapped gas is in the pure-2D regime if $k_{\mathrm{B}}T,\,\mu<\hbar\omega_{z}$. When $k_{\mathrm{B}}T$ and $\mu$ are comparable to $\hbar\omega_{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. In the experiment Dalibard2012 the trapped gas is stirred with a blue-detuned laser beam moving on a circular path around the trap center. For the circular motion of stirring we choose (x_{s},y_{s})=R\bigl{(}\cos(\omega_{m}t),\sin(\omega_{m}t)\bigr{)}, where $R$ and $\omega_{m}$ are the stirring radius and frequency, respectively. For the stirring potential we use the strength $V_{0}=k_{\mathrm{B}}\times 80\,\mathrm{nK}$ and the width $\sigma=1\,\mu\mathrm{m}$, in accordance with the experiment. The stirring sequence is the following: We linearly switch on the stirring potential over $5\,\mathrm{ms}$, let it stir the system for $200\,\mathrm{ms}$, and then switch it off over $5\,\mathrm{ms}$. This is again inspired by the experimental choices. We repeat this for various stirring velocities $v=R\omega_{m}$ by changing both $R$ and $\omega_{m}$. By choosing different $R$ we stir the different regimes of the trapped gas, the superfluid, the thermal, and the crossover regime.

After stirring we calculate the total energy $E={\langle}H_{0}{\rangle}$ using the unperturbed Hamiltonian in Eq. II, where we use $\psi$ instead of $\hat{\psi}$. From this energy we determine the equilibrium temperature $T_{\mathrm{eq}}$ of the stirred gas. We infer this temperature by numerically inverting the temperature dependence of the equilibrium state, $E_{\mathrm{eq}}=E_{\mathrm{eq}}(T)$. We elaborate on this in Appendix A. From the temperature difference between the stirred and initial system, the heating $\Delta T_{\mathrm{eq}}=T_{\mathrm{eq}}-T$ is determined. We also calculate the local energy, as well as the vortex and anti-vortex distribution. We define the local energy as $E_{i}=-\frac{J}{2}\sum_{j}(\psi_{i}^{\ast}\psi_{j}+\psi_{j}^{\ast}\psi_{i})+\frac{U}{2}n_{i}^{2}+V_{i}n_{i}$, where $j$ refers to the nearest neighbor sites. $\psi_{i}$, $n_{i}=|\psi_{i}|^{2}$, and $V_{i}$ are the complex-valued field, the density, and the trap potential at site $i$, respectively. $J$ and $U$ are the Bose-Hubbard parameters, see Appendix A. For the vortex distribution, we calculate the phase winding around the lattice plaquette of size $l\times l$, using $\sum_{\Box}\delta\phi(x,y)=\delta_{x}\phi(x,y)+\delta_{y}\phi(x+l,y)+\delta_{x}\phi(x+l,y+l)+\delta_{y}\phi(x,y+l)$, where the phase differences between sites is taken to be $\delta_{x/y}\phi(x,y)\in(-\pi,\pi]$. $\phi$ is the phase field of $\psi$. We identify a vortex and an antivortex by a phase winding of $2\pi$ and $-2\pi$, respectively. By counting all vortices and antivortices we determine the total number of vortices. We restrict this counting to the the superfluid region of the gas as we describe below.

III Superfluid response

To study the superfluid behavior we stir a 2D quasicondensate with a repulsive Gaussian potential. We prepare a trapped 2D quasicondensate of $N=93,267$ 87Rb atoms at temperature $T=85\,\mathrm{nK}$. We show the simulated density profile of the trapped gas in Fig. 1(a). We stir the gas with a circularly moving, repulsive stirring potential at stirring radius $R=14.4\,\mu\mathrm{m}$. As mentioned in Sec. II, we use the strength $V_{0}=k_{\mathrm{B}}\times 80\,\mathrm{nK}$ and the width $\sigma=1\,\mu\mathrm{m}$ for the stirring potential. This strength $V_{0}$ is well above the local mean-field energy $\mu_{\mathrm{loc}}\approx k_{\mathrm{B}}\times 21\,\mathrm{nK}$ at the stirring location. After stirring we determine the induced heating $\Delta T_{\mathrm{eq}}=T_{\mathrm{eq}}-T$ from the equilibrium temperature $T_{\mathrm{eq}}$ of the stirred gas, see Sec. II for details. By varying the stirring frequency $\omega_{m}$ we determine $\Delta T_{\mathrm{eq}}$ as a function of stirring velocity $v=R\omega_{m}$. We show $\Delta T_{\mathrm{eq}}$ determined for various $v$ in Fig. 1(b). The induced heating is almost negligible at low $v$, its onset occurs at a velocity $v_{c}$, and for $v>v_{c}$ it increases rapidly. We quantify the onset of heating using a fitting function,

[TABLE]

which is discussed in Ref. Pitaevskii2004 , with the free parameters $A$, $B$, and $v_{c}$. For the simulated heating shown in Fig. 1(b), this function gives a critical velocity of $v_{c}=0.58\pm 0.2\,\mathrm{mm/s}$. We compare this critical velocity to the Bogoliubov estimate of the phonon velocity $v_{\mathrm{B}}=1.4\,\mathrm{mm/s}$ at the stirrer location. The Bogoliubov velocity is determined by $v_{\mathrm{B}}=\hbar\sqrt{\tilde{g}n}/m$. The observed critical velocity is $v_{c}\approx 0.4\,v_{\mathrm{B}}$. This is notably different from the case of an attractive stirring potential, where $v_{c}\approx v_{\mathrm{B}}$ Singh2016 . We explain this reduction of $v_{c}$ for a repulsive stirring potential in Sec. IV.

By choosing different radii $R$ we explore the various regimes of the trapped gas. We use the same strength $V_{0}$ and the same width $\sigma$ as above. For each $R$, we first determine the induced heating $\Delta T_{\mathrm{eq}}$ as a function of $v$, and then by using the fitting function given in Eq. 4 we determine $v_{c}$. We show $v_{c}$ determined at various $R$ in Fig. 1(c). The stirring radii are in the range $R=10-20\,\mu\mathrm{m}$. For $R=10-16\,\mu\mathrm{m}$, there is no significant change of $v_{c}$. As $R$ reaches the crossover regime, $v_{c}$ is reduced sharply and for $R$ above the crossover regime, $v_{c}$ is zero. According to the BKT prediction in a uniform system Prokofev2001 with $\tilde{g}=0.093$ combined with local density approximation, the crossover regime should occur at $r_{\mathrm{BKT}}=16.8\,\mu\mathrm{m}$. This prediction is in good agreement with the crossover regime identified by the simulated $v_{c}$. Thus, we clearly identify the superfluid, the crossover, and the thermal regimes by the finite, the sharply decreasing, and the zero critical velocities $v_{c}$, respectively. We note that in the crossover region the decrease of $v_{c}$ is as sharp as the size of the stirrer allows. Furthermore, we note that the observed almost constant $v_{c}$ for $R<r_{\mathrm{BKT}}$ can be due to the accelerated circular motion and the large strength of the stirring potential Singh2016 .

IV Dissipation mechanism

The observed critical velocities are in the range $v_{c}=0.3-0.5\,v_{\mathrm{B}}$. To understand what leads to this reduction of the critical velocity with regard to the phonon velocity, we investigate the time evolution of the phase field $\phi({\bf r},t)$ of a single realization of the thermal ensemble. We obtain this phase field from the complex field $\psi({\bf r})$ via the phase-density representation $\psi=\sqrt{n}\exp(i\phi)$. In Fig. 2 we show the phase evolution of the trapped 2D quasicondensate stirred at $R=12\,\mu\mathrm{m}$. We use the velocity $v\approx 0.8\,v_{\mathrm{B}}$, which is above the steep onset of dissipation related to the breakdown of superfluidity. The phase evolution of the unperturbed gas shows rather weak phase gradients. As stirring is switched on, the phase field around the stirring potential starts to fluctuate. These fluctuations develop into strong phase gradients, which result in the creation of vortex-antivortex pairs. This can be confirmed by calculating the phase winding around each 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 and triangles, respectively. This indeed shows the creation of vortex-antivortex pairs during stirring. We recall that the stirring strength $V_{0}=k_{\mathrm{B}}\times 80\,\mathrm{nK}$ is much larger than the mean-field energy $\mu_{\mathrm{loc}}\approx k_{\mathrm{B}}\times 30\,\mathrm{nK}$ at the stirring location, which results in a strong reduction of the density at the stirrer location. This density reduction serves as a nucleation site for the creation of vortex pairs. We note that this mechanism of vortex-pair-induced dissipation is suppressed for an attractive stirring potential, as shown in Ref. Singh2016 . This scenario of dissipation induced by vortex pair creation is consistent with a recent experiment Shin2015 .

V Comparison to experiment

We now compare the results of our simulation with the experiment Dalibard2012 . We first show the comparison between the experiment and simulation for the heating as a function of $v$. In the superfluid regime, we stir the quasicondensate at the radius $R=14.4\,\mu\mathrm{m}$. The simulated density profile is shown in the inset of Fig. 3(a). After stirring we let the stirred gas relax for $100\,\mathrm{ms}$ of relaxation time and then determine the induced heating from the temperature of the wings of the cloud. We fit these wings to the Hartree-Fock prediction,

[TABLE]

with the fitting parameters $T$ and $\mu_{0}$. This method is adopted according to experiment, in which the temperature of the stirred gas is determined in the same way, following a relaxation of $100\,\mathrm{ms}$ as well. We denote this heating determined from the wing temperature $T_{\mathrm{w}}$ by $\Delta T_{\mathrm{w}}=T_{\mathrm{w}}-T$, with $T$ being the initial temperature. We show the simulated $\Delta T_{\mathrm{w}}$ and their comparison with the experimental ones for various $v$ in Fig. 3(a). The measured and simulated heating are found to be in good agreement if we base the comparison on $\Delta T_{\mathrm{w}}$. We also compare the measured $\Delta T_{\mathrm{w}}$ with the simulated $\Delta T_{\mathrm{eq}}$ determined from the equilibrium temperature $T_{\mathrm{eq}}$ of the stirred gas. We show $\Delta T_{\mathrm{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 difference at large $v$ is due to the absence of global thermal equilibrium of the stirred gas. As explained in Sec. VI, the stirred gas relaxes by transporting the excess energy between the superfluid in the central part and the thermal cloud in the periphery, which is a slow process. The absence of global thermal equilibrium leads to a smaller wing temperature than the equilibrium temperature.

The results shown in Fig. 3(a) indicate that the onset of heating occurs at a velocity $v_{c}$, and for $v>v_{c}$ heating increases rapidly. Both in experiment and simulation $v_{c}$ is determined using the fitting function in Eq. 4. In Fig. 3(b) we show the comparison between the experiment and simulation for stirring the thermal region of the trapped gas of $N=38,162$ atoms. The simulated density profile of the gas is given in the inset of Fig. 3(b). Both the measured and simulated heating $\Delta T_{\mathrm{w}}$ are in good agreement. The simulated $\Delta T_{\mathrm{eq}}$ determined from the equilibrium temperature of the system are below the measured $\Delta T_{\mathrm{w}}$ at large $v$. As we will explain in Sec. VI, this is again due to the absence of global thermal equilibrium. As the stirred thermal cloud has more excess energy than the condensate, the wing temperature is larger than the equilibrium temperature. The results shown 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 the experiment and simulation for $v_{c}$ that are determined by stirring the superfluid, the crossover, and the thermal regime. In the experiment Dalibard2012 $v_{c}$ is measured for different configurations of the total number of atoms $N$, the temperatures $T$, and the stirring radii $R$. We compare the measured $v_{c}$ with the simulated $v_{c}$ determined by stirring the 2D gas in Sec. III. We show both the measured and simulated $v_{c}$ as a function of the dimensionless parameter $\mu_{\mathrm{loc}}(r)/k_{\mathrm{B}}T$. The parameter $\mu_{\mathrm{loc}}/k_{\mathrm{B}}T$ characterizes the degree of degeneracy of the cloud and is the relavant parameter in the sense that the thermodynamic properties of the gas depend only on the ratio $\mu/k_{\mathrm{B}}T$ Yefsah2011 ; Prokofev2001 ; Hung2011 . We refer to $v_{c}$ determined from the wing temperature $T_{\mathrm{w}}$ and from the equilibrium temperature $T_{\mathrm{eq}}$ as $v_{c,\mathrm{w}}$ and $v_{c,\mathrm{eq}}$, respectively. Both the measured and simulated $v_{c,\mathrm{w}}$ show good agreement. The measured $v_{c,\mathrm{w}}$ and the simulated $v_{c,\mathrm{eq}}$ agree in the superfluid and thermal regime, while they differ in the crossover regime. For the measured and simulated $v_{c,\mathrm{w}}$, the crossover regime occurs at $\mu_{\mathrm{loc}}/k_{\mathrm{B}}T\approx 0.24$ and $0.22$, respectively. However, for the simulated $v_{c,\mathrm{eq}}$, it occurs at $\mu_{\mathrm{loc}}/k_{\mathrm{B}}T\approx 0.11$. The theoretical prediction for the BKT transition in a uniform gas Prokofev2001 with $\tilde{g}=0.093$ occurs at $(\mu/k_{\mathrm{B}}T)_{c}\approx 0.13$. This prediction is comparable to the simulated crossover regime identified by $v_{c,\mathrm{eq}}$, whereas its comparison with the crossover regimes identified by the measured and simulated $v_{c,\mathrm{w}}$ shows a shift. This shift was observed in Ref. Dalibard2012 , but could not be explained. We conclude that the experiments of Ref. Dalibard2012 can be reproduced quantitatively if the wing temperature is used, rather than $T_{\mathrm{eq}}$. This suggests that the system has not relaxed to thermal equilibrium after the waiting time of $100\,\mathrm{ms}$. We confirm and elaborate on this point and the underlying mechanism in the following sections.

VI Relaxation dynamics

We now investigate the relaxation of the system, following the stirring process in the superfluid regime. This includes a discussion of the influence of the confinement of the system in the $z$ direction. For strong confinement, the system approaches a purely 2D limit, while it is quasi-2D for intermediate confinement.

VI.1 Energy-flow dynamics

We first analyze the energy-flow dynamics of a stirred trapped gas in the purely 2D limit, and then compare this dynamics with a quasi-2D gas. For a pure-2D trapped gas, we consider a gas of $N=64,079$ 87Rb atoms, which is strongly confined in the transverse direction by the harmonic potential. The temperature $T=85\,\mathrm{nK}$ is smaller than the transverse trap energy $\hbar\omega_{z}/k_{\mathrm{B}}=144\,\mathrm{nK}$, so that the gas is in the ground state in this direction. As the width of the condensate in the $z$ direction is smaller than the lattice discretization length $l$, we simulate this system using a single $xy$-layer of lattice only, see Sec. II. We stir the gas at $R=12\,\mu\mathrm{m}$ for $200\,\mathrm{ms}$ at a velocity $v>v_{c}$. After that we switch off the stirring potential and let the gas relax. We calculate the local energy $E_{i}$ of the stirred gas and its final equilibrated local energy $E_{i}^{\mathrm{eq}}$, as described in Sec. II. We show the evolution of the excess energy \Delta\tilde{E}_{i}=\bigl{(}E_{i}(t)-E_{i}^{\mathrm{eq}}\bigr{)}/n_{\mathrm{max}} for various relaxation times $t_{\mathrm{relax}}$ in Fig. 4(a). $n_{\mathrm{max}}$ is the maximum density of the system. The evolution of $\Delta\tilde{E}_{i}$ after $t_{\mathrm{relax}}=0.1\,\mathrm{s}$ shows that most of the stirring-induced energy resides within the superfluid region. The system then relaxes by transporting this excess energy to the thermal cloud. This process occurs slowly and the system achieves fully equilibration only after about $2.5\,\mathrm{s}$ relaxation time, remarkably.

We now study this energy-flow dynamics for the case of a quasi-2D gas. We consider the quasicondensate that we use in Sec. III. The initial temperature of the gas and the harmonic potential in the transverse direction are equal and half of those in the pure-2D case, respectively. The resulting system is a quasi-2D gas. We simulate this system using five $xy$-layers of lattice in the $z$ direction. We stir the gas using the same stirring parameters as for the pure-2D case. We show the evolution of the excess energy $\Delta\tilde{E}_{i}$ of the stirred gas for various $t_{\mathrm{relax}}$ in Fig. 4(b). In this case, $E_{i}$, $E_{i}^{\mathrm{eq}}$ and $n_{\mathrm{max}}$ are the column (i.e. integrated along the $z$ axis) quantities. The evolution of $\Delta\tilde{E}_{i}$ after $t_{\mathrm{relax}}=0.1\,\mathrm{s}$ is similar to the pure-2D gas. Again, the system relaxes by transporting the excess energy to the thermal cloud. The equilibration process is slightly slower than for the pure-2D gas but again of the same order of $2.5\,\mathrm{s}$. We note that the measured relaxation times in Ref. Shin2012 are indeed of the order of the relaxation times that we find here.

In the experiment Dalibard2012 the waiting time after stirring and before measurement is $0.1\,\mathrm{s}$, which is shorter than the relaxation times that we observe here. This indicates that in the experiment thermal equilibrium between the superfluid and thermal cloud is not fully established, which influences the measured heating. It leads to respectively lower and higher measured heating for stirring the superfluid and thermal parts of the cloud.

VI.2 Vortex dynamics

To understand what causes this slow relaxation for the system, we now examine the evolution of the density and vortices of the system. We calculate the local density as $n_{i}=|\psi_{i}|^{2}$ and vortices as described in Sec. II. We show the density and vortex evolution of a single realization of the stirred pure- and quasi-2D gas after various $t_{\mathrm{relax}}$ in Fig. 5. For both systems, the density relaxation is hard to recognize, whereas the vortex evolution clearly exhibits decay of vortices. Thus, the system relaxes via decay of the induced vortices. Vortices can decay via both annihilation of a vortex with an antivortex, and drifting out to the thermal region of the cloud. For the pure-2D gas the number of vortices after $t_{\mathrm{relax}}=0.1\,\mathrm{s}$ is larger than the quasi-2D gas. Ref. Anderson2010 ; Rooney2011 reported that vortex annihilation in a pure-2D system is strongly suppressed as compared to a quasi-2D system because the vortex lines are impermeable to tilting Haljan2001 and bending Bretin2003 . So, the suppression of vortex annihilation can be a reason for the long-lived vortices for the pure-2D gas.

To make a quantitative comparison for the vortex relaxation between the two system, we count the total number of vortices within the superfluid region of the cloud at a detection radius $R_{\mathrm{det}}=14.5\,\mu\mathrm{m}$, and average it over $200$ realizations. We show the averaged vortex number $N_{v}$ as a function of time $t$ for both systems in Fig. 6. As stirring is switched on at time $t_{\mathrm{on}}$, $N_{v}$ starts to increase approximately linearly. It reaches its maximum $N_{\mathrm{max}}$ at $t_{\mathrm{off}}$. After the stirring is switched off, it decays approximately exponentially. For both systems the nature of vortex growth and decay are the same, but the rates with which they grow and decay are different. For the pure-2D gas the growth and decay rates are larger and smaller than those for the quasi-2D gas, respectively. The enhanced growth and the suppressed decay rate for the pure-2D gas can be due to the suppression of vortex annihilation, as mentioned above, and a slow vortex drift. We quantify the vortex decay rate using the function,

[TABLE]

with the free parameters $N_{1},\,N_{2},\,N_{0},\,\tau_{d_{1}},\,\tau_{d_{2}}$. From the fit, we determine $\tau_{d_{1}},\,\tau_{d_{2}}\approx 0.13,\,0.87\,\mathrm{s}$ and $0.05,\,0.87\,\mathrm{s}$ for pure- and quasi-2D gas, respectively. These decay times are similar to those determined from the mean excess energy in Appendix B.1. The fast decay $\tau_{d_{1}}$ and the slow decay $\tau_{d_{2}}$ are essentially connected to the vortex annihilation and drift lifetime, respectively. For the pure-2D gas $\tau_{d_{1}}$ and $\tau_{d_{2}}$ are larger than and equal to those for the quasi-2D gas, respectively.

We show in Table 1 $N_{\mathrm{max}}$ and the extracted $\tau_{d_{1}},\,\tau_{d_{2}}$ at varying $R_{\mathrm{det}}$, for both systems. $\tau_{d_{1}}$ and $\tau_{d_{2}}$ increase weakly as $R_{\mathrm{det}}$ is increased. However, the following conclusions are essentially independent of the choice of $R_{\mathrm{det}}$. Overall, $N_{\mathrm{max}}$ and $\tau_{d_{1}}$ are larger for pure-2D gas than those for quasi-2D gas, respectively, while $\tau_{d_{2}}$ are similar for both systems. We compare the simulated $\tau_{d_{1}},\,\tau_{d_{2}}$ of quasi-2D gas to the waiting time of $0.1\,\mathrm{s}$ in the experiment Dalibard2012 . This time is twice as large as the fast decay $\tau_{d_{1}}$, whereas it is smaller than the slow decay $\tau_{d_{2}}$. This suggests that most vortex recombination processes have occurred at the time of the measurement. However, the vortex drift to the thermal cloud has not occurred, and the system is in a metastable state, not in the equilibrated state. This is the mechanism that is responsible for the difference between the wing temperature and the equilibrium temperature.

We note that in Ref. Fedichev1999 an estimate for the time of 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. Fedichev1999 gives a timescale that is consistent with our simulation. We also note that the vortex lifetime is suppressed at high tempearures Ketterle2002 ; Jackson2009 ; Moon2015 .

VII Conclusions

We have studied the superfluid to thermal transition of a trapped 2D Bose gas of 87Rb atoms by stirring it with a repulsive stirring potential on a circular path around the trap center. The superfluid transition was probed by choosing different radii of the circular motion. We have identified the superfluid, the crossover, and the thermal regime by the finite, the sharply decreasing, and the zero critical velocity, respectively. The superfluid region of the gas yields critical velocities that are in the range $v_{c}=0.3-0.5\,v_{\mathrm{B}}$, where $v_{\mathrm{B}}$ is the phonon velocity. We have demonstrated that the onset of dissipation is due to the creation of vortex-antivortex pairs. The comparison of the simulation with the experiment shows good agreement if the temperature measurement of the experiment is imitated in the simulation, i.e. by extracting the wing temperature. However, we confirm the systematic shift that was observed in experiment, if thermal equilibrium is assumed. We have demonstrated that the absence of thermal equilibrium after the waiting time that was used in experiment is due to a remarkably slow relaxation mechanism: The energy transport across the superfluid to thermal interface occurs only on timescales of seconds. This slow transport mechanism is due to the slow drift of vortices out of the superfluid into the thermal wings of the system. We emphasize that this mechanism is relevant for many on-going experiments in the field of ultracold atoms, and their temperature measurements. Furthermore, this effect of suppressed transport across critical interfaces is in itself intriguing, and could be studied in a future cold atom experiment with clarity.

acknowledgements

We acknowledge support from the Deutsche Forschungsgemeinschaft through Grants No. MA 5900/1-1 and No. SFB 925, the Hamburg Centre for Ultrafast Imaging, and from the Landesexzellenzinitiative Hamburg, supported by the Joachim Herz Stiftung. JD thanks Jérôme Beugnon for many fruitful discussions and acknowledges support from ERC (Synergy grant UQUAM).

Appendix A Simulated heating

In this section we show how we determine the equilibrium temperature $T_{\mathrm{eq}}$ of a stirred gas using the c-field method described in Sec. II. We discretize the continuum Hamiltonian in Eq. II by the Bose-Hubbard Hamiltonian Jaksch1998 on a 2D square lattice,

[TABLE]

$\psi_{i}$ and $n_{i}=|\psi_{i}|^{2}$ are the complex-valued field and the density at site $i$, respectively. ${\langle}ij{\rangle}$ indicates nearest-neighbor bonds. For a lattice discretization length $l$, the Bose-Hubbard parameters are related to the continuum parameters via $J=\hbar^{2}/(2ml^{2})$ and $U=gl^{-2}$. The 2D coupling parameter $g$ is given by $g=\tilde{g}\hbar^{2}/m$, where $\tilde{g}=\sqrt{8\pi}a_{s}/l_{z}$ is the dimensionless interaction, $m$ is the atomic mass, $a_{s}$ is the 3D s-wave scattering length, and $l_{z}=\sqrt{\hbar/(m\omega_{z})}$ is the harmonic oscillator length of the confining potential $m\omega_{z}^{2}z^{2}/2$ in the $z$ direction. $\omega_{z}$ is the trap frequency in the $z$ direction. The harmonic trapping potential is $V_{i}=m\omega_{r}^{2}r^{2}/2$. $\omega_{r}$ is the trap frequency in the radial direction and $r=(x^{2}+y^{2})^{1/2}$ is the radial coordinate.

We first initialize the system in a thermal state at temperature $T$ via classical Monte Carlo, and then calculate its energy $E={\langle}H_{0}{\rangle}$ using the Hamiltonian in Eq. 7. By varying the temperature of the system $T$ while keeping the total number of atoms $N$ fixed, we calculate the energy $E$ as a function of $T$. In Fig. 7 we show the temperature dependence of the energy $E$ for the pure- and quasi-2D gas that are described in Sec. VI.

To determine the heating, we first stir the gas with the repulsive stirring potential as described in Sec. II and then after stirring calculate its energy $E$ using Eq. 7. We numerically invert this energy $E$ to the equilibrium temperature $T_{\mathrm{eq}}$ using the temperature dependence shown in Fig. 7. Finally, from the temperature difference between the stirred and initial system, the heating $\Delta T_{\mathrm{eq}}=T_{\mathrm{eq}}-T$ is determined.

Appendix B Relaxation dynamics

In this section we elaborate on the relaxation dynamics of the stirred trapped gas that we discuss in Sec. VI. We first discuss the energy flow dynamics in Sec. B.1 and then the vortex dynamics in Sec. B.2.

B.1 Energy-flow dynamics

Here we elaborate on the energy flow dynamics for the pure- and quasi-2D trapped system. We stir both systems with the stirring potential at velocity $v\approx 0.8\,v_{\mathrm{B}}$. After stirring we calculate the excess energy \Delta\tilde{E}_{i}=\bigl{(}E_{i}(t)-E_{i}^{\mathrm{eq}}\bigr{)}/n_{\mathrm{max}} as described in Sec. VI and by averaging this energy over the superfluid region of the gas, we calculate the mean energy $\Delta E_{\mathrm{mean}}$. We show the evolution of $\Delta E_{\mathrm{mean}}$ for both systems in Fig. 8(a). $\Delta E_{\mathrm{mean}}$ decays approximately exponentially as the excess energy $\Delta\tilde{E}_{i}$ within the superfluid region outflows to the thermal cloud. We quantify the energy decay time using the fitting function in Eq. 6. From the fit, we determine the decay times $\tau_{d_{1}},\,\tau_{d_{2}}\approx 0.13,\,0.98\,\mathrm{s}$ and $0.08,\,1.02\,\mathrm{s}$ for pure- and quasi-2D gas, respectively. These decay times are similar to the ones that we determine from the vortex decay in Sec. VI.2. In addition to the stirring at $v\approx 0.8\,v_{\mathrm{B}}$, we show $\Delta E_{\mathrm{mean}}$ corresponding to stirring at $v\approx 0.6\,v_{\mathrm{B}}$ for the pure- and quasi-2D gas in Fig. 8(b) and (c), respectively.

B.2 Vortex relaxation

Next, we turn to the vortex relaxation of the stirred gas. We stir the pure- and quasi-2D gas using different stirring radii $R$ and stirrer strengths $V_{0}$. We calculate the total number of vortices within the superfluid region of the cloud, as described in Sec. VI.2, and average it over $128$ realizations. We show the averaged vortex number $N_{v}$ normalized by its maximum vortex number $N_{\mathrm{max}}$ as a function of time $t$ for the pure- and quasi-2D gas in Figs. 9(a) and 9(b), respectively. For the pure-2D gas the relaxation of vortices is slower than that for the quasi-2D gas, as shown in Sec. VI.2.