Critical velocity of a finite-temperature Bose gas

We use classical field simulations of the homogeneous Bose gas to study the breakdown of superflow due to vortex nucleation past a cylindrical obstacle at finite temperature. Thermal fluctuations modify the vortex nucleation from the obstacle, turning anti-parallel vortex lines (which would be nucleated at zero temperature) into wiggly lines, vortex rings and even vortex tangles. We find that the critical velocity for vortex nucleation decreases with increasing temperature, and scales with the speed of sound of the condensate, becoming zero at the critical temperature for condensation.


I. INTRODUCTION
A defining feature of superfluids is the absence of excitations when the flow (relative to some obstacle or boundary) is slower than a critical velocity; above this velocity, the flow becomes dissipative. This can be understood in terms of the Landau criterion, which predicts excitations when the local fluid velocity exceeds v L = min[E(p)/p], where p is the momentum of elementary excitations and E(p) their energy [1]. In weakly-interacting atomic Bose-Einstein condensates, and for infinitesimally small perturbations, one obtains v L = c, the speed of sound. The breakdown of superfluidity has been experimentally probed by introducing a localized repulsive obstacle, engineered via the repulsive force generated by focussed blue-detuned laser beam, and moving the condensate relative to the obstacle [2][3][4][5][6][7][8][9]. This has enabled measurement of the critical velocity and the direct observation of the ensuing excitations, that is, pairs of quantized vortex lines with opposite polarity. In flattened condensates, this scenario currently provides a route to engineer states of two-dimensional quantum turbulence [2,3]; it also gives insight into the deep link between quantum fluids and their classical counterparts, where it has been predicted that the wake of quantized vortices produced downstream of the obstacle can collectively mimick the classical wakes, including the Bérnard-von Kármán vortex street [10][11][12].
The motion of an obstacle in the zero-temperature Bose gas, described by the Gross-Pitaevskii equation, is a well-studied problem, particularly for circular obstacles in 2D geometries. The pioneering simulations by Frisch et al. [13] of an impenetrable circular obstacle moving within the 2D nonlinear Schrödinger equation (NLSE), demonstrated the existence of a critical velocity of value v c ∼ 0.4c above which vortex-antivortex pairs are nucleated. For small obstacles, boundary effects tend to suppress vortex nucleation, and, as the obstacle's size in- * nick.parker@ncl.ac.uk creases, the critical velocity reduces towards an asymptotic value [14][15][16]. The critical velocity also depends on the shape of the obstacle, for example, obstacles with elliptical cross-section lead to reduced/heightened v c , depending on the orientation relative to the flow [11,17]. Similar behaviour holds for spherical obstacles, albeit with the emission of vortex rings and increased critical speeds of circa 0.7c [11,18,19]. In current condensate experiments [2][3][4][5][6][7][8][9], the obstacles are penetrable, corresponding to a Gaussian potential of finite amplitude, produced via an incident blue-detuned laser beam. The same qualitative behaviour emerges for impenetrable obstacles, although the critical velocity and vortex nucleation patterns become modified [10].
Very recently, Kwon et al. have undertaken a systematic experimental analysis of the critical velocity for vortex shedding, exploring the dependence of the nucleation on height and width of the penetrable obstacle and the crossover from penetrable to impenetrable obstacles [8]. Their results, obtained in a condensate with temperature much lower than the critical temperature for condensation, are in agreement with previous zerotemperature predictions based on the Gross-Pitaesvkii equation. Their work has made a significant step in consolidating our theoretical and experimental understanding of the critical velocity in a condensate in the zerotemperature limit. At the same time, it has highlighted the need to extend the study of the critical velocity to finite temperatures. While the role of finite temperature has been explored considerably for another vortex nucleation scenario, namely within deformed, rotating traps [20][21][22][23][24][25] (for which unstable surface modes underpin the vortex nucleation), there is a paucity of literature relating to the finite temperature behaviour of vortex nucleation by a translating obstacle. Indeed, to our knowledge, the only finite-temperature analysis of a moving obstacle in a three-dimensional condensate is that of Leadbeater et al. [26], who found that the critical velocity of a hard sphere decreases with temperature.
In this work we study the motion of a cylindrical Gaussian-shaped obstacle through a three-dimensional homogeneous Bose gas at finite temperature via classical field simulations. We find that the critical velocity decreases with temperature and increases with condensate fraction (ratio of condensate to total density). Indeed, the critical velocity is found to be closely proportional to the speed of sound of the condensate, which scales as the square root of the condensate fraction. Above the critical velocity, vortex nucleation occurs either through pairs of vortex lines, collections of vortex rings, or direct formation of a vortex tangle, and we indicate the occurrence of these structures in the parameter space of condensate fraction and flow speed.

II. CLASSICAL FIELD METHOD
We consider a weakly-interacting Bose gas with N atoms in a periodic box of volume 3 . The atoms have mass m and their interactions are approximated by a contact pseudo-potential V int (r − r ) = gδ(r − r ), where g is a coefficient which characterises the atomic interactions and δ is the Dirac delta function [27].
In order to theoretically model thermal excitations of the weakly-interacting Bose gas, one must progress beyond the standard mean-field approximation to include both the condensate and the thermal fraction atoms in the gas. Various methods have been proposed for this purpose, as reviewed elsewhere [28][29][30][31][32][33]. Among these methods, a popular one is the classical field method [33][34][35][36][37][38][39]. This method is based on the observation that, providing the modes of the gas are highly occupied (an a priori assumption in our work), then the gas can be approximated by a classical field ψ(r, t) whose equation of motion is the Gross-Pitaevskii equation (GPE). However, whereas the GPE conventionally describes the condensate only, ψ(r, t) now describes the entire multi-mode 'classical' gas [28,29]. The classical field method has been used to model phenomena beyondmean-field effects, including thermal equilibration dynamics [36,38,40,41], condensate fractions [39], critical temperatures [42], correlation functions [43], spontaneous production of vortex-antivortex pairs in quasi-2D gases [44], thermal dissipation of vortices [45], and related effects in binary condensates [40,46,47].
We parameterize the gas by the classical field ψ(r, t). The density distribution of atoms is then |ψ(r, t)| 2 . The evolution of ψ is governed by the GPE where V obj (r, t) is the externally applied potential. The GPE conserves the total number of particles, N = |ψ| 2 dV , and the total energy, In what follows we express all quantities in terms of the natural units of the homogeneous Bose gas: density in terms of a uniform value ρ, length in terms of the healing length ξ =h/ √ mgρ, speed in terms of the speed of sound c = ρg/m, energy in terms of the chemical potential of the homogeneous condensate µ = ρg, and time in terms of τ =h/gρ. We label the modes of the system through the wavevector k. To allow for occupation across all classical modes of the system, the initial condition is highly nonequilibrium, where the coefficients a k are uniform and the phases are distributed randomly [36]. The occupation of mode k is n k = |a k | 2 . The final temperature/condensate fraction of our simulations is varied through a rescaling of ψ, so as to fix the quantities N and H. The GPE is evolved numerically, in the absence of any potential V obj , using a fourth-order Runge-Kutta method on a 192 3 periodic grid with time step ∆t = 0.01τ and isotropic grid spacing ∆ = 0.75ξ. The spatial discretization of our numerical grid implies that high momenta are not described in our simulations. In effect, an ultraviolet cutoff is introduced, n k (t) = 0 for k > k max , where k = |k| and the maximum described wave vector amplitude is k max = √ 3π/∆. The ensuing evolution from the strongly nonequilibrium initial conditions has been outlined previously [36,40]. Initially the mode occupation numbers n k are uniformly distributed over wavenumber k, up to the cutoff. Self-ordering leads to the rapid growth in the occupation of low-k modes, which initially evolves in a state of weak turbulence. Then the distribution evolves to a bimodal form. The high-k part of the distribution is associated with the thermal excitations and low mode occupations. The low-k part of the field is the quasi-condensate, characterised by macroscopic mode populations and superfluid ordering.
From the bimodal distribution, a wavenumber k c can be chosen as the boundary in k-space between the quasicondensate and the thermal gas, as performed in [36].
Here we take k c ≈ 13 (2πN −1 ξ −1 ), although our qualitative results are insensitive to the precise definition of k c . The condensate density, ρ 0 , is then calculated as the density within the quasi-condensate, i.e. a coarse-grained averaging over the quasi-condensate modes. This is then used to define the condensate fraction, ρ 0 /ρ, where ρ is the total density of the gas.
While the raw wavefunction is too noisy to allow direct visualization of vortical structures, this can be overcome by defining a quasi-condensate wavefunctionψ, as established in [36]. This wavefunction is constructed by filtering out high-frequency spatial modes from the classical field wavefunction, by transforming the complex amplitudes viaâ k = a k × max{1 − k 2 /k 2 c , 0}.ψ represents the long-wavelength component of the classical field.
The quasi-condensate features a tangle of quantized vortices which relaxes over very long times, and the final equilibrium state is free of vortices. Its physical properties, e.g. temperature and condensate fraction, are uniquely determined by the number of particles N and the kinetic energy E = (h 2 /2m)|∇ψ| 2 dV of the system [38]. The equilibrium state of the non-condensed particles follows the Rayleigh-Jeans distribution, modified by the nonlinear interaction with the condensed particles [38]. It is interesting to note that the equilibrium condensate fraction is insensitive to the number of modes, providing that the number of modes is, or exceeds, 16 3 modes. This suggests that this number of modes is sufficient to model the thermodynamic limit of the system. For comparison, we employ 192 3 modes.
Here we parametrise the system in terms of its particle density ρ = N/ 3 and average energy density H / 3 . Note that the total energy H = E + E 0 , where E is the kinetic energy of the system and E 0 is the energy of the condensate [38]. The temperature is evaluated from the condensate fraction using the following empirical relationship established in Ref. [45]: where T λ is the critical temperature for condensation and α = 0.2275 is a fitting parameter. Table I lists the parameters chosen in our simulations and the resulting condensate fractions and temperatures of the ensuing equilibrated classical field states.

A. Critical velocity for vortex nucleation
Having obtained the equilibrated finite-temperature states of the Bose gas, we now move on to consider a laser-induced obstacle moving through the gas. The obstacle, uniform in z, is translated in the x-direction at speed v. Our simulations are conducted in the frame moving with the obstacle, modelled by the inclusion of a Galilean shift term ihv∂ x ψ to the righthand side of the GPE. In this frame the obstacle is imposed through the time-independent potential V obj (r) = V 0 exp −(x 2 + y 2 )/d 2 , where d and V 0 parameterize the width and amplitude of the potential. The amplitude is linearly increased from V 0 = 0 at first introduction to its maximal value V 0 = 5µ over a period of 200τ . The frame speed is increased adiabatically to the required value according to the temporal profile v tanh(t/200τ ), wheret denotes the time from introduction of the obstacle. Simulations are repeated (from identical initial conditions) with increasing terminal speeds (in steps of 0.057c) until vortices are detected. Vortex detection is by visual inspection of the filtered density, up to a maximum simulation timet = 500τ (which is long enough to ensure that the obstacle is fully introduced and at terminal speed, but otherwise arbitrary). This procedure defines the critical velocity v c . There is a systematic uncertainty in our measurement of v c , arising from the discrete terminal speeds employed. Note that we have repeated this process for multiple randomized initial conditions, and find no change in our measurement of v c ; that is, the systematic uncertainty due to using discretized speeds is larger than the statistical uncertainty arising from different random initial conditions. Figure 1 shows the variation of v c with both condensate fraction ρ 0 /ρ (lower abscissa) and temperature T /T λ (upper abscissa), for two example obstacles widths. The critical velocity has a maximum value at zero temperature (unit condensate fraction), and decreases nonlinearly as temperature increases (condensate fraction decreases), reaching zero at the critical point for condensation.
At zero temperature, the critical velocity is of the order of the condensate speed of sound c = ρg/m, with a general form v c (T = 0) = βc, where β is a parameter which depends solely on the shape of the obstacle (here d and V 0 ). The simulated v c data in Figure 1 closely follow the simple functional form v c (T ) = v c (0) ρ 0 /ρ, as shown by the dashed lines. An expression for the critical velocity valid at zero and non-zero temperatures is v c (T ) = β ρ 0 g/m.
In other words, for a given obstacle, the critical velocity is a fixed fraction of the speed of sound based on the condensate density rather than the total particle density [26].
The inset of Figure 1 shows the variation of v c with the obstacle width d at finite temperature, for the example of T /T λ = 0.56. The qualitative behaviour is consistent with that seen at zero temperature [11,15,48]: for small d the critical velocity is sensitive to d (due to the prominence of boundary layer effects) but as d increases v c decreases towards a limiting values (the Eulerian limit). However, the critical velocities are systematically reduced compared to the zero temperature case due to the reduced condensate speed of sound at finite temperature.

B. Vortex nucleation pattern
Finally we examine the manner in which vortices are nucleated from the obstacle. At zero temperature, one expects the nucleation of straight anti-parallel vortex lines from the obstacle, either released in unison or staggered in time [10,11,13], which move downstream relative to the obstacle. At finite temperature, we observe three general regimes of vortex nucleation, with representative examples shown in Fig. 2: Vortex lines: A pair of "wiggly" vortex lines is produced [ Fig. 2(a)]. The wiggles are driven by the thermal fluctuations, which cause the vortex elements to be nucleated at slightly different times along the obstacle; this is visible at intermediate times (snapshots (iii) and (iv)). These elements ultimately merge together along the axis of the obstacle to form a wiggly vortex/anti-vortex line. Similar vortex configurations in the form of lines which are partially attached to a thin wire were also observed in liquid helium [49].
Vortex rings: Here vortices predominately form vortex rings [ Figure 2(b)]. The vortex loops generated at the obstacle rapidly peel away from the obstacle, reconnecting with adjacent loops to form rings. Due to the way the vortex rings form initially along the obstacle, they are elliptical and polarised such that they are longer along the obstacle axis. While the vortex line regime is analogous to the zero temperature case, no analog occurs for the ring and tangle regimes. We note that even a small amount of thermal fluctuations is enough to vastly change the form of vortex nucleation, such as the vortex rings produced in Fig. 2(b) for a condensate fraction of 0.91.
To systematically map the occurrence of these regimes, we measure the vortex line-length density L (length of vortex line per unit volume) and vortex polarity R (described below) at a fixed observation time oft = 500τ , throughout the parameter space of flow velocity and condensate fraction. Our method to evaluate the vortex linelength density is described in Appendix A. The results are presented in Fig. 3. Below the critical velocity (solid black line) no vortices are produced, and thus L = 0. Above the critical velocity, L increases strongly with the flow velocity. This is to be expected since the frequency of vortex nucleation increases with flow velocity [13]. L also increases with decreasing condensate fraction (increasing temperature), indicating the significant role of thermal fluctuations in enhancing vortex production.
Just above the critical velocity, where the vortex linelength density is relatively small, vortex nucleation occurs through vortex lines and rings. The low flow velocity ensures that the vortex nucleation frequency is low, thereby suppressing strong interaction or reconnection between nucleated vortices. Here, whether lines or rings are produced is sensitive to the random initial conditions, and so it is not possible to further distinguish these nucleation regimes within this parameter space. In these cases a more consistent characterisation of the vortex form is given by R, described below. At higher flow velocities, where the vortex line-length density is relatively high, the nucleation frequency becomes sufficiently high that vortices immediately undergo strong interactions with each other, reconnecting and developing into a vortex tangle. The transition in the parameter space from vortex lines/rings to tangles is indicated approximately by the dashed line, although statistical effects blur the true boundary.
We further characterise the vortex distribution by its polarisation through the quantity R = A z /(A y + A z ), where A y and A z are the total area of vortices when the vortex distribution is projected along the y and z directions, respectively. A value R ≈ 0 corresponds to vortex lines aligned predominantly along the z axis, R ≈ 1 corresponds to lines aligned predominantly along y, and R ≈ 0.5 corresponds to an isotropic vortex distribution (in the yz plane). The parameter space of R has the same qualitative form as that for L, increasing with velocity and decreasing with condensate fraction. R typically lies in the range 0.1 < ∼ R < ∼ 0.4 for the lines/rings regime, consistent with the presence of lines which are predominantly aligned along z and rings which are elongated along z. It is worth noting that while the occurrence of lines or rings, for a given flow velocity and condensate fraction, is sensitive to the initial conditions, the value of R is highly reproducible (to within a few percent). For the vortex tangle regime, 0.4 < ∼ R < ∼ 0.5. It is worth noting that this shows that the produced tangle can be highly isotropic, despite two-dimensional nature of the obstacle that generates it.

IV. CONCLUSIONS
Using classical field simulations, we have analysed the nucleation of vortices past a moving cylindrical obstacle in a finite temperature homogeneous Bose gas. We have evolved the classical field from highly non-equilibrium initial conditions to thermalized equilibrium states with ranging temperatures and condensate fractions. We have then inserted a cylindrical obstacle with Gaussian profile into the system, and imposed a flow relative to the gas. We have found that, above the critical velocity, vortices are nucleated forming wiggly anti-parallel pairs of vortex lines, vortex rings, or as a vortex tangle. The critical velocity decreases with increasing temperature, becoming zero at the critical temperature, and scales with the speed of sound of the condensate, i.e. as the square root of the condensate fraction. While our work is based on a homogeneous system, in reality Bose-Einstein condensates are experimentally confined in traps, rendering the gas inhomogeneous. Then one can expect corrections to the critical velocity due to density gradients, as well as modifications to the vortex nucleation pattern. These higher order effects could be studied in future work. However, we note that recent advances have led to the formation of quasi-homogeneous condensates in box-like traps [50,51], where these corrections should have minimal effect.
Data supporting this publication is openly available under an Open Data Commons Open Database License [52]. For a given wavefunction, ψ, featuring a vortex distribution, the vortex volume V t (the total volume associated with the vortex cores) is evaluated by numerical integration of the inside of the vortex isosurface tubes obtained from the filtered density |ψ| 2 , with an integration region of the whole numerical box. Note that the isosurface level should be low enough to pick out vortex cores only (and not, e.g. fluctuations and waves), while large enough to contain sufficient grid points to allow a reasonable numerical evaluation of volume. In this work we use the isosurface level 0.04 |ψ| 2 (chosen so as to produce filtered vortex cores that are similar in radius to the true vortex core). The volume calculation can be written V t = Θ(0.04 |ψ| 2 − |ψ(r)| 2 ) dV , where Θ is the Heaviside step function. In practice the calculation of the vortex core volume can be efficiently performed by assigning a value of unity/zero to grid points located within/outside the isosurface tubes and directly integrating the result using the trapezium or Simpson's rule.
The total line length is then deduced by dividing V t by the cross-sectional area of a vortex core, A t (in effect, the average cross-sectional area of the isosurface tubes).
The measured values of V t and A t will depend on the chosen isosurface level but, providing the vortex tubes are well-separated, their ratio (and hence the evaluated line length) will remain constant. For closely-positioned vortex tubes, the isosurface level can affect whether the tubes appear as two separate tubes, or start to merge, and so will lead to deviations in this ratio. We have tested the effect of an alternative isosurface value. For twice the original isosurface value, the difference in the calculated line length is negligible for systems with low vortex density. For cases with the highest density of vortices, the difference remains less than 10%.