Exact solutions and stability of rotating dipolar Bose-Einstein
condensates in the Thomas-Fermi limit
Citation for published version (APA):
Bijnen, van, R. M. W., Dow, A. J., O'Dell, D. H. J., Parker, N. G., & Martin, A. M. (2009). Exact solutions and stability of rotating dipolar Bose-Einstein condensates in the Thomas-Fermi limit. Physical Review A : Atomic, Molecular and Optical Physics, 80(3), 033617. [033617]. https://doi.org/10.1103/PhysRevA.80.033617
DOI:
10.1103/PhysRevA.80.033617 Document status and date: Published: 01/01/2009
Document Version:
Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)
Please check the document version of this publication:
• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.
• The final author version and the galley proof are versions of the publication after peer review.
• The final published version features the final layout of the paper including the volume, issue and page numbers.
Link to publication
General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain
• You may freely distribute the URL identifying the publication in the public portal.
If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:
www.tue.nl/taverne Take down policy
If you believe that this document breaches copyright please contact us at: openaccess@tue.nl
Exact solutions and stability of rotating dipolar Bose-Einstein condensates
in the Thomas-Fermi limit
R. M. W. van Bijnen,1,2,3A. J. Dow,1D. H. J. O’Dell,3N. G. Parker,3,4and A. M. Martin1
1
School of Physics, University of Melbourne, Parkville, Victoria 3010, Australia 2
Department of Applied Physics, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands 3
Department of Physics and Astronomy, McMaster University, Hamilton, Ontario, Canada L8S 4M1 4
School of Food Science and Nutrition, University of Leeds, Leeds LS2 9JT, United Kingdom 共Received 14 May 2009; published 17 September 2009兲
We present a theoretical analysis of dilute gas Bose-Einstein condensates with dipolar atomic interactions under rotation in elliptical traps. Working in the Thomas-Fermi limit, we employ the classical hydrodynamic equations to first derive the rotating condensate solutions and then consider their response to perturbations. We thereby map out the regimes of stability and instability for rotating dipolar Bose-Einstein condensates and, in the latter case, discuss the possibility of vortex lattice formation. We employ our results to propose several routes to induce vortex lattice formation in a dipolar condensate.
DOI:10.1103/PhysRevA.80.033617 PACS number共s兲: 03.75.Kk, 34.20.Cf, 47.20.⫺k
I. INTRODUCTION
The successful Bose-Einstein condensation of52Cr atoms 关1–3兴 enables the realization of Bose-Einstein condensates
共BECs兲 with significant dipole-dipole interactions. These long-range and anisotropic interactions introduce rich physi-cal effects, as well as opportunities to control BECs. A basic example is how dipole-dipole interactions modify the shape of a trapped BEC. In a prolate共elongated兲 dipolar gas with the dipoles polarized along the long axis the net dipolar in-teraction is attractive, whereas for an oblate共flattened兲 con-figuration with the dipoles aligned along the short axis the net dipolar interaction is repulsive. As a result, in comparison to s-wave BECs共which we define as systems in which atom-atom scattering is dominated by the s-wave channel兲, a di-polar BEC elongates along the direction of an applied di- polar-izing field关4,5兴.
A full theoretical treatment of a trapped BEC involves solving the Gross-Pitaevskii equation共GPE兲 for the conden-sate wave function 关6,7兴. The nonlocal nature of the
mean-field potential describing dipole-dipole interactions means that this task is significantly harder for dipolar BECs than for s-wave ones. However, in the limit where the BEC contains a large number of atoms the problem of finding the ground-state density profile and low-energy dynamics simplifies. In a harmonic trap with oscillator length aho=
冑
ប/共m兲, a BEC containing N atoms of mass m which have repulsive s-wave interactions characterized by scattering length a enters the Thomas-Fermi共TF兲 regime for large values of the parameter Na/aho关6,7兴. In the TF regime the zero-point kinetic energycan be ignored in comparison to the interaction and trapping energies and the Gross-Pitaevskii equation reduces to the equations of superfluid hydrodynamics at T = 0 关6–8兴. When
applied to a trapped s-wave BEC these equations are known to admit a large class of exact analytical solutions 关9兴. The
TF approximation can also be applied to dipolar BECs关10兴.
Although the resulting superfluid hydrodynamic equations for a dipolar BEC contain the nonlocal dipolar potential, ex-act solutions can still be found关11,12兴 and we make an
ex-tensive use of them here. The calculations in this paper are all made within the TF regime.
Condensates are quantum fluids described by a macro-scopic wave function =
冑
exp关iS兴, whereis the conden-sate density and S is the condenconden-sate phase. This constrains the velocity fieldvគ =共ប/m兲ⵜគS to be curl free ⵜគ⫻vគ=0គ. In an experiment the rotation of the condensate can be accom-plished by applying a rotating elliptical deformation to the trapping potential 关13,14兴. At low rotation frequencies theelliptical deformation excites low-lying collective modes 共quadrupole, etc.兲 with a quantized angular momentum which may be viewed as surface waves 共and which obey ⵜគ ⫻vគ=0គ兲. Above a certain critical rotation frequency vortices are seen to enter the condensate and these satisfy theⵜគ⫻vគ = 0គ condition by having a quantized circulation. The hydro-dynamic equations for a BEC provide a simple and accurate description of the low-lying collective modes. Furthermore, they predict that these modes become unstable for certain ranges of rotation frequency 关15,16兴. Comparison with
ex-periments关13,14兴 and full numerical simulations of the GPE
关17–19兴 have clearly shown that the instabilities are the first
step in the entry of vortices into the condensate and the for-mation of a vortex lattice. Crucially, the hydrodynamic equa-tions give a clear explanation of why vortex lattice formation in s-wave BECs was only observed to occur at a much greater rotation frequency than that at which they become energetically favorable. It is only at these higher frequencies that the vortex-free condensate becomes dynamically unstable.
Individual vortices关20–23兴 and vortex lattices 关24–26兴 in
dipolar condensates have already been studied theoretically. However, a key question that remains is how to make such states in the first place. In this paper we extend the TF ap-proximation for rotating trapped condensates to include di-polar interactions, building on our previous work 关27,28兴.
Specifically, starting from the hydrodynamic equations of motion, we obtain the stationary solutions for a condensate in a rotating elliptical trap and find when they become dy-namically unstable to perturbations. This enables us to pre-dict the regimes of stable and unstable motions of a rotating dipolar condensate. For a nondipolar BEC 共in the TF limit兲 the transition between stable and unstable motions is
pendent of the interaction strength and depends only on the rotation frequency and trap ellipticity in the plane perpen-dicular to the rotation vector 关15,16兴. We show that for a
dipolar BEC it is additionally dependent on the strength of the dipolar interactions and also the axial trapping strength. All of these quantities are experimentally tunable and this extends the routes that can be employed to induce instability. Meanwhile, the critical rotation frequency at which vortices become energetically favorable ⍀v is also sensitive to the trap geometry and dipolar interactions 关21兴 and means that
the formation of a vortex lattice following the instability can-not be assumed. Using a simple prediction for this frequency, we indicate the regimes in which we expect vortex lattice formation to occur. By considering all of the key and experi-mentally tunable quantities in the system, we outline several accessible routes to generate instability and vortex lattices in dipolar condensates.
This paper is structured as follows. In Sec.IIwe introduce the mean-field theory and the TF approximation for dipolar BECs, in Sec.IIIwe derive the hydrodynamic equations for a trapped dipolar BEC in the rotating frame, and in Sec. IV
we obtain the corresponding stationary states and discuss their behavior. In Sec.Vwe show how to obtain the dynami-cal stability of these states to perturbations and in Sec.VIwe employ the results of the previous sections to discuss pos-sible pathways to induce instability in the motion of the BEC and discuss the possibility that such instability leads to the formation of a vortex lattice. Finally in Sec.VIIwe conclude our findings and suggest directions for future work.
II. MEAN-FIELD THEORY OF A DIPOLAR BEC
We consider a BEC with long-range dipolar atomic inter-actions, with the dipoles aligned in the z direction by an external field. The condensate wave function共mean-field or-der parameter兲 for the condensate ⬅共rគ,t兲 satisfies the GPE which is given by关4,29,30兴
iប t =
冋
− ប2 2mⵜ 2+ V共rគ,t兲 + ⌽ dd共rគ,t兲 + g兩兩2册
, 共1兲 where m is the atomic mass. Theⵜ2term arises from kinetic energy and V共rគ,t兲 is the external confining potential. BECs typically feature s-wave atomic interactions which gives rise to a local cubic nonlinearity with coefficient g = 4ប2a/m, where a is the s-wave scattering length. Note that a, and therefore g, can be experimentally tuned between positive 共repulsive interactions兲 and negative 共attractive interactions兲 values by means of a Feshbach resonance关2,3兴. The dipolarinteractions lead to a nonlocal mean-field potential ⌽dd共rគ,t兲 which is given by关5兴
⌽dd共rគ,t兲 =
冕
d3r Udd共rគ − rគ⬘
兲共rគ⬘
,t兲, 共2兲 where共rគ,t兲=兩共rគ,t兲兩2 is the condensate density andUdd共rគ兲 = Cdd
4
1 − 3 cos2
兩rគ兩3 共3兲
is the interaction potential of two dipoles separated by a vec-tor rគ, where is the angle between rគ and the polarization
direction, which we take to be the z axis. The dipolar BECs made to date have featured permanent magnetic dipoles. Then, assuming the dipoles to have moment dm and be aligned in an external magnetic field Bគ =kˆគB, the dipolar cou-pling is Cdd=0dm
2 关
29兴, where0is the permeability of free space. Alternatively, for dipoles induced by a static electric field Eគ =kˆគE, the coupling constant Cdd= E2␣2/⑀0 关30,31兴, where␣is the static polarizability and⑀0 is the permittivity of free space. In both cases, the sign and the magnitude of Cdd can be tuned through the application of a fast-rotating external field关32兴.
We will specify the interaction strengths through the parameter
dd= Cdd
3g , 共4兲
which is the ratio of the dipolar interactions to the s-wave interactions 关32兴. We take the s-wave interactions to be
re-pulsive, g⬎0, and so where we discuss negative values of dd, this corresponds to Cdd⬍0. We will also limit our analy-sis to the regime of −0.5⬍dd⬍1, where the Thomas-Fermi approach predicts that nonrotating stationary solutions are robustly stable关11兴. Outside of this regime the situation
comes more complicated since the nonrotating system be-comes prone to collapse关33兴.
We are concerned with a BEC confined by an elliptical harmonic trapping potential of the form
V共rគ兲 =12m⬜2关共1 −⑀兲x2+共1 +⑀兲y2+␥2z2兴. 共5兲 In the x-y plane the trap has mean trap frequency ⬜ and ellipticity ⑀. The trap strength in the axial direction, and in-deed the geometry of the trap itself, is specified by the trap ratio␥=z/⬜. When␥Ⰷ1 the BEC shape will typically be oblate共flattened兲 while for␥Ⰶ1 it will typically be prolate 共elongated兲, although for strong enough dipolar interactions the electrostrictive or magnetostrictive effect can cause a BEC in an oblate trap to become prolate itself.
The time-dependent GPE 共1兲 can be reduced to its
time-independent form by making the substitution 共rគ,t兲 =
冑
共rគ兲 exp共it/ប兲, whereis the chemical potential of the system. We employ the TF approximation whereby the ki-netic energy of static solutions is taken to be negligible in comparison to the potential and interaction energies. The va-lidity of this approximation in dipolar BECs has been dis-cussed elsewhere 关10兴. Then, the time-independent GPEre-duces to
V共rគ兲 + ⌽dd共rគ兲 + g共rគ兲 =. 共6兲 For ease of calculation the dipolar potential ⌽dd共rគ兲 can be expressed as ⌽dd共rគ兲 = − 3gdd
冉
2 z2共rគ兲 + 1 3共rគ兲冊
, 共7兲 where共rគ兲 is a fictitious “electrostatic” potential defined by 关11,12兴共rគ兲 = 1 4
冕
d3r
⬘
共rគ⬘
兲兩rគ − rគ
⬘
兩 . 共8兲 This effectively reduces the problem of calculating the dipo-lar potential共2兲 to the calculation of an electrostatic potentialof form 共8兲, for which a much larger theoretical body of
literature exists. Exact solutions of Eq.共6兲 for共rគ兲,共rគ兲, and
hence⌽dd共rគ兲 can be obtained for any general parabolic trap, as proven in Appendix A of Ref. 关12兴. In particular, the
solutions of共rគ兲 take the form 共rគ兲 =0
冉
1 − x2 Rx2 − y 2 R2y − z 2 Rz2冊
for 共rគ兲 ⱖ 0, 共9兲 where 0= 15N/共8RxRyRz兲 is the central density. Remark-ably, this is the general inverted parabola density profile fa-miliar from the TF limit of nondipolar BECs. An important distinction, however, is that for the dipolar BEC the aspect ratio of the parabolic solution differs from the trap aspect ratio.III. HYDRODYNAMIC EQUATIONS IN THE ROTATING FRAME
Having introduced the TF model of a dipolar BEC we now extend this to include rotation and derive hydrodynamic equations for the rotating system. We consider the rotation to act about the z axis, described by the rotation vector ⍀គ, where⍀=兩⍀គ兩 is the rotation frequency and the Hamiltonian in the rotating frame is given by
Heff= H0−⍀គ · Lˆគ, 共10兲 where H0 is the Hamiltonian in the absence of the rotation and Lˆគ =−iប共rគ⫻ⵜគ兲 is the quantum-mechanical angular mo-mentum operator. Using this result with the Hamiltonian H0 from Eq. 共1兲, we obtain 关34,35兴
iប⌿共rគ,t兲 t =
冋
− ប2 2mⵜ 2+ V共rគ兲 + ⌽ dd共rគ,t兲 + g兩⌿共rគ,t兲兩2 −⍀ប i冉
x y− y x冊
册
⌿共rគ,t兲. 共11兲 Note that all space coordinates rគ are those of the rotating frame and the time-independent trapping potential V共rគ兲, given by Eq.共5兲, is stationary in this frame. Momentumco-ordinates, however, are expressed in the laboratory frame 关34–36兴.
We can express the condensate mean field in terms of a density 共rគ,t兲 and phase S共rគ,t兲 as 共rគ,t兲=
冑
共rគ,t兲 exp关iS共rគ,t兲兴, so that the condensate velocity is vគ=共ប/m兲ⵜគS. Substituting into the time-dependent GPE共11兲 and equatingimaginary and real terms leads to the following equations of motion: t = −ⵜគ · 关共vគ − ⍀គ ⫻ rគ兲兴, 共12兲 mvគ t = −ⵜគ
冉
1 2mvគ · vគ + V共rគ兲 + ⌽dd共rគ兲 + g− mvគ · 关⍀គ ⫻ rគ兴冊
. 共13兲 In the absence of dipolar interactions共⌽dd= 0兲 Eqs. 共12兲 and 共13兲 are commonly known as the superfluid hydrodynamicequations关6–8兴 since they resemble the equation of
continu-ity and the Euler equation of motion from dissipationless fluid dynamics. Here, we have extended them to include dipolar interactions.
Note that the form of condensate velocity leads to the relation
ⵜគ ⫻ vគ =mបⵜគ ⫻ ⵜគS = 0គ, 共14兲 which immediately reveals that the condensate is irrotational. The exceptional case is when the velocity potential共ប/m兲S is singular, which arises when a quantized vortex occurs in the system.
IV. STATIONARY SOLUTION OF THE HYDRODYNAMIC EQUATIONS
We now search for stationary solutions of the hydrody-namic equations共12兲 and 共13兲. These states satisfy the
equi-librium conditions
t = 0,
vគ
t = 0. 共15兲
Following the approach of Recati et al.关15兴 we assume the
velocity field ansatz
vគ =␣ⵜគ共xy兲. 共16兲
Here,␣is a velocity field amplitude that will provide us with a key parameter to parametrize our rotating solutions. Note that this is the velocity field in the laboratory frame ex-pressed in the coordinates of the rotating frame, and also that it satisfies the irrotationality condition共14兲. Actually, the
ve-locity field amplitude ␣ can be given even more physical meaning by noting that, according to the continuity Eq.共12兲,
it can be written as关7兴
␣= −D⍀ 共17兲
whereD is the deformation of the BEC in the x-y plane D =具y2− x2典
具y2+ x2典=
y2−x2 y2+x2
. 共18兲
In the first term on the right-hand side of this expression 具¯典 signifies the expectation value in the stationary state, and in the second term on the right-hand side we have intro-duced the condensate aspect ratios defined asx= Rx/Rzand y= Ry/Rz.
Combining Eqs.共13兲 and 共16兲 we obtain the relation =m 2共˜x 2x2+˜ y 2y2+ z 2z2兲 + g共rគ兲 + ⌽ dd共rគ兲, 共19兲
˜x2=⬜2共1 −⑀兲 +␣2− 2␣⍀, 共20兲 ˜y2=⬜2共1 +⑀兲 +␣2+ 2␣⍀. 共21兲 The dipolar potential inside an inverted parabola density pro-file共9兲 has been found in Refs. 关12,27兴 to be
⌽dd 3gdd =0xy 2
冋
001− x2101+ y2011+ 3z2002 Rz 2册
− 3, 共22兲 where the coefficientsijkare given byijk=
冕
0 ⬁ ds 共x 2+ s兲i+1/2共 y 2+ s兲j+1/2共1 + s兲k+1/2, 共23兲 where i, j, and k are integers. Note that for the cylindrically symmetric case, where x=y=, the integralsijkevaluate to关37兴 ijk= 2 2F1冉
k + 1 2,1;i + j + k + 3 2;1 − 2冊
共1 + 2i + 2j + 2k兲2共i+j兲 , 共24兲 where 2F1denotes the Gauss hypergeometric function 关38兴. Thus, we can rearrange Eq.共19兲 to obtain an expression forthe density profile,
= − m 2共˜x 2 x2+˜y2y2+z 2 z2兲 g共1 − dd兲 + 3gdd n0xy 2Rz 2 关x2101+ y2011+ 3z2002− Rz2001兴 g共1 − dd兲 . 共25兲 Comparing the x2, y2, and z2 terms in Eqs.共9兲 and 共25兲 we
find three self-consistency relations that define the size and the shape of the condensate:
x 2 =
冉
z ˜x冊
21 +dd冉
3 2x 3 y101− 1冊
, 共26兲 y2=冉
z ˜y冊
21 +dd冉
3 2y 3 x011− 1冊
, 共27兲 Rz 2 =2g0 mz 2, 共28兲where = 1 −dd关1−共9xy/2兲002兴. Furthermore, by insert-ing Eq. 共25兲 into Eq. 共12兲 we find that stationary solutions
satisfy the condition
0 =共␣+⍀兲
冉
˜x2− 3 2dd ⬜2xy␥2 101冊
+共␣−⍀兲冉
˜y2−3 2dd ⬜2xy␥2 011冊
. 共29兲We can now solve Eq. 共29兲 to give the velocity field
ampli-tude ␣ for a given dd,⍀, and trap geometry. In the limit dd= 0 this amplitude is independent of the s-wave interac-tion strength g and the trap ratio␥. However, in the presence of dipolar interactions the velocity field amplitude becomes dependent on both g and␥. For fixedddand trap geometry, Eq. 共29兲 leads to branches of ␣ as a function of rotation frequency ⍀. These branches are significantly different be-tween traps that are circular共⑀= 0兲 or elliptical 共⑀⬎0兲 in the x-y plane, and so we will consider each case in turn. Note that we restrict our analysis to the range ⍀⬍⬜: for ⍀ ⬃⬜the static solutions can disappear, with the condensate becoming unstable to a center-of-mass instability关15兴.
A. Circular trapping in the x-y plane:⑀=0
We first consider the case of a trap with no ellipticity in the x-y plane共⑀= 0兲. In Fig.1共a兲we plot the solutions of Eq. 共29兲 as a function of rotation frequency ⍀ for a spherically
symmetric trap ␥= 1 and for various values of dd. Before discussing the specific cases, let us first point out that for eachddthe solutions have the same qualitative structure. Up to some critical rotation frequency, only one solution exists corresponding to␣= 0. At this critical point the solution bi-furcates, giving two additional solutions for ␣⬎0 and ␣
0.5 0.6 1.0 0.2 0.6 1.0 -0.2 -0.6 -1.0 (a) 0.7 0.8 0.9 x x y y Ω/ω α/ω εdd 1 γ Ω/ ω b (b) 102 104 10-2 10-4 0.55 0.6 0.65 0.7 0.75 0.8 εdd
FIG. 1. 共a兲 Irrotational velocity field amplitude␣ of the static condensate solutions as a function of the trap rotation frequency⍀ in a spherically symmetric trap共␥=1 and ⑀=0兲. Various values of ddare presented:dd= −0.49, 0, 0.5, and 0.99. Insets illustrate the
geometry of the condensate in the x-y plane. 共b兲 The bifurcation frequency⍀b关the point at which the solutions of␣ in 共a兲 bifurcate兴 according to Eq.共32兲 versus the trap ratio␥. Plotted are the results
fordd= −0.49, −0.4, −0.2, 0, 0.2, 0.4, 0.6, 0.8, 0.9, and 0.99. In共a兲 and共b兲 ddincreases in the direction of the arrow.
⬍0 on top of the original␣= 0 solution. We term this critical frequency the bifurcation frequency⍀b.
For dd= 0 we regain the results of Refs. 关15,16兴 with a bifurcation point at ⍀b=⬜/
冑
2 and, for ⍀⬎⍀b, nonzero solutions given by ␣=⫾冑
2⍀2−⬜2/⬜ 关15兴. The physicalsignificance of the bifurcation frequency has been estab-lished for the nondipolar case and is related to the fact that the system becomes energetically unstable to the spontane-ous excitation of quadrupole modes for ⍀ⱖ⬜/
冑
2. In the TF limit, a general surface excitation with angular momen-tum បl=បqlR, where R is the TF radius and qlis the quan-tized wave number, obeys the classical dispersion relation l2
=共ql/m兲ⵜRV involving the local harmonic potential V = m⬜2R2/2 evaluated at R 共see p. 183 of 关7兴兲. Consequently, for the nonrotating and nondipolar BEC, l=
冑
l⬜. Mean-while, inclusion of the rotational term in the Hamiltonian 共10兲 shifts the mode frequency by −l⍀. Then, in the rotatingframe, the frequency of the l = 2 quadrupole surface excita-tion becomes 2共⍀兲=
冑
2⬜− 2⍀ 关7兴. The bifurcation fre-quency thus coincides with the vanishing of the energy of the quadrupolar mode in the rotating frame, and the two addi-tional solutions arise from excitation of the quadrupole mode for ⍀ⱖ⬜/冑
2.For the nondipolar BEC it is noteworthy that⍀bdoes not depend on the interactions. This feature arises because the mode frequencieslthemselves are independent of g. How-ever, in the case of long-range dipolar interactions the poten-tial⌽ddof Eq.共7兲 gives nonlocal contributions, breaking the simple dependence of the force −ⵜគV upon R 关11兴. Thus, we
expect the resonant condition for exciting the quadrupolar mode, i.e., ⍀b=l/l 共with l=2兲, to change with dd. In Fig. 1共a兲we see that this is the case: as dipole interactions are introduced, our solutions change and the bifurcation point⍀b moves to lower 共higher兲 frequencies for dd⬎0 共dd⬍0兲. Note that the parabolic solution still satisfies the hydrody-namic equations providing −0.5⬍dd⬍1. Outside this range the parabolic solution may still exist but it is no longer guar-anteed to be stable against perturbations.
Density profiles for ␣= 0 have zero ellipticity in the x-y plane. By contrast, the 兩␣兩⬎0 solutions have an elliptical density profile, even though the trap itself has zero ellipticity. This remarkable feature arises due to a spontaneous breaking of the axial rotational symmetry at the bifurcation point. For ␣⬎0 the condensate is elongated in x while for ␣⬍0 it is elongated in y, as can be seen from Eq.共17兲 and as illustrated
in the insets in Fig. 1共a兲. In the absence of dipolar interac-tions the兩␣兩⬎0 solutions can be interpreted solely in terms of the effective trapping frequencies˜xand˜ygiven by Eqs. 共20兲 and 共21兲. The introduction of dipolar interactions
con-siderably complicates this picture since they also modify the shape of the solutions. Notably, fordd⬎0 the dipolar inter-actions make the BEC more prolate, i.e., reduce x andy, while for dd⬍0 they make the BEC more oblate, i.e., in-crease xandy.
In Fig. 1共a兲 we see that as the dipole interactions are increased the bifurcation point ⍀b moves to lower frequen-cies. The bifurcation point can be calculated analytically as follows. First, we note that for␣= 0 the condensate is cylin-drically symmetric and x=y=. In this case the aspect
ratio is determined by the transcendental equation 关11,12,30兴
冋
冉
␥2 2 + 1冊
f共兲 1 −2− 1册
+ 共dd− 1兲共2−␥2兲 32dd = 0, 共30兲 where f共兲 =2 + 2关4 − 3 000兴 2共1 −2兲 共31兲with 000=共1/
冑
1 −2兲ln关共1+冑
1 −2兲/共1−冑
1 −2兲兴 for the prolate case共⬍1兲, and000=共2/冑
2− 1兲arctan关冑
2− 1兴 for the oblate case 共⬎1兲. For small ␣→0+, we can calculate the first-order corrections toxandywith respect tofrom Eqs. 共26兲 and 共27兲. We can then insert these values in Eq.共29兲 and solve for ⍀, noting that in the limit␣→0 we have
⍀→⍀b. Thus, we find ⍀b ⬜=
冑
1 2+ 3 4 2 dd␥2 2 201−101 1 −dd冉
1 − 9 2 2 002冊
. 共32兲In Fig. 1共b兲 we plot ⍀b 关Eq. 共32兲兴 as a function of ␥ for various values ofdd. Fordd= 0 we find that the bifurcation point remains unaltered at ⍀b=x/
冑
2 as ␥=z/x is changed关15,16兴. As ddis increased the value of␥for which ⍀bis a minimum changes from a trap shape which is oblate 共␥⬎1兲 to prolate 共␥⬍1兲. Note that for dd= 0.99 the mini-mum bifurcation frequency occurs at ⍀b⬇0.55, which is over a 20% deviation from the nondipolar value. For more extreme values of dd we can expect ⍀b to deviate even further, although the validity of the inverted parabola TF solution does not necessary hold. For a fixed␥ we also find that as dd increases the bifurcation frequency decreases monotonically.B. Elliptical trapping in the x-y plane:⑀⬎0
Consider now the effect of finite ellipticity in the x-y plane 共⑀⬎0兲. Rotating elliptical traps have been created ex-perimentally with laser and magnetic fields关13,14兴.
Follow-ing the experiment of Madison et al.关13兴, we will employ a
weak trap ellipticity of⑀= 0.025. In Fig.2共a兲we have plotted the solutions to Eq.共29兲 for various values of ddin a ␥= 1 trap. As predicted for nondipolar interactions 关15,16兴, the
solutions become heavily modified for⑀⬎0. There exists an upper branch of␣⬎0 solutions which exists over the whole range of ⍀ and a lower branch of ␣⬍0 solutions which back-bends and is double valued. We term the frequency at which the lower branch back-bends to be the back-bending frequency ⍀b. The bifurcation frequency in nonelliptical traps can be regarded as the limiting case of the back-bending frequency, with the differing nomenclature em-ployed to emphasize the different structures of the solutions at this point. However, for convenience we will employ the same parameter for both, ⍀b. No ␣= 0 solution exists 共for any nonzero ⍀兲. In the absence of dipolar interactions the effect of increasing the trap ellipticity is to increase the back-bending frequency ⍀b. Turning on the dipolar interactions,
as in the case of⑀= 0, reduces⍀bfordd⬎0, and increases ⍀b fordd⬍0. This is more clearly seen in Fig. 2共b兲 where ⍀bis plotted versusddfor various values of the trap ratio␥. Importantly, the back-bending of the lower branch can introduce an instability. Consider the BEC to be on the lower branch at some fixed rotation frequency ⍀. Now consider a decreasing dd. The back-bending frequency ⍀b increases and at some point can exceed ⍀. In other words, the static solution of the BEC suddenly disappears and the BEC finds itself in an unstable state. We will see in Sec. VI that this type of instability can also be induced by variations in␥and ⑀.
As in the⑀= 0 case, increasingdddecreases bothxand y, i.e., the BEC becomes more prolate. As explained in the Introduction, this distortion is expected because of the aniso-tropy of dipolar interactions. However, because the dipolar interactions are isotropic in the x-y plane, it is perhaps sur-prising to find that they increase the deformation of the BEC in that plane. This can be clearly seen in Fig.2共a兲where we see that the magnitude of␣increases asddis increased共for any fixed value of ⍀兲. See Eq. 共17兲 for the relationship
be-tween␣and the deformation of the BEC in the x-y plane.
V. DYNAMICAL STABILITY OF STATIONARY SOLUTIONS
Although the solutions derived above are static solutions in the rotating frame, they are not necessarily stable, and so in this section we analyze their dynamical stability. Consider small perturbations in the BEC density and phase of the form =0+␦ and S = S0+␦S. Then, by linearizing the hydrody-namic equations共12兲 and 共13兲, the dynamics of such
pertur-bations can be described as
t
冋
␦S ␦册
= −冋
vគc·ⵜគ g共1 + ddK兲/m ⵜគ ·0ⵜគ 关共ⵜគ · vគ兲 + vគc·ⵜគ兴册冋
␦S ␦册
, 共33兲 wherevគc=vគ −⍀គ⫻rគ and the integral operator K is defined as共K␦兲共rគ兲 = − 32 z2
冕
␦共rគ
⬘
兲drគ⬘
4兩rគ − rគ
⬘
兩−␦共rគ兲. 共34兲 The integral in the above expression is carried out over the domain where0⬎0, that is, the general ellipsoidal domain with radii Rx, Ry, Rz of the unperturbed condensate. Ex-tending the integration domain to the region where 0+␦ ⬎0 adds higher-order effects since it is exactly in this do-main that0=O共␦兲. To investigate the stability of the BEC, we look for eigenfunctions and eigenvalues of operator共33兲:dynamical instability arises when one or more eigenvalues possess a positive real part. The size of the real eigenvalues dictates the rate at which the instability grows. Note that the imaginary eigenvalues of Eq.共33兲 relate to stable collective
modes of the system 关39兴, e.g., sloshing and breathing, and
have been analyzed elsewhere for dipolar BECs关40兴. In
or-der to find such eigenfunctions we follow Refs.关16,27兴 and
consider a polynomial ansatz for the perturbations in the co-ordinates x, y, and z of total degree N. All operators in Eq. 共33兲, acting on polynomials of degree N, result in
polynomi-als of 共at most兲 the same degree, including the operator K. This latter fact was known to 19th century astrophysicists who calculated the gravitational potential of a heterogeneous ellipsoid with a polynomial density关41,42兴. The integral
ap-pearing in Eq.共34兲 is exactly equivalent to such a potential.
A more recent paper by Levin and Muratov summarizes these results and presents a more manageable expression for the resulting potential 关43兴. Hence, using these results the
operator K can be evaluated for a general polynomial density perturbation␦= xpyqzr, with p, q, and r being non-negative integers and p + q + rⱕN. Therefore, the perturbation evolu-tion operator共33兲 can be rewritten as a scalar matrix
opera-tor, acting on vectors of polynomial coefficients, for which finding eigenvectors and eigenvalues is a trivial computa-tional task.
Using the above approach we determine the real positive eigenvalues of Eq. 共33兲 and thereby predict the regions of
dynamical instability of the static solutions. We focus on the case of an elliptical trap since this is the experimentally rel-evant case. Recall the general form of the branch diagram for this case, i.e., Fig. 2共a兲. In the ␣⬍0 half plane, the static solutions nearest the ␣= 0 axis never become dynamically unstable, except for a small region⍀⯝⬜, due to a center-of-mass instability of the condensate关44兴. The other
lower-branch solutions are always dynamically unstable and there-fore expected to be irrelevant to experiment. Thus, we only consider dynamical instability for the upper-branch solu-tions, i.e., the branch in the upper half plane共where␣⬎0兲. In Fig. 3 we plot the maximum positive real eigenvalues of the upper-branch solutions as a function of ⍀ for a fixed ellipticity ⑀= 0.025. The maximum polynomial perturbation was set at N = 3 since for this ellipticity it was found that higher-order perturbations did not alter the region of
insta-x y x y 0.2 0.6 1.0 -0.2 -0.6 -1.00 0.2 0.4 0.6 0.8 1.0
(a)
Ω/ω
α/ω
Ω/
ω
ε
dd 0.6 0.65 0.7 0.75 0.8 0.85 -0.5 0 0.5 1.0ε
dd b(b)
FIG. 2. 共a兲 Irrotational velocity field amplitude␣ as a function of the trap rotation frequency⍀ for a trap ratio␥=1 and ellipticity ⑀=0.025. Various values of dd are presented, dd= −0.49, 0, 0.5,
and 0.99, withdd increasing in the direction of the arrow. Insets
illustrate the geometry of the condensate in the x-y plane.共b兲 Back-bending point⍀bversusddfor⑀=0.025 and ␥=0.5 共solid curve兲,
bility, and such modes are therefore not displayed 共but see the end of this section for further comment on the higher-order perturbations兲.
For a givenddand␥there exists a dynamically unstable region in the⑀-⍀ plane. An illustrative example is shown in Fig. 3 共inset兲 for dd= 0 and ␥= 1. The instability region 共shaded兲 consists of a series of crescents 关16兴. Each crescent
corresponds to a single value of the polynomial degree N, where higher values of N add extra crescents from above. At the high-frequency end these crescents merge to form a main region of instability, characterized by large eigenvalues. At the low-frequency end the crescents become vanishingly thin and are characterized by very small eigenvalues which are at least one order of magnitude smaller than in the main insta-bility region 关19兴. As such these regions will only induce
instability in the condensate if they are traversed very slowly. This was confirmed by numerical simulations in Ref. 关19兴
where it was shown that the narrow instability regions have negligible effect when ramping ⍀ at rates greater than d⍀/dt=2⫻10−4
⬜
2. It is unlikely that an experiment could be sufficiently long lived for these narrow instability regions to play a role. For this reason we will subsequently ignore the narrow regions of instability and define our instability region to be the main region, as bounded by the dashed line in Fig. 3 共inset兲. For the experimentally relevant trap
ellip-ticities ⑀ⱗ0.1 the unstable region is defined solely by the N = 3 perturbations.
We define the lower bound of the instability region to be ⍀i共this corresponds to the dashed line in the inset兲. This is the key parameter to characterize the dynamical instability. As dd is increased ⍀i decreases and, accordingly, the un-stable range of ⍀ widens. Note that the upper bound of the instability region is defined by the end point of the upper branch at⍀⯝⬜.
Finally, we consider the higher perturbations N⬎3. We find that as we increase the size of matrix 共33兲 to N
= 3 , 4 , 5 , . . . the higher-lying modes that are thereby de-scribed also develop real eigenvalues as ⍀ is increased, but these lie within the region of instability already shown in Fig.3for N = 3 and so do not alter the region of instability, as
mentioned above. Significantly, we find that the modes be-come unstable in order, i.e., the perturbations contained in N = 4 that are not present in N = 3 become unstable at higher values of ⍀ than those in N=3, and similarly for N=5 in comparison to N = 4. We take this as circumstantial evidence that there is no “roton” minimum in the energy spectrum for the parameters we have considered. The possibility of a roton minimum in the Bogoliubov energy spectrum of a dipolar BEC has been widely discussed in recent literature共see, e.g., 关45–47兴兲. In analogy to the celebrated dispersion relation of
liquid 4He, the roton minimum refers to a minimum in the energy spectrum at a finite value of the eigenvalue 共e.g., momentum p兲 labeling the excitation. This means that, coun-terintuitvely, some higher-lying modes can have lower en-ergy than lower-lying modes and this causes important ef-fects in flowing systems. Pitaevskii 关48兴 discussed the case
of superfluid4He flowing through a pipe at velocityv. In the laboratory frame the energy spectrum is Galilean shifted such that E→E−pv and this leads, for large enough v, to the roton mode pr being brought down to zero energy first. Crudely speaking, this is expected to trigger an instability to the formation of a density wave with wavelength ⬀pr−1. In the present case we have rotational flow and the Galilean shifted energy E→E−L⍀ can presumably result in angular roton modes关46兴 becoming unstable as ⍀ is increased.
How-ever, our empirical observation that the modes become un-stable in order as⍀ is increased seems to rule out the pres-ence of an angular roton minimum at finite angular momentum Lrfor our parameters, at least up to N = 5. This is not surprising because roton minima in dipolar BECs have so far only been predicted to occur outside of the range −0.5 ⬍dd⬍1 where the system is stable against dipolar collapse. A proper treatment outside this stable range requires going beyond the TF approximation since the zero-point energy must be included. The interplay between rotational instabili-ties and dipolar collapse instabiliinstabili-ties remains a fascinating topic for future exploration, although very relevant theoreti-cal work has been performed关46兴.
VI. ROUTES TO INSTABILITY AND VORTEX LATTICE FORMATION
A. Procedures to induce instability
For a nondipolar BEC the static solutions and their stabil-ity in the rotating frame depend only on rotation frequency⍀ and trap ellipticity ⑀. Adiabatic changes in⑀ and⍀ can be employed to evolve the condensate through the static solu-tions and reach a point of instability. Indeed, this has been realized both experimentally 关13,14兴 and numerically
关17,18兴, with an excellent agreement with the hydrodynamic
predictions. For the case of a dipolar BEC we have shown in Secs. IVandVthat the static solutions and their instability depend additionally on the trap ratio ␥ and the interaction parameter dd. Since all of these parameters can be experi-mentally tuned in time, one can realistically consider each parameter as a distinct route to traverse the parameter space of solutions and induce instability in the system.
Examples of these routes are presented in Fig.4. Specifi-cally, Fig. 4 shows the static solutions ␣ of Eq. 共29兲 as a FIG. 3. The maximum positive real eigenvalues of Eq. 共33兲
共solid curves兲 for the upper-branch solutions of␣ as a function of ⍀. We assume⑀=0.025, ␥=1, and N=3 and present various dipolar strengths dd= −0.49, 0, 0.5, and 0.99, with dd increasing in the
direction of the arrow. The inset shows the full region of dynamical instability in the ⑀-⍀ plane for dd= 0. The narrow regions have
negligible effect and so we only consider the main instability region 共bounded by the dashed line兲.
function of⍀ 关Fig.4共a兲兴,⑀关Fig.4共b兲兴, dd关Fig.4共c兲兴, and␥ 关Fig.4共d兲兴. In each case the remaining three parameters are fixed at ⑀= 0.025, ␥= 1,⍀=0.7⬜, anddd= 0.99. Dynami-cally unstable solutions are indicated with red circles. Gray arrows mark routes toward instability共the point of onset of instability being marked by an asterisk兲, where the free pa-rameter⍀, ⑀,dd, or␥ is varied adiabatically until either a dynamical instability is reached or the solution branch back-bends and so ceases to exist. For solutions with ␣⬎0, the instability is always due to the system becoming dynamically unstable共dashed arrows兲, whereas for␣⬍0 the instability is always due to the solution branch back-bending on itself 共solid arrows兲 and so ceasing to exist. Numerical studies 关18兴
indicate that these two types of instability involve different dynamics and possibly have distinct experimental signatures. Below we describe the adiabatic variation of each param-eter in more general detail, beginning with the established routes toward instability in which共i兲 ⍀ and 共ii兲⑀are varied, and then routes, which are unique to dipolar BECs, based on adiabatic changes in 共iii兲 dd and 共iv兲␥. In each case it is crucial to consider the behavior of the points of instability, namely, the back-bending point⍀band the onset of dynami-cal instability of the upper branch⍀i.
共i兲 Adiabatic introduction of ⍀. The relevant parameter space of⑀and⍀ is presented in Fig.5共a兲, with the instability frequencies ⍀b共⑀兲 共solid curves兲 and ⍀i共⑀兲 共dashed curves兲 indicated. For a BEC initially confined to a nonrotating trap with finite ellipticity ⑀, as the rotation frequency ⍀ is in-creased adiabatically the BEC follows the upper-branch so-lution关Fig.4共a兲共dashed arrow兲兴. This particular route traces out a horizontal path in Fig.5共a兲until it reaches⍀i共⑀兲, where
the stationary solution becomes dynamically unstable. For the specific parameters of Fig. 4共a兲the system becomes un-stable at ⍀=⍀i共⑀兲⬇0.65⬜. More generally, Fig. 5共a兲 shows that asddis increased,⍀i共⑀兲 is decreased and as such instabilities in the stationary solutions will occur at lower rotation frequencies. At ⑀⯝0.1, the curve for ⍀idisplays a sharp kink, arising from the shape of the dynamically un-stable region, as shown in Fig.3 共inset兲.
共ii兲 Adiabatic introduction of ⑀. Here, we begin with a cylindrically symmetric 共⑀= 0兲 trap, rotating at a fixed fre-quency ⍀. The trap ellipticity ⑀ is then increased adiabati-cally, and in the phase diagram of Fig. 5共a兲the BEC traces out a vertical path starting at ⑀= 0. The ensuing dynamics depend on the trap rotation speed relative to⍀b共⑀= 0兲.
共a兲 For ⍀⬍⍀b共⑀= 0兲 the condensate follows the upper branch of the static solutions shown in Fig.4共a兲. This branch moves progressively to larger␣. For ⍀⬍⍀i共⑀兲 the BEC re-mains stable but as ⑀ is increased further the condensate eventually becomes dynamically unstable. Figure5共a兲shows that as dd is increased ⍀i共⑀兲 is decreased and as such the dynamical instability of the stationary solutions occurs at a lower trap ellipticity.
共b兲 For ⍀⬎⍀b共⑀= 0兲 the condensate accesses the lower-branch solutions nearest the ␣= 0 axis. These solutions are always dynamically stable and the criteria for instability are
0.4 α Ω/ω 0.6 0.8 1.0 0 0.5 1.0 -0.5 -1.0 (a) α0 0.5 -0.5 0 0.4 (c) 0.2 0.6 0.8 1.0 εdd α0 0.5 -0.5 0 0.05 0.1 0.15 ε (b) α0 0.5 -0.5 0.1 1.0 10 γ (d)
FIG. 4.共Color online兲 Stationary states in the rotating trap char-acterized by the velocity field amplitude ␣, determined from Eq. 共29兲. Dynamically unstable solutions are marked with red circles. In
each of the figures the共a兲 trap rotation frequency ⍀, 共b兲 trap ellip-ticity⑀, 共c兲 dipolar interaction strength dd, and共d兲 axial trapping strength␥ are varied adiabatically, while the remaining parameters remain fixed at⍀=0.7⬜,⑀=0.025, dd= 0.99, and␥=1. The adia-batic pathways to instability 共onset marked by red asterisk兲 are schematically shown by the dashed and solid arrows. Dashed ar-rows indicate a route toward dynamical instability, whereas solid arrows indicate an instability due to disappearance of the stationary state. ε 0 0.05 0.1 0.5 0.6 0.7 Ω/ω (i) (ii) (iii) (a) 0.8 0.9 Ω/ω 0.8 0.7 0.6 -0.5 0 0.5 1.0 0.1 1.0 10 εdd (i) (ii) (iii) (iv) (i) (iii) (ii) (b) (c) γ
FIG. 5. 共Color online兲 共a兲 Phase diagram of⑀ versus ⍀b共solid
curves兲 and ⍀i共dashed curves兲 for␥=1 and 共i兲 dd= −0.49,共ii兲 0.5,
and共iii兲 0.99. 共b兲 Phase diagram of ddversus⍀b共solid curves兲 and ⍀i共dashed curves兲 for⑀=0.025 and 共i兲 ␥=0.5, 共ii兲 1, and 共iii兲 2. 共c兲
Phase diagram of␥ versus ⍀b共solid curves兲 and ⍀i共dashed curves兲 for ⑀=0.025 and 共i兲 dd= −0.49,共ii兲 0, 共iii兲 0.5, and 共iv兲 0.99. In
each case the solid共dashed兲 arrows depict the routes to instability shown in Fig.4.
instead determined by whether the solution exists. As ⑀ is increased the back-bending frequency ⍀b共⑀兲 increases. Therefore, when ⑀ exceeds some critical value, the lower-branch solutions disappear for the chosen value of rotation frequency ⍀. This occurs when ⍀⬍⍀b共⑀兲. Figure 5共a兲 shows that asddis increased⍀b共⑀兲 is decreased and as such instabilities in the system will occur at a higher trap elliptic-ity. At this point the parabolic condensate density profile no longer represents a stable solution. The particular route indi-cated in Fig. 4共b兲 is included in Fig.5共a兲as a vertical solid gray arrow.
共iii兲 Adiabatic change in dd. The relevant parameter space ofddand⍀ is shown in Fig.5共b兲for several different trap ratios. Consider that we begin from an initial BEC in a trap with finite ellipticity⑀= 0.025 and rotation frequency⍀. 共This can be achieved, for example, by increasing ⑀ from zero at fixed⍀.兲 Then, by changing ddadiabatically an in-stability can be induced in two ways:
共a兲 For ⍀⬍⍀b共dd兲 the condensate follows the upper-branch solutions until they become unstable. This route to instability is shown in Fig.4共c兲by the dashed arrow, with the corresponding path in Fig.5共b兲shown by the vertical dashed arrow. Thus, for ⍀⬍⍀i共dd兲 the motion remains stable. However, for ⍀⬎⍀i共dd兲 the upper branch becomes dy-namically unstable. In Fig. 5共b兲,⍀i共dd兲 共dashed curves兲 is plotted for different trap ratios. As can be seen, the stable region of the upper branch becomes smaller as dd is increased.
共b兲 For ⍀⬎⍀b共dd兲 the condensate follows the lower-branch solutions nearest the ␣= 0 axis. These solutions are always stable and hence an instability can only be induced when this solution no longer exists, i.e.,⍀⬍⍀b共dd兲. Figure 5共b兲 shows ⍀b共dd兲 共solid curves兲 for various trap aspect ratios. As can be seen the back-bending frequency ⍀b de-creases asddis increased. Thus, ifddis increased the sys-tem will remain stable. However, ifddis decreased then the system will become unstable when⍀=⍀b共dd兲.
共iv兲 Adiabatic change in␥. Figure5共c兲shows the param-eter space of␥and⍀. Consider, again, an initial stable con-densate with finite trap rotation frequency ⍀ and ellipticity ⑀= 0.025. Then through adiabatic changes in ␥ the conden-sate can traverse the parameter space and, depending on the initial conditions, the instability can arise in two ways:
共a兲 For ⍀⬍⍀b共␥兲 the condensate exists on the upper branch. It is then relevant to consider the onset of dynamical instability ⍀i共␥兲 关dashed curves in Fig. 5共c兲兴. Providing ⍀ ⬍⍀i共␥兲 the solution remains dynamically stable. However, once⍀⬎⍀i共␥兲 the upper-branch solutions become unstable. 共b兲 For ⍀⬎⍀b共␥兲 the condensate exists on the lower branch nearest the␣= 0 axis. These solutions are always dy-namically stable and instability can only occur when the mo-tion of the back-bending point causes the solumo-tion to disap-pear. This occurs when⍀⬍⍀b共␥兲, with ⍀b共␥兲 shown in Fig. 5共c兲by solid curves for various dipolar interaction strengths. These two paths to instability are shown in Fig.4共d兲and are also indicated in Fig.5共c兲as vertical gray arrows, where the dashed共solid兲 arrow corresponds to the␣⬎0 共␣⬍0兲 path.
B. Is the final state of the system a vortex lattice?
Having revealed the points at which a rotating dipolar condensate becomes unstable, we will now address the
ques-tion of whether this instability leads to a vortex lattice. First, let us review the situation for a nondipolar BEC. The pres-ence of vortices in the system becomes energetically favor-able when the rotation frequency exceeds a critical frequency ⍀v. Working in the TF limit, with the background density taking the parabolic form 共9兲, ⍀v can be approximated as 关49兴 ⍀v= 5 2 ប mR2ln 0.67R s . 共35兲
Here, the condensate is assumed to be circularly symmetric with radius R and s=ប/
冑
2m0g is the healing length that characterizes the size of the vortex core. For typical conden-sate parameters ⍀v⬃0.4⬜. It is observed experimentally, however, that vortex lattice formation occurs at considerably higher frequencies, typically ⍀⬃0.7⬜. This difference arises because above ⍀v, the vortex-free solutions remain remarkably stable. It is only once a hydrodynamic instability occurs 共which occurs in the locality of ⍀⬇0.7⬜兲 that the condensate has a mechanism to deviate from the vortex-free solution and relax into a vortex lattice. Another way of visu-alizing this is as follows. Above⍀v the vortex-free conden-sate resides in some local energy minimum, while the global minimum represents a vortex or vortex lattice state. Since the vortex is a topological defect, there typically exists a consid-erable energy barrier for a vortex to enter the system. How-ever, the hydrodynamic instabilities offer a route to navigate the BEC out of the vortex-free local energy minimum toward the vortex lattice state.Note that vortex lattice formation occurs via nontrivial dynamics. The initial hydrodynamic instability in the vortex-free state that we have discussed in this paper is only the first step 关18兴. For example, if the condensate is on the upper
branch of hydrodynamic solutions 共e.g., under adiabatic in-troduction of ⍀兲 and undergoes a dynamical instability, this leads to the exponential growth of surface ripples in the con-densate 关13,18兴. Alternatively, if the condensate is on the
lower branch and the static solutions disappear共e.g., follow-ing the introduction of ⑀兲, the condensate undergoes large and dramatic shape oscillations. In both cases the destabili-zation of the vortex-free condensate leads to the nucleation of vortices into the system. A transient turbulent state of vortices and density perturbations then forms, which subse-quently relaxes into a vortex lattice configuration关18,50兴.
In the presence of dipolar interactions, however, the criti-cal frequency for a vortex depends crucially on the trap ge-ometry ␥ and the strength of the dipolar interactions dd. Following Ref. 关21兴 we will make a simple and
approxi-mated extension of Eq. 共35兲 to a dipolar BEC. We will
con-sider a circularly symmetric dipolar condensate with radius R = Rx= Ry that satisfies Eqs. 共26兲–共28兲 and insert this into Eq.共35兲 for the condensate radius. This method still assumes
that the size of the vortex is characterized by the s-wave healing length s. Although one does expect the dipolar in-teractions to modify the size of the vortex core, it should be noted that Eq. 共35兲 only has a logarithmic accuracy and is
relatively insensitive to the choice of vortex core length scale. The dominant effect of the dipolar interactions in Eq.
共35兲 comes from the radial size and is accounted for. Note
also that this expression is for a circularly symmetric system while we are largely concerned with elliptical traps. How-ever we will employ a very weak ellipticity ⑀= 0.025 for which we expect the correction to the critical frequency to be correspondingly small.
As an example, we take the parameter space of rotation frequency ⍀ and dipolar interactions dd. We first consider the behavior in a quite oblate trap with␥= 10. In Fig.6共a兲we plot the instability frequencies⍀iand⍀bfor this system as a function of the dipolar interactions dd. Depending on the specifics of how this parameter space is traversed, either by adiabatic changes in ⍀ 共vertical path兲 or in dd 共horizontal path兲, the condensate will become unstable when it reaches one of the instability lines 共short and long dashed lines兲. These points of instability decrease weakly with dipolar in-teractions and have the approximate value ⍀i⬇⍀b ⬇0.75⬜. On the same plot we present the critical rotation frequency⍀vaccording to Eq.共35兲. In order to calculate this
we have assumed a BEC of 150 000 52Cr atoms confined within a trap with ⬜= 2⫻200 Hz. In this oblate system we see that the dipolar interactions lead to a decrease in⍀v, as noted in关21兴. This dependence is very weak at this value
of ␥, and throughout the range ofddpresented it maintains the approximate value⍀v⬇0.1⬜. Importantly these results show that when the condensate becomes unstable a vortex or vortex lattice state is energetically favored. As such, we ex-pect that in an oblate dipolar BEC a vortex lattice will ulti-mately form when these instabilities are reached.
In Fig.6共b兲, we make a similar plot but for a prolate trap with ␥= 0.1. The instability frequencies show a somewhat similar behavior to the oblate case. However, ⍀v is drasti-cally different, increasing significantly withdd. We find that
this qualitative behavior occurs consistently in prolate sys-tems, as noted in关21兴. This introduces two regimes
depend-ing on the dipolar interactions. Forddⱗ0.8, ⍀i,b⬎⍀v, and so we expect a vortex or vortex lattice state to form follow-ing the instability. However, forddⲏ0.8 we find an intrigu-ing regime in which ⍀i,b⬍⍀v. In other words, while the instability in the vortex-free parabolic density profile still occurs, a vortex state is not energetically favorable. The final state of the system is therefore not clear. Given that a prolate dipolar BEC is dominated by attractive interactions 共since the dipoles lie predominantly in an attractive end-to-end con-figuration兲 one might expect similar behavior to the case of conventional BECs with attractive interactions共g⬍0兲 where the formation of a vortex lattice can also be energetically unfavorable. Suggestions for final state of the condensate in this case include center-of-mass motion and collective oscil-lations, such as quadrupole modes or higher angular-momentum-carrying shape excitations关51–53兴. However, the
nature of the true final state in this case is beyond the scope of this work and warrants further investigation.
VII. CONCLUSIONS
By calculating the static hydrodynamic solutions of a ro-tating dipolar BEC and studying their stability, we have pre-dicted the regimes of stable and unstable motions. In general we find that the back-bending or the bifurcation frequency ⍀b decreases with increasing dipolar interactions. In addi-tion, the onset of dynamical instability in the upper-branch solutions,⍀i, decreases with increasing dipolar interactions. Furthermore these frequencies depend on the aspect ratio of the trap.
By utilizing the features of dipolar condensates we detail several routes to traverse the parameter space of static solu-tions and reach a point of instability. This can be achieved through adiabatic changes in trap rotation frequency ⍀, trap ellipticity ⑀, dipolar interactionsdd, and trap aspect ratio␥, all of which are experimentally tunable quantities. While the former two methods have been employed for nondipolar BECs, the latter two methods are unique to dipolar BECs. In an experiment the latter instabilities would therefore demon-strate the special role played by dipolar interactions. Further-more, unlike for conventional BECs with repulsive interac-tions, the formation of a vortex lattice following a hydrodynamic instability is not always favored and depends sensitively on the shape of the system. For a prolate BEC with strong dipolar interactions, there exists a regime in which the rotating spheroidal parabolic Thomas-Fermi den-sity profile is unstable and yet it is energetically unfavorable to form a lattice. Other outcomes may then develop, such as a center-of-mass motion of the system or collective modes with angular momentum. However, for oblate dipolar con-densates, as well as prolate condensates with weak dipolar interactions, the presence of vortices is energetically favored at the point of instability and we expect the instability to lead to the formation of a vortex lattice.
We acknowledge financial support from the Australian Research Council 共A.M.M.兲, the Government of Canada 共N.G.P.兲, and the Natural Sciences and Engineering Research Council of Canada共D.H.J.O.D.兲.
Ω/ω
ε
dd(a)
(b)
0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1.0Ω/ω
FIG. 6. 共Color online兲 The relation between the instability fre-quencies,⍀b共long dashed red curve兲 and ⍀i共short dashed curve兲,
and the critical rotation frequency for vorticity⍀v共solid curve兲 for 共a兲 an oblate trap␥=10 and 共b兲 a prolate trap ␥=0.1. The instability frequencies are based on a trap with ellipticity⑀=0.025 while ⍀vis obtained from Eq.共35兲 under the assumption of a 52Cr BEC with 150 000 atoms and scattering length as= 5.1 nm in a circularly symmetric trap with⬜= 2⫻200 Hz.
关1兴 A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T. Pfau, Phys. Rev. Lett. 94, 160401共2005兲.
关2兴 T. Lahaye, T. Koch, B. Frohlich, M. Fattori, J. Metz, A. Gries-maier, S. Giovanazzi, and T. Pfau, Nature共London兲 448, 672 共2007兲.
关3兴 T. Koch, T. Lahaye, J. Metz, B. Frohlich, A. Griesmaier, and T. Pfau, Nat. Phys. 4, 218共2008兲.
关4兴 L. Santos, G. V. Shlyapnikov, P. Zoller, and M. Lewenstein, Phys. Rev. Lett. 85, 1791共2000兲.
关5兴 S. Yi and L. You, Phys. Rev. A 63, 053607 共2001兲.
关6兴 C. J. Pethick and H. Smith, Bose-Einstein Condensation in Dilute Gases 共Cambridge University Press, Cambridge, En-gland, 2002兲.
关7兴 L. Pitaevskii and S. Stringari, Bose-Einstein Condensation 共Oxford University Press, New York, 2003兲.
关8兴 P. Nozieres and D. Pines, Theory Of Quantum Liquids 共West-view, New York, 1999兲.
关9兴 S. Stringari, Phys. Rev. Lett. 77, 2360 共1996兲.
关10兴 N. G. Parker and D. H. J. O’Dell, Phys. Rev. A 78, 041601共R兲 共2008兲.
关11兴 D. H. J. O’Dell, S. Giovanazzi, and C. Eberlein, Phys. Rev. Lett. 92, 250401共2004兲.
关12兴 C. Eberlein, S. Giovanazzi, and D. H. J. O’Dell, Phys. Rev. A
71, 033618共2005兲.
关13兴 K. W. Madison, F. Chevy, W. Wohlleben, and J. Dalibard, Phys. Rev. Lett. 84, 806共2000兲; K. W. Madison, F. Chevy, V. Bretin and J. Dalibard, ibid. 86, 4443共2001兲.
关14兴 E. Hodby, G. Hechenblaikner, S. A. Hopkins, O. M. Marago, and C. J. Foot, Phys. Rev. Lett. 88, 010405共2001兲.
关15兴 A. Recati, F. Zambelli, and S. Stringari, Phys. Rev. Lett. 86, 377共2001兲.
关16兴 S. Sinha and Y. Castin, Phys. Rev. Lett. 87, 190402 共2001兲. 关17兴 E. Lundh, J.-P. Martikainen, and K.-A. Suominen, Phys. Rev.
A 67, 063604共2003兲.
关18兴 N. G. Parker, R. M. W. van Bijnen, and A. M. Martin, Phys. Rev. A 73, 061603共R兲 共2006兲.
关19兴 I. Corro, N. G. Parker, and A. M. Martin, J. Phys. B 40, 3615 共2007兲.
关20兴 S. Yi and H. Pu, Phys. Rev. A 73, 061602共R兲 共2006兲. 关21兴 D. H. J. O’Dell and C. Eberlein, Phys. Rev. A 75, 013604
共2007兲.
关22兴 R. M. Wilson, S. Ronen, J. L. Bohn, and H. Pu, Phys. Rev. Lett. 100, 245302共2008兲.
关23兴 M. Abad, M. Guilleumas, R. Mayol, M. Pi, and D. M. Jezek, Phys. Rev. A 79, 063622共2009兲.
关24兴 N. R. Cooper, E. H. Rezayi, and S. H. Simon, Phys. Rev. Lett.
95, 200402共2005兲.
关25兴 J. Zhang and H. Zhai, Phys. Rev. Lett. 95, 200403 共2005兲. 关26兴 S. Komineas and N. R. Cooper, Phys. Rev. A 75, 023623
共2007兲.
关27兴 R. M. W. van Bijnen, D. H. J. O’Dell, N. G. Parker, and A. M. Martin, Phys. Rev. Lett. 98, 150401共2007兲.
关28兴 A. M. Martin, N. G. Parker, R. M. W. van Bijnen, A. Dow, and D. H. J. O’Dell, Laser Phys. 18, 322共2008兲.
关29兴 K. Góral, K. Rzążewski, and T. Pfau, Phys. Rev. A 61,
051601共R兲 共2000兲.
关30兴 S. Yi and L. You, Phys. Rev. A 61, 041604共R兲 共2000兲. 关31兴 M. Marinescu and L. You, Phys. Rev. Lett. 81, 4596 共1998兲. 关32兴 S. Giovanazzi, A. Gorlitz, and T. Pfau, Phys. Rev. Lett. 89,
130401共2002兲.
关33兴 N. G. Parker, C. Ticknor, A. M. Martin, and D. H. J. O’Dell, Phys. Rev. A 79, 013617共2009兲.
关34兴 A. J. Leggett, Quantum Liquids: Bose-Condensation and Coo-per Pairing in Condensed Matter Systems共Oxford University Press, New York, 2006兲.
关35兴 A. J. Leggett, in Bose-Einstein Condensation: From Atomic Physics to Quantum Fluids, Proceedings of the Thirteenth Physics Summer School, edited by C. M. Savage and M. P. Das共World Scientific, Singapore, 2001兲, pp. 1–42.
关36兴 L. D. Landau and E. M. Lifshitz, Course of Theoretical Phys-ics: Mechanics, 3rd ed. 共Butterworth-Heinemann, Oxford, 1982兲.
关37兴 Table of Integrals, Series, and Products, 6th ed., edited by L. S. Gradshteyn and I. M. Ryzhik共Academic Press, San Diego, 2000兲.
关38兴 Handbook of Mathematical Functions, edited by M. Abramowitz and I. Stegun共Dover, New York, 1974兲. 关39兴 Y. Castin, in Coherent Matter Waves, Lecture Notes of Les
Houches Summer School, edited by R. Kaiser, C. Westbrook, and F. David共Springer-Verlag, Berlin, 2001兲, pp. 1–136. 关40兴 R. M. W. van Bijnen, N. G. Parker, A. M. Martin, and D. H. J.
O’Dell共unpublished兲.
关41兴 N. M. Ferrers, Quart. J. Pure and Appl. Math 14, 1 共1877兲. 关42兴 F. W. Dyson, Quart. J. Pure and Appl. Math 25, 259 共1891兲. 关43兴 M. L. Levin and R. Z. Muratov, Astrophys. J. 166, 441 共1971兲. 关44兴 P. Rosenbusch, D. S. Petrov, S. Sinha, F. Chevy, V. Bretin, Y. Castin, G. Shlyapnikov, and J. Dalibard, Phys. Rev. Lett. 88, 250403共2002兲.
关45兴 D. H. J. O’Dell, S. Giovanazzi, and G. Kurizki, Phys. Rev. Lett. 90, 110402共2003兲; L. Santos, G. V. Shlyapnikov, and M. Lewenstein, ibid. 90, 250403 共2003兲; S. Giovanazzi and D. O’Dell, Eur. Phys. J. D 31, 439共2004兲; U. R. Fischer, Phys. Rev. A 73, 031602共R兲 共2006兲.
关46兴 S. Ronen, D. C. E. Bortolotti, and J. L. Bohn, Phys. Rev. Lett.
98, 030406共2007兲; R. Wilson, S. Ronen, and J. Bohn, Phys.
Rev. A 80, 023614共2009兲. .
关47兴 O. Dutta, R. Kanamoto, and P. Meystre, Phys. Rev. Lett. 99, 110404 共2007兲; O. Dutta and P. Meystre, Phys. Rev. A 75, 053604共2007兲.
关48兴 L. P. Pitaevskii, Zh. Eksp. Teor. Fiz. 39, 423 共1984兲. 关49兴 E. Lundh, C. J. Pethick, and H. Smith, Phys. Rev. A 55, 2126
共1997兲.
关50兴 N. G. Parker and C. S. Adams, Phys. Rev. Lett. 95, 145301 共2005兲; J. Phys. B 39, 43 共2006兲.
关51兴 N. K. Wilkin, J. M. F. Gunn, and R. A. Smith, Phys. Rev. Lett.
80, 2265共1998兲.
关52兴 B. Mottelson, Phys. Rev. Lett. 83, 2695 共1999兲.
关53兴 C. J. Pethick and L. P. Pitaevskii, Phys. Rev. A 62, 033609 共2000兲.
