ISSN 1054-660X, Laser Physics, 2008, Vol. 18, No. 3, pp. 322–330.

PHYSICS OF COLD TRAPPED ATOMS

© MAIK “Nauka /Interperiodica” (Russia), 2008. Original Text © Astro, Ltd., 2008.

Instabilities and Vortex-Lattice Formation in Rotating Conventional and Dipolar Dilute-Gas Bose–Einstein Condensates A. M. Martina, N. G. Parkera, R. M. W. van Bijnena, b, A. Dowa, and D. H. J. O’Dellc a School

of Physics, University of Melbourne, Parkville, Victoria, 3010 Australia University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands c Department of Physics and Astronomy, McMaster University, Hamilton, Ontario, L8S 4M1 Canada b Eindhoven

e-mail: [email protected] Received October 15, 2007

Abstract—A theoretical study of vortex-lattice formation in atomic Bose–Einstein condensates confined by a rotating elliptical trap is presented. For the conventional case of purely s-wave interatomic interactions, this is done through a consideration of both hydrodynamic equations and time-dependent simulations of the Gross– Pitaevskii equation. We discriminate three distinct, experimentally testable regimes of instability: ripple, interbranch, and catastrophic. Additionally, we generalize the classical hydrodynamical approach to include longrange dipolar interactions, showing how the static solutions and their stability in the rotating frame are significantly altered. This enables us to examine the routes towards unstable dynamics, which, in analogy to conventional condensates, may lead to vortex-lattice formation. PACS numbers: 03.75.Kk, 34.20.Cf, 47.20.-k DOI: 10.1134/S1054660X08030225

1. INTRODUCTION Vortex lattices are a striking manifestation of superfluidity in rotating systems. Such states have been generated in atomic Bose–Einstein condensates (BECs) via rotation in an elliptical harmonic trap [1–3]. Vortex lattices form for trap rotation frequencies Ω ~ 0.7ω⊥, where ω⊥ is the mean harmonic-trap frequency in the rotating plane. Although a vortex lattice is thermodynamically favorable for much lower rotation frequencies, the instabilities necessary to seed vortex-lattice formation are only present in the region of Ω ~ 0.7ω⊥. Such instabilities have been predicted by hydrodynamic studies of condensate solutions in the rotating frame [4] and the dynamical perturbations of these states [5]. Although numerical simulations of the Gross–Pitaevskii equation (GPE) have also observed vortex-lattice formation in this region [6–13], the results are mixed: for Ω
0.7ω⊥, [6–8] require thermal effects to enable vortex nucleation, while [9–13] do not; for Ω  0.7ω⊥, [6, 9–12] observe a shape instability before nucleating vortices even in the absence of thermal effects. This is consistent with experiments [3, 14] which indicate that lattice formation is temperatureindependent. Breaking the rotational symmetry of such simulations has been shown to be crucial to enable the realistic nucleation of vortices [10–13]. The above results apply to conventional BECs composed of atoms of mass m with short-range s-wave interactions, parameterized via g = 4π2a/m, where a is the s-wave scattering length. However, a recent experiment has formed a BEC of chromium atoms with dipolar interactions [15]. Chromium has an anomalously

large magnetic dipole moment of 6 Bohr magnetons, which leads to magnetic dipole–dipole interactions that are 36 times stronger than those found in most alkalis. Theoretical work, using a modified GPE, has studied the effect of such long-range interactions on groundstate vortex-lattice solutions [16–18]. In a recent paper [19], we explored, for the first time, the route to generate such states in dipolar BECs. This was achieved by solving the hydrodynamic equations of motion for a dipolar BEC in a rotating anisotropic harmonic trap. In this paper, we present a theoretical analysis of rotating condensates, with and without dipolar interactions, and indicate the routes to vortex-lattice formation. In Section 2, we derive static solutions to the hydrodynamic equations for a dipolar BEC in the rotating frame. In Section 3, we consider the case of short-range interactions only and compare the static hydrodynamic solutions to the GPE. Then, in Section 4, we investigate the hydrodynamical solutions for dipolar gases. Finally, in Section 5, we discuss the results presented and comment on future directions of research. 2. HYDRODYNAMICAL EQUATIONS 2.1. Static Solutions Consider a BEC with long-range dipole–dipole interactions. The potential between dipoles, separated by r and aligned by an external electric or magnetic field along a unit vector eˆ , is given by [20]

322

C dd ( δ ij – 3rˆ i rˆ j ) - eˆ eˆ ---------------------------. U dd ( r ) = ------3 4π i j r

(1)

INSTABILITIES AND VORTEX-LATTICE FORMATION

For two atoms with dipoles induced by a static electric field E = Eeˆ , the coupling constant Cdd = E2α2/0 [21, 22]. Alternatively, if the atoms have permanent magnetic dipoles dm aligned in an external magnetic 2

field B = Beˆ , one has Cdd = µ 0 d m [23]. Denoting ρ as the condensate density, the dipolar interactions give rise to a mean-field potential

323

where v = (/m)∇S is the fluid velocity field in the laboratory frame, expressed in the coordinates in the rotating frame. Setting ∂ρ/∂t = ∂ v /∂t = 0, we look for stationary solutions with an irrotational velocity field of the form [4, 5, 19] v = α ( yiˆ + x ˆj ),

(9)

(2)

where α, the amplitude of the velocity field, is to be determined. Then, Eq. (8) leads to

which can be included in a generalized GPE [18, 23, 24] for the BEC

m ˜2 2 ˜2 2 2 2 2 µ = ---- ( ω x x + ω y y + γ ω ⊥ z ) + gρ ( r ) + Φ dd ( r ),(10) 2

Φ dd ( r ) =

∫

3

d r'U dd ( r – r' )ρ ( r' ),

where

∂ψ ( r ) i --------------∂t 2

 2 2 = – -------∇ + V ( r ) + Φ dd ( r ) + g ψ ( r ) ψ ( r ). 2m

(3)

2

(4)

where µ is the chemical potential and γ = ωz /ω⊥ is the trap-aspect ratio. For ease of calculation, the dipolar potential Φdd( r ) can be expressed in terms of a fictitious “electrostatic” potential φ( r ) [26] δ Φ dd ( r ) = – 3gε dd eˆ i eˆ j ⎛ ∇ i ∇ j φ ( r ) + ----ij-ρ ( r )⎞ , ⎝ ⎠ 3

(11)

˜ 2y = ω ⊥2 ( 1 +  ) + α 2 + 2αΩ ω

(12)

and

In the Thomas–Fermi regime [25], the GPE describing a static dipolar BEC in a harmonic trapping potential is mω 2 2 2 2 µ = -----------⊥ [ ( 1 –  )x + ( 1 +  )y + γ z ] 2 + gρ ( r ) + Φ dd ( r ),

˜ 2x = ω ⊥2 ( 1 –  ) + α 2 – 2αΩ ω

(5)

are effective trap frequencies. The form of Eq. (10) is identical to Eq. (4). Hence, we can use the methodology presented in [26] to calculate Φdd( r ). An exact solution to Eq. (10) is given by 2 2 2 ⎛ x y z ⎞ ρ = n 0 ⎜ 1 – ------2 – ------2 – ------2⎟ ⎝ R x R y Rz ⎠

for

ρ ≥ 0,

(13)

where n0 = 15N/(8πRx Ry Rz) is the central density. Following [26], the dipole potential for a polarizing field aligned along the z axis is n0 κ x κ y x β x + y β y + 3z β z ρ Φ dd - β 0 – --------------------------------------------------------- = --------------- – --- , (14) 2 3 2 3gε dd R 2

2

2

z

where φ(r) =

∫

where

3

d r'ρ ( r' ) --------------------- , 4π r – r'

(6)

and εdd = Cdd /3g parameterizes the relative strength of the dipolar and s-wave interactions. Self-consistent solutions for Eq. (4) for ρ( r ), φ( r ), and, hence, Φdd( r ) can be found for any general parabolic trap [26]. Consider atoms in a harmonic potential rotating at a frequency Ω about the z axis. In the mean-field approximation, the evolution of the condensate field ψ( r , t) is described by the time-dependent GPE. Writing the condensate field as ψ( r , t) = ρ ( r, t ) exp{is( r , t)} and neglecting the quantum pressure, we obtain the superfluid hydrodynamic equations ∂ρ ------ + ∇ ⋅ [ ρ ( v – Ω × r ) ] = 0, ∂t

(7)

∂v v ⋅ v V gρ Φ dd – v ⋅ [ Ω × r ]⎞ = 0, (8) ------- + ∇ ⎛ ------------ + ---- + ------ + -------⎝ ⎠ m ∂t 2 m m LASER PHYSICS

Vol. 18

No. 3

2008

∞

dσ β k = ------------------------------------------------------------------------------------------------------ (15) 2 2 2 0 ( 1 + σ )( κk + σ ) ( κ x + σ )( κ y + σ )( 1 + σ )

∫

R with k = x, y, z, κk = -----k , and Rz ∞

β0 =

dσ

∫ ---------------------------------------------------------------------------------. (1 + σ) (κ + σ)(κ + σ)(1 + σ) 0

2 x

2 y

(16)

Thus, we can rearrange Eq. (10) to obtain the density m ˜2 2 ˜2 2 2 2 2 µ – ---- ( ω x x + ω y y + γ ω⊥ z ) 2 ρ = -----------------------------------------------------------------------g ( 1 – ε dd ) n0 κ x κ y 2 2 2 2 - [ x β x + y β y + 3z β z – R z β 0 ] 3gε dd --------------2 2R z -. + ---------------------------------------------------------------------------------------------------g ( 1 – ε dd )

(17)

324

MARTIN et al.

Comparing the x2, y2, and z2 terms in Eq. (13) and Eq. (17), we find the three self-consistency relations 3 3 1 + ε dd ⎛ --- κ x κ y β x – 1⎞ ⎝ ⎠ 2 γ ω 2 κ x = ⎛ ---------⊥-⎞ ---------------------------------------------------, ⎝ ω ˜x ⎠ ζ 2

3 3 1 + ε dd ⎛ --- κ y κ x β y – 1⎞ 2 ⎝ ⎠ 2 γ ω 2 κ y = ⎛ ---------⊥-⎞ ---------------------------------------------------, ⎝ ω ⎠ ˜y ζ

(18)

(19)

and 2gn 0 2 - ζ, R z = ---------------2 2 mγ ω ⊥

(20)

9κ x κ y - βz . ζ = 1 – ε dd 1 – ------------2

(21)

with

Using Eq. (17), the stationary solutions to Eq. (7) are governed by 2

2

ω⊥ κ x κ y γ ⎞ ˜ 2 – 3--- ε -----------------------βx 0 = (α + Ω)⎛ω ⎝ x 2 dd ⎠ ζ 2

2

ω⊥ κ x κ y γ ⎞ ˜ 2y – 3--- ε dd ----------------------- βy . + (α – Ω)⎛ω ⎝ ⎠ 2 ζ

(22)

2.2. Stability of Static Solutions The stationary solutions of Eqs. (7) and (8) are not necessarily stable solutions. To investigate their stability, we follow the approach of [5, 12, 13, 19] by considering small perturbations δρ and δS of stationary solutions of density ρ0 and phase S0. Taking variational derivatives of Eqs. (7) and (8), we obtain the time-evolution equations for the perturbations ∂ ----- δS ∂t δρ =–

[v – Ω × r] ⋅ ∇ ∇ ⋅ ρ0 ∇ 2

(23)

g --- ( 1 + ε dd K ) δS ,  δρ [(∇ ⋅ v ) + [v – Ω × r] ⋅ ∇]

∂ where K = – 3 -------2 d x dy dz /(4π| r' – r |) – 1. As in [5, ∂z 19], we consider a polynomial ansatz of order n in the coordinates x, y, and z and evaluate the evolution operator for the perturbations. If one or more of the eigenvalues λ has a positive real component, the stationary solution is dynamically unstable. However, imaginary eigenvalues correspond to stable oscillatory modes of the system [27].

∫

3. VORTEX-LATTICE FORMATION IN CONVENTIONAL ROTATING BECS (εdd = 0) We initially consider the case where there are no dipolar interactions (εdd = 0) and show that vortex-lattice formation can be understood in the context of the stationary hydrodynamical solutions [Eq. (22)] and their stability [Eq. (23)]. Figure 1a shows the stationary solutions in the Ω–α space for various values of trap ellipticity . In the limit of  = 0 (solid line), a nonrotating (α = 0) solution occurs for all Ω, with two additional solutions bifurcating from the α = 0 axis at Ωb( = 0) = ω⊥/ 2 . For finite  (dotted and dashed lines), the α = 0 solution disappears and the plot consists of two distinct branches. The upper branch (positive α) is single valued and exists for all Ω, while the lower branch (negative α) is double valued and exists only when Ω is greater than the bifurcation frequency Ωb(). As  is increased from zero, the branches move away from the α = 0 axis, as can be seen in Fig. 1a. Furthermore, the bifurcation point of the lower branch Ωb() shifts to higher Ω, as shown in Fig. 1b (solid line). The branch diagrams, which have been probed experimentally [2], are key to understanding the response of the system to the adiabatic introduction of trap ellipticity  or rotation frequency Ω. Before any rotation/ellipticity is applied, the BEC has α = 0. When Ω is increased adiabatically (for fixed ), the BEC follows the upper branch, with an increasing α and ellipticity in the density profile. When  is introduced adiabatically (for fixed Ω), the BEC can follow two routes, depending on the value of Ω relative to the bifurcation point Ωb( = 0) [vertical dotted line in Fig. 1b]. For Ω < Ωb( = 0), the lower branch is nonexistent and the BEC follows the upper branch to increasing values of α. For Ω > Ωb( = 0), the lower branch moves from α = 0 to negative α, and the BEC follows this route. However, as  is increased, the edge of the lower branch Ωb() shifts to higher frequencies, and, when Ωb() > Ω, the lower branch no longer exists. In this manner, the evolution of the branches can induce instability, and has been linked to lattice formation [2, 12]. Using Eq. (23) to assess the stability of the upper branch solutions, it is possible to predict dynamical instability when Ω exceeds a critical value Ωi() [5], as indicated in Fig. 1a (circles). The unstable modes have length scales on the order of the condensate size, much greater than the healing length ξ = / mn 0 g , which characterizes the vortex size. In the limit of  = 0, Ωi ≈ 0.78ω⊥. As  is increased, Ωi() is reduced (dashed line in Fig. 1b). Note that the outlying point in Fig. 1a for  = 0.1 (dashed line) at Ω ≈ 0.61ω⊥ is not considered to be in the region of instability due to its narrow width and comparatively small eigenvalues [29]. For the parameter space of interest, the lower branch solution closest to the α = 0 axis is never dynamically unstable. LASER PHYSICS

Vol. 18

No. 3

2008

INSTABILITIES AND VORTEX-LATTICE FORMATION

α(ω⊥)

0.8

(a)

0.4

Ωi(ε)

0 Ωb()

–0.4

–0.8 0.15 (b)

I

II Ωi()

III Ωb()

0.10 ε

A key rotation frequency in our work is ΩX, which is the crossing point of Ωb() and Ωi(), and has the value ΩX ≈ 0.765ω⊥ [vertical dot-dashed line in Fig. 1b]. Based on the solutions to Eq. (22) and their dynamical instability Eq. (23), we can predict the stability of a BEC following the adiabatic introduction of Ω or . However, to determine how the instability manifests itself and whether it leads to lattice formation, we perform time-dependent simulations of the GPE (Eq. (3)) in the laboratory frame. Following previous approaches [7–12] and noting that the solutions to Eqs. (22) and (23) are independent of z [30], we consider a 2D system. The initial state is found by imaginary-time propagation. Then, in real time, either  is ramped up linearly at a rate of d/dt = 10–4ω⊥ (for fixed Ω) or Ω is

325

0.05

2

ramped up linearly at a rate dΩ/dt = 10–2 ω ⊥ (for fixed ). We derive α by fitting the velocity field in a central region [3 × 3] µm to the form of Eq. (9). In such idealized simulations, the trap and BEC can maintain a twofold rotational symmetry to unrealistic levels. This strongly inhibits vortex nucleation since the vortices must enter in pairs. Symmetry-breaking, rather than thermal effects, has been shown to be crucial in simulating lattice formation [10–13]. Indeed, Eq. (23) predicts that only odd modes of the system are dynamically unstable. Previous studies [6, 8] have broken this symmetry through the introduction of thermal noise. To overcome this problem, we shift the trap by 0.5ξ before running in it real time, thereby allowing the excitation of odd modes [28]. We first consider the adiabatic increase of trap ellipticity  for fixed Ω. We discriminate three cases of instability, which each occur in distinct frequency regimes, as indicated in Fig. 1b: (I) ripple instability, (II) interbranch instability, and (III) catastrophic instability. We will discuss each in turn. (I) Ripple instability, Ω < ω⊥ / 2 . The case for Ω = 0.7ω⊥ is shown in Fig. 2a (I). The velocity-field amplitude α (dots) follows the upper branch of Eq. (22) (solid line), for which the BEC and trap axes rotate in phase, as noted experimentally [2, 3]. As  is increased, the BEC moves along the branch to a higher α. However, when  exceeds a critical value (corresponding to when Ω > Ωi ), the solution becomes dynamically unstable according to Eq. (23). This occurs for  ≈ 0.09 in Fig. 2a (I) (vertical dashed line). Subsequently, α (dots) deviates from the upper branch (solid lines), consistent with the onset of dynamical instability. For low , α (dots) features small oscillations due to the centre-of-mass motion caused by the initial offset of the condensate. dev

At a critical ellipticity  cr , low density ripples form on the condensate edge [Fig. 2b (I)], each on the scale of the healing length and featuring a phase singularity (“ghost” vortices [6–8]). These ripples grow in amplitude as  is increased. If  becomes fixed when the ripples have a very low amplitude, they do not grow over LASER PHYSICS

Vol. 18

No. 3

2008

0 0.6

0.7

Ω(ω⊥)

0.8

0.9

Fig. 1. (a) Velocity-field amplitude α of the stationary solutions to Eq. (22) versus the rotation frequency Ω for  = 0 (solid line), 0.02 (dotted line), and 0.1 (dashed line). Regions of dynamical instability for  = 0.02 and 0.1 are shown (circles); (b) phase diagram of  versus Ω. Plotted are the bifurcation point Ωb() (solid line) from Eq. (22), the onset of dynamical instability in the upper branch Ωi() from Eq. (23) (dashed line), and experimental data of Hodby et al. [3] (circles). Furthermore, simulations of Eq. (3) (with µ = 10ω⊥) show the critical ellipticities beyond which the condensate deviates from an elliptical dev

shape  cr

(crosses) and beyond which lattice formation i

ultimately occurs  cr (points with error bars). The bifurcation point for  = 0 (dotted line) and crossing frequency ΩX (dot-dashed line) at which Ωb() = Ωi() are indicated.

the time scales considered. However, once  exceeds a i

second critical value  cr (corresponding to when the ripples have an amplitude on the order of 10%n0), the dynamical instability of Eq. (23) is triggered by the ripples. This instability generates large-scale shape oscillations [Fig. 2c (I)], disrupting the condensate, and enabling “ghost” vortices to nucleate into the condensate, which slowly crystallize into a lattice [Fig. 2d (I)] i

[10]. Once  cr is reached, lattice formation occurs independently of whether  becomes fixed or continuously increased. Surface ripples have been observed experimentally to precede vortex nucleation at this Ω [3]. The gradual growth of the ripples leads to a slow injection of energy/angular momentum into the system, as observed elsewhere at similar Ω [6, 9]. Note that, according to Eq. (23), the dynamical instability on the upper branch couples to the odd modes only. If the symmetry is preserved, we do not expect the instability to develop [6–8].

326

MARTIN et al. Ι(Ω = 0.7 ω⊥) (a)

ΙΙ(Ω = 0.7 ω⊥) (b)

ΙΙΙ(Ω = 0.7 ω⊥) (c)

α (ω⊥)

0.5 0

–0.5 0

0.05

0.10 0

0.05 ε

0.10 0

0.05

ε = 0.088

ε = 0.019

(b)

(b)

ε = 0.113

ε = 0.03

ε = 0.085

(c)

(c)

(c)

ε = 0.194

ε = 0.083

ε = 0.18

(d)

(d)

0.10

ε = 0.07

(b)

(d)

Fig. 2. Dynamics under a continuous increase of  (at a rate of d/dt = 10–4ω⊥) for (I) ripple instability (Ω/ω⊥ = 0.7), (II) interbranch instability (Ω/ω⊥ = 0.75), and (III) catastrophic instability (Ω/ω⊥ = 0.8). (a) Velocity-field amplitude α versus  according to Eq. (22) (solid lines) and GPE simulations (dots). To the right of the dashed line, Eq. (9) is no longer a valid fit to the simulatedvelocity field (the standard deviation of the fit becomes on the order of α). The regions of dynamical instability are indicated (circles). Density snapshots show the (b) onset of instability, (c) disrupted state seeded by the instability, and (d) vortex lattice. Once the instability occurs, the dynamics are qualitatively similar whether  is continuously increased or becomes fixed. Dark/light regions represent high/low density. Each box represents a region (12 × 12) µm (the simulation box is much greater than this). In Ib and IIb, the density scale is limited to 0.1%n0 to highlight low-density features.

ω⊥/ 2 , α (dots) initially follows the lower branch solutions, where the BEC and trap axes are π/2 out of phase. As  is increased, a point is reached when Ω < Ωb(). Here, the lower branch is no longer a solution to Eq. (22). For the example shown, this occurs for  ≈ 0.02 [dashed line in Fig. 2a (II)]. Due to the non-Thomas–Fermi nature of the numerical solutions to Eqs. (7) and (8), α does not perfectly fit the branch solutions to Eq. (22).

an almost circular and highly elliptical shape. If  becomes fixed close to this critical ellipticity, the shape oscillations are stable. However, if  is increased further, the shape oscillations destabilize, with the BEC shedding low-density material at its extrema in a spiral pattern [Fig. 2b (II)]. This situation is closely analogous to when Ω is suddenly turned on, with the fate of the condensate being qualitatively similar [10, 11]. The growth of the ejected material gradually destabilizes the condensate [Fig. 2c (II)], leading to vortex nucleation and, ultimately, a lattice [Fig. 2c (II)]. This is fully consistent with the observations in [2].

When α (dots) reaches the cusp of the lower branch, it deviates nonadiabatically, as observed experimentally [2]. Since Ω < ΩX, the upper branch is dynamically stable at this . The BEC tries to transform to the upper branch, but, without dissipation, it cannot relax to this state. Instead, α (dots) oscillates between positive and negative values, and the density oscillates between

(III) Catastrophic instability, Ω > ΩX. The case for Ω = 0.8ω⊥ is shown in Fig. 2a (III). Again, α (dotted line) follows the lower branch, which ceases to be a solution at some critical  (dashed line). However, since Ω > ΩX, the upper branch is dynamically unstable at this point and no stable solutions exist. The BEC undergoes a quick and catastrophic instability, with α (dots)

(II) Interbranch instability, ω⊥/ 2 < Ω < ΩX. The case for Ω = 0.75ω⊥ is shown in Fig. 2a (II). Since Ω >

LASER PHYSICS

Vol. 18

No. 3

2008

INSTABILITIES AND VORTEX-LATTICE FORMATION α/ω⊥ Ωb/ω⊥ (a) (b) 1.0 0.72 0.6 0.68 0.2 0.64 εdd –0.2 0.60 –0.6 0.58 –1.0 0.5 0.6 0.7 0.8 0.9 1.0 10–2 10–1 100 101 102 Ω/ω⊥ γ

Fig. 3. (a) Irrotational fluid-velocity amplitude α as a function of the trap-rotational frequency Ω, with γ = 1 and  = 0 for εdd = 0 (solid curve), εdd = 0.5 (dashed curve), and εdd = 0.99 (dotted curve); (b) bifurcation point Ωb versus γ for εdd = 0, 0.2, 0.4. 0.6, 0.8, and 0.99; εdd increases in the direction of the arrow.

deviating rapidly from the rotating solutions to Eq. (22). The BEC density becomes strongly contorted into a spiral shape [Fig. 2b (III)]. The arms of the spiral collapse inwards and trap phase singularities to form vortices. Energy and angular momentum are very rapidly injected into the BEC (in contrast to the ripple and interbranch instabilities), as observed in [6, 9] in this frequency regime. The BEC becomes highly excited and turbulent [Fig. 2c (III)], with the structure on the length scales less than the healing length. Although we observe this state to ultimately settle into a lattice [Fig. 2c (II)], one may question the validity of our zerotemperature simulations for such a “heated” state. Indeed, for Ω  0.78ω⊥, turbulent states, rather than vortex lattices, were observed experimentally [1]. In the time-dependent GPE simulations, we have dev measured two distinct critical ellipticities:  cr [crosses in Fig. 1b] is when, for a continuously increasing , we observe the density deviating from a smooth ellipse (on i the level of 0.1%n0);  cr [points with error bars in Fig. 1b] is when, for  ramped up to some final value, instability and lattice formation ultimately occur. In dev i regime I, typically,  cr ≤  cr , since surface ripples are generated which are stable for a narrow range of , above which they induce instability and lattice formadev i tion. In regimes II and III,  cr ≈  cr , indicating the relatively sudden onset of this instability once the density deviates from a smooth ellipse. According to the GPE simulations, the region above the points in Fig. 1b is unstable, leading to lattice formation. The prediction of Eq. (23) (dashed line) gives reasonable agreement with the GPE results in region I, while Eq. (22) gives excellent agreement in region III. Also plotted in Fig. 1b are the experimental results of Hodby et al. [3] (circles). The GPE results give good LASER PHYSICS

Vol. 18

No. 3

2008

α/ω⊥ 1.0 0.6

327 Re(λ/ω⊥) 0.25 (b) 0.20

(a)

0.2

0.15

–0.2

0.10

–0.6

0.05

–1.0 0.4

0.6

0.8

1.0

0 0.65 0.75

Ω/ω⊥

εdd

0.85

0.95

Fig. 4. (a) α vs. Ω, as in Fig. 1a but for  = 0.02. (b) The maximum positive real eigenvalues of Eq. (17) (solid curves) as a function of Ω, for  = 0.02, γ = 1, n = 3, and εdd = 0, 0.2, 0.4, 0.6, 0.8, 0.95, and 0.98; εdd increases in the direction of the arrow. The short and long dashed curves are additional positive eigenvalue solutions for εdd = 0.95 and 0.98, respectively.

agreement with the experimental data throughout, with the agreement being particularly good in region III. So far, we have considered the adiabatic introduction of ellipticity for fixed . However, the results in Fig. 1b also apply to the case when Ω is introduced adiabatically for fixed . Here, the condensate follows the upper branch until it becomes dynamically unstable to the ripple instability. According to Eq. (23), dynamical instability occurs when Ω > Ωi() [dashed line in Fig. 1b]. In [2], for  = 0.025, vortex-lattice formation was observed to occur when Ω  0.75ω⊥, which agrees well with our GPE simulations and Fig. 1b. 4. HYDRODYNAMICAL SOLUTIONS IN THE PRESENCE OF DIPOLAR INTERACTIONS (εdd > 0) We now show how the introduction of long-range dipolar interactions can significantly alter the hydrodynamical solutions and their stability. 4.1. Circular Trap ( = 0) Figure 3a shows the solutions to Eq. (22) for various values of εdd with γ = 1 and  = 0. For εdd = 0 (solid curve), where the static solutions are independent of both the strength of the s-wave interactions (g) and aspect ratio of the trap (γ), we reproduce the results in Section 2. It is a remarkable feature of the pure s-wave case that these solutions do not depend upon g. This is because, in the Thomas–Fermi limit, surface excitations with angular momentum l = ql R, where R is the Thomas–Fermi radius and ql is the quantized wave 2

number, obey the classical dispersion relation ω l = (ql /m)∇RV involving the local harmonic potential V = 2

m ω ⊥ R2/2 evaluated at R [31]. Consequently, ωl =

328

MARTIN et al.

lω ⊥ , which is independent of g. However, in the case of long-range dipolar interactions, the potential Φdd of Eq. (5) gives nonlocal contributions, breaking the simple dependence of the force –∇V upon R [20]. Thus, we expect the resonant condition for exciting the quadrupolar mode to change with εdd, resulting in a shift of Ωb. In Fig. 3a, we see that, as the dipole interactions are introduced, the bifurcation point Ωb moves to lower frequencies. In contrast to the s-wave case, the shape of the BEC determines the potential Φdd. For an oblate (κx, y > 1) BEC, more dipoles lie side-by-side, giving a nonrepulsive interaction, in comparison to the prolate (κx, y < 1) case, where a majority sit end-to-end, inducing a net-attractive interaction. To investigate how this affects the static solutions, we consider a small perturbation to the shape of the dipolar BEC at the bifurcation point of the form κx = κ(1 + δ) and κy = κ(1 – δ). From this, we find that the frequency of the bifurcation point is given by 2

Ω ------b = ω⊥

κ β3 – β3 ---, 3 ---, 2 2 1 3 2 2 2 -, --- + --- κ ε dd γ -----------------------------------------------2 4 9 2 ⎛ ⎞ 1 – ε dd 1 – --- κ β 5 ⎝ ---, 1⎠ 2 2

(24)

where ∞

β a, b =

dσ

∫ (----------------------------------------1 + σ) (κ + σ) a

(25)

b

2

0

and κ is defined by [20] f (κ) 2 γ 3κ ε dd ⎛ ----- + 1⎞ --------------2 – 1 ⎝2 ⎠1 – κ 2

2

(26)

2

+ ( ε dd – 1 ) ( κ – γ ) = 0, where 2

2 + κ [ 4 – 3Ξ ( κ ) ] -, f ( κ ) = ------------------------------------------2 2(1 – κ )

(27)

with Ξ(κ) being defined by the shape of the BEC [26] 2

1 1+ 1–κ Ξ ( κ ) = ------------------ ln --------------------------2 2 1–κ 1– 1–κ 2 2 Ξ ( κ ) = ------------------ arctan κ – 1 2 κ –1

for

for

κ < 1 , (28)

κ > 1.

(29)

In Fig. 3b, we plot Ωb [Eq. (24)] as a function of γ for various values of εdd. For εdd = 0, we find that the bifurcation point remains unaltered at Ωb = ωx / 2 as γ =

ωz /ωx is changed [4, 5]. As εdd is increased, the value of γ for which Ωb is a minimum changes from a trap shape which is oblate (γ > 1) to prolate (γ < 1). 4.2. Elliptical Trap ( > 0) In Fig. 4a, we have plotted the solutions to Eq. (22) for various values of εdd with γ = 1 and  = 0.02. As in the case without dipolar interactions [4, 5, 12], the solution α = 0 is no longer a solution for all Ω. The effect of introducing the anisotropy, in the absence of dipolar interactions, is to increase the bifurcation frequency Ωb. Turning on the dipolar interactions, as in the case of  = 0, reduces the bifurcation frequency. In Fig. 4b, we have plotted the real positive eigenvalues Re(λ) of Eq. (23) for the upper branch solution to Eq. (22) as a function of Ω for various values of εdd with n = 3. For higher values of εdd [0.95 and 0.98 in Fig. 4b], there can be more than one real positive eigenvalue; thus, we define the region of instability as the range over which max[Re(λ) > 0], as shown by the solid curves in Fig. 2b [32]. As the dipolar interaction strength is increased, the lower bound in Ω for the unstable region is reduced. For example, for εdd = 0.6, the range of rotation frequencies where the upper branch solution is unstable is [0.75ωx, ωx]; this increases to [0.67ωx, ωx] for εdd = 0.98. In analogy to the case for conventional BECs (Section 2), for an adiabatic increase in the trap ellipticity with a fixed trap-rotation frequency, we identify three regimes for inducing an instability in the dipolar BEC. (I) Ripple instability, Ω < Ωb ( = 0). As for the case of conventional BECs, this instability occurs when the upper branch solutions become unstable. The introduction of dipolar interactions, for a fixed Ω, reduces the value of  at which this instability occurs. Thus, the introduction of dipolar interactions will have the effect of reducing the trap anisotropy at which a ripple instability will occur. (II) Interbranch and (III) Catastrophic instabilities, Ω > Ωb ( = 0). As for conventional BECs, as  is increased adiabatically from zero, the α = 0 solution moves to negative values of α and the BEC follows this route. However, as  is increased further, the edge of the lower branch Ωb() shifts to higher frequencies. At some critical value of , Ωb() = Ω, the lower branch ceases to be a solution for this value of Ω. As the dipole interactions are increased, the bifurcation frequency is reduced and the range of Ω for which this type of instability can occur increases from [ω⊥/ 2 , ω⊥] to [0.5ω⊥, ω⊥]. Thus, dipolar interactions increase the value of  for which lower branch solutions exist. LASER PHYSICS

Vol. 18

No. 3

2008

INSTABILITIES AND VORTEX-LATTICE FORMATION

5. DISCUSSION 5.1. Conventional BECs We have shown that vortex-lattice formation is inherently a two-dimensional and zero-temperature effect. We have theoretically mapped out the condensate stability as a function of rotation frequency Ω and trap ellipticity , with our interpretation being consistent with previous experimental and theoretical results. Specifically, for fixed Ω and the adiabatic introduction of , we find three distinct regimes of instability—ripple (Ω < ω⊥/ 2 ), interbranch (ω⊥/ 2 < Ω < ΩX), and catastrophic (Ω > ΩX). Each instability manifests itself in a characteristic manner, which could be observed experimentally. Ultimately, in each case, the instability seeds vortex-lattice formation. This formation process not only relies on the presence of an instability, but is crucially dependent on the breaking of the twofold rotational symmetry of the system, as inevitably occurs experimentally. 5.2. Dipolar BECs By calculating the static hydrodynamic solutions of a rotating dipolar BEC and studying their dynamical stability, we have predicted the regimes of instability of the condensate. By relating to analogous work for conventional BECs, we identify how the regimes of ripple, interbranch, and catastrophic instability change as dipolar interactions are introduced. In general, we find that the bifurcation frequency Ωb decreases in the presence of dipolar interactions and the upper branch solutions become unstable at lower rotation frequencies. Thus, for a fixed Ω [at Ω < Ωb( = 0)] and an adiabatic increase in , the critical anisotropy at which we expect a ripple instability to occur will be lower than for a conventional BEC. For a fixed Ω [at Ω > Ωb( = 0)] and an adiabatic increase in , the critical anisotropy at which we expect interbranch and catastrophic instabilities to occur will be higher than for a conventional BEC. Furthermore, we find that the size of this shift not only depends on the strength of the dipolar interactions, but also on the aspect ratio of the trap. These shifts can be understood in the context that the introduction of dipolar interactions shifts both the dashed line [Ωi()] and the solid line [Ωb()] in Fig. 1b to lower frequencies, with the size of the shift being controlled by the strength of the dipolar interactions and the aspect ratio of the trap. In conventional BECs, these instabilities relate to vortex-lattice formation. This occurs primarily because the instability disrupts the BEC at an Ω which is greater than the rotation frequency at which it is energetically favorable to have a vortex state [33]. However, in a prolate trap, the rotational frequency at which it is energetically favorable to form a vortex in a dipolar BEC grows rapidly as εdd is increased [34] and can exceed the frequency at which we expect an instability to LASER PHYSICS

Vol. 18

No. 3

2008

329

occur. The final state under these circumstances warrants further investigation. ACKNOWLEDGMENTS We thank Y. Castin, C.S. Adams, and C.J. Foot for helpful discussions and C.J. Foot for the use of experimental data in Fig. 1b. We acknowledge funding from the ARC. REFERENCES 1. K. W. Madison, F. Chevy, W. Wohlleben, and J. Dalibard, Phys. Rev. Lett. 84, 806 (2000). 2. K. W. Madison, F. Chevy, V. Bretin, and J. Dalibard, Phys. Rev. Lett. 86, 4443 (2001). 3. E. Hodby et al., Phys. Rev. Lett. 88, 010405 (2001). 4. A. Recati, F. Zambelli, and S. Stringari, Phys. Rev. Lett. 86, 377 (2001). 5. S. Sinha and Y. Castin, Phys. Rev. Lett. 87, 190402 (2001). 6. C. Lobo, A. Sinatra, and Y. Castin, Phys. Rev. Lett. 92, 020403 (2004). 7. M. Tsubota, K. Kasamatsu, and M. Ueda, Phys. Rev. A 65, 023603 (2002); K. Kasamatsu, M. Tsubota, and M. Ueda, Phys. Rev. A 67, 033610 (2003). 8. A. A. Penckwitt, R. J. Ballagh, and C. W. Gardiner, Phys. Rev. Lett. 89, 260402 (2002). 9. E. Lundh, J.-P. Martikainen, and K.-A. Suominen, Phys. Rev. A 67, 063604 (2003). 10. N. G. Parker and C. S. Adams, Phys. Rev. Lett. 95, 145301 (2005). 11. N. G. Parker and C. S. Adams, J. Phys. B 39, 43 (2006). 12. N.G. Parker, R. M. W. van Bijnen, and A. M. Martin, Phys. Rev. A 73, 061603(R) (2006). 13. I. Corro, N. G. Parker, and A. M. Martin, J. Phys. B 40, 3615 (2007). 14. J. R. Abo-Shaeer, C. Raman, and W. Ketterle, Phys. Rev. Lett. 88, 070409 (2002). 15. A. Griesmaier et al., Phys. Rev. Lett. 94, 160401 (2005). 16. N. R. Cooper, E. H. Rezayi, and S. H. Simon, Phys. Rev. Lett. 95, 200402 (2005). 17. J. Zhang and H. Zhai, Phys. Rev. Lett. 95, 200403 (2005). 18. S. Yi and H. Pu, Phys. Rev. A 73, 061602(R) (2006). 19. R. M. W. van Bijnen, D. H. J. O’Dell, N. G. Parker, and A. M. Martin, Phys. Rev. Lett. 98, 150401 (2007). 20. D. H. J. O’Dell, S. Giovanazzi, and C. Eberlein, Phys. Rev. Lett. 92, 250401 (2004). 21. M. Marinescu and L. You, Phys. Rev. Lett. 81, 4596 (1998). 22. S. Yi and L. You, Phys. Rev. A 61, 041604(R) (2000). . 23. K. Góral, K. Rza¸z ewski, and T. Pfau, Phys. Rev. A 61, 051601(R) (2000). 24. L. Santos, G.V. Shlyapnikov, P. Zoller, and M. Lewenstein, Phys. Rev. Lett. 85, 1791 (2000). 25. The limit where the zero-point kinetic energy (quantum pressure) is negligible compared to the potential and interaction energies.

330

MARTIN et al.

26. C. Eberlein, S. Giovanazzi, and D. H. J. O’Dell, Phys. Rev. A 71, 033618 (2005). 27. Y. Castin, in Coherent Matter Waves, Lecture Notes of Les Houches Summer School, Ed. by R. Kaiser, C. Westbrook, and F. David (Springer-Verlag, Berlin, 2001), p. 1–136. 28. Over long times, symmetry can break intrinsically through the stimulated growth of numerical noise [10, 11]. 29. Numerical simulations do not indicate any dynamical instability within these narrow, low eigenvalue regions.

30. As in [5] we find that the regions of instability according to Eq. (23) are independent of ωz. 31. L. Pitaevskii and S. Stringari, Bose–Einstein Condensation (Oxford Univ. Press, London, 2003). 32. For n > 3 we find that although there are more eigenvalues the region of instability, defined by max[Re(λ) > 0], remains the same. 33. E. Lundh, C.J. Pethick and H. Smith, Phys. Rev. A 55, 2126 (1997). 34. D. H. J. O’Dell and C. Eberlein, Phys. Rev. A. 75, 013604 (2007).

LASER PHYSICS

Vol. 18

No. 3

2008