Generating vortex rings in Bose-Einstein condensates in the line-source approximation

We present a numerical method for generating vortex rings in Bose-Einstein condensates confined in axially symmetric traps. The vortex ring is generated using the line-source approximation for the vorticity, i.e., the rotational of the superfluid velocity field is different from zero only on a circumference of given radius located on a plane perpendicular to the symmetry axis and coaxial with it. The particle density is obtained by solving a modified Gross-Pitaevskii equation that incorporates the effect of the velocity field. We discuss the appearance of density profiles, the vortex core structure and the vortex nucleation energy, i.e., the energy difference between vortical and ground-state configurations. This is used to present a qualitative description of the vortex dynamics.


I. INTRODUCTION
Since 1999, when vortex lines in a trapped Bose-Einstein condensate (BEC) were first experimentally obtained [1], their study has received great experimental and theoretical interest as it constitutes a clear signature of superfluidity effects in these confined systems. Some remarkable experimental achievements are, among others, the study of the dynamics of single vortex lines [2] and the formation of small [3] and large [4] vortex arrays. A review of the research done in this field is presented in Ref. [5]. It is worth to note that the experimental situation in BEC regarding the formation and detection of vortices is at variance with that in helium II droplets. For these drops, a paradigm of superfluid finite systems, the problem of nucleating vortices and their stability has only been addressed very recently from the theoretical point of view (see Refs. [6][7][8]), and their experimental detection is still an open question.
Vortex rings are vortices whose core is a closed loop with quantized circulation around it [9]. They are complex topological structures that have attracted and will continue to attract some experimental and theoretical interest. Based on numerical simulations, different methods have been proposed to generate vortex rings in BEC. As in bulk liquid helium, vortex rings may be produced introducing an impurity in the condensate with a definite velocity, whose displacement causes a vortex ring [10]. Another mechanism [11] consists in using dynamical instabilities in the condensate to cause dark solitons to decay into vortex rings. A well controlled method to produce vortex rings by electromagnetically induced atomic transitions in two-component condensates has also been put forward [12]. The method proposed in Ref. [11] has been successfully applied in Ref. [13] to generate vortex rings experimentally.
Rather than proposing a method that could be implemented to create a vortex ring experimentally, our aim here is to set up a numerical method simple yet accurate enough to generate quantized vortex rings with a definite radius R in one-component condensates. We restrict our study to vortex states such that the divergence of their velocity is vanishingly small. This assumption yields analytical expressions for the velocity field, and provides a fair approximation to the superfluid flow around the vortex core for rings in the bulk of large condensates (Thomas-Fermi limit). We have considered large condensates at zero temperature and therefore, dissipation has not been taken into account. They are axially symmetric about the z axis, and have the z = 0 plane as symmetry plane, and may host a vortex ring of radius R, coaxial with the symmetry axis of the trap, and placed on a plane at a distance Z ≥ 0 from the symmetry plane. Although generalization to situations with more than one such vortex rings is straightforward, we have not considered this possibility. This paper is organized as follows. In Sect. II we present the method used to generate vortices. The way to obtain the velocity field of a vortex ring is described in Sect. III, and the explicit expressions of the velocity components are given in an Appendix. Sect. IV is devoted to the analysis of the particle density profiles and vortex nucleation energies, which allows a qualitative description of their dynamics. We also present results obtained from a completely different, more involved method we have set up to generate ring vortices that permits to test the approximation of zero divergence of the velocity. Finally, a brief summary of the results is presented in Sect. V.

II. ENERGY FUNCTIONAL
We consider a weakly interacting Bose-condensed gas confined in a harmonic trap V ext (r) at zero temperature. In the Gross-Pitaevskii (GP) theory, the ground state (g.s.) energy of the condensate is given by the functional [14] where Ψ(r) is the condensate wave function. The first term in Eq. (1) is the kinetic energy of the condensate, the second term is the harmonic oscillator energy arising from the trapping potential, and the third term is the mean-field interaction energy. The coupling constant is g = 4πh 2 a/m, where a is the s-wave scattering length, and m is the atomic mass. The number of atoms in the condensate is dr|Ψ| 2 = N. The g.s. wave function is determined by solving the GP equation obtained minimizing the energy functional. The wave function Ψ can be written in terms of the particle density ρ(r) = |Ψ(r)| 2 and phase S(r) as and the superfluid velocity is given by v = (h/m)∇S. In this work we will use the equivalent quantum hydrodynamic description of the condensate in terms of the density and the superfluid velocity [14], since it allows a straightforward generalization of the energy functional Eq. (1) to include vortex states. Using Eq. (2) it follows that The first term is only density-dependent: The first term in E 0 [ρ] is the quantum kinetic energy. The second term of the energy functional Eq. (3) corresponds to the kinetic energy associated to the flow velocity v, and is given by The g.s. wave function in the absence of vortices has a spatially constant phase and therefore zero velocity. It can be obtained by minimizing E 0 [ ρ ]. We consider condensates with positive scattering length and axially symmetric traps V ext (r) = m[ω 2 ⊥ (x 2 + y 2 ) + ω 2 z z 2 ]/2 [= m(ω 2 ⊥ r 2 +ω 2 z z 2 )/2 in cylindrical (r, z) coordinates], with different values of the asymmetry parameter λ = ω z /ω ⊥ . The trap harmonic frequency ω ⊥ provides a length scale for the system, a ⊥ = (h/mω ⊥ ) 1/2 . We will use a ⊥ ,hω ⊥ , and N/a 3 ⊥ as units of length, energy, and density, respectively.
If we consider large condensates in which the Thomas-Fermi (TF) approximation holds [14], the quantum kinetic energy can be neglected compared to the interaction energy in Eq. (4), and the g.s. density of the condensate in the absence of vortices is given by ρ 0 (r, z) = µ(1 − r 2 /R 2 TF − z 2 /Z 2 TF )/g in the region where this expression is positive, and zero elsewhere. The TF extents of the condensate in the radial and axial directions are R TF = (2µ/mω 2 ⊥ ) 1/2 and Z TF = R TF /λ. The chemical potential µ is fixed by normalization µ =hω ⊥ (15λaN/a ⊥ ) 2/5 /2. We recall that the validity of the Thomas-Fermi approximation is guaranteed if Na/a ⊥ ≫ 1.
We turn now our attention to the case of condensates with vortex states characterized by a given irrotational velocity field associated with a non-vanishing quantized circulation. The total energy of the system is given by the energy functional Eq. (3). If the velocity field is known, for a given number of particles the density profile can be obtained minimizing Eq. (3). This yields the equation where ψ = √ ρ is the modulus of the complex wave function Eq. (2). This is the GP equation expressed in terms of the hydrodynamic variables.
In the presence of a quantized vortex, the density of the system drops to zero at its core, whose size is characterized by a healing length ξ. For very large condensates the healing length can be approximated by ξ = (8πρ 0 a) −1/2 , where ρ 0 is the density of the condensate before creating the vortex. For a centered vortex line in the TF approximation ρ 0 = µ/g, and the corresponding healing length ξ 0 is with ξ 0 ≪ a ⊥ ≪ R T F . In this approximation a local healing length can be defined [15] as Note that the size of the core is larger for a vortex in the low-density region.

III. SUPERFLUID VELOCITY FIELD
The velocity field around a straight vortex line has an analytical expression when the vortex is along the symmetry axis [5,16]. Approximate analytical expressions can be found in the case of vortex lines off the symmetry axis in large and very elongated, quasi-two dimensional condensates (see Ref. [17] and references therein). However, vortex lines generally bend in three-dimensional condensates [5,18] which renders impracticable an analytical treatment of the velocity field.
In the case of quantized vortex rings we proceed as Schwarz and Jang have done in the case of helium II [19] (see also Ref. [20]). For a vortex ring characterized by the values of (R, Z) already defined, circulation number n = 1, 2, . . . and quantum circulation k 0 = n h/m, we write the vorticity in the line-source approximation: where (r, z, φ) are the cylindrical coordinates, andφ is the unit vector in the azimuthal directionφ = (− sin φ, cos φ, 0). The superfluid velocity field that arises from this distributed vorticity fulfills ω = ∇ × v. Hence, the velocity field around the vortex is irrotational except on the vorticity line, where the density of the condensate is zero. If ∇ · v = 0 to a good approximation [21], a velocity vector potential A(r) can be introduced such that v = ∇×A. If the vorticity is specified, A(r) is determined by the equation whose integral solution is [19] The radial and z-components of the velocity are obtained as whereas the azimuthal component of the velocity field around the vortex ring Eq. (8) is zero. We give in the Appendix the general expressions of A 0 , v r and v z written in terms of hypergeometric functions [22]. For a confined system, the existence of a boundary and the fact that the density is inhomogeneous may have an effect on the actual velocity field [23]. In our case, since the condensate extends up to 'infinite distances', there is no need to introduce any image vortex to ensure that there is no particle flow across the boundary due to the superfluid motion. Moreover, it has been shown [15] that in the TF limit the corrections to the velocity due to density inhomogeneities can be safely neglected. Thus, we approximate the velocity field by Eqs. (10) and (11) (see Appendix), but have restricted ourselves to study vortex ring configurations (R, Z) whose vorticity line is inside the domain where the density is positive in the TF approximation, i.e., only vortex rings (R, Z) satisfying the condition (R/R TF ) 2 + (Z/Z TF ) 2 < 1 are considered. We call the (R, Z) line defined by the condition (R/R TF ) 2 + (Z/Z TF ) 2 = 1 the TF boundary. Once the the superfluid velocity field has been fixed, the density profile of the condensate is obtained solving Eq. (6).

IV. RESULTS
We present numerical results for large condensates in the Thomas-Fermi limit hosting a quantized vortex ring with circulation number equal to one (nucleating vortex rings with n > 1 is energetically less favorable [9,24]). For a singly quantized vortex ring configuration (R, Z) with vorticity given by Eq. (8) and circulation number n = 1, we have computed the velocity field, Eq. (11), and have obtained the density profile of the confined condensate by solving Eq. (6) in cilyndrical coordinates using the imaginary time method [25]. Note that although we are in the TF limit, we do solve the complete GP equation to obtain the density profiles. Figure 1 shows density profiles in the z = 0 plane as function of r for the experimental parameters of Ref. [13], that is, N = 3 × 10 5 atoms of 87 Rb (scattering length a = 5.82 × 10 −9 cm) in a spherical trap with ω ⊥ /2π = 7.8 Hz. We have plotted two configurations with a vortex ring located in Z = 0 having a radius R = 2 a ⊥ (dotted line), and R = 4 a ⊥ (dot-dashed line), respectively. The density profile of the condensate without vortex is also shown (solid line).
The density is zero on the vorticity line, as seen in the density profiles. As expected, the size of the core is of the order of the healing length. Indeed, for this condensate one has R TF ≃ 5.8 a ⊥ and ξ 0 ∼ 0.17 a ⊥ . Assuming that the vortex diameter is twice the local healing length, and using the TF expressions we get as core diameters the values 2×ξ(2, 0) = 0.37 a ⊥ , and 2 × ξ(4, 0) = 0.47 a ⊥ . One can see from Fig. 1 that the core sizes are in agreement with these estimates. When the radius of the vortex ring increases, the core lies at a lower density region and therefore its size increases.
To study the vortex ring energetics we have chosen a larger condensate made of N = 10 6 atoms of 87 Rb confined in an axially symmetric trap with axial frequency ω z /2π = 220 Hz and three different geometries, namely spherical, disk-shaped and cigar-shaped, with λ = 1, profile. Yet, for large condensates it is a rather local effect, as can be seen from Fig. 2, and also from Fig. 1 for the smaller condensate.
The nucleation of a vortex has an energy cost, since the energy of a condensate with a vortex is always larger than the energy of the condensate without it, E GS . The vortex nucleation energy, E − E GS , corresponding to the condensate of the spherical trap in Fig.  2 is plotted in Fig. 3 as a function of Z for different values of R. Fixed R, the nucleation energy is maximum when the ring is in the z = 0 plane, and it decreases as Z increases, since displacing the ring outwards implies that the superfluid flow affects less atoms in the condensate. Note that to have the vorticity line within the TF boundary, the R-curves end at different values of Z/Z TF given by 1 − (R/R TF ) 2 . The nucleation energy does not vanish when the ring is located on the TF boundary as there still exists a superfluid velocity field which produces some effect. Since we treat the superfluid velocity as an external field, the nucleation energy of very superficial vortex configurations may be somewhat overestimated. This drawback is apparent beyond the TF boundary, and shows up as a long-tailed nucleation energy as a function of either Z/Z T F or R/R T F , although it goes eventually to zero.
The nucleation energy is plotted in Fig. 4 as a function of R for different Z values. The smaller ring we have considered has R/R TF ∼ 0.1. The middle panel corresponds to the same spherical condensate as in Fig. 3. Besides, the top panel corresponds to the cigar-shaped trap, and the bottom panel to the disk-shaped trap. The top curve in all panels corresponds to vortex rings with Z = 0, and for a given trap geometry (fixed panel), the higher the curve, the lower the Z. As in Fig. 3, we have only described vortex ring configurations inside the TF boundary. For this reason, fixed a Z value the corresponding nucleation energy curve stops at the value R/R TF = 1 − (Z/Z TF ) 2 . Figure 4 shows that for each asymmetry parameter λ, given a Z there exists a vortex configuration corresponding to a certain radius R that maximizes the nucleation energy. This radius increases as Z decreases. In the (R, Z) two-dimensional configuration space, the nucleation energy E(R, Z) − E GS has only one maximum at (R eq /R TF ≃ 0.6, Z = 0). It is in agreement with the results obtained in Ref. [26] for a spherically symmetric trap using a simplified model of a vortex ring.
This maximum corresponds to the only stationary, although unstable, vortex ring configuration (equivalent to a vortex line nucleated along the z-axis), whose location is almost independent of λ but its value depends on it, decreasing as λ does. This means that the energy to nucleate a vortex ring characterized by the values of (R, Z) is lower the more elongated the trap is.
For each asymmetry parameter, the 'iso-nucleation energy' lines in the (R, Z) configuration space corresponding to a value not too far from the maximum value are closed curves around the maximum, and the allowed ring radii are in the range defined by the intersection of a straight line representing the given nucleation energy and the Z = 0 curve in Fig. 4. If we relax the condition that the vorticity line must be inside the TF boundary, this should be always like this [26]. Otherwise, the isoenergy lines eventually do not close. Yet, the allowed ring radii are determined in a similar way.
We have considered the vortex ring as if it were a static object. However, in a trapped condensate a vortex ring will move with a velocity that results from the interplay between the effect caused by the inhomogeneity of the condensate and self-induced effects arising from its own local curvature [5,9,26]. A qualitative analysis of the dynamics of coaxial vortex rings can be carried out from Figs. 3 and 4, as proposed in Ref. [26], and tested in a real time dynamics for two-component condensates [12,26]. It can be seen that a coaxial vortex ring oscillates inside the condensate, moving towards the surface and instead of being annihilated at the boundary, the ring generates a background flow in the rest of the cloud that draws it back towards the center. The ring dynamics in a trapped condensate results in an oscillation along the symmetry axis of the trap. The radius of the vortex ring grows or shrinks depending on its Z-position, in such a way that the ring motion almost follows a trajectory with constant energy in the (R, Z) configuration space.
Finally, we have used a completely different method, numerically more involved, to generate quantized vortex rings. It does not make any assumption on the superfluid velocity field, nor on the radius or position of the vortex ring either, and can be used to test the results we have obtained based on the validity of the approximation ∇ · v = 0.
We have proceeded by analogy with the case of a vortex line, and have added a constraint to the energy functional that forces to form a toroidal-like hole in the condensate with a flux of atoms around its core with zero azimuthal velocity component. We have used as constraining operator the angular momentum about theφ axis, i.e.,φ·L ≡ L φ = ih(r∂ z −z∂ r ) [27].
Introducing the functional where the Lagrange multiplier Ω φ can be understood as the local angular velocity about the azimuthal direction, we have obtained the associated GP equation. For a given Ω φ , we have solved the partial differential equations obeyed by the real and imaginary part of the condensate wave function, which are coupled by the presence of the constraint. These equations have been discretized on a (r, z) mesh using seven-point formulas to represent the differential operators, and have been solved using again the imaginary time method [25].
Taking as an example Ω φ = 0.2 ω ⊥ , we have obtained a coaxial vortex ring configuration with its vorticity located on the z = 0 plane whose equidensity lines cannot be distinguished from these represented in the middle panel of Fig. 2 at the scale of this figure (whose parameters were chosen before to allow us to make this statement). Figure 5 presents a comparison between the two methods. It shows the particle density profile as a function of r at z = 0 for the spherical condensate of Fig. 3. The solid line is the configuration that minimizes Eq. (12). The density goes to zero in R = 3.1 a ⊥ , and the vortex ring has an energy equal to 37.164hω ⊥ . The dashed line is the configuration obtained by the method of Sect. 2 with the ring vorticity Eq. (8) placed in (R = 3.1 a ⊥ , Z = 0) which has an energy of 37.167hω ⊥ . Thus, the energy and density profile of both configurations are almost identical, and the velocity fields are also in good agreement.
We have also checked that minimizing the functional Eq. (12) yields quantized ring vortices with n = 1. To this end, we have proceeded to a direct numerical integration of the circulation v · dl using several arbitrary closed paths in the (r, z) 'plane'. We have found that the circulation for paths enclosing the vortex core is h/m with a good accuracy, and is zero otherwise. The velocity field has been calculated from the wave function recalling that the current field is j The current field j(r) in the y = 0 plane is plotted in arbitrary units in Fig. 6. The axes are in units of a ⊥ . It can be clearly distinguished in this figure the position of the vortex ring core at Z = 0 and R ≃ 3.1 a ⊥ .

V. SUMMARY
We have developed a method to study vortex rings hosted in confined condensates in the TF limit. It provides analytical expressions for the superfluid velocity field that satisfy the irrotational condition and the quantization of the circulation of the superfluid flow. The velocity around the vortex core is introduced as an external field in the Gross-Pitaevskii energy functional, which is minimized to obtain the wave function of the condensate.
Using this method, we have computed the density profile and nucleation energy of vortex rings in one-component condensates. We have found that the size of the vortex ring core is of the order of the healing length as in the case of vortex lines. The presence of the vortex causes a sizeable distortion of the density, as atoms are pushed off the core, but the effect is rather local.
The analysis of the dependence of the nucleation energy on the position and radius of the ring allows one to conclude that a Hamiltonian dynamics will lead to oscillations of the vortex ring along the symmetry axis of the condensate, increasing and decreasing its radius in accordance with previous predictions. A dissipative dynamics would cause the vortex to decay, for example disappearing across the border of the condensate, as occurs for vortex lines [28].
The method provides a handleable way to generate vortex rings in confined condensates that can be used as a starting point to study the dynamics. A detailed calculation of the vortex ring dynamics in a confined condensate, either Hamiltonian or dissipative, is beyond the scope of the present work and will be developed elsewhere.
To evaluate E and K we have used polynomial approximations [29].