Instability and Vortex Rings Dynamics in a Three-Dimensional Superfluid Flow Through a Constriction

We study the instability of a superfluid flow through a constriction in three spatial dimensions. We consider a Bose-Einstein condensate at zero temperature in two different geometries: a straight waveguide and a torus. The constriction consists of a broad, repulsive penetrable barrier. In the hydrodynamic regime, we find that the flow becomes unstable as soon as the velocity at the classical (Thomas-Fermi) surface equals the sound speed inside the constriction. At this critical point, vortex rings enter inside the bulk region of the cloud. The nucleation and dynamics scenario is strongly affected by the presence of asymmetries in the velocity and density of the background condensate flow.


Introduction.
Even when flowing past an obstacle, a superfluid can be stationary, since its excitations appear only under certain conditions [1]. For instance, a homogeneous Bose-Einstein condensate (BEC) of density n running through a weak obstacle, whose presence only slightly changes the cloud chemical potential, is energetically unstable towards phononic excitations only above a critical velocity v c equal to the sound speed c ∝ √ n, according to the Landau criterion. On the other hand, the critical velocity and the nature of the excitations change in the presence of a strong obstacle that considerably perturbs the homogeneous flow. Helium flowing through narrow channels becomes unstable at velocities much lower than predicted by the Landau criterion, which this time would correspond to rotonic excitations, since, as first predicted by Feynman, vortices are nucleated inside the channel and carry energy away from the superfluid through the phase-slip mechanism [2][3][4][5]. Feynman's energetics argument is analogous to the Landau criterion, and is based on the choice of vortex rings, instead of rotons, as the excitations coming at the critical velocity. His conjecture has been refined by Varoquaux [5], who used half-rings as the excitations to be nucleated inside the channel. In general, the type of excitation and the critical velocity depend on the flow configuration.
In contrast to nearly-incompressible fluids like helium, atomic BECs, due to an exquisite control over the confining potential, offer a wide choice of flow geometries and obstacles, and also permit the modification of the effective dimensionality. Imaging techniques also provide a means of addressing the system on the healing length scale, at which vortex dynamics takes place. The superfluid critical velocity in an elongated BEC can be studied for instance by sweeping a laser beam through the cloud [6]. At the critical velocity, both vortices in an elongated trap [7] and solitons for tight confinement transverse to the flow [8] have been experimentally observed. Another way to study the stability of the superflow would be to perturb a persistent current, which has been experimentally demonstrated with a BEC in a torus [9]. Apart from being an hallmark of superfluidity, BECs in toroidal geometries and the study of their critical velocities have very important technological applications for sensing devices, like SQUIDs already realized with superfluid helium and superconductors. BECs are often very well described by the mean-field Gross-Pitaevskii equation (GPE), which provides, unlike for helium, a reliable theoretical model to study the instability mechanism and its dynamics. The nucleation of a vortex-antivortex pair at the critical velocity has been studied numerically in two and three dimensions with an obstacle moving through a uniform condensate [10], and in a two-dimensional waveguide with an orifice [11]. Flow dissipation caused by the formation of solitons in one dimension has been investigated in [12]. In a previous work [13], we studied the scenario of superfluid flow dissipation in an effective two-dimensional toroidal geometry.
In this manuscript, we study the critical velocity and superfluid flow dissipation mechanism in a three-dimensional constricted flow configuration, where the size of the cloud along all the three directions is much larger than the healing length ξ = / √ 2mgn. We consider a subsonic flow of a zero-temperature BEC in two different geometries: a wave guide with periodic boundary conditions, which can mimic an elongated cloud along the flow direction, and a torus. The unstable regime is reached by raising a repulsive penetrable barrier perpendicular to the flow. This barrier is broader than the cloud dimension and extends over a few (typically 5 to 10) healing lengths along the flow direction. In this way, we create a constriction for the flow in the barrier region, as shown in Fig. 1. Starting from a stationary flow with a given velocity, the barrier is adiabatically raised during the dynamical evolution until the instability sets in, and we observe vortex rings penetrating the cloud. In this way, the superflow energy is dissipated through the phase-slip mechanism [2]: vortices are able to take energy away from the flow since the velocity drops by a quantized amount in a region crossed by a vortex core [13]. The two geometries under study present significative differences in the vortex ring nucleation and dynamics. We find that vortex rings, which can find a stationary configuration after entering an axially symmetric waveguide, are instead always transient in the torus, and, more generally, as soon as the axial symmetry in the direction of flow is broken.

Formulation and Computation
In order to evolve these various cases, we solve a GPE (1) for three spatial dimensions in scaled form as where length, time, and energy are given in units of d o = [ mωo ] 1 2 , 1/ω o , and ω o respectively for a representative harmonic frequency ω o that characterizes the trap. V (r, t) is the external potential and g = 4πa do with a the inter-particle s-wave scattering length and m the atomic mass. ψ(r, t) is the condensate wavefunction normalized to the total number of particles N. The external potential has components associated with the trapping and the barrier potentials of the form In both cases, we find the ground state of the system with V b = 0.
The waveguide geometry is implemented by choosing the trapping potential with ω o ≡ ω x = 30 × 2πHz, γ = ω z /ω x , and periodic boundary conditions along the flow direction y. We considered three different values of the γ, namely 1, 1.05, and 1.2, which correspond, respectively, to a ground state chemical potential µ = 11.7, 12.0, 16.5 with N = 3 × 10 5 87 Rb atoms and a nonlinear scaling value g N = 10134 .
The toroidal geometry is implemented by choosing the trapping potential where α = ω x /ω z and β = ω y /ω z . We take α = β = 0.5 with ω o ≡ ω z = 25 × 2πHz and form the torus by including a core potential with parameters V c = 144, and σ c = 1.88 with The ground state chemical potential is µ = 7.6 for N = 2.5x10 5 23 Na atoms, corresponding to a nonlinear scaling value g N = 2028.
During the dynamical evolution, we use a time-dependent barrier potential of the form: with f (t) = t/t r (f (t) = 1 for t > t r ) and V bx = tanh( x−Rx+x 0 bs ) + tanh( −x+Rx+x 0 bs ). Here R x is the x-shift of the center of the barrier while its width is w x ∼ 2x 0 . V by (y) and V bz (z) have the same form as V bx . The final height of the barrier is V s as long as x 0 , y 0 , z 0 ≫ b s .
Before starting the dynamics, we put the condensate in motion by imprinting an appropriate spatially dependent phase θ(r) on the wavefunction. In the waveguide case, θ(r) = mvy/ generates a uniform flow of velocity v along the y direction. In the torus, θ(r) = mlφ/ , with l integer, generates a tangential flow of speed v(r) = l/ρ, where The results presented in this manuscript are obtained by solving numerically the GPE.
After finding the ground state with an imaginary-time propagation, we perform the dynamical evolution in real time. For all the simulations performed we used a finite-element discrete variable representation (DVR) [22]. In the waveguide geometry, we employed a splitoperator method, with a Fast-Fourier-Transform algorithm used for the kinetic part of the evolution, while for the simulations in the torus we used a Real Space Product Formula [22]. The spatial grid for the waveguide calculations was n x = 60, n y = 120, n z = 60 (or

Criterion for instability
In the hydrodynamic regime, when the healing length ξ is much smaller than all other length scales (the smallest of which is the barrier width along the flow direction), we find that the onset of instability coincides with a simple condition: inside the barrier region, the local fluid velocity at the classical (Thomas-Fermi) surface of the cloud equals the sound speed, as depicted in the inset of Fig. 2. It is important to precisely define the sound speed which sets the threshold for the critical velocity. The latter is the Bogoliubov sound speed, calculated inside the barrier region, for the low-lying modes propagating along the flow, taken as if the system was homogeneous in this direction. In the case of a waveguide with a harmonic transverse confinement, and within the Thomas-Fermi approximation, this sound speed is simply the average of the local sound speed on a plane perpendicular to the flow is the local sound speed at the center of the transverse harmonic trap [17].
We verified this numerically in the waveguide case, as shown in Fig. 2. In a waveguide with cylindrical symmetry, this happens simultaneously at points on a circle perpendicular to the flow, from where a vortex ring will then enter. In the torus, due to a higher flow speed at the inner edge of the annulus, the critical condition is reached first on the interior of the cloud, with the consequence that vortex rings, if ever formed, must be asymmetric, as we discuss below.
Slightly below the critical point, we observe vortices getting closer to the edges of the cloud but failing to enter. The presence of these "ghost" vortices [13,14] suggests that a pre-instability is triggered at the edges of the condensate, but, since the vortices do not enter the bulk region, it is not sufficient to dissipate the superflow.
In a previous work [13], we verified the same condition for instability to hold for a twodimensional toroidal BEC. The results reported here represent an extension to a threedimensional case. As discussed above, this criterion allows for a simple understanding of the details of vortex penetration dynamics, based on the observation that vortex cores enter first where the critical condition if first reached.

Instability dynamics in the waveguide
In analogy with the standard situation where a superfluid flows through a channel, we dynamics. Such a setup has been used in [16], where, moreover, the possibility of creating a penetrable repulsive barrier with a control over a length comparable to the healing length has been demonstrated.
In this configuration, as soon as the critical barrier height is reached, a vortex ring detaches from the system boundaries and starts shrinking into the cloud inside the barrier region, as shown in Fig. 3, a). In the figure, black points indicate the position of vortex cores while the gray surface corresponds to the Thomas-Fermi surface of the cloud. Detection of vortex cores is performed using a plaquette method, described in [21].
Depending on the initial flow velocity, waveguide transverse section, and barrier height at the critical point, the ring attains a certain radius and velocity with which it propagates in the flow direction in a stable fashion, as long as axial symmetry is preserved, as depicted in Fig. 3, c) and d). Here the vortex ring eventually propagates at the speed u r = 0.65, and a radius R = 2.9.
For sufficiently strong barriers, the ring shrinks to a point crossing the transverse section, thereby annihilating and completing a full single phase-slip, as shown in Fig. 4. After this process, the velocity has dropped everywhere by the same quantized amount. We observe that such an event is the three-dimensional analogue of vortex-anti-vortex annihilation in two dimensions [13]. It is interesting to analyze the ring annihilation event in more detail.
In Fig. 5 a) and b), respectively, we show the position of the vortex cores toghether with the density on a plane parallel to the flow direction, at a time just before and just after the ring has shrunk to a point. In Fig. 5 c), the density on the same plane is plotted at four subsequent times after the annihilation. The details of self-annihilation process which we observe are consistent with the previous studies of axisymmetric vortex ring solutions [23][24][25][26][27]. Namely, if we consider the position of the vortex cores, we see that the points of phase singularity form a loop which shrinks inside the constriction, Fig. 5 a), and whose radius eventually becomes zero, Fig. 5 b). At this moment, we see the zero-density core being filled with atoms, thereby transforming into a density depression, which further propagates out of the constriction as a rarefaction pulse, as shown in Fig. 5 c). Finally, this rarefaction pulse decays into sound. Figs. 3 and 4, we see that, at least initially, the vortex cores are subjected to transformed into a rarefaction pulse, whose propagation and decay into sound appears in c).

Both in
an essentially radial motion, giving rise to the shrinking of the vortex ring. This behavior can be understood by considering what contributes to the motion of vortex cores in nonhomogeneous systems [19]. Let us consider a single vortex core located at some position x c . Its velocity would be the sum of two terms: i) the background flow velocity at x c , and ii) a term perpendicular to the gradient of the density at x c . Both contributions must be calculated as if the vortex was not present. Since the relative weight of term ii) is proportional to the value of the healing length calculated at x c , in high density regions, a vortex core will move mainly with the background superfluid velocity, while in low density regions, it will move mainly due to the gradient of the density. Therefore, when the vortex ring is inside the constriction, where the density is low, and especially when it is close to the Thomas-Fermi surface, it will mainly move perpendicular to the gradient of the density, whose main contribution comes from the density modulation induced by the barrier along the flow direction. This results in an essentially radial motion of the cores, and thus into the shrinking of the ring.

Instability dynamics in the torus
In the toroidal geometry, the scenario is richer since the tangential flow velocity at a given total angular momentum, due to the conservation of quantized circulation, decreases like 1/r, where r is the distance from the center of the torus. This introduces an asymmetry between the inner and the outer edges of the cloud, which is not present in the waveguide case. Consequences of this asymmetry on vortex nucleation in two dimensions are discussed in [13]. In three-dimensional configurations, the result is that vortex rings are transient features in a toroidal geometry.
The ring vortex either breaks or shrinks to a point and annihilates for sufficiently strong barriers. While the annihilation process, shown in Fig. 6, is very similar to the one taking place in a waveguide, the vortex ring breaking mechanism reflects instead more evidently the asymmetry in the velocity field. As shown in Fig. 7, the ring is strongly deformed since the inner part moves faster (the velocity of a vortex core, when the healing length is much smaller than the length scale of density variation, is essentially given by the background flow velocity [19]). The deformation increases in time up to the point at which the ring bends in on itself, forming a right angle at whose vertex a kink is present, as shown in Fig. 7 b). Meanwhile, a vortex line coming from the inner core of the torus approaches the vertex, also forming a kink in correpondence to the latter. Eventually, the line and the ring connect, joining each other at the position of the kinks. Such an event produces a vortex line plus a new ring vortex, appearing in  Fig. 7 c). Thus, the breaking is actually a vortex reconnection. Vortex reconnections play an important role in turbulent scenarios, being a very efficient mechanism for increasing the number of vortices in the system [20]. When the barrier is not strong enough to make both the edges of the annulus unstable, we also observed cases in which the ring is not formed at all, and a strongly bent vortex line enters the cloud from its inner edge, to circulate around the torus. The fact that the vortex line enters the inner edge of the annulus is due to the above mentioned asymmetry in the velocity field, decaying as the inverse of the distance from the center of the torus, which makes the instability set in there first, according to the criterion discussed in section .
The formation of a vortex ring in the torus is very interesting, since it can be seen as a reconnection of a vortex line with itself. As shown in Fig. 8, a bended vortex line is always present at first. While the bending increases, the vortex line develops two sharp kinks whose tips get closer to each other, up to when they join, thereby cutting the original line into a vortex ring plus a yet another line. The latter is then reabsorbed at the system's boundary.
Connection between effective two-and three-dimensional geometries When the condensate is effectively in two dimensions, meaning that the cloud size along the third direction is comparable with the healing lentgth ξ, vortex lines oriented along the direction of tighter confinement can become the preferred excitations with respect to vortex rings [18]. In studying a condensate inside an effective two-dimensional toridal trap (the system's size along the axis of the torus was about one healing length) [13], we indeed observed that superfluid flow dissipation took place through the formation of vortex lines entering the cloud. For the trap configurations considered in this manuscript, we have observed vortex ring formation in the very low density regions outside the condensate classical Thomas-Fermi surface. As suggested above, we could call these ghost [14] vortices since they are not visible from the condensate density profile. In some cases, as discussed above, these ghost vortex rings are able to enter completely the classical surface of the cloud, transforming into what we should then call "real" ring vortices.
In the crossover between effective two-and three-dimensional regimes, moving from a scenario in which only vortex lines are present to one in which real vortex rings come into play, it is reasonable to expect an intermediate regime in which ghost vortex rings are formed, but the condensate is sufficiently squashed along the third dimension that the full vortex loop is not able to enter the cloud's classical surface. In this regime, superfluid would be dissipated by vortex rings which are partly real and partly ghost, appearing as simple bent vortex lines in the density profile.
An example of such situation is given in Fig. 9, where the waveguide is non-axially symmetric about the flow direction (γ = 1.2). After the instability sets in, a ghost vortex ring forms and shrinks around the cloud in the barrier region up to when we observe a fullfledged ring vortex, which is only partially inside the Thomas-Fermi surface of the condensate (Fig. 9a)). However, since the part of the ring which is inside the Thomas-Fermi surface moves with a larger speed along the flow direction with respect to the part which remains outside, the ring vortex is soon deformed (Fig. 9b)), and eventually breaks up (Fig. 9c)).
The vortex cores located inside the Thomas-Fermi surface of the cloud move, to a good approximation, along with the background velocity field. On the other hand, cores in the low density region outside the surface of the cloud essentially do not feel the background velocity and move along the flow direction only because of the presence of transverse density gradients. After the vortex ring breaks up, the two lines move downstream and continue to deform, and eventually re-join to form a vortex ring (Fig. 9d)). The latter undergoes the same deformation described above, leading to another break up.
We also verified that, even in a very slightly non-axially symmetric waveguide, ring vortices eventually break up. With γ = 1.05, the deformation is created more slowly with respect to the γ = 1.2 case, but the ring breaks up anyway, though at a later time.

Conclusions
We studied the superfluid instability condition and dynamics in the presence of a threedimensional constriction across the flow. We solved the GPE, which would accurately describe an ultracold dilute Bose gas, but would anyway provide a useful insight for the under-standing of the instability in strongly correlated superfluids like Helium II. We considered two different setups: a waveguide, which could be experimentally realized by trapping an ultracold Bose gas inside a cigar-shaped trap, and a torus. We found that the instability sets in when the fluid velocity at the Thomas-Fermi surface of the cloud equals the sound speed inside the constriction region. This critical condition explains the strong dependence of the instability dynamics on the geometry considered, and in particular on the asymmetries present in the flow velocity and density. In an axially-symmetric waveguide, the superfluid dissipates through the nucleation of circular ring vortices that shrink across the flow within the barrier regions. As soon as the axially symmetry is broken, the rings are strongly deformed and can break up into vortex lines or not form at all. The phase slip instability observed here can be considered as a superfluid flow dissipation mechanism, since it subtracts energy from the potential flow of the superfluid. However, GPE does not include the dissipative processes which bring the system to thermal equilibrium. In a realistic situation, the vortices, carrying the energy subtracted from the superfluid, are supposed to eventually thermalize, with the result of heating the system. The effect of dissipation on vortex nucleation has been studied by adding a dissipation term to the GPE [15]. The role of finite temperature on the superfluid instability scenario presented here deserves further study.