# Vortex stability of interacting Bose-Einstein condensates confined in anisotropic harmonic traps

\[

## Abstract

Vortex states of weakly-interacting Bose-Einstein condensates confined in three-dimensional rotating harmonic traps are investigated numerically at zero temperature. The ground state in the rotating frame is obtained by propagating the Gross-Pitaevskii (GP) equation for the condensate in imaginary time. The total energies between states with and without a vortex are compared, yielding critical rotation frequencies that depend on the anisotropy of the trap and the number of atoms. Vortices displaced from the center of nonrotating traps are found to have long lifetimes for sufficiently large numbers of atoms. The relationship between vortex stability and bound core states is explored.

\]

The recent experimental achievement of Bose-Einstein condensation (BEC) in trapped ultracold atomic vapors  has provided a unique opportunity to investigate the superfluid properties of weakly-interacting dilute Bose gases. Mean-field theories, which are usually based on the Bogoliubov approximation  or its finite-temperature extensions , yield an excellent description of both the static and dynamic properties of the confined gases . These theories also predict that a continuum Bose condensate with repulsive interactions should be a superfluid, which can exhibit second sound, quantized vortices, and persistent currents. While there exists some evidence for second sound in trapped condensates , vortices in these systems have never been observed despite considerable experimental effort . Numerous techniques for the generation of vortices have been suggested, including stirring the condensate with a blue detuned laser , adiabatic population transfer via a Raman transition into an angular momentum state , spontaneous vortex formation during evaporative cooling , and rotation of anisotropic traps .

Several studies of vortex stability recently have been carried out . The free energy of a singly-quantized vortex attains a local maximum when the vortex is centered in a stationary trap . In the presence of dissipation, such a vortex would migrate to the edge of the trap and eventually disappear . It has been suggested  that this mechanical instability may be related to a bound state in the vortex core, corresponding to a negative-energy ‘anomalous’ dipole mode found numerically in the vortex state at low densities . It is possible to stabilize singly-quantized vortices by rotating the trap at an angular frequency $`\mathrm{\Omega }`$. When $`\mathrm{\Omega }`$ is larger than the ‘metastability’ frequency $`\mathrm{\Omega }_0`$, infinitesimal displacements of the vortex no longer decrease the system’s free energy, and the vortex becomes locally stable; above the critical frequency $`\mathrm{\Omega }_c\mathrm{\Omega }_0`$, the vortex state becomes the ground state of the condensate .

While rotation of the confining potential has been proposed as a method to both nucleate and stabilize vortices in trapped Bose gases, the relevant critical angular frequencies are not presently known. In the present work, numerical results are obtained for Bose-condensed atoms confined in three-dimensional rotating anisotropic traps at zero temperature. The critical frequency $`\mathrm{\Omega }_c`$ is found to increase with the degree of anisotropy in the plane of rotation. In order to nucleate vortices, however, the trapped gas must be rotated either more rapidly than $`\mathrm{\Omega }_c`$, or at temperatures above the BEC transition . While vortices in nonrotating traps are found to be always unstable, for large numbers of atoms their lifetimes can be very long compared with a trap period.

The trapped Bose condensate, comprised of $`N_0`$ repulsively-interacting Rb atoms with mass $`M=1.44\times 10^{25}`$ kg and scattering length $`a100a_0=5.29`$ nm , obeys the time-dependent Gross-Pitaevskii (GP) equation in the rotating reference frame :

$$i_\tau \psi (𝐫,\tau )=\left[\frac{1}{/}2\stackrel{}{}^2+V_\mathrm{t}+V_\mathrm{H}\mathrm{\Omega }L_z\right]\psi (𝐫,\tau ),$$
(1)   
where the trap potential is $`V_\mathrm{t}=\frac{1}{/}2\left(x^2+\alpha ^2y^2+\beta ^2z^2\right)`$, the Hartree term is $`V_\mathrm{H}=4\pi \eta |\psi |^2`$, and the condensate is rotated about the $`z`$-axis at the trap center. The effects of gravity (along $`\widehat{z}`$) are presumed negligible. The constant rotation at frequency $`\mathrm{\Omega }`$ induces angular momentum per particle given by the expectation value of $`L_z=i\left(y_xx_y\right)`$. The trapping frequencies are $`(\omega _x,\omega _y,\omega _z)=\omega _x(1,\alpha ,\beta )`$ with $`\omega _x=2\pi \times 132`$ rad/s, $`\alpha 1`$, and $`\beta =\sqrt{8}`$ . Choosing the condensate to be normalized to unity yields the scaling parameter $`\eta =N_0a/d_x`$. Note that energy, length, and time are given throughout in scaled harmonic oscillator units $`\mathrm{}\omega _x`$, $`d_x=\sqrt{\mathrm{}/M\omega _x}0.94\mu `$m, and $`\mathrm{T}=2\pi /\omega _x7.6`$ ms, respectively.

The ground-state of the GP equation is found within a discrete-variable representation (DVR)  by imaginary time propagation using an adaptive stepsize Runge-Kutta integrator. A total of between $`\mathrm{40\hspace{0.17em}000}`$ and $`\mathrm{130\hspace{0.17em}000}`$ DVR points of a Gauss-Hermite quadrature are used, and all calculations are performed on a standard workstation. The stationary ground state in the rotating frame is found by setting $`\stackrel{~}{\tau }i\tau `$ and solving the diffusion equation:

$$_{\stackrel{~}{\tau }}\psi (𝐫,\stackrel{~}{\tau })=\left(H\mu \right)\psi (𝐫,\stackrel{~}{\tau }),$$
(2)   
where $`H`$ is the GP operator appearing on the right side of Eq. (1) and $`\mu `$ is the chemical potential. The condensate wavefunction is assumed to be even under inversion of $`z`$, and is initially taken to be the vortex-free Thomas-Fermi result, which is the time-independent solution of Eq. (2), neglecting the kinetic energy operator and $`L_z`$. A vortex is generated by imposing one quantum of circulation $`\kappa `$, a $`2\pi `$-winding of the phase around the $`z`$-axis, on the condensate wavefunction at $`\stackrel{~}{\tau }=0`$. At each imaginary timestep, the chemical potential $`\mu `$ is readjusted in order to preserve the norm of the wavefunction (i.e. the value of $`N_0`$). The propagation continues until the right side of Eq. (2) is equal to a tolerance $`\delta 10^{10}`$ defining the error in the dimensionless chemical potential. Stationary solutions are verified by subsequently integrating in real time; any deviations from self-consistency would be made manifest by collective motion.

The solution to the GP equation for a vortex located at the center of a nonrotating anisotropic trap containing $`N_0=10^5`$ atoms is shown in Fig. 1. In general, the vortex core is found to become decreasingly anisotropic as $`N_0`$ increases. Furthermore, the condensate density preserves its overall vortex-free shape far from the origin except for a slight overall bulge in order to preserve the norm. The structure of the vortex indicates that the healing length $`\xi `$ is governed largely by the local density and is only weakly dependent on trap geometry. In the Thomas-Fermi (TF) approximation, which is valid for large $`N_0`$, the healing length scales with the TF $`\widehat{x}`$-axis radius $`R=(15\alpha \beta \eta )^{1/5}d_x`$ as $`\xi (d_x/R)d_x`$. Indeed, the numerics clearly indicate that the mean vortex core radius (approximately $`d_x`$ at low densities) shrinks very slowly as the TF limit is approached.

A superfluid subjected to a torque will remain purely irrotational until the critical frequency $`\mathrm{\Omega }_c`$ is reached, at which point it becomes globally favorable for the system to contain a vortex with a single quantum of circulation $`\kappa `$. In cylindrically-symmetric systems where the Hamiltonian commutes with $`L_z`$, the circulation and angular momentum (with quantum $`m`$) are identical; in the rotating frame, the free energies of the $`m0`$-states are shifted by $`m\mathrm{\Omega }`$, and $`\mathrm{\Omega }_c`$ is simply the difference in energy between the $`m=1`$ and $`m=0`$ states (divided by $`\mathrm{}`$). In fully anisotropic traps, however, even the $`\kappa =0`$ state is shifted, so the applied $`\mathrm{\Omega }`$ in Eq. (1) must be increased until the free energy curves cross. It is straightforward to extend the TF estimate of $`\mathrm{\Omega }_c`$  to include a small deviation from cylindrical symmetry ; neglecting the shift of the vortex-free chemical potential (valid for $`\alpha 1`$), one obtains:

$$\mathrm{\Omega }_c\frac{5\alpha }{2}\left(\frac{d_x^2}{R^2}\right)\mathrm{ln}\left(\frac{R}{\xi }\right)\omega _x.$$
(3)   
Fig. 2 shows the critical frequencies for the global stability of a vortex with $`\kappa =1`$ at the trap center. For all geometries, the critical frequency drops monotonically as $`N_0`$ is increased. For a given number of atoms, the value of $`\mathrm{\Omega }_c`$ increases with in-plane anisotropy, similar to the behavior found for liquid helium in rotating elliptical containers . The energy of vortex formation must exceed that of the irrotational velocity field, which is finite for a vortex-free condensate in a rotating anisotropic trap. The TF result (3) agrees well with the numerical data in its regime of validity $`\alpha 1`$, though it tends to slightly overestimate the value of $`\mathrm{\Omega }_c`$.

While $`\mathrm{\Omega }_c`$ provides the criterion for the global stability of a vortex, it does not necessarily indicate the critical frequency for vortex nucleation. When initially vortex-free condensates are placed in anisotropic traps rotating at a frequency $`\mathrm{\Omega }<\omega _x`$, the velocity field of the stationary solution is found to be irrotational even for $`\mathrm{\Omega }\mathrm{\Omega }_c`$. In a harmonic trap with smooth edges, it is not clear if there exists any suitable locus for vortex formation. The vortices are most likely to originate at the condensate surfaces normal to the axis of weak confinement, where the local critical velocity is small  but the tangential superfluid velocity in the laboratory frame is largest . While these issues are beyond the scope of the present issue of vortex stability, there is evidence that multiple vortices appear at higher frequencies . For smaller $`N_0`$, it would likely be easier to generate a vortex experimentally by rotating the anisotropic trap before the condensate is cooled below the BEC transition .

When $`\alpha >1`$, the angular momentum per particle $`l_\kappa `$ is a nontrivial function of $`N_0`$, $`\alpha `$, and $`\mathrm{\Omega }`$. In a nonrotating system with unit vorticity, $`l_\kappa `$ increases with $`N_0`$. In the absence of a vortex, $`l_\kappa `$ is finite for a given $`\mathrm{\Omega }`$, and increases with $`\alpha `$; the superfluid velocity $`𝐯_s`$ can be locally appreciable but still remain irrotational $`\times 𝐯_s=0`$. At the critical frequency, the difference between $`l_1`$ and $`l_0`$ is always less than unity; for the most extreme case considered here, a system with $`N_0=10^6`$ and $`\alpha =3`$ rotating at $`\mathrm{\Omega }_c=0.14\omega _x`$, one obtains $`l_1=2.63\mathrm{}`$ and $`l_0=1.77\mathrm{}`$. As $`\alpha \mathrm{}`$, the angular momentum approaches that of a non-superfluid TF cloud $`l_0I_{\mathrm{sb}}\mathrm{\Omega }`$ with ‘solid-body’ moment of inertia $`I_{\mathrm{sb}}=\frac{1}{/}7MR^2`$.

An anisotropic harmonic oscillator potential becomes unconfining when it is rotated at a frequency between the smallest and largest trapping frequencies. Since $`\mathrm{\Omega }_c`$ exceeds $`\omega _x`$ for sufficiently large $`\alpha `$, there exists a critical minimum number of condensed atoms $`N_c`$ able to support a vortex. The value of $`N_c`$ increases with $`\alpha `$ and is given by the intercept of the $`\mathrm{\Omega }/\omega _x=1`$ line in Fig. 2. In cylindrically-symmetric systems $`N_c=1`$, since in the rotating frame the free energies for all the $`m`$-states become degenerate at $`\mathrm{\Omega }=\omega _x`$. In the limit of extreme anisotropy $`\alpha \mathrm{}`$ vortices can never be stabilized.

It should be noted that states with vortices at the center of anisotropic harmonic traps are found to be stationary solutions of Eq. (2) for all values of $`N_0`$ and $`\mathrm{\Omega }0`$ considered; such configurations do not appear to decay in either real or imaginary time. A vortex at the center of a nonrotating cylindrical trap increases the system’s free energy, but is stationary because both the vorticity and angular momentum commute with the Hamiltonian; in principle, the angular momentum can be eliminated, and the free energy reduced, only if this symmetry is broken by displacing the vortex from the center. Since angular momentum is not conserved in anisotropic traps, the apparent vortex stability is likely due to the free energy maximum at the trap center . In the absence of an external pinning mechanism, any such configuration should be unstable against infinitesimal displacements.

In order to further explore the issue of vortex stability in nonrotating traps, the initial condensate phase is wound by $`2\pi `$ a small distance $`x_00.2d_x`$ from the origin of a trap with $`\alpha =1`$. For all values of $`N_010^6`$, the condensate wavefunction rapidly (by $`\stackrel{~}{\tau }`$ T) converges to a metastable solution with a vortex, where the fluctuations in $`\mu `$ become smaller than $`\delta 10^7`$ per timestep $`\mathrm{\Delta }\stackrel{~}{\tau }10^3`$T. This wavefunction subsequently decays to the true ground state, but both the real and imaginary time required to do so is found to increase with $`N_0`$ . To an excellent approximation, the total time diverges like $`\stackrel{~}{\tau }N_0^{2/5}`$T; for $`N_0\genfrac{}{}{0pt}{}{>}{}10^5`$, the time required $`30`$T becomes computationally inaccessible and the vortex state becomes numerically indistinguishable from stationary. The numerics suggest that while vortices in nonrotating traps are always unstable against off-center displacements, they may be very long-lived.

The observed $`x_0>0`$ instability of the vortex state is likely due to the existence of an ‘anomalous’ collective mode $`\omega _a`$ at low densities . This dipole mode, which has positive norm but negative energy (or vice versa), is associated with a zero angular momentum bound state in the vortex core ; its value corresponds to the precession frequency of the vortex relative to the cloud . Previous numerical calculations  found $`|\omega _a|>0`$ for all $`N_010^4`$. As the core radius shrinks with larger $`N_0`$, however, the anomalous energy might be pushed to zero, yielding long-lived or even stable vortices in the TF limit.

The low-lying excitation frequencies of a nonrotating condensate in the vortex state are calculated using the Bogoliubov equations . For completely anisotropic geometries, however, the Bogoliubov operator is too large to diagonalize explicitly. Calculations are therefore restricted to the cylindrical case $`\alpha =1`$, where the vortex condensate is $`\psi \psi _1(\rho ,z)e^{i\varphi }`$ and the quasiparticle amplitudes $`u`$ and $`v`$ are labeled by $`m`$, the projection of the angular momentum operator $`L_z`$. The Bogoliubov equations are then

$$\left(\begin{array}{cc}\widehat{O}& V_\mathrm{H}\\ V_\mathrm{H}& \widehat{O}^{}\end{array}\right)\left(\begin{array}{c}u_m\\ v_{2m}\end{array}\right)=ϵ_m\left(\begin{array}{c}u_m\\ v_{2m}\end{array}\right),$$
(4)   
where $`\widehat{O}\frac{1}{/}2\rho \frac{}{/}\rho \rho \frac{}{/}\rho \frac{1}{/}2\frac{^2}{/}z^2+\frac{m^2}{/}2\rho ^2+V_\mathrm{t}+2V_\mathrm{H}`$, $`\widehat{O}^{}\widehat{O}+\frac{2(1m)}{/}\rho ^2`$, and $`u_1=v_1=\psi _1`$ when $`ϵ_1=0`$. In the $`\widehat{\rho }`$-direction, the points of the DVR grid correspond to those of Gauss-Laguerre quadrature, and the kinetic energy matrix elements are obtained using the prescription of Baye and Heenen .

The anomalous mode $`\omega _a`$, which is labeled by $`m=2`$, is shown as a function of $`N_0`$ in Fig. 2. The results indicate that $`0<|\omega _a|\mathrm{\Omega }_c`$ for all $`N_010^6`$ considered. Our calculations suggest that for large numbers of atoms, $`\omega _a`$ coincides with the metastability rotation frequency $`\mathrm{\Omega }_0`$ discussed above; the numerical value of $`\omega _a`$ is consistent with the TF result $`\mathrm{\Omega }_0=\frac{3}{/}5\mathrm{\Omega }_c`$ . Indeed, in the frame of a condensate rotating at $`\mathrm{\Omega }=\omega _a`$, the frequency of the vortex oscillation would be Doppler shifted to zero. Alternatively, it can be shown in both the weakly-interacting and TF limits that $`\mathrm{\Omega }_0`$ is also the frequency at which the chemical potentials $`\mu `$ for the vortex and vortex-free states become equal; in the TF limit, $`\omega _a`$ vanishes when the vacua (or energy zero) for quasiparticle excitations for both states coincide.

Since the anomalous mode corresponds to the precession of the vortex about the trap origin, one may make a crude estimate of the vortex lifetime $`\tau `$. In the presence of dissipation, the vortex will spiral out of the condensate after a few orbit periods $`\omega _a^1`$. Assuming that $`\omega _a=\frac{5}{/}3\mathrm{\Omega }_c`$, then with Eq. (3) one obtains $`\tau N_0^{2/5}`$T in the TF limit neglecting logarithmic factors. This result is consistent with the imaginary time $`\stackrel{~}{\tau }`$ required to yield the vortex-free ground state in the fully three-dimensional numerical calculations discussed above. Similar decay times have been obtained for solitons and vortices in the presence of a small noncondensate component .

In summary, we have obtained numerically the critical frequencies $`\mathrm{\Omega }_c`$ for the stabilization of a vortex at the center of a rotating anisotropically-trapped Bose condensate. Since $`\mathrm{\Omega }_c`$ increases with the in-plane anisotropy $`\alpha =\omega _y/\omega _x`$ and the condensate becomes unconfined for $`\mathrm{\Omega }>\omega _x`$, there is a minimum number of atoms able to support a vortex state. Vortices in nonrotating traps are found to be unstable against small off-center displacements, but their decay time diverges with the total number of atoms.

###### Acknowledgements.

The authors are grateful for many stimulating discussions with M. Brewczyk, J. Denschlag, M. Edwards, A. L. Fetter, D. Guery-Odelin, D. A. W. Hutchinson, M. Matthews, W. D. Phillips, W. Reinhardt, and C. J. Williams. This work was supported by the U.S. Office of Naval Research.