Vortex line and ring dynamics in trapped Bose-Einstein condensates
B. Jackson,1 J. F. McCann,2 and C. S. Adams1
arXiv:cond-mat/9907325v1 [cond-mat.stat-mech] 21 Jul 1999
Department of Physics, Rochester Building, University of Durham, South Road, Durham, DH1 3LE, UK. Department of Applied Mathematics and Theoretical Physics, Queen’s University, Belfast, BT7 1NN, UK. (February 7, 2008)
Vortex dynamics in inhomogeneous Bose-Einstein condensates are studied numerically in two and three dimensions. We simulate the precession of a single vortex around the center of a trapped condensate, and use the Magnus force to estimate the precession frequency. Vortex ring dynamics in a spherical trap are also simulated, and we discover that a ring undergoes oscillatory motion around a circle of maximum energy. The position of this locus is calculated as a function of the number of condensed atoms. In the presence of dissipation, the amplitude of the oscillation will increase, eventually resulting in self-annihilation of the ring. PACS numbers: 03.75.Fi, 05.30.Jp, 67.40.Vs
Among the most important phenomena associated with Bose-Einstein condensation (BEC) is the quantization of vorticity, which is intimately connected with the existence of persistent currents and super?uidity in quantum ?uids. Study of quantized vortices has been con?ned mainly to liquid HeII , where detailed comparison to mean-?eld theory is complicated by strong interactions between atoms. However, such considerations are much less important for the recently achieved BEC in atomic vapors [2–5]. In this case the condensate can be accurately described by the Gross-Pitaevskii (GP) equation, an example of a nonlinear Schr¨ odinger equation, whose properties are well-known. This equation admits vortex solutions, where a non-zero circulation is accompanied by a zero in the condensate density. The density variation de?nes the vortex core, with a size of ? 1 ?m (cf. ? 1 ? A in HeII). Thus, vortices may be directly observable by absorption imaging [2,6] or other detection schemes [7–10]. Apart from their intrinsic interest, vortices play an important role in the breakdown of super?ow in Bose ?uids [1,11,12]. Despite the enormous progress of experiments on atomic condensates , a vortex state has yet to be observed. One possible procedure involves rotating the condensate [14,15] at a rate exceeding a critical angular velocity [16–18], creating a singly-quantized vortex line along the axis of rotation [19,20]. Larger angular velocities might be expected to produce vortices with higher quanta of circulation. However, previous studies [21,22] suggest that such vortices are unstable, leading to arrays of singly-charged vortices [23,24] similar to those found in HeII . 1
A separate but complementary idea is to use a tightlyfocussed far-o? resonant laser beam , which creates vortex pairs when dragged through the condensate. Alternatively, by stirring the condensate angular momentum is transferred, and a single vortex may be produced . A vortex ring may be formed by translating one condensate through another , or by 3D soliton decay [28,29]. Experimentally, the Bose-condensed gas is usually con?ned in a magnetic trap, often modelled by a harmonic potential. This profoundly alters the condensate properties in many ways [16,30], the most apparent of which is its spatial inhomogeneity. This results in a variation of vortex energy as a function of position, with a maximum at the center of a non-rotating trap. As a consequence, a single vortex precesses around the condensate center [8,31]. In addition, the vortex is thermodynamically unstable, and dissipation at ?nite temperatures leads to its expulsion from the cloud [32,33]. The instability is also apparent in the excitation spectrum, with the existence of a mode possessing negative energy with respect to the ground state [22,34,35]. In this paper, we study the motion of vortices in trapped Bose condensates, by numerical solution of the GP equation. As discussed in Section II, this equation is valid in the limit of low temperatures, and describes the conservative motion of a vortex. Section III presents measurements of the precession frequency of a single vortex in two and three dimensions, and compares the results to analytical expressions. We use a Magnus force argument to estimate the precession frequency in 2D. In Section IV we study vortex rings, and ?nd that they perform a cyclical motion. A stationary state can be found which corresponds to a ring with maximum energy. This point of unstable equilibrium is equivalent to that of a single vortex in the center of a non-rotating condensate. Finally, we summarize in Section V, and brie?y discuss ?nite temperature e?ects.
II. THEORY A. The Gross-Pitaevskii equation
In experiments on atomic vapours [2–4], evaporative cooling can be extended to very low e?ective temperatures, such that the non-condensate fraction is very small. The densities are su?ciently low that interactions can be represented by a pseudopotential of the form U0 δ (r ? r ′ ),
where U0 = 4π 2 a/m, and a is the s-wave scattering length. This leads to the Gross-Pitaevskii (GP) equation [36,37] for the condensate wavefunction, Ψ(r ,t), given by: i ?Ψ = ?t
B. The Vortex State
?2 + V + N U0 |Ψ|2 Ψ,
where N is the number of atoms, each of mass m. Consequently, N |Ψ|2 is the condensate number density, which satis?es the normalization condition: |Ψ(r ,t)|2 d3 r = 1. (2)
Using the Madelung transformation , the condensate wavefunction can be represented in terms of its den√ sity ρ(r , t) = |Ψ(r, t)|2 and phase S (r ,t) by Ψ = ρ eiS . ? The super?uid velocity is given by vs = (Ψ ?Ψ ? Ψ?Ψ? )/imρ = ( /m)?S . So, the circulation around an arbitrary closed loop is: κ= m ?S · dl = nh , m (6)
The harmonic trapping potential is denoted by V (r) = 2 2 2 2 2 2 m(ωx x + ωy y + ωz z )/2. It is convenient to scale (1) in dimensionless units. We use harmonic oscillator units (h.o.u.) , where the units ?1 of length, time and energy are ( /2mωx)1/2 , ωx and ωx respectively. We also consider a frame rotating about the z -axis with angular velocity ?, where the angular momentum operator is given by Lz = i(y?x ? x?y ). Retaining the normalization condition (2), Eq. (1) becomes: ? + C |Ψ|2 ? ?Lz Ψ, i?t Ψ = ??2 + V (3)
where the property of a single-valued wavefunction restricts n to an integer value. From (6), ? × v s = i 2πni δ (r ? r 0i ), and the super?uid is irrotational everywhere except for vortices at r = r0i . We de?ne a circulation vector κ on the vortex axis, which has magnitude |κ| = κ and direction ? × v s . A stationary state exists for a single vortex line at r0 = 0. If, for simplicity, we assume a non-rotating isotropic trap in the x ? y plane, then this state can be represented as Ψ(r, φ, z ) = ρ(r, z )einφ . Substitution into (1) gives:
2 2 1 ?f ?2f n2 f ωx 2 ?2f + + 2 ? 2 ? (? r + ?z ?2 )f 2 ?r ? r ? ?r ? ?z ? r ? 4? +f ? f 3 = 0,
? = 1 (x2 + ηy 2 + ?z 2 ) and the anisotropy where V 4 2 2 parameters are de?ned as η = ωy /ωx and ? = 2 2 ωz /ωx . The interaction parameter is given by C = (N U0 / ωx )(2mωx / )α/2 , where α is the number of dimensions. In 2D, where atoms are con?ned in the x ? y plane, N represents the number of atoms per unit length along z . Eq. (3) can be integrated by various numerical methods . As in our previous work [25,27], we utilize a Fast Fourier Transform (FFT) technique. Stationary solutions of (3), can be represented by Ψ(r, t) = Ψ(r)e?i?t where ? is the chemical potential. They are found by propagating a trial wavefunction (e.g. a Gaussian so?. In lution for C = 0) in imaginary time, t → ?it contrast to real time, the resulting evolution operator is non-unitary, so the wavefunction must be renormalized after each time-step. The ratio of norms provides a convenient estimate for the chemical potential, ? = (2?t)?1 ln[ |Ψ(t)|2 / |Ψ(t + ?t)|2 ]. In imaginary time, excitations are exponentially damped, and both Ψ and ? rapidly converge to a stationary solution, providing accurate initial conditions for time-dependent simulations. The free energy per particle, E , is given by: E= ? |Ψ|2 + |?Ψ|2 + V C E =?? 2 C |Ψ|4 ? ?Ψ? Lz Ψ d3 r , (4) 2
where r ? = (2m?)1/2 r/ , z ? = (2m?)1/2 z/ , and f = 1/2 1/2 (N U0 ρ/?) = (ρ/ρ0 ) . Eq. (7) yields the following asymptotic forms: ρ ? ρ0 nξ r r ξ
r ? ξ,
ρ ? ρ0
ξ ? r < R⊥ ,
where for a stationary solution it can be shown that: |Ψ| d r .
2 2 and ρ = 0 for r ≥ R⊥ and |z | ≥ Rz , where R⊥ = 2?/mωx 2 2 1/2 and Rz = 2?/mωz . The parameter ξ = /(2m?) = (8πρ0 a)?1/2 is the coherence length in the center of the condensate , and determines the size of the vortex core. Eq. (9) is presented in [17,19,32] and is analogous to the well-known Thomas-Fermi (TF) approximation for the ground state, valid for large N U0 . In this limit, ρ0 corresponds to the density in the center of the ground state condensate. For the special case n = 1, Eqs. (8) and (9) can be interpolated by a simple function, which generalized to a straight, o?-axis vortex line at r 0 is:
ρ ≈ ρ0
|r ? r 0 |2 |r ? r 0 |2 + ξ 2
Imaginary time propagation minimizes the chemical potential, ?. 2
where ρ = 0 when the right-hand side is negative. We compare Eq. (10) to the numerical solution of the GP equation in Fig. 1. The latter is found by imaginary time propagation, imposing a phase variation of 2π around the vortex line at each time-step. Eq. (10) provides a good estimate especially in the high-N limit, in a similar way to the TF approximation to the ground state.
Erot = En ? n? . Appearance of a vortex becomes energetically favorable when Erot < E0 , so that the critical angular velocity is simply: ?c = En ? E0 . n (11)
In the TF limit [17,31,40]: ?c = 5 2 2 2 (ω + ωy ) ln 8? x R⊥ ξ . (12)
|Ψ| (10 )
For |?| > ?c , an on-axis vortex attains global stability; however, there remains an energy barrier for vortices entering from the edge  (see Fig. 2). Inspection of Fig. 2 also reveals that that above an angular velocity ?m , the vortex attains a local minimum, so a vortex is metastable when ?m < |?| < ?c . In the TF limit, ?m = 3?c /5 , while in the non-interacting limit ?m → ?c . In both the weak and strong coupling limits, ?m ? ?ωa , where ωa is the frequency of the socalled ‘anomalous’ mode obtained from solution of the Bogoliubov equations [31,34,35,41]. It is thought that ωa corresponds to the frequency of vortex precession, thus linking vortex dynamics and instability.
FIG. 1. Cross-section through a singly-quantized vortex line, showing condensate density as a function of position. The solid line plots the exact 3D wavefunction as calculated from imaginary time propagation of Eq. (3), for C = 4000, ? = 0 and ? = η = 1. The dashed line represents the analytic approximation (10), with ρ0 ? 2.443 × 10?3 , ξ ? 0.3199, r0 = 0, and R⊥ = Rz ? 6.252. These parameters are calculated from the chemical potential ?TF = (15C/64π )2/5 , given by the normalization condition of the ground-state Thomas-Fermi wavefunction.
III. SINGLE VORTEX MOTION A. Two Dimensions
In this section we study the dynamics of a singlyquantized (n = 1) vortex line in a non-rotating condensate. Imaginary time propagation of the GP equation (3) is used to provide an initial condition for the real-time simulation. If we consider a vortex at position r0 relative to the axis of a non-rotating trap, then the free energy of the system, Eq. (4), is found to attain a maximum when r0 = 0 (see Fig. 2). So, a vortex initially at the origin will remain stationary, but is unstable to in?nitesimal displacements. The GP equation implies that the system is Hamiltonian, and therefore an o?-axis vortex will follow a path of constant energy corresponding to precession around the trap center. The presence of dissipation will lead to drift towards lower energies, causing the vortex to spiral out of the condensate [32,33]. If the trap is rotated with angular velocity, ?, then the energy of a central vortex, Erot , decreases such that
FIG. 2. Energy, E , as a function of n = 1 vortex displacement, in a 2D condensate rotating with angular frequency, ? (C = 1000, η = 1). The top solid curve corresponds to ? = 0, where ? increases in steps of 0.05 as one moves towards the lowest curve. The dashed line marks the energy of the condensate without a vortex.
First, we consider the simplest case of a vortex in two dimensions, which corresponds experimentally to a condensate con?ned in an axisymmetric cylindrical trap (where ? → 0, η = 1). For ? = 0, simulations show that an o?-center vortex accelerates from its initial condition, soon attaining a near-constant angular velocity ω around the trap center, such that the instantaneous velocity is v L = ω κ ? × r . The angular velocity, ω , is plotted as a function of interaction strength and initial position in Fig. 3. For small C , ω is averaged over a few cycles 3
(e.g. 3 revolutions for C = 200). However, for higher C , numerical instabilities can restrict the simulations to less than a half-cycle: the error bars re?ect the resulting uncertainty.
may be shown by using the decomposition Ψ(x, y ) = Φ(x, y )Θ(x, y ), where an o?-set vortex core (Φ) is imprinted on the TF ground-state (Θ), as implied by Eq. (10). Starting from the GP equation (1), and noting that in the TF limit Θ is slowly varying, then it is simple to show that:
?2 Φ + ?TF |Φ|2 Φ = ?Φ, (13) 2m where ?TF is the chemical potential of Θ. As the Laplacian is spatially invariant, it follows that ? is independent of the o?set at small r0 . This can also be justi?ed numerically, though the approximation only becomes valid at high C . Using Eq. (5) it follows that: ?
0 2 r0 4 6
?E C ? ≈? ?r0 2 ?r0
ρ2 d2 r .
Substituting Eq. (10) for ρ then gives an estimate for the precession frequency (to logarithmic accuracy):
100 1000 C 10000
|ω | ≈
2 ωx ln ?
5 . 4
FIG. 3. Vortex precession frequency, ω , in a 2D condensate at r0 = 0.5, as a function of interaction parameter C = 8πN a. Filled circles show the results of numerical simulations. The triangles indicate the Magnus force estimates, obtained from the gradient of the numerical values of E0 . The analytical Magnus force estimate, Eq. (15), and Eq. (16) are plotted with dot-dashed and dashed lines, respectively. The TF vortex metastability frequency, 3 ?E , is plotted as a solid line. Inset: 5 ω as a function of vortex position r0 , for C = 1000.
This result may be compared to the expression obtained by Svidzinsky and Fetter , using a timedependent variational analysis: |ω | =
2 3 ωx ln 4?
The precession frequency decreases as a function of increasing interaction strength, C . An intuitive semianalytical argument for this behavior can be formulated in terms of the Magnus e?ect, familiar from classical hydrodynamics, and in super?uids and superconductors [42,43]. When the background ?uid ?ows past the circulating ?uid connected with the vortex, a pressure imbalance is created perpendicular to the direction of the background ?ow. The resulting Magnus force must balance the force due to the variation of energy E with position, i.e. ?E/?r0 = mρκ × vL , where v L is the velocity of the vortex line relative to the ambient condensate. Note that the Magnus force can also be produced by ?ow of the condensate around a stationary vortex, as in a rotating trap. So, one expects that ω ? ?m . We ?nd E by evaluating the functional, Eq. (4), with a wavefunction grown in imaginary time with ? = 0. Using the Madelung transformation, the ?rst term splits into a ‘quantum pressure’ and a ‘kinetic energy’ term: √ |?Ψ|2 = (? ρ)2 + ρ(?S )2 . A numerical di?erentiation of E with respect to r0 gives an estimate for ω using the Magnus force argument (plotted as triangles in Fig. 3). However, these estimates are sensitive to small numerical errors in the energy. To obtain an analytical estimate, we observe that ? is approximately constant at small r0 and high C . This
again valid for small r0 . These expressions are plotted together with the numerical results in Fig. 3. In addition, we plot ?m = 3 5 ?E , where ?E is the energy di?erence between the ground and n = 1 vortex states. Recall that one expects that ω ? ?m in the TF limit. The numerical results lie between the results of (16) and the metastability curve. All of the curves reproduce the correct functional dependence at high C . Note that these expressions are only valid for small r0 ; the vortex precesses faster as it nears the edge of the condensate, as shown in Fig. 3 (inset). Compressibility e?ects become important when the vortex is accelerating or when the velocity is an appreciable fraction of the speed of the sound, cs = (ρN U0 /m)1/2 . In an in?nite compressible ?uid, phonons may be emitted by a moving vortex, leading to a drift of the vortex to lower energies . However in a ?nite condensate, where excitations remain con?ned in the region of the vortex, no net drift is expected. At the beginning of the simulations, we observe an increase in radius of precession, together with excitation of an elliptical center of mass mode at the trap frequency. In addition, surface waves are created when the vortex is near the condensate edge. However, as expected we do not observe a sustained vortex drift (for T up to 110, corresponding to ? 6 full cycles for C = 200). A drift to lower energies would be expected where a thermal cloud damps the motion (i.e. at ?nite temperatures). Nevertheless, the vortex lifetime is expected to be long, especially for large numbers of atoms . 4
B. Three Dimensions
Vortex dynamics become more complex in 3D, as the vortex line can deform along its length. In classical and quantum ?uids, this can result in the propagation of helical waves along the line—so-called Kelvin modes [1,36]. In simulations of 3D vortex motion, we have observed line deformation and oscillations. However, the inhomogeneity of the condensate complicates matters, and the motion is di?cult to resolve into simple Kelvin waves characteristic of the bulk condensate. The amplitude of the oscillations are typically small, and as a consequence helical waves are likely to be di?cult to detect experimentally. Moreover, it is worth noting that the energy of the vortex increases as it lengthens. Hence, in the presence of the dissipation the line will tend to straighten, e?ectively damping the Kelvin modes. Fig. 4 compares numerically measured values of the precessional frequency with the TF result of Svidzinsky and Fetter Eq. (16)  . It can be seen that the frequency dependence is well described; however, the numerical results are signi?cantly higher (? 20%), but converge slowly to the analytical expression towards high C . This disparity may be due to e?ects resulting from the curvature of the line, which will modify predictions that assume a rigid line motion. Numerical values of the TF metastability frequency, 3 5 ?E , are also found to be lower than the observed precession frequencies.
precession due to the inhomogeneity of the condensate, as discussed for a single vortex in Sec. III, and second, the velocity induced by the remainder of the ring, vin , which is directed along its axis (de?ned as the z -axis). For a spherical condensate, the total velocity on each element is given by: v = vin z ? + ωκ ? × r, (17)
where κ ? de?nes the direction of the circulation at the element, and ω is the precession frequency. In a homogeneous Bose ?uid vin = ( /2mRr )[ln(8Rr /ξ ) ? 0.615] , where Rr is the ring radius. Consider a ring at z = 0, r = Rr . If the radius is small, the induced velocity dominates and the ring moves in the +z ? direction, while if Rr is large the precession dominates and it travels backwards. In addition, the precessional term leads to ring expansion for z > 0 and contraction for z < 0. Thus, the two terms produce an oscillatory motion of the ring.
6 4 2 0 ?2 ?4 ?6 ?6
FIG. 5. Vortex ring energy in a cylindrically symmetric condensate (η = ? = 1, C = 2000) as a function of radius, r , and z -position. The energy contours are equally spaced between 5.6179 and 6.1159. The bold line shows the motion of one element of the vortex ring, where the circles represent the position at equally-spaced times (every T = 1). The ring begins at (1.5, 0), marked by an arrow, and cycles around the energy maximum in a clockwise direction.
0 1000 2000
FIG. 4. Vortex line frequency, ω , in a 3D condensate (η = 1, ? = 9) at r0 = 0.5, plotted as a function of interaction parameter, C = 8πN a(2mω/ )1/2 . Squares display numerical results, while the dashed line plots Eq. (16). The 3 solid line shows the TF metastability frequency, 5 ?E .
IV. VORTEX RING MOTION
The motion of a vortex ring in a trapped BEC may be understood in terms of a sum of two contributions to the velocity of each element in the ring. First, the
One can also understand the ring motion as a trajectory around an energy maximum, in analogy with the single line vortex. To demonstrate this, in Fig. 5 we plot the energy of an on-axis ring as a function of its radius, r, and z -position. Without dissipation one would expect the ring motion to follow an energy contour; however as is apparent in Fig. 5 this is not exactly true. Acceleration of the ring at the beginning of its motion results in a back-action on the condensate, exciting a center-of-mass mode, and the subsequent ring dynamics are complicated by the underlying motion of the condensate. In addition, for small C we observe a decay of the ring to lower energies over the ?rst cycle of its motion. As for a single vortex, this e?ect is associated with the compressibility, 5
which results in acoustic emission from the moving ring . The energy maximum at r = Req , z = 0, corresponds to the point where the two velocity contributions in (17) are equal and opposite, leading to a ring in unstable equilibrium. To obtain an analytical estimate for this position, we approximate the ring energy by taking the energy of a single 2D vortex, Ev , and integrating around a circle of radius Rr , such that Er = 2πRr Ev . The dominant contribution to Ev is given by the kinetic energy, so 2 3 Ev = (m/2) ρvs d r. Taking v s = κ × r /2πr2 , where we translate the cylindrical coordinate system so that the origin lies on the vortex axis, and using (10) gives: Er ≈ 2π 2 ρ0 Rr m
FIG. 6. The equilibrium ring radius, Req , as a function of nonlinear parameter, C = 8πN a(2mω/ )1/2 , in a spherically symmetric condensate. The lower plot shows the ratio of Req to the Thomas-Fermi radius, RTF ≡ R⊥ = (2?T F /mω 2 )1/2 , where the chemical potential is given by the normalization condition of the TF wavefunction, such that ?TF = (15C/64π )2/5 .
2 r0 2 R⊥
2 ? r2 R⊥ 0 ξ
An experimental technique for ring production was proposed in Ref. . The method utilizes a twocomponent BEC, such that when the smaller component, |2 , is translated with respect to the other, vortex rings are created in the larger condensate |1 . To model the two condensates, a pair of coupled GP equations are solved for the wavefunctions Ψ1 (r , t) and Ψ2 (r, t): i?t Ψi = ??2 + Vi + C |Ψi |2 + |Ψj |2 Ψi , (20)
1 r2 + 0 2 ? 2 , R⊥
neglecting terms of order ξ 2 and higher. For a ring of 2 2 2 radius Rr at z0 , r0 = Rr + z0 . This expression describes the qualitative features of Fig. 5; however, it tends to over-estimate the energy near to the peak (? 10% at C = 2000) and is a poor approximation as Rr → 0. Nevertheless, we can obtain an estimate for the equilibrium position from Eq. (18), which yields z0 = 0 and: Req ≈ R⊥
4)?1 ln(βR⊥ 4 ) ? 4, 3 ln(βR⊥
2 where √ β = (mω/ ) . In the TF limit (R⊥ → ∞), Req ≈ R⊥ / 3 ? 0.577R⊥, which is close to the results from numerical solution, where Req ? 0.54R⊥ at high C (Fig. 6).
0.52 0.51 0.49
where i, j = 1, 2 (i = j ) . The initial states of the simulations are created by imaginary time propagation as previously, where the normalizations are set so that a variable fraction of atoms are in each condensate, χ = |Ψ1 |2 / |Ψ2 |2 , and |Ψ1 |2 + |Ψ2 |2 = 1. Our simulations then follow the creation and subsequent dynamics of the vortex ring, an example of which is shown in Fig. 7. The trajectories roughly follow a contour of constant energy (see Fig. 5). The ring is created from zero radius at t ≈ 1.6 and y ≈ 1.8. It then expands and travels forward, before turning and progressing backwards along the edge of the condensate. Finally, it turns again and collapses to a point, where the ring is annihilated. The annihilation produces a sound wave, which decreases in amplitude as it propagates along the z -direction. The sound wave then disappears at the edge of the cloud. In the presence of dissipation, annihilation will eventually occur for any ring as a culmination of a decay to lower energies. This is equivalent to the instability mechanism for a single vortex line (Sec. III). Vortex ring detection is likely to present considerable challenges to experimentalists. The simplest and most widely used method of condensate imaging is by measurement of probe laser absorption, after release of the condensate from the trap. The subsequent ballistic expansion results in an e?ective magni?cation of a vortex line, which should then be visible as a density ‘hole’ . However, for rings the density minima are obscured by the rest of the condensate along any line-of-sight. One solution is to view slices of the condensate after ballistic expansion, using light sheets. An alternative method is to study the center-of-mass motion of coupled condensates, yielding details of the mutual drag that reveal vortex ring formation . Collective excitations are also utilized in another scheme, proposed in [9,10], where a vortex line splits the degeneracy of the quadrupole mode of the condensate. However, further work is needed to extend this analysis to vortex rings.
2 6 4 2 0 -2 -4 2 0 -2 -4
expulsion of the vortex from the condensate. A full microscopic model of vortex dynamics is needed to describe this behaviour and will form the focus of future work.
Financial support for this work was provided by the EPSRC.
FIG. 7. Vortex ring motion in a condensate (η = ? = 1) after creation from an object. The trajectory of the ring is determined by solving Eq. (20), with v = 1.75, C = 1100 and χ = 10/11. The upper plot shows the ring radius (solid line) and z -coordinate (dashed) as a function of time, while a parametric plot (bottom) displays the ring radius, r , against position z .
We have studied the motion of vortex lines and rings in Bose-Einstein condensates in harmonic traps, by numerical solution of the Gross-Pitaevskii equation. We considered a single vortex in two and three dimensions. At the center of a non-rotating condensate the vortex state possesses maximum energy, corresponding to unstable equilibrium, while an o?-center vortex undergoes precession around this maximum. The precession frequency was measured and compared to theoretical models. Vortex rings were also found to undergo an oscillatory motion. For a particular radius, the ring energy is a maximum, corresponding to state of unstable equilibrium. A lower energy ring precesses around the locus of maximum energy. Also, a ring which is created at a point will eventually collapse to a point, resulting in self-annihilation. The study of dissipative vortex dynamics at ?nite temperatures is particularly interesting. In this case, interaction of the vortex with the thermal cloud leads to a transfer of energy, and consequently a decay of the vortex state. One can consider this to be due to anisotropic scattering of thermal excitations , where the momentum transfer is apparent as a frictional force on the vortex. The transverse component of this force results in 7
 R. J. Donnelly, Quantized Vortices in Helium II (Cambridge University Press, Cambridge, 1991).  M. H. Anderson et al., Science 269, 198 (1995).  K. B. Davis et al., Phys. Rev. Lett. 75, 3969 (1995).  C. C. Bradley, C. A. Sackett, and R. A. Hulet, Phys. Rev. Lett. 78, 985 (1997); C. C. Bradley, et al., Phys. Rev. Lett. 75, 1687 (1995).  D. G. Fried et al., Phys. Rev. Lett. 81, 3811 (1998).  E. Lundh, C. J. Pethick, and H. Smith, Phys. Rev. Lett. 58, 4816 (1998).  E. V. Goldstein, E. M. Wright, and P. Meystre, Phys. Rev. A 58, 576 (1998).  E. L. Bolda and D. F. Walls, Phys. Rev. Lett. 81, 5477 (1998).  F. Zambelli and S. Stringari, Phys. Rev. Lett. 81, 1754 (1998).  A. A. Svidzinsky and A. L. Fetter, Phys. Rev. A 58, 3168 (1998).  T. Frisch, Y. Pomeau, and S. Rica, Phys. Rev. Lett. 69, 1644 (1992).  T. Winiecki, J. F. McCann, and C. S. Adams, Phys. Rev. Lett. 82, 5186 (1999).  see e.g. W. Ketterle, D. S. Durfee, and D. M. StamperKurn, preprint cond-mat/9904034.  K. -P. Marzlin and W. Zhang, Phys. Rev. A 57, 4761 (1998).  S. Stringari, Phys. Rev. Lett. 82, 4371 (1999).  G. Baym and C. J. Pethick, Phys. Rev. Lett. 76, 6 (1996).  S. Sinha, Phys. Rev. A 55, 4325 (1997).  E. Lundh, C. J. Pethick, and H. Smith, Phys. Rev. A 55, 2126 (1997).  S. Stringari, Phys. Rev. Lett. 76, 1405 (1996).  F. Dalfovo and S. Stringari, Phys. Rev. A 53, 2477 (1996).  H. Pu, C. K. Law, J. H. Eberly, and N. P. Bigelow, Phys. Rev. A 59, 1533 (1999).  J. J. Garcia-Ripoll, V. M. Perez-Garcia, preprint condmat/9903353.  D. A. Butts and D. S. Rokshar, Nature 397, 327 (1999).  D. L. Feder, C. W. Clark, and B. I. Schneider, preprint (1999).  B. Jackson, J. F. McCann, and C. S. Adams, Phys. Rev. Lett. 80, 3903 (1998).  B. M. Caradoc-Davies, R. J. Ballagh, and K. Burnett, preprint cond-mat/9902092.
 B. Jackson, J. F. McCann, and C. S. Adams, preprint cond-mat/9901087.  C. A. Jones, S. J. Putterman and P. H. Roberts, J. Phys. A 19, 2991 (1986).  C. Josserand and Y. Pomeau, Europhys. Lett. 30, 43 (1995).  F. Dalfovo, L. Pitaevskii, S. Giorgini, and S. Stringari, Rev. Mod. Phys. 71, 463 (1999).  A. A. Svidzinsky and A. L. Fetter, preprint condmat/9811348.  D. S. Rokhsar, Phys. Rev. Lett. 79, 2164 (1997).  P. O. Fedichev and G. V. Shlyapnikov, preprint condmat/9902204.  R. J. Dodd, K. Burnett, M. Edwards, and C. W. Clark, Phys. Rev. A 56, 587 (1997).  T. Isoshima and K. Machida, J. Phys. Soc. Jpn. 67, 1840 (1998).  L. P. Pitaevskii, Zh. Eksp. Teor. Fiz. 40, 646 (1961) [Sov. Phys. JETP 13, 451 (1961)].  E. P. Gross, Nuovo Cimento 20, 454 (1961).  P. A. Ruprecht, M. J. Holland, K. Burnett, and M. Edwards, Phys. Rev. A 51, 4704 (1995).  See e.g. M. D. Feit, J. A. Fleck, and A. Steiger, J. Comp. Phys. 47, 412 (1982).  D. L. Feder, C. W. Clark, and B. I. Schneider, Phys. Rev. Lett. 82, 4956 (1999).  M. Linn and A. L. Fetter, preprint cond-mat/9906139.  E. B. Sonin, Phys. Rev. B 55, 485 (1997), and references therein.  P. Ao and D. J. Thouless, Phys. Rev. Lett. 70, 2158 (1993).  Yu. N. Ovchinnikov and I. M. Sigal, Nonlinearity 11, 1295 (1998).  P. H. Roberts and J. Grant, J. Phys. A 4, 55 (1971).  L. S. Pismen and A. A. Nepomnyashchy, Physica D 69, 163 (1993).  The trap potentials are given by V1 = (x2 + y 2 )/4 and V2 = V1 + (x2 + (y ? vt)2 ), where v is a constant translational velocity. The scattering lengths between hyper?ne states are assumed to be equal. Simulation of the two-component dynamics involves propagating (20) for Ψ1 and Ψ2 consecutively at each time step. We restrict attention to the cylindrically symmetric case (η = 1), so that ?2 = (1/r )?r + ?rr + ?zz . To propagate (20), the resulting Hamiltonian is split into operators involving the ?rst derivative term, which is dealt with by ?nite di?erences, and the remainder which is treated by the FFT approach mentioned previously. Care must be taken with the singularity at r = 0, where a forward di?erencing approximation is utilized (1/r )?r Ψ|r=0 = 2[Ψ(?r ) ? Ψ(0)]/(?r )2 + O(?r 2 ), whilst elsewhere central di?erencing is applied.