Superfluid flow past an obstacle in annular Bose--Einstein condensates

We investigate the flow of a one-dimensional nonlinear Schrodinger model with periodic boundary conditions past an obstacle, motivated by recent experiments with Bose--Einstein condensates in ring traps. Above certain rotation velocities, localized solutions with a nontrivial phase profile appear. In striking difference from the infinite domain, in this case there are many critical velocities. At each critical velocity, the steady flow solutions disappear in a saddle-center bifurcation. These interconnected branches of the bifurcation diagram lead to additions of circulation quanta to the phase of the associated solution. This, in turn, relates to the manifestation of persistent current in numerous recent experimental and theoretical works, the connections to which we touch upon. The complex dynamics of the identified waveforms and the instability of unstable solution branches are demonstrated.

A characteristic feature associated with superfluidity is the existence of a critical velocity above which its breakdown leads to the creation of excitations. In experiments with BECs, evidence for a critical velocity was obtained by moving an obstacle, i.e. a tightly focused laser beam, through a BEC [12,13]. This setting has been demonstrated to be prototypical for dark soliton formation in one-dimensional (1D) [14,15] (see [16] for experiments, although the latter were only quasi-1D), and for vortex formation in 2D [17], which can be thought of as a type of nonlinear Cherenkov radiation. In the case of obstacles in a supersonic flow of the BEC, the formed Cherenkov cone [18][19][20] transforms into a spatial shock wave consisting of a chain of dark solitons [21]. The appearance of such radiation in photonic crystals [22] is yet another illustration of the importance of the fundamental study of critical velocity. The formation of vortex dipoles in a similar setting was also directly observed experimentally in the work of [23]. The case of a heavy impurity and the associated drag force were studied in [24].
In a homogeneous weakly interacting Bose gas the critical velocity is the same as the speed of sound, as per the associated Landau criterion [25]. Moving inhomogeneities can alter this critical value. For a ring geometry it has been shown in [26] that the instability of the superfluid is caused by outer and inner edge surface modes, in a similar fashion as in an infinite cylindrically symmetric tube with transverse harmonic confinement [27,28]. The different mechanism is due to the presence of a centrifugal force arising from the nature of the rotation. The effect of potential barriers in BECs confined in a ring trap has been studied experimentally and theoretically [7,[29][30][31][32][33]. The weak link due to the barrier, which affects the current around the loop, has a promising application as a closed-loop atomic circuit (atomtronics), e.g. as analogs of superconducting quantum interference devices [34,35]. The current-phase relation of a BEC flowing through a weak link was explored for a repulsive square barrier in [36]. The existence of a critical velocity above which superfluid flow stops in the ring is connected to Cherenkov radiation through the excitations of vortex-antivortex pairs [7,29], in analogy to the rectilinear case [23]. Nevertheless, the relevant instability remains somewhat inconclusive (in connection to corresponding experiments [7]) with different mechanisms proposed to account for discrepancies between theory and observation including thermal fluctuations [37] and imaging system resolution [38].
In the present study we consider a BEC confined in an effectively 1D annular trap with a moving potential barrier, which is equivalent to a stationary barrier and a moving condensate as realized in [7,30]. Using the mean-field, i.e. Gross-Pitaevskii (GP) approximation in the infinite domain, it was shown that the critical velocity v c corresponds to a saddle-center bifurcation of two branches of solutions [14], a stable (center) and an unstable (saddle) one. Using a 1D approximation (for a narrow ring geometry) in an effectively periodic domain, we reveal in an analytical and corroborate in a numerical fashion the existence of a sequence of saddlecenter bifurcations and associated critical velocities. These, in turn, correspond to different topological charges that are all connected within the same bifurcation diagram. We present numerical simulations as well as analytical calculations, where it is shown that the critical points can be obtained from solving two coupled nonlinear equations. We also observe the presence of a critical strength of the inhomogeneity (or length of the domain) above (respectively, below) which there is no critical velocity, i.e. the inhomogeneity can move with any velocity while preserving the superfluidity. This occurs when the ring circumference becomes shorter than the healing length (a setting that may thus be less relevant from a physical perspective) or when the obstacle is strong enough. Our examination reveals a series of unstable branches in the relevant dynamics; we explore the dynamical evolution of the solutions associated with these branches by means of direct numerical simulations.
The results presented here are intimately connected with recent experimental and theoretical observations. One of the early attempts to identify topological winding (and unwinding) in atomic BECs resulted in the seminal findings of [39][40][41]. In these works, rather than a defect rotating inside a BEC, a setting where a rotation was imposed on the entire quasi-1D ring BEC was examined. This has similarities but also substantial differences from our setup. A similarity is that the system is analytically tractable; in fact, it is a genuinely homogeneous system (1D in the co-rotating frame) where the effective 1D GP equation associated with the dynamics (including the rotational term) can be solved analytically by means of elliptic function cnoidal wave solutions which account for the phase slip events also identified here. On the other hand, there are nontrivial differences from that case. In particular, in our setting (and in recent experiments such as [7,30,31,34,35,42], the phase slips do not arise in a 'distributed' manner, associated with these periodic cnoidal solutions, but rather in a localized manner being co-located with the defect. Hence, the analytical considerations presented herein are expected to be more closely connected to recent experiments. Among the latter, a few [34,35] have been more directly related to the case with a pair of weak links or Josephson junctions, aiming at least in part to potential superconducting quantum interference related applications, while here we will focus solely on the realm of a single weak link. Arguably, the recent experimental settings most clearly related to our own work are those of [30,31]. The former one measures experimentally the emergence of the phase slips and uses a qualitative model based on the Bohr-Sommerfeld quantization condition and an approximate current-phase relation at the weak link to theoretically trace a structure similar to the one that we analytically identify in the present work (see their figure 4 and our bifurcation diagrams of figures 2 and 3 below). The latter work of [31], in fact, explicitly identifies the hysteretic dynamics that has been proposed to be a key characteristic of this system in the above figures. However, it also recognizes the disparity of the experimental observations from the GP findings (a feature that retraces discussions of earlier work mentioned above [37,38]). The analytical tractability of our findings in this system (in a sense, adapting to it the spirit of calculations performed earlier in the homogeneous rotated system of [39][40][41]) may offer further insight in the relevant comparison.
As a final step in this theme of comparisons, we would like to mention recent work, which has explored the case of a rotating weak link as a function not of the potential/domain parameters considered here (such as the barrier strength or the domain length), but rather as a function of the interaction strength [43]. This elaborate task requires different approaches in the weakly interacting limit (treated by means of a GP equation) and in the strongly interacting limit (treated by means of a Luttinger liquid approach and in the case of a Tonks gas by a Bose-Fermi mapping to the case of noninteracting fermions). Intermediate regimes were treated by density-matrix renormalization group computations which, in fact, revealed an unexpected optimality in the observed persistent currents at some intermediate interaction strengths between the above limits.

Theoretical setup
The 3D GP equation is given by [44]   y y y y y y , y t x, ( ) is the mean-field wave function, m is the atomic mass, D = ¶ + ¶ + ¶ xx yy zz is the Laplacian, μ is the chemical potential, 4 2 s the atomic interaction strength, which is proportional to the atomic scattering length a s , n is the number of atoms, V is the external trapping potential, and U is a short range potential representing the moving obstacle with an angular velocity ω. Here, the wave function is scaled by the integral BECs on a ring, with radius R, can be described by the GP equation (1) with a trapping potential that can be written in cylindrical coordinates q r z , , ( ) as where w w  , 1 r z , such that the dynamics of the BECs would be confined at r=R and z=0, which is the minimum of the confining potential V r z , ( ). Typical parameters used in experiments [30]  In this case, D = ¶ + ¶ + ¶ + ¶ qq rr r r r zz 1 1 2 and the moving potential U can be treated as a δ potential of strength α, i.e., q ad where y r 1 ( ) and y z 2 ( ) are the ground states satisfying the equations ( ) ( ) ( ) ( ) into (1) evaluated at r = R and z = 0 under the assumption that the atoms are concentrated at the minimum of the potential, one obtains approximately   y q y q y y y q y q m m y q (2) can be reduced into the scaled equation (after dropping all the tildes) y y k y y y a d q w y = - where the periodic boundary conditions along the azimuthal direction are , , , , .
A static BEC with a moving potential is equivalent to the flow of a nonlinear Schrödinger fluid past an immobile obstacle. Writing , , and considering the traveling frame (i.e., x x vt), equation (3) for the effectively 1D problem can be written as  y y y ky y y ad y = - In order to study the existence of persistent superflow, we search for a steady state solution (4), where * m is the chemical potential, which leads to (6) is solved simultaneously with the scaling equation, without loss of generality, = N L 2 L , using a Newton-Raphson method by discretizing the Laplacian with a central finite difference method [45].
We then examine the (linear) stability of a solution u(x) for which we introduce the linearization ansatz , where  is a formal small parameter, w  an eigenfrequency and (r, s) an eigenvector. Substituting it into equation (4) and keeping the linear terms in r and s, one obtains the linear eigenvalue problem ⎡ Note that eigenvalues in Hamiltonian systems come in complex quartets [46]. In particular, in the realm of the Schrödinger systems (3), if w  is an eigenfrequency of (7) with a corresponding eigenfunction r x s x T [ ( ) ( )] and T represents the transpose, then ) and * w  are also eigenfrequencies with corresponding eigenfunctions for all eigenfrequencies w .

Numerical results and connections to theory
In figure 1 we show one of the principal results of the present work, namely the bifurcation diagram of superfluid flow for varying velocity v starting from the static solution v=0 of (3) for a system with = = N L 2 10 L and a = 0.5. As we fix the norm, the bifurcation diagram is depicted in μ as a function of v. From figure 1 we observe that the solution experiences many saddle-center bifurcations (turning points) as indicated by black dots. Since steady flows do not exist beyond the turning points (i.e. for either larger or smaller v for the respective turning point), the abscissa of the bifurcation points corresponds to critical velocities v c . Note that hysteresis in 1D rings and in optical lattices has been reported before in [47], where it was argued that superfluidity can be naturally viewed as a hysteretic response to rotation.
In a ring trap geometry the phase of the BEC circulates around the center by an integer multiple of p 2 (see figure 1). The so-called topological charge q corresponds to how many times the phase winds along the ring. Macroscopic states with different q have distinct energies and the effect of q on persistent flow has been recently studied experimentally [8].
Considering the phase of the solutions along the branch in figure 1, it is interesting to note that the topological charge jumps along the branch segments that correspond to decreasing velocity v. More precisely, q increases at the points where the density at the obstacle vanishes. Hence, the solutions for all values of q are smoothly connected along the diagram. In figure 1, the upper insets show the density and the phase profile right before the charge jump (phase slip).
To provide a better understanding of the relation between the bifurcation diagram in ring systems (figure 1) and that of the infinite domain, which only has one saddle-center bifurcation [14], we now study bifurcation diagrams for different values of the domain length L. For this calculation it is preferable to fix the chemical potential μ and let the solution norm vary. Using a = 0.5 and m = 1, the results are shown in figure 2. Plotted is the square-root density of the stationary solution u 0 | ( )| against the velocity. It is known that in the infinite domain, the critical velocity corresponds to a saddle-center bifurcation between a dark soliton pinned to the obstacle and the uniform solution that is modified due to the inhomogeneity U(x) [14]. The continuation diagram of the solution in this case forms a loop with its symmetric counterpart (note that equation (4) is invariant under the transformation v v and x x) shown as black curve in figure 2. As v varies further, one will go around in the closed loop. In a ring system, when L is finite, we do not obtain a closed loop, but have connected 'loops' instead. The situation, when the curves touch the horizontal axis for finite L, corresponds to the creation of dark-soliton-like states at the position of the impurity. Exactly at these points the topological charge increases. While in the infinite domain there is only one velocity point, where dark soliton-like-states can be created by the impurity, there are several of these points in the ring systems. The implications of this feature lead to complex dynamics as will be discussed further below.  With decreasing length L the initial 'loops' become smaller and the distance between two consecutive 'loops' increases. There will be critical values of L when the shrinking 'loops' become points, i.e. pairs of saddle-node bifurcations collide. In that case, when one decreases the length further, the corresponding 'loops' will disappear as is the case in the green curve in figure 2.
We have also studied the effect of varying the potential strength α on the existence of steady state solutions. From our computations shown in figure 3, we obtain that as α increases, the bifurcation 'loops' are getting smaller and two consecutive saddle-center bifurcation points get closer to each other and then finally disappear at some value α. This implies that there is a critical potential strength parameter a cr above which there is no saddle-node bifurcation, i.e. the obstacle can move with any velocity. This may be interpreted as a condition when the obstacle is strong enough to pin the ring BEC such that moving the obstacle means moving the BEC as a whole and hence there is no relative velocity between the two. This informs us that there is a limiting α above which it is unfavorable for the BEC to have a non-vanishing density at the position of the barrier, i.e. the BEC moves together with the barrier. The green curve ( = L 1 2) in figure 2 also corresponds to such a situation. This is a feature that is not present in the infinite domain.
The bifurcation diagrams can be analyzed as follows (a similar derivation was presented in the supplemental material of [43] without connection with bifurcations that occur in the system). Using the Madelung's transformation ( ) ( ) ( ) , one obtains from the static version of (4)

Multiplying (8) with A and integrating yields
where C 1 is a constant of integration, which can be directly taken to be any number in the infinite domain [14]. Notice that equation (10) is suggestive (as discussed above in connection with figure 2) of the fact that where  A 0 0 ( ) , sharp gradients in the phase may arise; see the top profiles in figure 1. For a δ-potential U(x), one then obtains from (9) the first integral with C 2 being a constant of integration, and boundary conditions Using the equations in (2), which are equivalent to  which can be expressed in terms of the incomplete elliptic integral of the first kind [48]. Then, the solution of (11) is Hence, A(L) can be written in terms of A (0). Finally, using (10), one obtains the nonlinear algebraic equation that will yield the diagrams in figure 2, i.e.
is the topological charge. Figure 1 can be obtained similarly from solving (16) simultaneously with the )) for μ and A(0). It is then interesting to investigate the dynamics of solutions of equation (4) in two specific regions of the For the two unstable solutions depicted in the insets of figure 1, we show the time evolution dynamics in the laboratory frame in figure 4 by plotting the density distribution along the ring. In both situations a dark soliton is released from the inhomogeneity after initially traveling along the obstacle. The detachment of the dark soliton from the 'pinning' potential is the only dynamics of instability that we observed in all our simulations. Whether the soliton travels ahead or behind the obstacle depends sensitively on the perturbation. Note as well that the density of the cloud shows that after detachment the dark soliton in panel (b) moves faster than that in panel (a).
We then analyze the dynamics beyond a critical value, i.e. we take the steady state solution at a critical velocity v c as the initial condition and then compute the evolution for = + D v v v c . In figure 5(a) we show the dynamics in moving coordinate frame for clarity with = + v v 0.02 c for the second bifurcation point in figure 1 at v c = 1.2. We observe that a dark soliton is emitted from the impurity along the evolution.  In striking difference to the infinite domain [14], the released dark soliton interacts periodically with the impurity, which then traps the soliton for some time (e.g. between t = 40 and t = 50) before releasing it again. Note that the trapping time after each round trip is not constant. For a wide interval of D > v 0 ( < < v 1.2 1.3) we obtain similar dynamics with only one soliton present in the annular BEC. In the 'laboratory frame' (where the obstacle is rotating), this situation implies that the emitted dark soliton will tend to be standing, even though it also slowly drifts due to the temporary entrapment that it suffers from the moving obstacle.
When D = v 0.6, one or two dark solitons can be created within the ring, as is shown in figure 5(b). For the first few time units, there is only one soliton present within the ring. After some time an additional dark soliton can be created, and as a result at the time instant indicated with arrows, two solitons exist in the trap. The presence of several solitons in the trap and repeated interaction with the impurity can lead to complex dynamics including collisions between the solitons or annihilation of the soliton by the impurity.

Discussion and future challenges
In current state of the art ring traps that have been created e.g. by a spatial light modulator (SLM) [8] or by magnetic potentials [49], the impurity can easily be added through a focused blue detuned laser beam or by the SLM as well. In the setup of [7,34] the impurity is present and excitations have indeed been observed [30]. The predictions of this paper are therefore directly relevant for the effects observed in current experiments.
In summary, we have studied the creation of dark solitary waves (and more generally the generation of phase slips/ persistent currents) in a system with periodic boundary conditions. In contrast to the infinite domain case, for our bounded domain setting we find the existence of several critical velocities corresponding to different charges q of the stable solution. A somewhat unexpected feature was also the existence of sufficiently narrow domains or sufficiently strong obstacles for which no critical velocity could be identified. The ability to create coherent structures by increasing the velocity or to annihilate them through the impurity allows the creation-via the bifurcation diagram presented herein-of a controllable number of p 2 phase windings within the ring trap. The analytical tractability of this formation through the quasi-1D theoretical formulation proposed herein is a feature adding to the controllability of the process. The exact dynamics of the resulting structures can be highly complex including possible collisions and interactions and will be an interesting object for further study. It is especially relevant to systematically extend such considerations to higher dimensional contexts not only in 2D but also in 3D. In higher dimensions, the situations are distinct as the superfluid does not necessarily rotate like a solid body, see e.g., [29] for numerical studies of flow dissipation in 2D in the mean field regime in the presence of a static barrier, [33] for superfluid with q-topological charge in the presence of rotating barrier and [32] for persistent currents in the 3D case.