# Nucleation of vortex arrays in rotating anisotropic Bose-Einstein condensates

\[

## Abstract

The nucleation of vortices and the resulting structures of vortex arrays in dilute, trapped, zero-temperature Bose-Einstein condensates are investigated numerically. Vortices are generated by rotating a three-dimensional, anisotropic harmonic atom trap. The condensate ground state is obtained by propagating the Gross-Pitaevskii equation in imaginary time. Vortices first appear at a rotation frequency significantly larger than the critical frequency for vortex stabilization. This is consistent with a critical velocity mechanism for vortex nucleation. At higher frequencies, the structures of the vortex arrays are strongly influenced by trap geometry.

\]

Since the experimental achievement of Bose-Einstein condensation (BEC) in confined alkali gases , the possibility of generating vortices in confined weakly-interacting dilute Bose gases has been intensively studied . While theoretical investigations of stability have generally been restricted to the case of a single vortex , the proposed experimental techniques may induce several vortices simultaneously . Under appropriate stabilizing conditions, such as a continuously applied torque, these vortices would form an array akin to those obtained in rotating superfluid helium .

A standard approach used to ‘spin-up’ superfluid helium is to rotate the container at an angular frequency $`\mathrm{\Omega }`$. Aside from significant hysteresis effects , vortices tend to first appear at a frequency $`\mathrm{\Omega }_\nu `$, whose value is comparable to the critical frequency $`\mathrm{\Omega }_c`$ at which the presence of vortices lowers the free energy of the interacting system . Energy minimization arguments have also yielded vortex arrays that are very similar to those observed experimentally . Despite these successes, the mechanisms for vortex nucleation by rotation remain poorly understood; important factors are thought to include presence of a normal fluid, impurities, and surface roughness .

It has been suggested that vortices may be similarly generated in the dilute Bose gases by rotating the trap about its center . Evidently, a harmonic potential can transfer angular momentum to the gas only if it is anisotropic in the plane of rotation. While vortices in such a system at zero temperature have been shown to become energetically stable for $`\mathrm{\Omega }>\mathrm{\Omega }_c`$ , the particle flow could remain irrotational at these angular frequencies since there exists an energy barrier to vortex formation . Suppression of this barrier could be induced by application of a perturbing potential near the edge of the confined gas, as has been simulated in the low-density limit . One of the primary motivations for the present work, however, is to determine if there exists any intrinsic mechanism for vortex nucleation in a dilute quantum fluid that is free of impurities, surface effects, and thermal atoms. We find that vortices can indeed be generated by rotating Bose condensates confined in an anisotropic harmonic trap. The value of $`\mathrm{\Omega }_\nu `$ at which vortices are spontaneously nucleated is somewhat larger than $`\mathrm{\Omega }_c`$. For $`\mathrm{\Omega }>\mathrm{\Omega }_\nu `$ multiple vortices appear simultaneously, in patterns that depend upon the geometry of the trap.

The dynamics of a dilute Bose condensate at zero temperature are governed by the time-dependent Gross-Pitaevskii (GP) equation . Previous simulations of the GP equation have demonstrated that vortex-antivortex pairs or vortex half-rings can be generated by superflow around a stationary obstacle  or through a small aperture . In , the vortex pairs form when the magnitude of the superfluid velocity exceeds a critical value which is proportional to the local sound velocity; recent experimental results support this conclusion . To our knowledge, no numerical investigation of vortex nucleation in three-dimensional inhomogeneous rotating superfluids has hitherto been attempted.

The numerical calculations presented here model the experimental apparatus of Kozuma et al. , where <sup>23</sup>Na atoms are confined in a completely anisotropic three-dimensional harmonic oscillator potential. In the presence of a constant external torque, the condensate obeys the time-dependent GP equation in the rotating reference frame:

$$i_\tau \psi (𝐫,\tau )=\left[T+V_{\mathrm{trap}}+V_\mathrm{H}\mathrm{\Omega }L_z\right]\psi (𝐫,\tau ),$$
(1)   
where the kinetic energy is $`T=\frac{1}{/}2\stackrel{}{}^2`$, the trap potential is $`V_{\mathrm{trap}}=\frac{1}{/}2\left(x^2+\alpha ^2y^2+\beta ^2z^2\right)`$, and the Hartree term is $`V_\mathrm{H}=4\pi \eta |\psi |^2`$. The angular momentum operator $`L_z=i\left(y_xx_y\right)`$ rotates the system about the $`z`$-axis at the trap center at a constant angular frequency $`\mathrm{\Omega }`$. The trapping frequencies are $`(\omega _x,\omega _y,\omega _z)=\omega _x(1,\alpha ,\beta )`$ with $`\omega _x=2\pi \times 26.87`$ rad/s, $`\alpha =\sqrt{2}`$, and $`\beta =1/\sqrt{2}`$. Normalizing the condensate $`𝑑𝐫|\psi (𝐫,\tau )|^2=1`$ yields the scaling parameter $`\eta =N_0a/d_x`$, where $`a=2.75`$ nm is the s-wave scattering length for Na and $`N_0`$ is the number of condensate atoms. Unless explicitly written, energy, length, and time are given throughout in scaled harmonic oscillator units $`\mathrm{}\omega _x`$, $`d_x=\sqrt{\mathrm{}/M\omega _x}4.0\mu `$m, and $`\mathrm{T}=\omega _x^16`$ ms, respectively.

The stationary ground-state solution of the GP equation, defined as that which minimizes the value of the chemical potential, is found by norm-preserving imaginary time propagation (the method of steepest descents) using an adaptive stepsize Runge-Kutta integrator. The complex condensate wavefunction is expressed within a discrete-variable representation (DVR)  based on Gauss-Hermite quadrature, and is assumed to be even under inversion of $`z`$. The numerical techniques are described in greater detail elsewhere . The initial state (at zero imaginary time $`\stackrel{~}{\tau }i\tau =0`$) is taken to be the vortex-free Thomas-Fermi (TF) wavefunction $`\psi _{\mathrm{TF}}=\sqrt{(\mu _{\mathrm{TF}}V_{\mathrm{trap}})/4\pi \eta }`$, which is the time-independent solution of Eq. (1), neglecting $`T`$ and $`L_z`$, with chemical potential $`\mu _{\mathrm{TF}}=\frac{1}{/}2(15\alpha \beta \eta )^{2/5}`$. The GP equation for a given value of $`\mathrm{\Omega }`$ and $`N_0`$ is propagated in imaginary time until the fluctuations in both the chemical potential and the norm become smaller than $`10^{11}`$. It should be emphasized that the equilibrium configuration is found not to depend on the choice of purely real initial state. Since the final state is unconstrained except for $`z`$-parity, the lowest-lying eigenfunction of the GP equation corresponds to a local minimum of the free energy functional.

In Fig. 1 are depicted the condensate density, which is stationary in the rotating frame, as well as the condensate phase and the velocity field in the laboratory and rotating frames, for $`\mathrm{\Omega }=0.45\omega _x`$ and $`N_0=10^6`$. The density profile at this angular frequency contains no vortices but is slightly extended from that of a non-rotating condensate due to the centrifugal forces. The velocity field in the laboratory frame is given by $`𝐯_s^l\stackrel{}{}\phi `$ in units of $`\omega _xd_x`$, where $`\phi `$ is the condensate phase. In the rotating frame, $`𝐯_s^r=𝐯_s^l\mathrm{\Omega }\widehat{z}\times 𝐫`$. There are no closed velocity streamlines found in Fig. 1(a). Such an irrotational flow $`\stackrel{}{}\times 𝐯_s=0`$ is characteristic of a superfluid, distinct from the related properties of vortex quantization and stability. The only solution of the GP equation satisfying irrotational flow in a cylindrically-symmetric trap is $`𝐯_s=0`$: rotating the trap is equivalent to doing nothing. The irrotational velocity field for an anisotropic trap is nontrivial, however. Since the density profile is independent of orientation, mass flow must accompany the rotation even though the superfluid prefers to remain at rest .

The condensate is found to remain vortex-free for angular velocities significantly larger than the expected critical frequency for the stability of a single vortex $`\mathrm{\Omega }_c^{(1)}`$ . In order to determine if irrotational configurations correspond to the global free energy minima of the system, vortex states are investigated by artificially imposing total circulation $`n\kappa `$ on the condensate wavefunction. By winding the phase at $`\stackrel{~}{\tau }=0`$ by $`2\pi n`$ about the trap center, imaginary-time propagation of the GP equation yields the minimum energy configuration with $`n`$ vortices if such a solution is stationary or metastable . The results for $`N_0=10^6`$ and $`\mathrm{\Omega }=0.45\omega _x`$ are summarized in Table I. At this angular frequency, states with $`n=1,2,3`$ are all energetically favored over the vortex-free solution. The vortices in these cases are predominantly oriented along the ($`\widehat{z}`$) axis of rotation, and are located symmetrically about the origin on the (loose) $`x`$-axis. The frequency chosen is too low to support the four vortex case $`\mathrm{\Omega }<\mathrm{\Omega }_c^{(4)}`$, but is larger than the frequency (which may correspond to metastability) at which the chemical potentials for $`n=3`$ and $`n=4`$ cross.

As shown in Fig. 2, vortices with the same circulation $`\kappa `$ (as opposed to vortex-antivortex pairs) begin to penetrate the cloud above a critical angular velocity for vortex nucleation $`\mathrm{\Omega }_\nu `$. The value of $`\mathrm{\Omega }_\nu `$ is found to not depend strongly on trap geometry and to decrease very slowly with $`N_0`$; for $`N_0=10^q`$ with $`q=\{5,6,7\}`$, we obtain $`\mathrm{\Omega }_\nu =\{0.65,\mathrm{\hspace{0.17em}0.50},\mathrm{\hspace{0.17em}0.36}\}\omega _x\pm 0.01\omega _x`$, respectively. In contrast, the critical frequency for the stabilization of a single vortex in an anisotropic trap is approximately given by the TF expression $`\mathrm{\Omega }_c^{(1)}(5\alpha /2R^2)\mathrm{ln}(R/\xi )\omega _x`$, where $`R=\sqrt{2\mu _{\mathrm{TF}}}`$ and $`\xi =\sqrt{\alpha }/R`$ are the dimensionless condensate radius along $`\widehat{x}`$ and the healing length, respectively . For the parameters considered here, the values are predicted to be $`\mathrm{\Omega }_c^{(1)}=\{0.61,\mathrm{\hspace{0.17em}0.33},\mathrm{\hspace{0.17em}0.16}\}\omega _x`$ and are found numerically to be $`\mathrm{\Omega }_c^{(1)}=\{0.54,\mathrm{\hspace{0.17em}0.29},\mathrm{\hspace{0.17em}0.14}\}\omega _x`$. The number of vortices $`n_\nu `$ present just above $`\mathrm{\Omega }_\nu `$ is found to increase with $`N_0`$; $`n_\nu =4`$ and $`8`$ for $`N_0=10^6`$ and $`10^7`$, respectively. The value of $`\mathrm{\Omega }_\nu `$ may be interpreted as the critical frequency $`\mathrm{\Omega }_c^{(n_\nu )}`$ for the stabilization of $`n_\nu `$ vortices. If $`n_\nu =n`$ for all $`N_0`$, then $`\mathrm{\Omega }_\nu \mathrm{\Omega }_c^{(n)}N_0^{2/5}`$. That $`\mathrm{\Omega }_\nu `$ decreases more slowly with $`N_0`$ implies that $`n_\nu `$ must increase with $`N_0`$. The small difference between $`\mathrm{\Omega }_c^{(1)}`$ and $`\mathrm{\Omega }_\nu `$ for $`N_0=10^5`$ reflects the instability of vortex arrays in the low-density limit. As $`N_0`$ decreases, the spacing between successive $`\mathrm{\Omega }_c^{(n)}`$ diminishes, and vanishes for $`N_0=0`$ in cylindrically-symmetric traps; for very large $`N_0`$, the spacing approaches a constant as the vortex-vortex interactions become negligible.

The numerical results for $`\mathrm{\Omega }_\nu `$ suggest that the criteria for vortex stabilization and nucleation are different. Superflow through microapertures, or the motion of an object or ion through a superfluid, can give rise to vortex half-ring  or vortex-pair  production through the accumulation of phase-slip. One might expect similar excitations in a rotating condensate : vortex half-rings would be nucleated at the condensate surface when the local tangential velocity exceeds a critical value. Indeed, the distinction between a half-ring and vortex becomes blurred in a trapped gas with curved surfaces, as discussed further below.

A crude estimate of $`\mathrm{\Omega }_\nu `$ may be obtained by invoking the Landau criterion for the critical velocity $`v_{\mathrm{cr}}=\mathrm{min}(\omega _q/q)`$, where $`\omega _q`$ is the frequency of the mode at wavevector $`q`$. Such a minimum corresponds to values of $`q_c`$ at which the hydrodynamic description of the collective excitations begins to fail . For a spherical trapped Bose gas, the crossover to a single-particle behavior occurs in a boundary-layer region at the cloud surface whose thickness is several $`\delta =(2R)^{1/3}d_x`$ . Minimizing $`\omega _q/q`$ using the dispersion relation for the planar surface modes  of such a system $`\omega _q^2\omega _x^2[qRd_x^4q^4\left(\mathrm{ln}(q\delta )0.15\right)]`$, one obtains $`q_c=(R/0.3)^{1/3}d_x^1\delta ^1`$ and $`\mathrm{\Omega }_\nu =v_{\mathrm{cr}}/RR^{2/3}`$. Since $`RN_0^{1/5}`$, the critical frequency $`\mathrm{\Omega }_\nu N_0^{2/15}`$ decreases far more slowly than does the TF estimate for $`\mathrm{\Omega }_cN_0^{2/5}`$. The number-dependence of $`\mathrm{\Omega }_\nu `$ is in reasonable agreement with the numerical data. Real-time simulations further confirm that high-frequency oscillations of the condensate are required for vortex production at the same $`\mathrm{\Omega }_\nu `$ found using the imaginary-time approach.

The above analysis does not clearly identify the instability of the surface modes with the penetration of vortices into the cloud, however. Further insight may be gained by considering the free energy $`F`$ of a single vortex in a cylindrical trap, relative to that of the vortex-free state, as a function of the vortex displacement $`\rho `$ from the trap center . In the TF limit, $`F`$ vanishes for $`\rho ^2=R^2`$ and $`\rho ^2=R^2(5/2\mathrm{\Omega })\mathrm{ln}(R/\xi )`$, corresponding to the right and left roots of the free energy barrier to vortex generation, respectively. As $`\mathrm{\Omega }`$ increases, the energy barrier at the surface narrows but remains finite. Yet, as discussed above, the hydrodynamic excitations begin to break down at a radius $`\stackrel{~}{\rho }R\delta `$. Vortices will therefore spontaneously penetrate the cloud when the angular frequency exceeds

$$\mathrm{\Omega }_\nu =\frac{5\sqrt[3]{2}}{4R^{2/3}}\mathrm{ln}\left(\frac{R}{\xi }\right)\omega _x,$$
(2)   
since the barrier effectively disappears when $`\stackrel{~}{\rho }\rho `$. Thus, the frequencies of nucleation and penetration have the same number-dependence and are defined by a single critical wavelength. Once the condensate contains vortices at a given $`\mathrm{\Omega }>\mathrm{\Omega }_\nu `$, the functional $`F`$ will again include a barrier to vortex penetration from the surface, reflecting the hydrodynamic stability of the vortex state. One may thus envisage a succession of multiple-vortex nucleation events at well-defined angular frequencies.

The stationary configurations of vortex arrays are shown as a function of applied rotation in Fig. 2. The condensate density is shown integrated down the axis of rotation $`\widehat{z}`$, in order to mimic an in situ image of the cloud. While the vortices near the origin appear to have virtually isotropic cores, those in the vicinity of the surface are generally wider and are noticeably distorted. The anisotropy is due in part to the divergence of the coherence length as the density decreases, but is mostly the result of vortex curvature. Off-center vortices are not fully aligned with the axis of rotation $`\widehat{z}`$, since they terminate at normals to the ellipsoidal condensate surface. Far from the origin, the vortex structure approaches that of a half-ring pinned to the condensate surface.

The symmetries of the confining potential impose constraints on the vortex arrays that may be produced by rotating anisotropic traps. Stationary configurations are found to always have the inversion symmetry $`(x,y,z)(x,y,z)`$. As shown in Fig. 2, the number of vortices is at least four and is even for each array; real-time simulations demonstrate that vortices with the same circulation are nucleated in pairs at inversion-related points on the surface. No vortex is found at the origin, since the odd number of remaining vortices cannot be distributed symmetrically. At low angular velocities, therefore, the array tends to approximate a regular tetragonal lattice. As the total number of vortices increases with $`\mathrm{\Omega }`$, however, a different pattern begins to emerge. While a triangular array is inconsistent with the twofold trap symmetries, it is more efficient for close packing; this geometry is favored for vortices near the rotation axis of rapidly rotating vessels of superfluid helium . If vortices in trapped condensates could be made sufficiently numerous, they would likely form a near-regular triangular array but with the central vortex absent.

In summary, the critical frequencies for the zero-temperature nucleation of vortices $`\mathrm{\Omega }_\nu `$ in rotating anisotropic traps are obtained numerically, and are found to be larger than the vortex stability frequencies $`\mathrm{\Omega }_c`$. The number-dependence of $`\mathrm{\Omega }_\nu `$ is consistent with a critical-velocity mechanism for vortex production. The structures of vortex arrays are strongly affected by trap geometry, but approach triangular at large densities.

###### Acknowledgements.

The authors are grateful to A. L. Fetter, R. L. Pego, S. L. Rolston, J. Simsarian, and S. Stringari for numerous fruitful discussions, and to P. Ketcham for assistance in generating the figures. This work was supported by the U.S. office of Naval Research.