Table Of Content

Vortex structure in spinor F=2 Bose-Einstein condensates W. V. Pogosov, R. Kawate, T. Mizushima, and K. Machida Department of Physics, Okayama University, Okayama 700-8530, Japan (Dated: 6th February 2008) Extended Gross-Pitaevskii equations for the rotating F = 2 condensate in a harmonic trap are solved both numerically and variationally using trial functions for each component of the wave 6 function. Axially-symmetricvortexsolutionsareanalyzedandenergiesofpolarandcyclicstatesare 0 calculated. Theequilibriumtransitionsbetweendifferentphaseswithchangingofthemagnetization 0 are studied. We show that at high magnetization the ground state of the system is determined by 2 interaction in ”density” channel, and at low magnetization spin interactions play a dominant role. n Althoughthereare fivehyperfinestates, all theparticles are always condensed in one, twoor three a states. Two noveltypesof vortex structuresare also discussed. J 4 PACSnumbers: 73.21.La,73.30.+y,79.60.Dp,68.65.La,36.40.Qv 2 ] I. INTRODUCTION mated spin interactions are situated in the close vicinity r e of the phase boundary between polar and cyclic phases h [22];Thedetailedstudyontherotatinggroundstatewith Properties of Bose-Einstein condensates (BEC) of al- t the cyclic spin interaction is an unexplored region. o kali atom gases attract a considerable current interest. . Recently quantized vortices and lattices of vortices have The aim of the present work is to study vortex struc- t a beenobtainedexperimentallyinBECcloudsconfinedby ture in rotating spinor F = 2 condensate having finite m magnetization. The condensate wave function has five magnetic traps [1, 2, 3, 4]. BEC’s can have internal - degrees of freedom associated with the hyperfine spin. components. We solved the extended Gross-Pitaevskii d equations both numerically and variationally using trial Such condensates are usually called spinor BEC’s. First n functionsforeachcomponentofthewavefunction. There examples of these systems with hyperfine spin F = 1 o c were found in optically trapped 23Na [5]. In zero mag- is a good agreement between the results of both meth- [ netic field,spinF =1condensatecanbe intwodifferent ods. We restricted our consideration only to the case states, which are called ferromagnetic and polar [6, 7]. of axially-symmetric solutions. Energies of polar, ferro- 1 magnetic and cyclic states with various sets of winding v Dependingonthevaluesofinteractionparameters,which numbers for different components of the order parame- 9 determine coupling between different hyperfine states, ter were evaluated. The equilibrium transitions between 3 one of these states has a lowest energy. Vortex matter 5 inspinorBEC’sis representedbya richvarietyofrather differentphaseswithchangingofthemagnetizationwere 1 studied. exotictopologicalexcitations. Thesevorticeswereinves- 0 tigatedinalargenumberoftheoreticalworksforthecase 6 F =1(see,e.g.,Refs. [8,9,10,11,12,13,14,15,16,17]). 0 II. THEORETICAL FORMALISM / In the most recent experiments, F = 2 spinor Bose- t a Einstein condensates have been created and studied [18, m 19, 20, 21]. However, superfluid phases in the F = 2 Consider two-dimensional F = 2 condensate with N particles confined by the harmonic trapping potential - BEC were analyzed only for the case of the absence of d magnetic field and rotation [22] (see also Refs. [23, 24]). mω2r2 n Spinor F = 2 BEC has one more interaction parameter U(r)= ⊥ , (1) o 2 as compared to the F = 1 BEC because of the larger c where ω is a trapping frequency, m is the mass of the v: spin value. Therefore, there are three possible phases in atom, an⊥d r is the radial coordinate. The system is ro- absenceofmagneticfield: polar,ferromagneticandcyclic i tated with the angular velocity Ω. The energy of the X states. system depends on three interaction parameters α, β, r Due to the internal degree of freedom, a rich variety and γ, which can be defined as [22] a of exotic vortices have been proposed in F = 1 spinor BEC’s by a large number of authors [25]. For instance, 1 F=1 spinor BEC’s with the ferromagnetic spin interac- α= 7(4g2+3g4), (2) tionexhibitSO(3)symmetryinspinspace,whichmeans 1 thatthelocalspinsmaysweepthe wholeorhalfthe unit β =−7(g2−g4), (3) sphere. It has been found that this topological excita- 1 2 tion, called the Mermin-Ho vortex, can be stabilized in γ = (g0 g4) (g2 g4), (4) 5 − − 7 − the rotating system [8]. In the case of F = 2 spinor where (q =0,2,4) BEC’s, the possibility of such the coreless vortex state has been predicted by Zhang et al. [26]. However,in the 4π~2 possiblekinds ofatoms,suchas87Rband23Na, the esti- gq = aq (5) m 2 abnetdwaeqenisatthome sscwatittheritnhgelteontgatlhsspicnha0r,a2c,tearnizdin4g. collisions Ψon−d2c=ase√12Ψeiϑ1, Ψ=−11,Ψei0ϑ,,ΨΨ1 =2,0Ψ,0Ψ,Ψ22==√120e,iυΨ,1in=the1 seeiυc-, TheorderparameterinF =2casehasfivecomponents − √2 − √2 and in the third case Ψ0 = 1 and all other components Ψ (i= 2, 1,0,1,2). Thefreeenergyofthesystemcan i − − are equal to zero. Here values of ϑ and υ are arbitrary be written as [6, 7] and the energy is degenerate with respects to them. De- pendingonvaluesofscatteringlengthsa ferromagnetic, α q F = dS Ψ∗jhΨj + 2Ψ∗jΨ∗kΨjΨk cyclic or polar phase has the lowest energy [22]. In fer- Z h romagnetic phase, Θ=0, f = 2; in cyclic phase, Θ=0, +βΨ Ψ (bF ) (F ) Ψ Ψ f =0; in polar phase, Θ|h =i|1, f =0. 2 ∗j ∗l a jk a lm k m h i | | h i γ + Ψ Ψ Ψ Ψ ( 1)j( 1)k Extended Gross-Pitaevskii equations can be obtained 2 ∗j ∗k −j −k − − inastandardwayfromtheconditionofminimumoffree −BzM −i~Ω·Ψ∗j(∇×r)Ψj , (6) energy of the system Eq. (6): where integration is performed over the syst(cid:3)em area, re- peated indices are summed, Bz is the magnetic field, {h−µ+αΨ∗kΨk}Ψj +β{(Fα)lm(Fα)jkΨ∗lΨmΨk} which is treated as a Lagrange multiplier, h and M are +bγ(−1)j(−1)kΨ∗kΨkΨj the one-body Hamiltonianandmagnetization,whichare i~Ω rΨ B jΨ =0, (12) j z j − ·∇× − given by b where a chemical potential µ is interpreted as the La- ~2 2 grange multiplier. We use the total number of particles h= ∇ +U(r), (7) − 2m N = dSΨiΨ∗i andthemagnetizationM asindependent variables. 2 bM = dS Ψ i. (8) R i | | Z Here F (a=x,y,z) is the angular momentum operator a and it can be expressed in a matrix form as III. VORTEX PHASES AND ENERGY 0 2 0 0 0 A. Classification of phases 2 0 √6 0 0 1  F = 0 √6 0 √6 0 , (9) x 2 Five nonlinearGross-Pitaevskiiequations Eq.(12)are 0 0 √6 0 2   coupled and they can be solved numerically. However, 0 0 0 2 0   some important consequences can be derived from the 0 2 0 0 0 preliminaryanalysisoftheseequations. Inthispaper,we − 2 0 √6 0 0 consider only axially-symmetric solutions of the Gross- F = i 0 √6 −0 √6 0 , (10) Pitaevskiiequations. Inthiscase,eachcomponentofthe y 20 0 √6 −0 2 order parameter Ψj can be represented as  −  0 0 0 2 0    Ψ (r,ϕ)=f (r)exp( L ϕ), (13) 2 0 0 0 0  j j − j 0 1 0 0 0 whereϕisapolarangleandL isawindingnumber. Ax- j Fz =0 0 0 0 0 . (11) ial symmetry of the solution implies that there are some 0 0 0 −1 0  constraints for the possible sets of Lj. It can be shown 0 0 0 0 2 from Eq. (12) that L obeys the following equations  −  j   It is convenient to introduce two additional or- L2+L1+L 1+L 2 4L0 =0, (14) der parameters f = Ψ F Ψ / Ψ2 and Θ = − − − h i ∗j jk k | | L2+L 2 L 1 L1 =0, (15) ( 1)jΨ Ψ / Ψ2 characterizingferromagneticordering − − − − a−nd form∗j a−tijon|of| singlet pairs, respectively [22]. In the L2−2L−2−2L1+2L−1 =0. (16) absence of magnetic field and rotation, BEC can be Eqs.(14)-(16)wereobtainedundertheconditionthatall in three different states [22] that is easily seen from the five components of the order parameter are nonzero. Eq. (6). These states are called ferromagnetic, cyclic Ifsome ofthese componentsare equalto zero identically and polar [22]. In ferromagnetic phase, only one com- then other possibilities appear for the L values. We j ponent of the order parameter is nonzero: Ψ 2 = 1. In list here all the possible phases different from ordi- cΨy2cli=c p21hea−sieθ,,Ψw−h1e,reΨ1θ =is 0ananadrbΨit−ra2ry=ph12eaisθe,−(Ψen0e=rgy√1o2f, n(1a,r×y,0v,o×rt,e−x-1f)r,ee stat(e×:,0,(×1,,11,,1×,)1,,1), (−(×1,,1×,,×0,,0×,,×1)),, the system is degenerate with respect to θ). In po- (0, ,1, ,2), (2, ,1, ,0), ( ,2, ,1, ), × × × × × × × lar phase, there are three possibilities: in the first case ( ,1, ,2, ), ( 2, 1,0,1,2), (2,1,0, 1, 2), × × × − − − − 3 ( 1,0,1,2,3), (3,2,1,0, 1), (4,3,2,1,0), (0,1,2,3,4). C. Results and discussion − − Here numbers denote values of L , ” ” denotes zero j × value of the corresponding component of the order We consider the situation, when the concentration of parameter. We restricted ourselves only to the cases, atomsinz-directionisequalto2000µm 1,andthescat- − when the largest winding number does not exceed 4, tering length a0 equals 5.5 nm. According to the esti- vorticeswith higher winding numbers are assumedto be mates made in Ref. [22], 23Na, 83Rb, 87Rb, 85Rb corre- nonstable. spondtopointsonthephasediagramina2 a4vs. a0 a4 − − plane (in absence of magnetic field and rotation), which aresituatedin the close vicinity to the phase boundaries B. Method betweenferromagnetic,polarandcyclicstates(seeFig.1 in Ref. [22]). There are even some uncertainties in po- Now we can find the solutions of the Gross-Pitaevskii sitions of these points on the phase diagram because of equationsforeachphaselistedabove. Forthesolutionof the error bars in a . It was assumed in Ref. [22] that q Eq. (12), we apply a numerical method, which was used 23Na,83Rb,and85RbBEC’sareinpolar,ferromagnetic, before in Ref. [27]. and cyclic states, respectively, and 87Rb corresponds to the phase boundary between the polar and cyclic states. Besides the numerical solution, we also use a varia- In this paper, we perform all the calculations for the po- tionalansatz basedontrialfunctionsforeachcomponent lar state at β = α, γ = α. For the cyclic state we of the order parameter. It follows from Eq. (12) that put β = α, γ = 5α0. And−fo5r0the state situated on the 50 50 each component of the order parameter has an asymp- phase boundary between cyclic and polar state, which totic fj(r) rLj at r 0 and that in the expansion of we call cyclic+polar, we use β = α, γ = 0. We calcu- ∼ → 50 fj(r) in powers of r there are only terms proportionalto latedthe dependencesoftheenergyofthe systemonthe rLj+2n, where n 0 is an integer number. Superfluid magnetization for these three phases for different vortex ≥ density vanishes at infinite distances from the center of states, when the system is rotated with the frequency the potential well and fj(r) 0 at r . Therefore, Ω= 0.4ω . Our results obtained by the numerical solu- → → ∞ we have chosen the following trial function tion to th⊥e Gross-Pitaevskii equations are presented on Fig.1forcyclic(a),polar(b),andcyclic+polar(c)states. r Lj r Lj+2 r2 In Table 1 we show numerically and variationally calcu- f (r)=C +a exp , j j"(cid:18)Rj(cid:19) j(cid:18)Rj(cid:19) # −2Rj2! lated values of the magnetization, at which the transi- (17) tionsoccurbetweendifferentphasesforthecaseofcyclic where C , a andR are variationalparameters. Param- state (Fig. 1 (a)). In this table, ”a” denotes the transi- j j j eters C are not completely independent since they are tionbetween( 1,0,1,2,3)and( ,0, ,1, )states,”b” j − × × × relatedby one equationthat is a condition of equality of between ( ,0, ,1, ) and ( 2, 1,0,1,2) states, and × × × − − number of atoms to the given number. ”c” between ( 2, 1,0,1,2) and ( 1, ,0, ,1) states. − − − × × In the limit of noniteracting gas, Gross-Pitaevskii There is a good agreement between the numerical and equation becomes linear with respect to the order pa- variational results. rameter. Inthis case,itis easyto see thata =0.Inthe Ingeneralcase,therecanbephasedifferencesbetween j Thomas-Fermiregimethisisnolongervalidandthesys- functions fj(r). Axial symmetry of the solution implies tem tries to minimize its energy by changing parameters thattherearesomeconstraintsonphases,whicharesim- a and R as compared to the limit of an ideal gas. Ac- ilar for the constraints on winding numbers and can be j j cording to our estimates, for the case of one-component also obtained from the GP equation (12). The energy of orderparameter,trialfunction(17)isabletogiverather the system is then degenerate with respect to remaining accurate results for the energy and for the rotation fre- phases, as in the nonrotating case, which was discussed quency of transition from the vortex-free to the vortex above. stateeveninthe’moderate’Thomas-Fermilimit. There- In all states under study, it also turns out that for the fore, in this paper we apply the method to the case of condensate it is energetically favorable to be distributed five-component order parameter. Note that variational between one, two of three hyperfine states and not be- approaches were applied before for the analysis of vortex structures in mesoscopic superconductors within the Ginzburg-Landau theory, see e.g. [28], and for vortices Table I: Values of the magnetization corresponding to the in scalar and spinor BEC [29, 30]. transitions between different ground states for the case of Using Eq. (6) and the normalizationcondition for the cyclic phase (Fig. 1(a)), which were calculated numerically order parameter one can find the energy of the system and variationally. analytically as a function of variational parameters for eachsetofLj. However,finalexpressionfortheenergyis a b c rathercumbersomeandwedonotpresentithere. Values numerical 0.34 0.79 1.36 of variational parameters can be calculated by a numer- variational 0.38 0.76 1.40 ical minimization of the energy. 4 tween four or five. zation, depends on many factors. For instance, at high It can be seen from Fig. 1 that, in all cases, at zero magnetization, close to M/N = 2, condensate has to be magnetization,stateswith nonzerowinding numbersare concentratedmostlyinthestatewithmF =2. Sinceour energeticallyfavorabledue tothe rotationofthe system. system is rotated with the frequency enough to create a Bychangingthemagnetizationitispossibletojumpfrom vortex,for the condensate it is favorableenergetically to onevortexphasetoanotherone. Thetransitionsbetween havewinding number 1in this state. For otherparticles, different phases are discontinuous. Note that in this pa- which are not in mF = 2 state, it is energetically favor- per we consider only the thermodynamical transitions able to be in a superfluid phase with winding number 0 between different states. Actual position of transition in order to occupy the inner part of the trap, where the linebetweendifferentvortexphasesiscontrolledbylocal trappingpotentialissmall. Thatiswhyinallthreephase stability of states and, therefore, the prehistory of the diagramspresentedinFig.1vortexphase( 1, ,0, ,1) − × × system. hasthelowestenergyathighmagnetization. Thedepen- dences of the density of particles in each hyperfine state Which state has the lowest energy at given magneti- on the distance from the potential well center is shown on Fig. 2(a) for the case of polar phase at M/N =1.87. Spatial variations of the order parameters Θ and f |h i| arepresented onFig.2(b). ItcanbeseenfromFig.2(b) that Θ is maximum in the center of the potential well and vanishes at the infinity. At the same time, f has |h i| a minimum at r = 0 and tends to 2 at r . This is → ∞ because close to r = 0 a component with m = 0 has F a maximum and all other components are small. There- fore, Θ, which characterizes a formation of singlet pairs, has a maximum at r = 0, and f = 0. Far from the |h i| center of the cloud, the density of particles with m =2 F is much largerthan that for m =0 and the densities in F allotherhyperfine statesareverysmall. Forthis reason, Θ=0 and f =2 at r . |h i| →∞ For lower magnetization, condensate has to be distributed between the states with m = 0,1,2. It turns F outthatagaininallphasediagramsinFig.1theground stateisrepresentedbythephase( 2, 1,0,1,2)inrather − − broad range of M. This is due to the fact that, in this case, it is favorable to put most of the particles in the state with winding number 0 in order to occupy the inner part of the trap. Most of remaining particles are condensed to the state with winding number 2 occupy- ing the outer part of the trap and thus decreasing the energy of interaction (in ”density” channel) of particles withm =0and2. Typicalprofilesoftheparticlesden- F sitiesindifferenthyperfinestatesforthisvortexphaseare shown in Fig. 2(c) for the cyclic state at M/N = 1.21. Fig. 2(d) indicates r-dependences of the order parameters Θ and f , which physically are similar to that in |h i| the case of ( 1, ,0, ,1)state shown in Fig. 2 (a) and − × × (b). At lower M particles have to be distributed between several states with different m . In this case, spin inter- F actions become important and, therefore, the sequences ofphasetransitionsfordifferentphasediagramsinFig.1 are different. We can conclude that at high magnetiza- tions, M & 1, the state with lowest energy is mostly determined by interactions in the ”density” channel, whereas at low magnetization M 0, spin interactions ∼ Figure 1: (color online) Dependencesof energies of different play an important role. In Fig. 2(e) and (f) we present vortex phases on the magnetization for cyclic (a), polar (b), the r-dependences of the superfluid density in different and polar+cyclic (c) states. The energy is measured in units hyperfinestatesandorderparametersΘand f forthe of ~ω⊥. cyclic+polar ( 1,0,1,2,3) state at M/N = 0|h.1i8|, where − 5 Figure 2: The spatial variation of the density of particles in different hyperfine states normalized by the total density at the system axis and the order parameters Θ and |hfi| for different vortex phases. Fig. 2 (a) and (b) correspond to the (−1,×,0,×,1) polar phase at M/N =1.87; Fig. 2 (c) and (d) to the(−2,−1,0,1,2) phase at M/N =1.21; Fig. 2 (e) and (f) tothecyclic+polar(−1,0,1,2,3)stateatM/N =0.18. Dotlinesshowthetotaldensityofparticles. Totalnumberofparticles is 10000. this phase has the lowest energy. In this case, particles thepurepolarstateforms(Θ=1),becausethespintex- are distributed between states m = 1 and m = 1. ture exhibits the 2-dimensional radial disgyration in the F F − The state with m = 1 has a winding number 0, and coreregion. Note that( 2, 1,0,1,2)and( 1,0,1,2,3) F − − − − theseparticlesoccupyaspacewithminimaltrappingpo- vortices were not described before in a literature. tential. All other particles are condensed in the state with winding number 2 thus decreasing the interaction energy in the ”density” channel. As a result, the order parameter Θ is nonzero everywhere except of the point r =0, since at any r >0 there are particles in the states IV. CONCLUSIONS with m = 1,andthe formationofsingletpairsis pos- F sible. The o±rder parameter f is nonzero at r =0 and |h i| r , since inboth casesthere areparticles condensed Insummary,weanalyzedthevortexstructureinspinor →∞ in the states with nonzero mF. At the same time, at F =2 condensateusing extendedGross-Pitaevskiiequa- small values of r most of particles are condensed in the tions. We considered only axially-symmetric vortices. state with mF = 1 and at larger r in the state with Based on symmetric configurations, all possible vortex − mF = 1. Therefore, there is an abrupt change in the states were classified. The Gross-Pitaevskii equations spindirectionatintermediate values ofr. This resultsin were solved both numerically and by the variational the vanishing of f near the vortex core. method using trial functions for the order parameter. |h i| In Fig. 3 we show the spin texture, f and f , for Spatialdistributionoftheparticlesdensityandtheorder x y the phases ( 1,0,1,2,3) and ( 2, 1h,0,i1,2) ihn ciyclic parameters Θ and f were obtained. Energies of dif- − − − |h i| state. In the first case, in the vortex-core region, the ferent vortex phases were found as a function of magne- spinsarepolarizedalongthe z-direction. Intheouterre- tization, when the system is rotated with the frequency gion the spin vanishes and Θ grows, where the spin am- Ω = 0.4ω . We found that at high magnetization the plitudehasanodeandthepurepolarstateforms. Inthe energy of⊥the system is mostly determined by the inter- second case, in the outer region, the spins are polarized action in the ”density” channel, whereas at low magne- along the z-direction. In the core-region, the spins lean tization spin interaction plays an important role. Also, toward the origin. At the origin, the spin vanishes and two new types of vortices were described. 6 Figure 3: Spin textures, hfxi and hfyi, for the cyclic (−1,0,1,2,3) phase at M/N = 0.16 (a) and the cyclic (−2,−1,0,1,2) phase at M/N =1.21 (b). Acknowledgments W. V. Pogosov is supported by the Japan Society for the Promotion of Science. [1] M.R. Matthews, B. P. Anderson, P. C. Haljan, D. S. [17] J. Ruostekoski and J. R. Anglin, Phys. Rev. Lett. 91, Hall, C.E. Wieman and E.A. Cornell, Phys. Rev. Lett. 190402 (2003). 83, 2498 (1999). [18] H. Schmaljohann, M. Erhard, J. Kronjäger, M. Kottke, [2] K. W. Madison, F. Chevy, W. Wohlleben and J. Dal- S.vanStaa,L.Cacciapuoti, J.J.Arlt,K.Bongs,andK. ibard, Phys.Rev.Lett, 84, 806 (2000). Sengstock, Phys. Rev.Lett. 92, 040402 (2004). [3] J. R. Abo-Shaeer, C, Raman, J. M. Vogels and W. Ket- [19] M.-S. Chang, C. D. Hamley,M. D.Barrett, J. A.Sauer, terle, Science 292, 476 (2001). K. M. Fortier, W. Zhang, L. You, and M. S. Chapman, [4] P.C.Haljan,I.Coddington,P.EngelsandE.A.Cornell, Phys. Rev.Lett. 92, 140403 (2004). Phys.Rev.Lett. 87, 210403 (2001). [20] A. Görlitz, T. L. Gustavson, A. E. Leanhardt, R. Löw, [5] J. Stenger, D. M. Stamper-Kurn, H. J. Miesner, A.P. A. P. Chikkatur, S. Gupta, S. Inouye, D. E. Pritchard, Chikkatur,W. Ketterle, Nature, 396, 345 (1999). and W. Ketterle, Phys.Rev. Lett. 90, 090401 (2003). [6] T. Ohmi and K. Machida, J. Phys. Soc. Jpn., 67, 1822 [21] T. Kuwamoto, K. Araki, T. Eno, and T. Hirano, Phys. (1998). Rev. A 69, 063604 (2004). [7] T.- L. Ho,Phys. Rev.Lett., 81, 742 (1998). [22] C. V. Ciobanu, S.-K. Yip, and Tin-Lun Ho, Phys. Rev. [8] T. Mizushima, K. Machida, and T. Kita, Phys. Rev. A 61, 033607 (2000). Lett., 89,, 030401 (2002). [23] M. Koashi and M. Ueda, Phys. Rev. Lett. 84, 1066 [9] S.-K.Yip,Phys. Rev.Lett. 83, 4677 (1999). (2000). [10] T.Isoshima,K.MachidaandT.Ohmi,J.Phys.Soc.Jpn. [24] M.UedaandM.Koashi,Phys.Rev.A65,063602(2002). 70, 1604 (2001). [25] Seeforreview,K.Kasamatsu,M.Tsubota,andM.Ueda, [11] T. Isoshima and K. Machida, Phys. Rev. A 66, 023602 Int.J. Mod. Phys. B 19, 1835 (2005). (2002). [26] H. Mäkelä, Y. Zhang, and K.-A. Suominen, J. Phys. A: [12] T.Mizushima,N.KobayashiandK.Machida,Phys.Rev. Math. Gen. 36, 8553 (2003). A 70, 043613 (2004). [27] T. Mizushima, K. Machida, and T. Kita, Phys. Rev. A [13] J. W. Reijnders, F. J. van Lankvelt, K. Schoutens, and 66, 053610 (2002). N.Read, Phys.Rev. Lett. 89, 120401 (2002). [28] W. V. Pogosov, Phys. Rev.B 65, 224511 (2002). [14] J.-P.Martikainen,A.Collin, andK.A.Suominen,Phys. [29] A. L. Fetter, J. Low Temp. Phys. 65, 224511 (1996). Rev.A 66, 061601 (2002). [30] K.Kasamatsu,M.Tsubota, andM.Ueda,Phys.Rev.A [15] T. Kita, T. Mizushima, and K. Machida, Phys. Rev. A 71, 043611 (2005). 66, 053610 (2002). [16] E. J. Mueller, Phys. Rev.A 69, 033606 (2004).