Australian Institute of Physics 17th National Congress 2006 – Brisbane, 3-8 December 2006 RiverPhys

INSTABILITIES LEADING TO VORTEX LATTICE FORMATION IN ROTATING BOSEEINSTEIN CONDENSATES A.M. MartinA, R.M.W. van BijnenA,B, A.J. DowA N.G. ParkerA and D.H.J. O’Dell,C A School of Physics, University of Melbourne, Parkville, Victoria 3010, Australia. B Eindhoven University of Technology, PO Box 513, 5600 MB Eindhoven, The Netherlands. C Centre for Cold Matter, Imperial College, Prince Consort Road, London, SW7 2BW, United Kingdom.

Abstract We calculate the hydrodynamic solutions for a dilute Bose-Einstein condensate with and without long-range dipolar interactions in a rotating, elliptical harmonic trap, and analyse their dynamical stability. When the trap anisotropy is introduced adiabatically at fixed rotation frequency, Ω, the upper bound for the unstable region in Ω is decreased by the dipolar interactions. Furthermore, for fixed trap anisotropy and adiabatically increasing rotational frequency the dynamically unstable region in Ω widens non-trivially as the strength of the dipole interactions is increased. We show, through numerical simulations of the Gross-Pitaevskii equation, that in the absence of dipolar interactions these unstable regions correspond to rotational frequencies at which vortex lattices are formed. Introduction In recent years a considerable amount of experimental [Madison00, Madison01, Hodby02] and theoretical [Recati01, Madison01, Sinha01, Tsubota02, Kasamatsu03, Lundh03, Lobo04, Parker05, Parker06a, Parker06b] work has been carried out on dilute Bose-Einstein condensates (BECs) in rotating anisotropic traps. Where short-range interactions dominate, a vortex lattice forms when the rotational frequency (Ω) of the system is ≈0.7ωr (where ωr is the trapping frequency perpendicular to the axis of rotation). Insight into the mechanism of vortex formation can be gained by noting that 0.7ωr closely coincides with the frequency at which certain hydrodynamic surface excitations become unstable [Sinha01, Parker06b]. Through comparison with experimental results [Madison00, Madison01, Hodby02] and numerical solutions of the Gross-Pitaevskii equation (GPE) [Lobo04, Parker05, Parker06a, Parker06b] such instability has been directly related to the formation of a vortex lattice. The above results apply to conventional BECs composed of atoms of mass m with short-range s-wave interactions, parameterised via g= h2 a/(πm), where a is the s-wave scattering length. However, a recent experiment has formed a BEC of chromium atoms with dipolar interactions [Griesmaier05]. This opens the door to experimentally study the effect of dipolar interactions in BECs. Parallel theoretical work, using a modified GPE, has studied the effect of such long-range interactions on the ground state vortex lattice solutions [Cooper05, Zhang05, Yi06]. However, the route to generating such states has not been explored. For this purpose we solve the hydrodynamic equations of motion for a dipolar BEC in rotating anisotropic harmonic traps [Bijnen06]. We show that the solutions depend on both the strength of the dipolar interactions, εdd, and the aspect ratio of the trap, γ = ωz /ωr, in stark contrast to conventional BECs where the instability is independent of both the strength of the interactions, g, and γ [Recati01, Sinha01]. In addition we evaluate the dynamical stability of our solutions, showing that the unstable region of Ω can be controlled via εdd and γ. Below we briefly present the methodology of the hydrodynamical model in the presence dipolar interactions. We then consider the case of only short range interactions, specifically comparing the solutions from a hydrodynamical approach to solutions of the GPE [Parker06b]. We then examine the solutions and their stability of the hydrodynamical model in the presence of long-range dipolar interactions. Hydrodynamical Model 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 ê is given by, in the notation of [O’Dell04],

U dd (r ) =

(δ ij − 3rˆi rˆj ) C dd eˆi eˆ j . 4π r3

(1)

Denoting ρ as the condensate density, the dipolar interactions give rise to a mean-field potential

Φ dd (r ) = ∫ d 3 r U dd (r − r ')ρ (r ') Paper No. WC0091

(2)

Australian Institute of Physics 17th National Congress 2006 – Brisbane, 3-8 December 2006 RiverPhys

which can be included in a generalized GPE [Góral00, Yi00, Santos00] for the BEC. In the Thomas-Fermi limit for a dipolar BEC in a harmonic trapping potential V(r)=m(ωx2x2 + ωy2y2 + ωz2z2)/2 rotating at Ω about the z-axis, it has been shown [Bijnen06] that for an irrotational quadrupolar velocity field, v = α∇xy , that the stationary solutions of the system are given by

⎛ ⎛ ω x2κ xκ y γ 2 ⎞ ω x2κ x κ y γ 2 ⎞ 3 3 0 = (α + Ω )⎜ ω~ x2 − ε dd β x ⎟ + (α − Ω )⎜ ω~ y2 − ε dd βy ⎟ ⎟ ⎟ ⎜ ⎜ 2 ς 2 ς ⎠ ⎠ ⎝ ⎝

(3)

where ε dd = C dd / 3 g characterises the strength of the dipolar interactions, κ x , κ y , ς , β x , β y are defined in [Bijnen06],

ω~x2 = ω r2 + α 2 − 2αΩ and ω~ y2 = ω r2 + α 2 + 2αΩ .

To investigate the stability of the stationary solutions of α [Sinha01, Bijnen06] we consider the dynamics of small perturbations to the density ( δρ ) and phase ( δS ) of the BEC

⎡(v − Ω × r ) ⋅ ∇ ∂ ⎡δS ⎤ = −⎢ ⎢ ⎥ ∂t ⎣δρ ⎦ ⎣ ∇ ⋅ ρ 0∇

g (1 − ε dd ) / m

⎤ ⎡δS ⎤ . [(∇ ⋅ v ) + (v − Ω × r ) ⋅ ∇]⎥⎦ ⎢⎣δρ ⎥⎦

(4)

The eigenfunctions of Eq. (4) grow in time as exp(λ t), where λ is the corresponding eigenvalue. As in Refs. [Sinha01, Parker06b, Bijnen06] we consider a polynomial ansatz, of order n in the coordinates x,y,z, for the perturbations, and evaluate the evolution operator for the perturbations: dynamical instability arises when one or more eigenvalues has a positive real part with the imaginary values relating to stable oscillatory modes. Conventional BECs (εdd = 0) Based on the stationary solutions of Eq. (3) and their dynamical stability, Eq. (4), we can predict the stability of a BEC following the adiabatic introduction of trap ellipticity ε (where ωx = (1- ε)1/2 ωr and ωy = (1+ ε)1/2 ωr) for a fixed trap rotation frequency Ω. However, to determine how the instability manifests itself and whether it leads to lattice formation we perform time dependent simulations of the GPE in two dimensions. The initial state is found by imaginary-time propagation. Then, in real time, ε is ramped up adiabatically. We monitor the evolution of α by fitting the velocity field in a central region to v = α∇xy . In doing this we find three distinct regimes of instability, (I) ripple, (II) intermediate and (III) catastrophic, all leading ultimately to vortex lattice formation [Parker06b]. (I) Ripple instability, Ω < ωr / 21/2. The case for Ω ≈ 0.7ωr is shown in Fig. 1[I(a)]. The velocity field amplitude α (dots) follows the stationary solutions (in the rotating frame) of Eq. (3) (solid curve). However, when the ellipticity exceeds a critical value the solution becomes dynamically unstable, according to Eq. (4) (dynamically unstable solutions denoted by red circles). Subsequently α (dots) deviates from the rotating solutions of Eq. (3) (solid curve), consistent with the onset of dynamical instability. At the critical ellipticity low density ripples form on the condensate edge, Fig. 1[I(b)], each on the scale of the healing length and featuring a phase singularity (ghost vortices [Tsubota02, Lobo04]). These ripples grow in amplitude as ε is increased. Once ε exceeds a second critical value, the dynamical instability of Eq. (4) is triggered by the ripples. This instability generates largescale shape oscillations, Fig. 1[I(c)], disrupting the condensate, and enabling ghost vortices to nucleate into the condensate, which slowly crystallise into a lattice Fig. 1[I(d)] [Parker06b]. (II) Interbranch instability, ωr / 21/2 < Ω < 0.765 ωr. The case for Ω ≈ 0.75ωr is shown in Fig. 1[II(a)]. Instead of only having one solution to Eq. (3) for α we now have three, at ε=0. From the numerical solutions of the GPE we find that for Ω > ωr / 21/2 our solution follows the lower branch (dots). As ε is increased a point is reached where the lower branch solutions cease to exist, ε ≈ 0.02 (vertical dashed line in Fig. 1[II(a)]). When our numerical solution reaches the cusp of the lower branch it deviates non-adiabatically, as observed experimentally [Madison01]. At this point the upper branch is dynamically stable and the condensate tries to transform to this branch, but without dissipation it cannot relax to this state. As ε is increased beyond this point the condensate sheds low density material at its extrema in a spiral pattern, Fig. 1[II(b)]. The growth of the ejected material gradually destabilises the condensate, Fig. 1[II(c)], leading to vortex nucleation and ultimately a lattice, Fig. 1[II(c)], consistent with the observations in [Madison01]. (III) Catastrophic instability, Ω > 0.765 ωr. The case for Ω ≈ 0.8ωr is shown in Fig. 1[III(a)]. From the numerical solutions to the GPE we find α (dotted line) follows the lower branch, which ceases to be a solution at some critical ε (vertical dashed line). However, the upper branch is dynamically unstable at this point. The BEC undergoes a quick and Paper No. WC0091

Australian Institute of Physics 17th National Congress 2006 – Brisbane, 3-8 December 2006 RiverPhys

catastrophic instability, with α (dots) deviating rapidly from the rotating solutions of Eq. (3). The BEC density profile becomes strongly contorted into a spiral shape, Fig. 1[III(b)]. 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). The BEC becomes highly excited and turbulent, Fig. 1[III(c)], with structure on length-scales less than the healing length. Figure 1. Dynamics under an adiabatic increase of ε for (I) ripple (Ω ≈ 0.7ωr), (II) interbranch (Ω ≈ 0.75ωr) and (III) catastrophic (Ω ≈ 0.8ωr) instability. (a) α versus ε according to Eq. (3) (green solid curves), with the blue circles denoting unstable solutions for α as determined from Eq. (4) and time dependent GPE simulations (black dots). Density snap shots, showing the (b) onset of instability (c) disrupted state and (d) vortex lattice.

Dipolar BECs (εdd > 0) In Fig. 2(a) we plot the solutions for α with ε = 0, γ = 1 and a range of values of εdd as a function of the rotational frequency of the trap, Ω [Bijnen06]. For εdd = 0 (solid curve) we regain the results of [Recati01, Sinha01] with a bifurcation point at Ωb = ωx / 21/2 (independent of the strength of the s-wave interactions and the aspect ratio of the trap, γ) which exactly coincides with the vanishing of the energy of the quadrupolar mode in the rotating frame. The two additional solutions, for Ω > Ωb are a consequence of the quadrupole mode being excited. In Fig. 1(a) we see that as the dipole interactions are introduced the bifurcation point Ωb moves to lower frequencies. Considering small a perturbation to the shape of the dipolar BEC at the bifurcation point of the form κx = κ(1 + δ) and κy = κ(1 – δ)

Ωb

ωx where β a , b =

∞

⎡ ⎤⎡ ⎛ 9 ⎞⎤ 1 3 2 + κ ε dd γ 2 ⎢κ 2 β 3 − β 3 ⎥ ⎢1 − ε dd ⎜⎜1 − κ 2 β 5 ⎟⎟⎥ ,3 ,2 ,1 2 4 2 2 ⎦⎢ 2 ⎠⎥ ⎣ ⎝ 2 ⎣ ⎦

=

dσ

∫ (1 + σ ) (κ a

0

2

+σ

)

b

(5)

, and κ is defined by [O’Dell04]

⎡⎛ γ 2 ⎤ ⎞ f (κ ) 2 2 + 1⎟⎟ − 3κ 2 ε dd ⎢⎜⎜ 1 ⎥ + (ε dd − 1) κ − γ = 0 2 ⎠1−κ ⎣⎝ 2 ⎦

(

where f (κ ) ≡

−1

)

(6)

2 + κ 2 [4 − 3Ξ(κ )] , with Ξ(κ) being defined by the shape of the BEC [Eberlein05]. 2 1−κ 2

(

)

Consider now the effect of finite trap anisotropies (ε > 0). In Fig. 2(b) we have plotted the solutions to Eq. (3) for various values of εdd = 0 with γ = 1 and ε = 0.02. As in the case without dipolar interactions [Recati01, Sinha01] α=0 is no longer Paper No. WC0091

Australian Institute of Physics 17th National Congress 2006 – Brisbane, 3-8 December 2006 RiverPhys

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. Finally, to quantify the regions of dynamical instability we have plotted in Fig. 2(c) the maximal real positive eigenvalues of Eq. (4) for the upper branch (α > 0) stationary solutions of Eq. (3). We observe that the introduction of dipolar interactions can significantly increase the range of dynamical instability [Bijnen06].

1.0

1.0

0.6

0.6

0.8

0.2

0.6

0.2

-0.2 -0.6 -1.0 0

εdd

1.0

α/ωx

α/ωx

Figure 2. (a,b) α as a function of Ω, as obtained from Eq. (3), for γ = 1, ε = 0 (a), ε = 0.02 (b), and εdd = 0 (solid curve), εdd = 0.25 (short-dashed curve), εdd = 0.5 (long-dashed curve), εdd = 0.75 (dashed-dotted curve) and εdd = 0.99 (dotted curve). (c) Density plot of the maximal real positive eigenvalues as obtained from Eq. (4) as a function of Ω and εdd for ε = 0.02 and γ = 1 (white = 0 and black = high).

-0.2 -0.6

(a) 0.6

0.8

Ω/ωx

1.0

-1.0

0.4 0.2

(b) 0 0.2 0.4 0.6 0.8 1.0

Ω/ωx

(c)

0 0.4

0.6

0.8

Ω/ωx

1.0

Conclusions We have shown that in the absence of dipolar interactions, through the comparison of numerical solutions of the GPE and the stationary hydrodynamic solutions of an irrotational quadrupolar velocity field, that vortex lattice formation is related to the stability of our stationary solutions in the rotating frame. It is important to note that to gain these results it was imperative [Parker05, Parker06a, Parker06b] that we broke the two fold rotational symmetry in our numerical simulations of the GPE. In addition, we have shown that it is possible to find the equivalent solutions and their stability for a dipolar BEC, noting in general that the introduction of dipolar interactions reduces the critical rotational frequency at which we expect instability to occur. It is an as yet an unanswered question as to whether this instability will lead to vortex lattice formation, since it has recently been shown [O’Dell06], for γ < 1, that Ω at which the instability occurs can be less than the critical rotational frequency at which it becomes energetically favourable for a vortex to enter the BEC. Acknowledgments We acknowledge the financial support of the ARC, the EPSRC and QIP IRC (GR/S82176/01). References Bijnen, R.M.W. van et al. (2006). Dynamical instability of a rotating dipolar Bose-Einstein condensate. cond-mat/0602572, accepted for publication in Physical Review Letters. Cooper, N.R. et al. (2005). Vortex lattices in rotating atomic Bose gases with dipolar interactions. Physical Review Letters 95, 200402. Eberlein, C. et al. (2005). Exact solution of the Thomas-Fermi equation for a trapped Bose-Einstein condensate with dipole-dipole interactions. Physical Review A 71, 033618. Góral, K. et al. (2000). Bose-Einstein condensation with magnetic dipole-dipole forces. Physical Review A 61, 051601. Griesmaier, A. et al. (2005). Bose-Einstein condensation of chromium. Physical Review Letters 94, 160401. Hodby, E. et al. (2002). Vortex nucleation in Bose-Einstein condensates in an oblate purely magnetic potential. Physical Review Letters 88, 010405. Kasamatsu, K. et al. (2003). Nonlinear dynamics of vortex lattice formation in a rotating Bose-Einstein condensate. Physical Review A 67, 033610. Lobo, C. et al. (2004). Vortex lattice formation in Bose-Einstein condensates. Physical Review Letters 92, 020403. Lundh, E. et al. (2003). Vortex nucleation in Bose-Einstein condensates in time-dependent traps. Physical Review A 67, 063604. Madison, K.W. et al. (2000). Vortex formation in a stirred Bose-Einstein condensate. Physical Review Letters 84, 806. Madison, K.W. et al. (2001). Stationary states of a rotating Bose-Einstein condensate: routes to vortex nucleation. Physical Review Letters 86, 4443. O'Dell, D.H.J. et al. (2004). Exact hydrodynamics of a trapped dipolar Bose-Einstein condensate. Physical Review Letters 92, 250401. O'Dell, D.H.J. and Eberlein, C. (2006). Vortex in a trapped Bose-Einstein condensate with dipole-dipole interactions. cond-mat/0608316, accepted for publication in Physical Review A. Parker, N.G. and Adams, C.S. (2005). Emergence and decay of turbulence in stirred atomic Bose-Einstein condensates. Physical Review Letters 95, 145301 Parker, N.G. and Adams, C.S. (2006a). Response of an atomic Bose-Einstein condensate to a rotating elliptical trap. Journal of Physics B 39, 43. Parker, N.G. et al. (2006b). Instabilities leading to vortex lattice formation in rotating Bose-Einstein condensates. Physical Review A 73, 061603(R). Recati, A. et al. (2001). Overcritical rotation of a trapped Bose-Einstein condensate. Physical Review Letters 86, 377. Santos, L. et al. (2000). Bose-Einstein condensation in trapped dipolar gasses. Physical Review Letters 85, 1791. Sinha, S. and Castin, Y. (2001). Dynamic instability of a rotating Bose-Einstein condensate. Physical Review Letters 87, 190402. Tsubota, M. et al. (2002). Vortex lattice formation in a rotating Bose-Einstein condensate. Physical Review A 65, 023603. Yi, S. and Pu, H. (2006). Vortex structures in dipolar condensates. Physical Review A 73, 061602(R). Yi, S. and You, L. (2000). Trapped atomic condensates with anisotropic interactions. Physical Review A 61, 041604. Zhang, J. and Zhai, H. (2005). Vortex lattices in planar Bose-Einstein condensates with dipolar interactions. Physical Review Letters 95, 200403. Paper No. WC0091