# Suppression of Kelvon-induced decay of quantized vortices in oblate Bose-Einstein Condensates

## Abstract

We study the Kelvin mode excitations on a vortex line in a three-dimensional trapped Bose-Einstein condensate at finite temperature. Our stochastic Gross-Pitaevskii simulations show that the activation of these modes can be suppressed by tightening the confinement along the direction of the vortex line, leading to a strong suppression in the vortex decay rate as the system enters a regime of two-dimensional vortex dynamics. As the system approaches the condensation transition temperature we find that the vortex decay rate is strongly sensitive to dimensionality and temperature, observing a large enhancement for quasi-two-dimensional traps. Three-dimensional simulations of the recent vortex dipole decay experiment of Neely et al. [Phys. Rev. Lett. 104, 160401 (2010)] confirm two-dimensional vortex dynamics, and predict a dipole lifetime consistent with experimental observations and suppression of Kelvon-induced vortex decay in highly oblate condensates.

## I Introduction

Quantum vortices in Bose-Einstein condensates (BECs) are topologically stable excitations of the superfluid matter wave field that carry angular momentum (1); ?. Unless the system is rotating sufficiently fast, vortices are thermodynamically unstable (3). At finite temperatures the thermal cloud can provide a source of dissipation and noise, and hence allow vortex decay. However, there is currently little understanding of the decay mechanism and its timescale. The dynamics of vortices in such systems is therefore of interest as a robust test of non-equilibrium field theory. A number of different methods have been developed to describe the dynamics of vortices in finite-temperature systems. Single vortex decay has been simulated using the Zaremba-Nikuni-Griffin formalism (4), the classical field method (5), and within a Gross-Pitaevskii theory with phenomenological damping (6). The nature of the thermal instability has also been studied within a variational approach (7) to the stochastic Gross-Pitaevskii equation (8); ?.

Classical vortices can support long-wavelength helical traveling waves known as Kelvin waves (10). Superfluid Kelvin waves, or Kelvons (11), have long be studied in Helium (12), and more recently have become a topic of interest in BECs. The Kelvon excitation mechanism has been theoretically investigated (13), along with the microscopic dynamics of coherently excited Kelvin waves (14); ?, with results consistent with the experimental evidence for Kelvons (16). Similar bending waves have also been shown to play an important role in the instability of vortex rings in trapped BECs (17).

While the coherent excitation of Kelvin modes in a BEC is well understood, their role in the finite-temperature system is largely unexplored. Here we show that the thermal activation of Kelvin waves is a dominant factor in the decay of a vortex for a finite-temperature three-dimensional (3D) BEC. The underlying reason is that Kelvin waves cause the vortex to wobble and emit acoustic radiation, thus allowing the vortex to more effectively dissipate energy, and hence decay by moving out to the boundary of the BEC. As Kelvin waves are fundamentally 3D excitations of the vortex line, they may be suppressed by flattening the system (see Fig. 1).

Indeed, beyond a certain oblateness we find that the Kelvon mode contribution is essentially frozen out and the vortex decay rate reduces to a geometrically-invariant value, characterizing a regime of two-dimensional vortex (2DV) dynamics. We also find that for tighter traps and higher temperatures the matter wave itself crosses over to being quasi-2D and we find a dynamical signature of the phase-fluctuating condensate: an anomalous increase in the rate of vortex decay.

As a further application of these concepts, we model a recent experiment (18) studying vortex dipoles in the 2DV regime, confirming the 2DV behavior in finite temperature dynamics. The decay characteristics are found to be qualitatively very different from 3D dipole decay which exhibits vortex ring and loop formation, accelerating the damping process.

Our theoretical framework for studying this system is the stochastic projected Gross-Pitaevskii equation (SPGPE) (19); ?, which is a grand-canonical c-field method (21) (also see (22); ?; ?; (25) for related theories). This formalism has been used to study vortex lattice formation from a rotating thermal cloud (26), a quantitative model of spontaneous vortex formation in experiments (26); (27), and the decay of a single vortex (28) (also see (29); ?; ?).

## Ii System and equation of motion

### ii.1 Physical system

We begin by describing the system of bosonic atoms in the cold-collision regime with the second-quantized many body Hamiltonian

(1) | |||||

where the single-particle Hamiltonian is

(2) |

the cylindrically symmetric harmonic trap is

(3) |

and , with the atomic mass and the s-wave scattering length. The field operators obey Bose commutation relations

(4) |

We define as the degree of trap oblateness, with enforcing a constant geometric mean frequency for arbitrary oblateness. This choice ensures the ideal gas transition temperature is invariant under scaling of at fixed total atom number .

### ii.2 Stochastic PGPE Formalism

In this work we use the SPGPE formalism to describe the process of vortex decay in a finite temperature BEC. The approach we follow in choosing parameters and preparing initial states for the SPGPE is extensively detailed in Ref. (28), hence we only briefly summarize it here.

The SPGPE is a c-field method where the system is decomposed in terms of the single-particle modes of the system satisfying , where represents all quantum numbers required to specify the modes. The system is divided into the coherent region (C) consisting of all states with , and the incoherent region (I) containing the remaining high energy modes (21). This division is made using a projection operator, which defines the C-region field operator as

(5) |

Equivalently, we have the expansion

(6) |

where . The I-region acts as a thermal reservoir for the C-region, and is assumed to be in thermal equilibrium at a temperature and chemical potential . The C-region dynamics are described by a stochastic differential equation for the c-field corresponding to (6):

(7) |

where are ordinary c-numbers. The equation of motion is known as the simple-growth SPGPE (26)

(8) |

where is the Hamiltonian evolution operator for the C region , and the complex Gaussian noise, , satisfies . For a full derivation of the SPGPE we refer the reader to (19); ?.

The first term of the SPGPE (8) describes Gross-Pitaevskii evolution over the C region, while the second and third terms account for scattering of two I-region atoms, resulting in growth of the C-region, and the corresponding time-reversed process. The growth rate, , can be calculated in near-equilibrium systems (26) and is independent of position over a large region of the trap with the value

(9) |

where and is the Lerch transcendent.

To obtain initial states for evolution according to (8), we follow the procedure used in (28), which we briefly summarize. We generate (non-vortex) finite-temperature equilibrium states by evolving the SPGPE with an appropriate choice of , , and ^{1}^{2}

## Iii Single Vortex Decay

In our simulations we use initial vortex states prepared for a system containing Rb atoms in a trap of (fixed) geometric mean trap frequency [for which nK], and a range of geometries from spherical to highly oblate. For each parameter set we calculate vortex decay properties by evolving 10-30 trajectories of the SPGPE ^{3}

Fig. 2 shows density isosurfaces of the vortex in a variety of condensate geometries after of evolution. The phase imprinting technique leads to a straight vortex line along the -axis, but under evolution the vortex rapidly thermalizes, typically within (see discussion below). We observe a significant amount of vortex bending in the spherical geometry and as the condensate becomes more two-dimensional, the magnitude of bending lessens and the vortex behaves essentially as a straight line. Although vortex bending on a short length scale is evident in spherical geometries, the vortex stays approximately aligned with the -axis throughout its evolution.

### iii.1 Kelvon power spectrum

To obtain the vortex line coordinates during evolution we locate the radial position of the phase singularity in each plane (i.e. value) and for all times . Fluctuations in the phase at the condensate boundary restrict the spatial range over which the vortex can be detected. This reduces the range of -values for which we can obtain the vortex-line coordinates, and, as time evolves and the vortex precesses toward the radial condensate boundary, this range of -values further reduces. For this reason we limit vortex detection to the range , where , with the Thomas-Fermi radius along the direction. We have found this to be the largest range which gives a well defined and unique curl signal over the vortex lifetime for our simulations.

We define to be the (complex) vortex line displacement from its instantaneous mean, where the barred quantities are the vortex positions averaged over . To quantify the extent of vortex bending at a given time, we take the Fourier transform with respect to (denoted by ) of , and calculate the power spectrum

(10) |

Note that this measure of the power spectrum is insensitive to varicose waves (15) and is thus entirely due to vortex bending.

Fig. 3 shows the power spectra vs Kelvon wavenumber for a range of geometries at s. We rescale the axis by , which is the wavenumber corresponding to the largest visible wavelength within our choice of window length. In all geometries we see the power spectra take a thermalized form, with maximum occupation of the long-wavelength bending modes (low values of ). The thermalization of the vortex line from the straight initial conditions is most directly monitored using the Kelvon power spectrum. We find that typically within ms the Kelvon spectrum reaches a stable state. This is much shorter than the vortex decay time (see Sec. III.2), so we can rule out the influence of any initial state artifacts on our predictions for the vortex decay process.

### iii.2 Effect of geometry: Kelvon power and vortex decay

The effect of geometry on the importance of the Kelvon modes is clearly visible comparing Figs. 3 (a)-(c): The total power in the spectrum decreases significantly as the system becomes more oblate.
It is useful to calculate this total power in the bending modes, , which by Parseval’s theorem is related to the line-averaged deviation of the vortex from its mean position ^{4}

Fig. 4(b) shows the effect of condensate geometry on the vortex lifetime, for a range of system temperatures. The lifetime is quantified here in terms of the *mean first exit time* (28), denoted by .
This is calculated as the ensemble average of the time it takes for the vortex to become indistinguishable from thermal fluctuations at the BEC edge (in the plane ). The vortex decay rate, , decreases with increasing , following the same trend as the total bending mode power [Fig. 4(a)]. However, rather than decreasing to zero, the decay rate approaches a constant value of order , being almost independent of for .
This constant is determined by purely 2D vortex dissipation processes which do not involve bending. Interestingly, for the highest temperatures considered, the decay rate starts to increase again at . This increase is not associated with an increase of bending mode power [Fig. 4(a)]. Aside from this anomaly, discussed below, we infer a critical oblateness of for the onset of 2DV dynamics in our system.

### iii.3 Effect of temperature: Anomalous decay in quasi-2D systems

The results of the previous section show that the vortex enters a 2D dynamics regime at , associated with strong suppression of bending modes, and a diminished vortex decay rate. However, for the highest temperatures the decay rates exhibit a revival as increases. We interpret this as arising from the system entering a phase fluctuating quasi-2D regime.

First, to more accurately characterize the critical temperature we include the downward shifts from due to interaction and finite-size effects (36). This gives a corrected critical temperature

(11) |

where is the first-order correction due to finite-size effects, and is the first order correction due to interactions, with and .

In Fig. 5 we show the decay rates versus verifying that these systems at the highest temperatures remain below the 3D critical temperature. However, for oblateness we see anomalously rapid increase in the vortex decay rate for temperatures exceeding . For these cases the system crosses over to being quasi-2D as the chemical potential becomes of the order of the energy to activate degrees of freedom in the tightly confined direction, i.e. (where is the single particle ground state energy). In particular, the highest temperature results for have (c.f. ). In the quasi-2D regime the dynamics of vortices are not well-understood, but strong phase fluctuations will likely facilitate the rate of vortex damping as we observe in our results ^{5}

## Iv Dipole decay

Having systematically studied the role of geometry and temperature in single vortex decay, we now consider two particular cases of the decay of a vortex dipole at finite temperature.

### iv.1 Spherical geometry: decay via ring and loop formation

We first consider dipole decay in a spherical trap. We use SPGPE simulation parameters , which give a total atom number and dimensionless damping rate (9) of . These parameters give a system with the same atom number and temperature as the oblate configuration considered in the experiment of Ref. (18) and in the next section, so the system serves as a benchmark for the effect of geometry on vortex dipole dynamics in this system. We create a finite-temperature vortex dipole state by imposing a vortex dipole phase pattern on a SPGPE equilibrium state via , where . This creates an axially aligned vortex dipole with positive and negative vortices located at , with . Figure 6 shows the decay sequence for the vortex dipole. The initially axially aligned vortices rapidly curve as the vortices precess, such that the vortex cores are normal to the surface. This long-wavelength bending leads to formation of a vortex ring, which then subsequently dissociates into two vortices at . As the vortices closely approach again at , both loops and rings form, reducing the total vortex length, and leading to total decay in less than 1s. We thus observe rapid vortex dipole decay via intermediate vortex ring and loop formation.

### iv.2 2DV regime: comparison of experiment and theory

In this section we use (8) to model the dynamics of a charge-one vortex dipole in a highly oblate finite-temperature Bose-Einstein condensate, as in the experiment by Neely *et al.* (18). The system consists of atoms at temperature nK, held in a harmonic trap with frequencies Hz. The trap aspect ratio exceeds identified in Sec. III.2 for closely related parameters. We can therefore expect that the experimental regime is also approximately that of 2DV. Experimentally, a single vortex dipole is created by pushing the BEC around an past an obstacle above a critical velocity. Ideally, we would like to compare the distance between the two vortices upon the creation of the dipole in both experiment and simulations, and to relate to the vortex dipole lifetime in each case. However, the experimental spatial resolution was , which is also the approximate vortex separation found immediately after dipole creation in Gross-Pitaevskii simulations of the expriment (18). Since the vortices are not individually resolvable immediately after creation, the experimental data can only place an upper bound on . We are thus unable to carry out quantitative comparison at this point. In this work we instead map out the dependence of dipole lifetime on by computing the decay dynamics for a range of initial vortex separations.

As in the spherical geometry of Sec. IV.1 we use SPGPE simulation parameters , which again gives a total atom number and dimensionless damping rate (9) of .
As in the spherical case, we phase imprint a vortex dipole into the BEC. We choose the initial separation between vortices as and . Note that the condensate radius is , so these separations are small relative to the system size, but larger than the healing length .
Qualitatively, our simulations differ from the experiment in several ways: in our simulations we have two vortices initially, a high degree of cylindrical and mirror symmetry, and a well defined total atom number. In the experiment there is occasionally a vortex present from the BEC formation process (27), an irreducible (albeit very small) cylindrical asymmetry of the harmonic potential, a mirror asymmetry of the stirring potential, and uncertainty in the total atom number. We do not attempt to model all of this complexity here, confining our study to the basic decay process of internal annihilation in the symmetric system, with fixed total atom number ^{6}

Figure 7 shows the mean vortex number over time from the experimental data (18), compared with the results of our simulations. In our simulations we always find that the vortices mutually annihilate each other in the center of the BEC, rather than damping at the boundary. The vortices move rapidly when closely approaching each other, nearing the speed of sound . Consequently, most of the orbital time is spent at the outer edge of the BEC, and the time interval in which annihilation may occur is a small fraction of the orbital period (). Thus the vortex number has an abrupt time dependence for a given initial separation, with characteristic timescale given by the orbital period of the vortex dipole, found numerically and experimentally to be and independent of over a wide range (18). In general, the value selects a given orbit of annihilation, within which there is very small variation of lifetime.

In detail, we find that for the mean lifetime is , corresponding to dipole orbits before annihilation. For example, with , the dipole self annihilates immediately after creation. The separations and lead to almost identical lifetimes of about , corresponding to the significant drop in vortex number seen experimentally for which we interpret as an experimental indication of the dipole lifetime.

In Figure 8 we plot a time series of simulation data for showing internal annihilation at the end of the second orbit. Interestingly, we observe a tilt in the axis of dipole precession caused by thermal fluctuations, which is also observed in experimental absorption images.

## V Conclusions and outlook

Here we summarize our results on single vortex and vortex dipole decay.

*Sub-critical Regime.—*For spherical systems at , a lone vortex is freely deformed from its initial straight-line configuration by the activation of vortex bending modes. This effect diminishes as the condensate is flattened, and bending modes become energetically inaccessible. The total bending mode power decreases exponentially as the condensate geometry becomes increasingly two-dimensional. For geometries with oblateness , the total power asymptotes to zero, indicating bending is strongly inhibited. The effective vortex decay rate also decreases as increases, approaching a constant two-dimensional value for . can therefore be identified as the critical oblateness required for the system to enter the 2D vortex regime.

*Critical Regime.—*The foregoing comments hold for temperatures outside an observed critical regime. When critical fluctuations associated with the phase transition become significant, we find that the vortex decay rate is highly sensitive to temperature and geometry. In particular, we see an anomalously large increase in the decay rate for geometries with , and that we attribute to phase-fluctuations associated with the system entering the quasi-2D regime.

*Vortex Dipole Decay.—*
We have modeled the experiment of Neely et al (18) which is in the 2DV regime (), and found that the approximate time for decay of a vortex dipole in the experiment (s) is consistent with numerical simulations with an initial vortex separation of m. In the simulations with these initial conditions mutual vortex annihilation occurs after two complete orbits, at s. In this work we have used as the fitting parameter; improved experimental resolution would be required for a quantitative comparison without any fitting parameters. In similar computations of spherical () BECs, the vortices decay much more quickly to rings and other bent geometries that would not be readily visible after only using the experimental methods of Ref. (18). We thus infer from the experiments of Ref. (18) and the calculations presented here that the experimental system consisting of a BEC with oblateness is within the 2D vortex dynamics regime, consistent with our calculations suggesting a critical oblateness of .

Acknowledgements: We thank Simon Gardiner for discussions regarding this work.

### Footnotes

- We require to be fixed and the highest energy mode in C to be appreciably occupied ( atoms).
- We do not find it necessary to imprint a vortex density profile because at the temperatures under consideration the extra energy of the core with phase-only imprinting is negligible compared with that from thermal fluctuations
- The number of trajectories required to obtain good statistics increases with the system temperature
- To obtain good statistics we integrate an instantaneous curve (as in Fig. 3), and then time-average the result for . This time is chosen to be much shorter than the vortex decay-time (which exceeds for all parameters), to give a short-time ensemble-average of the power data.
- The system still has a phase space density sufficiently high that spontaneous vortices are not expected to be significant (e.g. see (39); ?; ?).
- Our simulation initial conditions correspond to a grand canonical ensemble for the C-field region, with well defined average total atom number . We made use of a Hartree-Fock scheme for SPGPE parameter estimation developed in Ref. (28).

### References

- A. L. Fetter and A. A. Svidzinsky, J. Phys.: Condens. Matter 13, R135 (2001).
- A. L. Fetter, Rev. Mod. Phys. 81, 647 (2009).
- D. S. Rokhsar, Phys. Rev. Lett. 79, 2164 (1997).
- B. Jackson, N. P. Proukakis, C. F. Barenghi, and E. Zaremba, Phys. Rev. A 79, 053615 (2009).
- H. Schmidt, K. Góral, F. Floegel, M. Gajda, and K. Rza̧żewski, J. Opt. B 5, S96 (2003).
- E. J. Madarassy and C. F. Barenghi, J. Low. Temp. Phys. 152, 122 (2008).
- R. A. Duine, B. W. A. Leurs, and H. T. C. Stoof, Phys. Rev. A 69, 053623 (2004).
- H. T. C. Stoof, J. Low Temp. Phys. 114, 11 (1999).
- H. T. C. Stoof and M. J. Bilsma, J. Low Temp. Phys. 124, 431 (2001).
- J. Thompson, Phil. Mag. 10, 155 (1880).
- A. L. Fetter, Phys. Rev. A 69, 043617 (2004).
- M. Tsubota, J. Phys-Cond. Mat. 21, 164207 (2009).
- T. Mizushima, M. Ichioka, and K. Machida, Phys. Rev. Lett. 90, 180401 (2003).
- T. P. Simula, T. Mizushima, and K. Machida, Phys. Rev. Lett. 101, 020402 (2008).
- T. P. Simula, T. Mizushima, and K. Machida, Phys. Rev. A 78, 053604 (2008).
- V. Bretin, P. Rosenbusch, F. Chevy, G. V. Shlyapnikov, and J. Dalibard, Phys. Rev. Lett. 90, 100403 (2003).
- T.-L. Horng, S.-C. Gou, and T.-C. Lin, Phys. Rev. A 74, 041603 (2006).
- T. W. Neely, E. C. Samson, A. S. Bradley, M. J. Davis, and B. P. Anderson, Phys. Rev. Lett. 104, 160401 (2010).
- C. W. Gardiner, J. R. Anglin, and T. I. A. Fudge, J. Phys. B: At. Mol. Opt. Phys. 35, 1555 (2002).
- C. W. Gardiner and M. J. Davis, J. Phys. B: At. Mol. Opt. Phys. 36, 4731 (2003).
- P. B. Blakie, A. S. Bradley, M. J. Davis, R. J. Ballagh, and C. W. Gardiner, Adv. in Phys. 57, 363 (2008).
- N. P. Proukakis, Phys. Rev. A 74, 053617 (2006).
- N. P. Proukakis, J. Schmiedmayer, and H. T. C. Stoof, Phys. Rev. A 73, 053603 (2006).
- N. P. Proukakis and B. Jackson, J. Phys. B 41, 203002 (2008).
- N. Sepúlveda, Phys. Rev. A 83, 043603 (2011).
- A. S. Bradley, C. W. Gardiner, and M. J. Davis, Phys. Rev. A 77, 033616 (2008).
- C. N. Weiler, T. W. Neely, D. R. Scherer, A. S. Bradley, M. J. Davis, and B. P. Anderson, Nature 455, 948 (2008).
- S. J. Rooney, A. S. Bradley, and P. B. Blakie, Phys. Rev. A 81, 023630 (2010).
- T. M. Wright, A. S. Bradley, and R. J. Ballagh, Phys. Rev. A 80, 053624 (2009).
- T. M. Wright, R. J. Ballagh, A. S. Bradley, P. B. Blakie, and C. W. Gardiner, Phys. Rev. A 78, 063601 (2008).
- T. M. Wright, A. S. Bradley, and R. J. Ballagh, Phys. Rev. A 81, 013610 (2010).
- We require to be fixed and the highest energy mode in C to be appreciably occupied ( atoms).
- We do not find it necessary to imprint a vortex density profile because at the temperatures under consideration the extra energy of the core with phase-only imprinting is negligible compared with that from thermal fluctuations.
- The number of trajectories required to obtain good statistics increases with the system temperature.
- To obtain good statistics we integrate an instantaneous curve (as in Fig. 3), and then time-average the result for . This time is chosen to be much shorter than the vortex decay-time (which exceeds for all parameters), to give a short-time ensemble-average of the power data.
- S. Giorgini, L. P. Pitaevskii, and S. Stringari, Phys. Rev. A 54, R4633 (1996).
- The system still has a phase space density sufficiently high that spontaneous vortices are not expected to be significant (e.g. see (39); ?; ?).
- Our simulation initial conditions correspond to a grand canonical ensemble for the C-field region, with well defined average total atom number . We made use of a Hartree-Fock scheme for SPGPE parameter estimation developed in Ref. (28).
- Z. Hadzibabic, P. Kruger, M. Cheneau, B. Battelier, and J. Dalibard, Nature 441, 1118 (2006).
- T. P. Simula and P. B. Blakie, Phys. Rev. Lett. 96, 020404 (2006).
- R. N. Bisset, M. J. Davis, T. P. Simula, and P. B. Blakie, Phys. Rev. A 79, 033626 (2009).