Table Of Content

Structure of vortices in rotating Bose-Einstein condensates Gentaro Watanabe,a,b S. Andrew Gifford,c Gordon Baym,a,c and C. J. Pethicka aNORDITA, Blegdamsvej 17, DK-2100 Copenhagen Ø, Denmark bThe Institute of Chemical and Physical Research (RIKEN), 2-1 Hirosawa, Wako, Saitama 351-0198, Japan cDepartment of Physics, University of Illinois, 1110 W. Green Street, Urbana, IL 61801 WecalculatethestructureofindividualvorticesinrotatingBose-Einsteincondensatesinatrans- 7 verse harmonic trap. Making a Wigner-Seitz approximation for the unit cell of the vortex lattice, 0 0 we derive the Gross-Pitaevskii equation for the condensate wave function in each cell of the lat- 2 tice, includingeffectsof varyingcoarse grained density. Wecalculate theAbrikosov parameter, the fractional core area, and theenergy of individual cells. n a PACSnumbers: 03.75.Hh,03.75.Kk,05.30.Jp,67.40.Vs J 2 1 I. INTRODUCTION cal density (an effect taken into account only globally in Ref. [8]). Numerical calculations of the local structure ] r of vortices are given in Sec. III. In Sec. IV we calcu- e Rapid rotationof an atomic Bose-Einsteincondensate late the Abrikosov parameter, fractional core area, and h hasmadeitpossibletocreatevortexlatticesinwhichthe the energyofa Wigner-Seitz cell, anddiscuss the spatial ot vortex core size is comparable to the intervortexspacing dependence of the vortex structure. . [1, 2, 3]. This has opened up a new regime inaccessi- t a ble in superfluid He-II. At low rotation rates, the vortex m coresize is smallcomparedwith the intervortexspacing, II. BASIC FORMALISM - and it is given by the Gross-Pitaevskii healing length d [4,5, 6]. However,asthe angularvelocityΩ isincreased, n In this section, we derive the GP equation for a sin- the vortex core radii grow and eventually become com- o gle Wigner-Seitz cell of the vortex lattice. We assume a parable to the separation of vortices [7, 8, 9, 10]. In c separable trapping potential, [ a harmonic trap, when Ω is close to the transverse trap- 3 pingfrequencyω⊥,thecondensatewavefunctionisdom- V(r)=V⊥(r⊥)+Vz(z), (1) inated by the lowest Landau level (LLL) of the Coriolis v 1 force [11]. The criterion to be in the LLL regime is that where r⊥ is the coordinate perpendicular to the rota- 5 ¯hΩ gn, where gn is the interaction energy, n is the tion axis, which we take to be in the z direction. We 1 part≫icle density, g 4π¯h2as/m is the two-body interac- also assume that the confinement in the z direction ≡ 6 tionstrength,mistheparticlemass,andas isthes-wave is sufficiently tight that the system is essentially two- 0 scattering length. If the system at rest can be described dimensional,andthatthe condensatewavefunctionΨis 6 bytheGross-Pitaevskii(GP)equationforthecondensate separable, 0 wavefunction, then the GP approachis validfor any ro- / at tation rate until the number of vortices Nv becomes of Ψ=ψ(r⊥)H(z). (2) m order 10% of the number of particles N [12, 13, 14, 15]. We normalize Ψ to the number of particles by taking d- strTuhcteupreuropfothseeovfortthiicsespianpaerlaisttticoecfaolrcaurlabtiterathryerionttaertnioanl dr⊥|ψ|2 = N and dz|H|2 = 1. In the remainder of n this paper we consider only a two-dimensional system velocities and position in the lattice. Such calculations R R o and drop the subscript. c are needed to analyze current experiments in which gn The structu⊥re of the condensate wave function in : is comparable to h¯Ω, and one is between the limits of v the transverse direction is determined by the two- slow rotation and the LLL regime. This work general- i dimensional GP equation, X izes an earlier account [8], which assumed that the vor- r tex structure was the same throughout the lattice. It ¯h2 a is also complementary to the work of Cozzini, Stringari, 2+V +g ψ 2 ψ =µψ. (3) 2D −2m∇ | | and Tozzo [16], who investigated numerically the vortex (cid:18) (cid:19) structure in the case when the rotation rate equals the Here µ is the chemicalpotential, excluding the contribu- transversetrappingfrequency,andthustheparticleden- tions from the motion in the z direction. The effective sity is uniform. Our basic approach, as in Ref. [8], is two-dimensional interaction strength is to treat a unit cell of the vortex lattice in the Wigner- Seitz approximation, in which one replaces the actual g =g dz H 4. (4) unit cell, hexagonal for a triangular lattice, by a circu- 2D | | lar one. In Sec. II, we present the basic formalism. We Z find that the structure depends not only on the rotation For a system uniform in the z direction, g = g/Z, 2D rate and local density, but also on derivatives of the lo- whereZ istheextentofthecloudinthezdirection,while 2 if H is the ground-state wave function for an oscillator potentialwithfrequencyω ,g =g/√2πd ,whered = z 2D z z (h¯/mω )1/2 [17]. z We motivate our derivation of the basic equation for ρ thevortexstructurebyfocusingonthesimpleexampleof avortexatthecenterofthe trap. TosolvetheGPequa- l tion, subject to the boundary conditions imposed by the r WS cell presence of the other vortices, we make a Wigner-Seitz approximation to the central cell of the vortex lattice, R replacing the hexagonal unit cell by a circle of the same j area. The radius of the Wigner-Seitz cell is 1 ℓ= , (5) 0 √πnv where nv is the local vortex density. The wave function FIG.1: AWigner-Seitzcellofthevortexlattice. Theposition in the central cell is then a p-wave state, of the center of the cell j relative to the center of the trap is Rj,and ρ is the transverse coordinate measured from Rj. ψ =√n f(r)eiφ, (6) 0 where φ is the azimuthal angle measured relative to the center of the trap, and n is the average of the particle 0 axis, and that the system is in equilibrium. We then density over the centralcell. We take the averageof f 2 over the cell to be unity, d2r f2 = πℓ2. Without |los|s specialize, in calculating the vortex structure, to a uni- j form triangular lattice. In separate publications we will of generality, we may take f to be real; then f satisfies R investigate the equilibrium distortions of the lattice and calculate the elastic constants [18]. ¯h2 1 ∂ ∂f f r + V(r)f The local fluid velocity is given by h¯ Φ/m, which we − 2m r∂r ∂r − r2 ∇ (cid:20) (cid:18) (cid:19) (cid:21) write as a sum ofthe backgroundflow, v , plus a locally + g n f3 =µf. (7) b 2D 0 circulating flow: In order that the wave function in the central cell con- nect smoothly with that in neighboring cells, we require ¯h ¯h Φ=v + φ , (9) b j ∂f/∂r = 0 at the boundary of the Wigner-Seitz cell, m∇ m∇ r = ℓ. Note that the structure of the vortex depends on the local trapping potential. In the Wigner-Seitz ap- where φ is the azimuthal angle measured with re- j proximation, the structure of the vortex on the axis of spect to the center of the vortex in cell j. This thetrapdependsontheangularvelocityofthetraponly latter term accounts for the vorticity in cell j. We throughthe length ℓ that determines the size of the unit can write the background flow in general as v = b cell. If the symmetry ofthe lattice is takeninto account, (h¯/m)[ Φ(R )+ φ (r)],whereR isthecenterofcell j sj j there will be other terms that depend on Ω due to the jandφ∇isregular∇inthecell[19]. Hereweneglecteffects sj flow created by the rotating lattice of vortices outside due to φ . sj the central cell, but these will be small because the true We estimate v by calculating the circulation around hexagonal unit cell is well approximated by a circle. b a circle of radius R in terms of the number of vor- We now turn to the generalvortex lattice. As in Refs. j tices N within the circle: v = Ω zˆ R , where [7,8],weassumethatthenumberofvorticesislargeand Ω =h¯Nvj /mR2. Thus b v × j writethecondensatewavefunctionastheproductofthe v vj j square root of a slowly varying coarse-graineddensity n, a rapidly varying real function f, which describes the ¯h ¯h local structure of vortices, and a phase factor eiΦ: Φ=Ωv Rj + φj . (10) m∇ × m∇ ψ = n(r)f(r)eiΦ(r) . (8) Ifthelatticeisuniformwithvortexdensityn ,thenΩ = ˆ ˆv v As in the case of thepcentral cell, we normalize f2 by ¯hveπcntov/rmin. Itnheadadziitmiounth∇alφjdi=recφtjio/nρ,awbhoeurteRφj,isatnhdeρun=it setting its average over each unit cell of the lattice to j r R is the coordinate in the x-y plane measured from unity. j − the center of the cell (see Fig. 1). We assume, in setting up the general expressions for the energy and angular momentum of the vortex lattice, We now derive the Gross-Pitaevskii equation in a thatthelatticeisclosetouniform,withpossibledisplace- Wigner-Seitzcellforarotatingcondensateinaharmonic ments depending only on the distance from the central trap. Thetotalenergyinthelaboratoryframeis[Eq.(8) 3 of Ref. [8]]: Thus 1 E I(ω2+Ω2) ≃ 2 v E = d2r ¯h2 ψ 2+ mω2r2 nf2+ g2D n2f4 + n(R ) d2r ¯h2 ∂f 2+ f2 2m|∇ | 2 2 j 2m ∂ρ ρ2 Z (cid:26) (cid:27) j Zj ( "(cid:18) (cid:19) # 1 X Iω2 ≃ 2 1mΩ2ρ2f2+ g2Dn(R )f4 ¯h2 ¯h2 −2 v 2 j ) + d2r n ( f)2+ ( φ )2 Xj Zj (2m ∇ (cid:18)2m ∇ j −12 Rj∂∂R n(Rj) d2r (mΩ2vρ2−¯hΩv)f2 . +12mΩ2vR2j +h¯Ωv·(Rj ×∇φj)(cid:19)f2+ g22Dnf4) Xj j Zj (13) ¯h2 1 Here we have retained in the kinetic and trap energies + d2r ( √n)2f2+ f2 n , (11) termsoforderh¯Ωandh¯ω2/Ω perparticle,andthusthe 2m ∇ 2∇ ·∇ v Z (cid:26) (cid:27) last two terms in Eq. (11), which are smaller by a factor 1/N ,areneglected[24]. ToobtainEq.(13),weusedthe v fact that where I = d2r mn(r)f2r2 is the moment of inertia. In general, n(r) varies in space, and we continue to allow d2r nf2ρ r for the possRibility that Ω depends on position. j Zj · v X 1 ∂ The angularmomentumL [Eq.(11)ofRef. [8]] is sim- n(Rj)+ Rj n(Rj) d2r f2ρ2, ≃ 2 ∂R ilarly j (cid:18) j (cid:19)Zj X (14) wherethefactor1/2comesfromaveragingintwodimen- L = d2r nf2[r ¯h Φ] sions [25]. z × ∇ The angular momentum L [Eq. (12)] is Z d2r nf2[mΩ (R2+ρ R )+h¯ρ R /ρ2] ≃ j Zj v j · j · j L = d2r nf2[r×¯h∇Φ]z X Z +h¯N. (12) IΩ +h¯N n(R ) d2r mΩ f2ρ2 v j v ≃ − j Zj X 1 ∂ R n(R ) d2r mΩ ρ2 ¯h f2. (15) j j v −2 ∂R − In the remainder of this paper we neglect the weak j j Zj X (cid:0) (cid:1) spatial dependence of Ω arising from the small devi- v We write the moment of inertia as ations of the vortex lattice from regular triangularity [18, 20, 21, 22, 23], and assume that the lattice is uni- I = I¯+ d2r mn(f2 1)r2 form. Wealsoassumethatf isaxiallysymmetricwithin − Z each cell. Our objective is to calculate the energy in a ∂ frame rotatingwith angularvelocityΩ as a functionalof I¯+m n(R )+R n(R ) j j j f, and then to determine f by minimizing the energy. ≃ ∂Rj j X In Eq. (11), the energy per particle due to the trapping 1 ∂2 potential(Iω2/2)andthekineticenergyofbulkrotation + R2 n(R ) d2r (f2 1)ρ2 o(mtfhΩte2vhReΩ22vclo∼teurdm¯h.Ω)vTaNhreev,olerfaeodspridneegcrticvmoenlωyt2r,iRbw2uht∼eioren¯hs(Rωto2i/stΩhtvehs)eNervtaedraimnudss = I¯4 ¯hjN∂R+j2m j !n(ZRjj) d2r−f2ρ2 − 2Ω are independent of f, and therefore to determine f it v j Zj X is necessary to include the leading f-dependent contri- ∂ 1 ∂2 butions, which are smaller by a factor 1/Nv. The en- +m Rj∂R n(Rj)+ 4Rj2∂R2n(Rj) ergy per particle from interparticle interactions is of or- j j j ! X der Ng /R2, and since this depends on f it is unnec- 2D d2r (f2 1)ρ2 , (16) essary to include corrections of relative order 1/N . To v × − calculate the energy in a uniform lattice, we expand the Zj coarse-grained density about the center of the Wigner- where I¯ d2r mn(r)r2 is the momentof inertia calcu- Seitz cell, n(r) n(R )+ρ n(R )+ 1(ρ )2n(R ). lated for≡the coarse-grained density. Note that here the ≃ j ·∇ j 2 ·∇ j R 4 ′ secondderivativetermgivescontributionsoforderh¯Ωto minimizing E with respect to Ω keeping n and f fixed, v ′ the energy per particle. In the limit of a small core, i.e., ∂E /∂Ω =0. This condition has the form v for h¯Ω g n, I reduces to I¯. The angular momentum 2D written≪in terms of I¯is ∂E′ =I¯(Ω Ω)+ ∂Er′es =0, (23) v ∂Ω − ∂Ω 1 m v v L = I¯Ω + ¯hN + n(R )αΩ d2r ρ2f2 v 2 2 j v where the last term is given by E′ = E′ Xj Zj (I¯/2) (ω2 Ω2)+(Ω Ω )2 and is ofresorder N¯hΩ−. m − − v + j 4 n(Rj)βΩvZjd2r (f2−1)ρ2, (17) mFreonmt(cid:2)wthiitshiteafrollileorwwsotrhkast[(2Ω0,−(cid:3)2Ω1,v)2/2Ω]. ∝Si1n/cNevw,einaassgurmeee- X that N 1, we shall henceforth put Ω =Ω. Thus v v ≫ where β R ∂ ∂lnn(R ) ω˜2 =ω2+ α+ (ω2 Ω2). (24) α j n(R )= j , (18) 4 − j (cid:18) (cid:19) ≡ n(R )∂R ∂lnR j j j In the fast rotation regime, in which Ω ω, the cor- and rection to the trap frequency due to the α≃and β terms is small and ω˜ ω for all cells. Thus the vortex struc- R2 ∂2 ≃ j ture does notchange significantlyfromcell to cell. How- β n(R ). (19) ≡ n(Rj)∂Rj2 j ever,at lowerrotationvelocities, the position-dependent αandβ correctionsandthevariationofthevortexstruc- ′ From Eqs. (13) and (17), the total energy E in the turefromcelltocellarenon-negligible. Thefrequencyω˜ frame rotating with angular velocity Ω is equalsω inthe centralcellanddecreaseswith increasing R . j ′ E =E ΩL As explained in earlier works, for Na /Z 1 Ω/ω − s ≫ − 1 ¯hN [8,21,22,26](whichisequivalenttotheconditionNv = I¯ (ω2 Ω2)+(Ω Ωv)2 1) the global density profile is to a good approximatio≫n 2 − 2Ω − − (cid:18) v(cid:19) a Thomas-Fermi parabola, n(r) = n(0)(1 r2/R2) [27]. π¯h2(cid:2) ω2 ∂ (cid:3) Here R is the radius of the cloud [8, 21, 22−, 26], and the ¯hNΩ 1 R n(R ) − − 4m Ω2 − j∂R j central density is [20, 28] (cid:18) v (cid:19) j j X π¯h2 (ω2 Ω2)+(Ω Ω )2 R2 ∂2 n(R ) Nmω2 1/2 Ω2 1/2 − 16mΩ2 − − v j∂R2 j n(0)= 1 . (25) v (cid:2) (cid:3)Xj j (cid:18)πbg2D (cid:19) (cid:18) − ω2(cid:19) + E , (20) j Here b is the Abrikosov parameter, j X f4 d2rf4 where we have collected all terms dependent on f in b h i = j . (26) ≡ f2 2 d2r h i R j ¯h2 ∂f 2 f2 Ej = n(Rj) d2r 2m ∂ρ + ρ2 In deriving the Thomas-FermifRorm it has been assumed Zj ( "(cid:18) (cid:19) # that variations of b over the cloud, which amount to at the very most 16%, may be neglected. The radius of the 1 g + mω˜2ρ2f2+ 2D n(Rj)f4 . (21) cloud is given by 2 2 ) ω2 −1/2 The effective trap frequency ω˜ is given by R2 =κ′b1/2ℓ2 1 , (27) 0 Ω2 − (cid:18) (cid:19) β ω˜2 ≡ω2+α(ω2−ΩΩv)+ 4(ω2+Ω2v−2ΩΩv). (22) with 1/2 [The coefficient of the term ρ2f2 differs from that in Eq. ′ Ng2Dm κ 2 . (28) (15) of Ref. [8] because of the assumption made there 0 ≡ π¯h2 (cid:18) (cid:19) that f is independent of position.] The explicit form of E is found by setting Ω = 0 in Eq. (20). The term pro- In the experiments of Ref. [10], ωz =2π 5.3 Hz, and portionaltoI¯(ω2 Ω2)istheenergydue tothetrapping as =5.6 nm for the tripletstate of87Rb. T×he number of and centrifugal p−otentials, while the term proportional atomsvariesfromN 107 forΩ<0.95ω toN 105 for to I¯(Ω Ω )2 is the extra kinetic energy due to the fact Ω 0.99ω, and the c∼orrespondin∼g values of κ ∼ κ′/√b − v ≈ 0 ≡ 0 that when Ω = Ω there is bulk motion in the rotat- areoforder100and10,respectively. Wecalculatebelow v 6 ing frame. The equilibrium value of Ω is determined by for κ =10 and 100. v 0 5 For a Thomas-Fermi profile, α = β = 2[1 n(0)/n(r)]= 2r2/(R2 r2), and [29] − − − 5r2 ω˜2 =ω2 (ω2 Ω2). (29) − 2(R2 r2) − − Varying E with respect to f, subject to the normaliza- j tion condition d2r f2 = πℓ2, we obtain the GP equa- j tion in a Wigner-Seitz cell: R ¯h2 1 ∂ ∂f f 1 ρ + mω˜2ρ2f+g nf3 =µf . −2m ρ∂ρ ∂ρ − ρ2 2 2D (cid:20) (cid:18) (cid:19) (cid:21) (30) This equation, with boundary conditions f(0) = 0 and FIG. 2: The first five basis functions f(0) of the free Hamil- i ∂f(ℓ)/∂ρ = 0, determines the optimized form of f. It tonian H , for Ω=ω for thecentral cell, where ω˜ =ω. 0 is a generalization of Eq. (7) to vortices away from the center of the trap. The equation fails for vortices close to the edge of the cloud where the change of the coarse- where Ai is the normalization constant ensuring graineddensitywithinthecellislargeandtheexpansion d2η [(f(0))2 1] = 0, and F is Kummer’s conflu- ofn(r)aboutRj isinadequate;theexpansionisjustified enjt hyperigeom−etric function (s1ee1, e.g., Ref. [31]). The R provided eigenvalue is determined by the boundary condition 1 ∂n(r) 1 r 1 ∂ f(0)(η )=0 , (36) ℓ 1 , (31) η i cell n(r) ∂r ∼ √N R1 r2/R2 ≪ (cid:12)(cid:12) (cid:12)(cid:12) v − where η ℓ/d˜= ω˜/Ω. Figure 2 shows an example a condition(cid:12)(cid:12) that (cid:12)(cid:12)breaks down for the outermost single cell ≡ (0) of the first five f plotted as functions of ρ/ℓ. The layer of vortices, where r R. i p → ith eigenfunction f(0) has i 1 nodes. At Ω = ω and i − thus ω˜ = ω, the eigenvalue Λ(0) of the first basis state 1 III. VORTEX STRUCTURE f(0), determined by Eq. (36), equals 4, and f(0) reduces 1 1 to the p-wave function for a harmonic oscillator: f(0) Wenowinvestigatethevortexcorestructureusingthe (ρ/d)exp( ρ2/2d2). Therefore, in this limit, f(0) i1s th∼e framework derived in the previous section. It is conve- − 1 LLLcomponentofthewavefunctionintheWigner-Seitz nienttomeasurelengthsinunitsoftheeffectiveoscillator length d˜ ¯h/mω˜ [cf. Eq. (24)]; the GP equation (30) cell, and the fi(0), for i ≥ 2, describe the higher Landau ≡ level components. becomes p We expand f in Eq. (32) as 1 ∂ ∂f f − η∂η η∂η + η2 +η2f +κf3 =Λf , (32) f = cifi(0) , (37) (cid:18) (cid:19) i X where η ρ/d˜, Λ 2mµd˜2/¯h2, and κ 2mg nd˜2/¯h2. 2D and obtain the expansion coefficients c for given κ by ≡ ≡ ≡ i Thescaledstructureofavortexdependsonlyonthelocal propagating the solution numerically in imaginary time, parametersκand(viatheboundaryconditions)ℓ/d˜. For thusdampingthecontributionofhigherenergystates. In a Thomas-Fermi density profile, this analysis we use up to 15 basis functions for κ =10 0 andup to20basisfunctions forκ =100. Ingeneralthe ω2 Ω2 1/2 Rj2 number of basis functions require0d increases as the core κ=κ − 1 . (33) 0(cid:18) ω˜2 (cid:19) − R2! size decreases relative to the cell size. Before showing results, we comment on the limits of We solve Eq.(32) by constructing a basis set of eigen- the present analysis. The convergence of the expansion (0) (37)becomespoorintheslowrotationregimeatlargerκ, functions f of the free (κ = 0) Hamiltonian H = η−1∂ (η∂{ )i+}1/η2+η2 (see, e.g., Ref. [30]), satisfy0ing where higher energy states have significant amplitudes. η η − More importantly, we have implicitly assumed [e.g., in (0) (0) (0) Eq. (8)] that the cloud contains a large number of vor- H f =Λ f , (34) 0 i i i tices, which is not the case for slow rotation. From Eq. (27) one finds (0) (0) with eigenvalue Λ . The eigenfunction f of H is i i 0 R2 ω2 −1/2 fi(0) =Ai 1F1 1−Λ(i0)/4,2;η2 ηe−η2/2 , (35) Nv ≃ ℓ2 =κ0b(cid:18)Ω2 −1(cid:19) . (38) (cid:16) (cid:17) 6 FIG.3: (Coloronline) VortexcorestructureatvariousΩfor ρ thecentralcellobtainedfromtheGross-Pitaevskiiequationin c a Wigner-Seitz cell, with κ =10 (black lines) and κ =100 0 0 (bluelines). , ρ L L L which is large compared with unity for Ω/ω ≫ 1/ κ2+1. , 0 In the following, we illustrate the procedure by cal- ρ p culating the vortex structure close to the center of the trap,whereαandβ canbe neglectedandω˜ =ω. Figure 3 shows the function f obtained from the GP equation (32) in a Wigner-Seitz cell for 0.5 Ω/ω 0.99. With ≤ ≤ increasing Ω, the vortex core radius increases, becoming comparable to the Wigner-Seitz cell radius at Ω> 0.9ω. In Fig. 4 we compare f at κ = 10 obtained ∼numeri- 0 cally with analytic expressionsin the fast and slow rotation regimes. In this figure κ =10 0 (central cell) fLLL(ρ)= 1 ρe−ρ2/2ℓ2 (39) f 1 2/eℓ f − LLL fc isthe p-wavefunctionofpatwo-dimensionalharmonicos- cillatorwithd=ℓ,the lowestenergysolutionofEq.(32) without the interaction term in the limit Ω ω. ρ l → For slow rotation, the core structure is well approximated by the single vortex form [6], FIG.4: (Coloronline) Comparisonofthenumericalsolution ρ ofthelocalwavefunctionf withtheLLLstructurefLLL and f (ρ)=A , (40) c c 2ξ2+ρ2 with thesingle vortex form fc, calculated in thecentral cell. where the coherence lengthpξ is given by ¯h2 and, at Ω = 0.9ω, it is better approximated by fc than ξ2 , (41) by f . For slow rotation the single vortex form f is a ≡ 2mg n(0) LLL c 2D good approximation[Figs. 4(c) and 4(d)]. and Next we investigate the amplitude ci of the ith basis state f(0) in f, Eq. (37). As Fig. 5, a plot of c for each −1/2 i i 2ζ of the cases in Figs. 3 and 4 vs i, shows, c decreases A = 1+2ζln , (42) i c 1+2ζ with i almost exponentially, exp[ c(i 1)], where c (cid:20) (cid:18) (cid:19)(cid:21) ∼ − − is determined by κ and Ω. At Ω = 0.9ω, where the 0 with ζ ξ2/ℓ2, normalizes fc according to system is not yet in the LLL regime, the amplitudes for ≡ jd2r fc(ρ)2 =πℓ2. i 2, corresponding to higher Landau levels (LL) in the ≥ As Fig.4(a)shows,the numericalsolutionoff atΩ= rapidlyrotatinglimit,arelessthan5%(c 0.044);even R 2 ≃ 0.99ω is well described by f . However,the numerical atΩ=0.5ω,theyarelessthan10%(c 0.090),sothat LLL 2 ≃ solution starts to deviate from f with decreasing Ω the population in higher LL is less than 1%. A small LLL 7 FIG.5: Theamplitudesci ofthestatesfi(0) intheexpansion FIG. 6: Abrikosov parameter b vs Γ−1 2h¯Ω/g n(0), for of f, for the cases of κ0 = 10 shown in Figs. 3 and 4. In all κ =10 (for Ω 0.425ω) and κ =1L0L0L(f≡or Ω 02.D5ω). cases here, c1 >0.995 and ci <0.1 for i 2. 0 ≥ 0 ≥ ≥ TABLE I: Abrikosov parameter b and two measures of the fractional area of the core, r2/ℓ2 and r′2/ℓ2, at κ = 10, admixture of order 1% of the higher LL components is c c 0 calculated for thewave functions f in Figs. 3 and 4. sufficienttochangetheinternalvortexstructurefromthe LLL to the form for an isolated vortex. Ω/ω Γ−1 b r2/ℓ2 r′2/ℓ2 LLL c c 1.00 1.158 0.225 0.368 IV. EQUILIBRIUM PROPERTIES OF ∞ 0.99 2.807 1.141 0.217 0.316 VORTICES 0.90 0.826 1.115 0.199 0.231 0.75 0.454 1.087 0.177 0.171 We nowcalculatetheAbrikosovparameterb,thefrac- 0.50 0.231 1.056 0.117 0.103 tionalareaofavortexcore,andtheenergyofacellofthe vortex lattice, focusing on the central cell. In the final part of this section we discuss the spatial dependence of the vortex structure. The factor ρ2 makes this quantity sensitive to the be- havior of f(ρ) in the outer parts of the cell; in the slow rotationregime,the slowvariationof f inthe outer part A. Abrikosov parameter ofthe cell contributes significantlyto r2/ℓ2. Inaddition, c forslowrotation(forΩ<0.4ωinthepresentcalculations) f develops a local max∼imum for ρ < ℓ, which leads to a Figure 6 shows the Abrikosovparameter as a function uofesΓ−LoLf1Lb≡for2h¯tΩhe/gΩ2Dvna(l0u)e,swinhiFleigTsa.b3leaIndgiv4e.sIsnpethcieficsmvaall-l rfoarpiκd0d=ec1r0eaasnedofartc2Γ/−Lℓ2L1Lw<ith0.d0e5crfoearsκin0g=Ω1a0t0.ΓT−LL1hLe<∼me0a.2- core limit, where Γ−1 0, b = 1. Extrapolation of sure(43)offractionalarea∼isalsorelativelyinsensitiveto the curve in Fig. 6LtLoLΓ→−1 = 0, indicates a linear in- the variation of the vortex wave function f close to the LLL crease in b at small Γ−1 , in agreement with Ref. [7]. center of the cell. The curve increases mLoLnLotonically with Γ−1 towards An alternative measure of the core size which more LLL accuratelycharacterizesthe structureclosetothe vortex b = d2r f4 / d2r = (e2 5)/[4(e 2)2] 1.158 j LLL j − − ≃ is as Ω ω and thus Γ−1 . This value of b is re- R→ R LLL → ∞ −1 markably close to that for the exact LLL wave function ′ ∂f(0) r f(ℓ) . (44) for a triangular lattice, b=1.1596 [32]. c ≡ ∂ρ (cid:18) (cid:19) Inaddition,this quantityisnotinfluencedbythe behav- ior of f in the outer parts of a lattice cell. B. Fractional core area Figure 7 shows these two measures of the fractional area of the core for κ = 10 and 100 vs increasing ro- 0 Reference[8]introducedthe measureofthe fractionof tation rate. For slow rotation, the measure r′ more ac- c the area of a Wigner-Seitz cell taken up by the vortex curately describes the core size extracted from the data core: [9, 10]. In the rapidly rotating regime, where the vortex coreradiusiscomparabletotheintervortexdistance,the d2r [f(ℓ)2 f(ρ)2]ρ2 r2 j − . (43) rms radius of the density deficit rc is readily accessible c ≡ R jd2r [f(ℓ)2−f(ρ)2] to experiment. For Γ−LL1L > 0.3, the calculated values of ∼ R 8 100 E 0 a m bd 2−nh/ 10 E0,0 E0,int c of 2π s nit u n s i 1 E0,b e gi ner κ 0 = 100 E FIG. 7: (Color online) Fractional area of the core for the (central cell) centralcell,measuredbyr2/ℓ2(solidlines)andr′2/ℓ2(dashed 0.1 c c lines) (see text), as functions of Γ−1 . The lines a and c are LLL for κ = 100 and the lines b and d for κ = 10. The dotted 0 0 lines show the empirical expression for the fractional area, 0.1 1 1r2.3/4ℓΓ2−LfLo1Lr,laorbgteaiΓn−ed1 fo[8r].low Ω [9]and theasymptotic valueof ΓL-L1L=2−hΩ/g2Dn c LLL FIG. 8: (Color online) The energy of the central Wigner- Seitz cell E (solid line) in units of n¯h2/m. The kinetic plus 0 rc are consistent with the measurements of Refs. [9, 10] trapping potential energy E0,0 (dashed line) and interaction (see also Table I). As Γ−1 , r2/ℓ2 approaches energyE0,int (dash-dottedline)arealsoplotted. (Inthenon- LLL →∞ c interacting case, E =2π.) The “binding”energy E isthe 0 0,b r2 11 4e energy gain compared to the energy for the first basis state ℓc2 = 6 −2e ≃0.225 , (45) f1(0). − (cf. r′2/ℓ2 =e−1 0.368 as Γ−1 ). c ≃ LLL →∞ The theoreticalestimates of fractionalcore areaarein qualitative agreement with the published experimental data [9, 10]. However, it is important to observe that in the analysis of the experiments, core areas were determined from a Gaussian fit to the profiles of the density deficit in the core. In addition, in our calculations we have not allowed for the three-dimensionalnature of the clouds and the effects of expansion on vortex structure. For these reasons, and because of uncertainties in the number of particles in the experiments, it is premature to perform a detailed comparison between theory and experiment. C. Energy of the central cell In Fig. 8, we show for the central Wigner-Seitz cell, FIG.9: (Coloronline) Spatialdependenceofthevortexcore the energy Ej=0, the contributions from the kinetic structurefor κ0 =100 and Ω=0.9ω. and trap energies, E d2r n (h¯2/2m)[(∂ f)2 + 0,0 ≡ j=0 { ρ f2/ρ2]+mω2ρ2f2/2 ,andtheinteractionenergyE } R 0,int ≡ (g /2) d2r n2f4, calculated with the numerical so- 2D j=0 lutions for f as a function of Γ−1 . The dependence of R LLL E on Γ−1 is small, while E Γ . Thus in the sΓeln−o0e,w10rgryotda0ot.Lm5iLo4Lin,naEltimesi,t=E(Γ0E−LaL1lLso≪aronud10g,)iEh,ntlwy∼hsecbraeLelcLeotsLmhaeessinΓdtLoeLmrLai.cntaiAonntt EAΓ−t01aΓn−Ld=L1Lt1h=e(ae0nt.1eΩ,rgE=y0,cb0a.li9csωu∼laat7ne%ddκfoofrE=th0e,10ffia0rl,sltiΓnbg−a1tsois=∼sta00t..0e78%f21(60a))t.. ΓwL−LiTtLL1hLLoi→≃necxrh∞eiab,siiEtn0gth=0Γ,e0−LEL1e0Lff,.0ec(≃t0I,ni2notπtfhneh¯hin2g/ohmner-0.i,)n0LteLraccotimngpolnimenittswiftohr FFniooLfirrLcLκsamn0ta=lylle1fr0ro,ΓmE−10t,h,Eebi0ro,0tv,haElEu0e,isnt,faoanr0ndκd0EE=0,b/1E0d0oLLdnaLteoctΓrde−LaiL1ffsLeer>∼wsiitg1h-. LLL 0,b 0,b 0 slowerrotation,wealsoshowE ,thedifferencebetween decreasing κ . 0,b 0 9 D. Spatial variation of vortex structure dimensional nature of the system and the expansion of the cloud prior to the measurement of vortex structure. Although we have looked so far only at the central In the numerical calculation in this paper we have as- cell, two effects give rise to spatial dependence of the sumed a uniform triangular lattice. The next step is to vortex structure. The first is that the coarse grained extend the methods developed here to consider distor- density n depends on position; also the effective trap- tions of the vortex lattice, to calculate the elastic prop- ping frequency ω˜ depends on spatial derivatives of erties of the lattice, including the elastic constants C1 n(r). Both n and ω˜ decrease with increasing distance and C2 [13, 16, 28, 33, 34], for general rotation rates in from the center of the cloud, and this increases the thepresenceofaspatiallyvaryingcoarsegraineddensity. core radius in the outer region of the cloud. Figure 9 shows the variation of the vortex core structure with increasing distance r from the center. At r = 0.75R the Acknowledgments core radius is 30% greater than for the central cell. ∼ WethankH.M.Nilsen,L.M.Jensen,andH.Smithfor stimulating and helpful discussions. This work was supported in part by Grants-in-Aid for Scientific Research V. OUTLOOK provided by the Ministry of Education, Culture, Sports, ScienceandTechnologythroughResearchGrantNo. 14- Inthispaperwehaveconsideredtheequilibriumstruc- 7939, by the Nishina Memorial Foundation, and by the ture of vortices in two-dimensional situations. In order JSPS Postdoctoral Program for Research Abroad. This to be able to make detailed comparisonwith experimen- research was also supported in part by NSF Grants No. tal data, it is necessary to take into account the three- PHY03-55014and No. PHY05-00914. [1] M. R. Matthews, B. P. Anderson, P. C. Haljan, D. S. Phys. Rev.A 75, 013602 (2007). Hall,C.E.Wieman,andE.A.Cornell, Phys.Rev.Lett. [19] Tomakecontactbetweenthepresentdecompositionand 83, 2498 (1999). thelong-wavelength average in theelastohydrodynamics [2] K. W. Madison, F. Chevy, W. Wohlleben, and J. Dal- of Refs. [15, 28, 33], we note that the long-wavelength ibard, Phys.Rev.Lett. 84, 806 (2000). average,denotedby ... ,of Φismv(r)/¯h,wherev(r)is h i ∇ [3] J.R.Abo-Shaeer,C.Raman,J.M.Vogels, andW.Ket- thelongwavelengthflowvelocityinthelab.Thevelocity terle, Science 292, 476 (2001). vR(r) in a frame rotating at angular velocity Ω is given [4] E. P. Gross, NuovoCimento 20, 454 (1961). by vR(r)=v(r) Ω r. This is related to thefunction − × [5] L.P.Pitaevskii,Zh.Eksp.Teor.Fiz.40,646(1961)[Sov. ǫ(r) which gives the smoothed displacement of a vortex Phys.JETP 13, 451 (1961)]. fromitspositionforauniformlatticewithvortexdensity [6] A.L.Fetter,inLecturesinTheoreticalPhysics,editedby n0v =mΩ/π¯h byvR+2Ω×ǫ(r)=(h¯/m)∇φs,where φs K. Mahanthappa and W. E. Brittin (New York, 1969), describes vorticity-free superfluid flow. Thus Φ(r) = h∇ i Vol. XIB,p. 351. (m/¯h)Ω [r 2ǫ(r)]+ φs(r), from which we identify × − ∇ [7] U.R.FischerandG.Baym,Phys.Rev.Lett.90,140402 the long-wavelength limit of φsj(r) with the potential (2003). φs(r). The curl of this equation implies Φ(r) = ∇×h∇ i [8] G. Baym and C. J. Pethick, Phys. Rev. A 69, 043619 2(m/¯h)Ω(1 ǫ); since the mean density of vortices (2004). is given by n−v(∇r)·= n0v−∇·ˆn0vǫ(r)˜, we see that ∇× [9] I. Coddington, P. C. Haljan, P. Engels, V. Schweikhard, Φ(r) = 2πn (r), as expected. Note that this result v h∇ i S. Tung, and E. A. Cornell, Phys. Rev. A 70, 063607 may also be obtained by working in terms of a function (2004). giving the smoothed displacement of a vortex from the [10] V.Schweikhard,I.Coddington, P.Engels, V.P. Mogen- equilibrium lattice, which has a non-uniform density of dorff, and E. A. Cornell, Phys. Rev. Lett. 92, 040404 vortices. (2004). [20] D. E. Sheehy and L. Radzihovsky, Phys. Rev. A 70, [11] T.-L. Ho, Phys.Rev. Lett. 87, 060403 (2001). 051602(R) (2004); 70, 063620 (2004). [12] N. K. Wilkin, J. M. F. Gunn, and R. A. Smith, Phys. [21] G. Watanabe, G. Baym, and C. J. Pethick, Phys. Rev. Rev.Lett. 80, 2265 (1998). Lett. 93, 190401 (2004). [13] J. Sinova, C. B. Hanna, and A. H. MacDonald, Phys. [22] N. R. Cooper, S. Komineas, and N. Read, Phys. Rev. A Rev.Lett. 89, 030403 (2002). 70, 033604 (2004). [14] N. R. Cooper, N. K. Wilkin, and J. M. F. Gunn, Phys. [23] A.Aftalion,X.Blanc,andJ.Dalibard,Phys.Rev.A71, Rev.Lett. 87, 120405 (2001). 023611 (2005). [15] G. Baym, Phys. Rev.A 69, 043618 (2004). [24] Each of these terms is of order ¯h2/2mR2, and thus they [16] M.Cozzini, S.Stringari,andC.Tozzo,Phys.Rev.A73, scale in the same way with R as does the interaction 023615 (2006). energy g N/R2. They are thus negligible compared 2D ∼ [17] Y.Castin and R. Dum,Eur. Phys.J. D 7, 399 (1999). with the interaction energy provided Na /Z 1. For s [18] G.Baym,S.A.Gifford,C.J.Pethick,andG.Watanabe, 1 > Na /Z 1 Ω/ω, the term (¯h2/2m)(≫√n)2f2 s ∼ ≫ − ∇ 10 ′ wouldleadtoaGaussiandensityprofileifthevortexlat- const with E given by Eq. (20). To establish this result tice were forced to be perfect [8]. However, as has been withinthepresentframeworkfor1>Na /Z 1 Ω/ω, s demonstrated earlier, when distortion of the vortex lat- itisnecessarytoallowfordistortion∼ofthevor≫tex−lattice. tice is allowed for, these contributions are largely can- [28] G. Baym, Phys.Rev.Lett. 91, 110402 (2003). celedbycontributionstotheelastic energyofthevortex [29] Note that the effective oscillator potential within a cell latticeandthedensityprofilehastheThomas-Fermiform becomes inverted when n(r)/n(0) < 5(1 χ)/(7 5χ), provided N 1 [21]. Within the contex of the present where χ = Ω2/ω2; however, this inversion−, which−at the v ≫ paper, the Thomas-Fermi form for Na /Z 1 follows rotationsofinterestoccursonlyintheouterlayersofthe s from the condition δE′/δn=µ, where µ is t≫he chemical lattice, haslittle practical effect onthevortexstructure. potential in therotating frame. [30] T.Busch,B.-G.Englert,K.Rza¸z˙ewski,andM.Wilkens, [25] Bycontrast,thecalculation ofRef.[8]tookintoaccount Found.Phys. 28, 549 (1998). onlytheglobalaverageofthisterm,whichvanisheswhen [31] Handbook of Mathematical Functions with Formu- thedensityvanishesat large distances; asaconsequence las, Graphs, and Mathematical Tables, edited by M. the present equations differ from those of [8], e.g., the Abramowitz andI.A.Stegun(National Bureauof Stan- signs of the second and third terms of Eq. (16) are dif- dards, Washington, DC, 1972), Chap. 13, p.504. ferent from those of Eq. (13) of Ref. [8] for the global [32] W.H.Kleiner,L.M.Roth,andS.H.Autler,Phys.Rev. average. 133, A1226 (1964). [26] J. Sinova, C. B. Hanna, and A. H. MacDonald, Phys. [33] G. Baym and E. Chandler, J. Low Temp. Phys. 50, 57 Rev.Lett. 90, 120401 (2003). (1983). [27] In the present formalism the Thomas-Fermi form for [34] E. B. Sonin, Phys.Rev.A 72, 021606(R) (2005). Na /Z 1 follows from the condition δE′/δn(r) = s ≫