Phys.Rev. A74,065602(2006) Extension of theThomas-Fermi approximation for trappedBose-Einstein condensateswith an arbitrary number of atoms A. Mun˜oz Mateo∗ and V. Delgado† DepartamentodeF´ısicaFundamental II,UniversidaddeLaLaguna, LaLaguna, Tenerife,CanaryIslands,Spain By incorporating the zero-point energy contribution we derive simple and accurate extensions of the 7 usual Thomas-Fermi (TF) expressions for the ground-state properties of trapped Bose-Einsteinconden- 0 satesthatremainvalidforanarbitrarynumberofatomsinthemean-fieldregime.Specifically,weobtain 0 approximate analytical expressions for the ground-state properties of spherical, cigar-shaped, and disk- 2 shaped condensates that reduce to the correct analytical formulas in both the TF and the perturbative n regimes,andremainvalidandaccurateinbetweenthesetwolimitingcases.Mean-fieldquasi-1Dand-2D a condensatesappearassimpleparticularcasesofourformulation. Thevalidityofourresultsiscorrobo- J ratedbyanindependentnumericalcomputationbasedonthe3DGross-Pitaevskiiequation. 9 2 PACSnumbers:03.75.Hh,05.30.Jp,32.80.Pj ] r e The experimental realization of Bose-Einstein conden- ground state, ψ(r) = (πa2)−3/4exp( r2/2a2), and the r − r h sates(BECs)ofdiluteatomicgasesconfinedinopticaland chemicalpotentialsatisfies ot magnetictraps[1,2,3]hasstimulatedgreatactivityinthe . characterization of these quantum systems. Of particular (3/2)~ω +gn¯ = µ (3) t a interest are the ground-state properties of trapped BECs m wheren¯ = N/(√2πa )3 isthemeanatomdensity. Away with repulsive interatomic interactions [4]. These proper- r - tiesderivefromthecondensatewavefunctionψ(r)which, fromthesetwolimitingcases,inprinciple,onehastosolve d theGPEnumerically. Veryfewtheoreticalworkshavead- inthezero-temperaturelimit,satisfiesthestationaryGross- n dressed the question of looking for approximate analyti- o Pitaevskiiequation(GPE)[5] calsolutionsvalidinbetweenthetwoanalyticallysolvable c [ regimes. Themostrelevantproposalsarebasedonavaria- ~2 1 2+V(r)+gN ψ 2 ψ = µψ, (1) tionaltrialwavefunction[6],oronthesemiclassicallimit −2m∇ | | oftheWignerphase-spacedistributionfunctionofthecon- v (cid:18) (cid:19) 5 where N is the number of atoms, g = 4π~2a/m is densate[7]. However,thepracticalusefulnessoftheseap- 9 proaches turns out to be somewhat limited in comparison the interaction strength, a is the s-wave scattering length, 6 withthesimpleTFapproximation. 1 V(r) = 1m(ω2r2 +ω2z2) is the harmonic potential of 2 ⊥ ⊥ z In this work we address the above question from a dif- 0 theconfiningtrap,andµisthechemicalpotential. ferent point of view. We start from the usual TF approxi- 7 OnlyintwolimitingcasescanEq. (1)besolvedanalyti- 0 mationandmodify itconvenientlyto account,in a simple cally: intheThomas-Fermi(TF)andperturbativeregimes. / manner, for the zero-point energy contribution. This en- t When N is sufficiently large that µ ~ω , ~ω , one a ≫ ⊥ z ables us to derive simple and accurate extensions of the m enters the TF regime. In this case the kinetic energy can TF expressions that remain valid for an arbitrary number beneglectedincomparisonwiththeinteractionenergyand - ofatomsinthemean-fieldregime. Specifically,weobtain d the GPE reduces to a simple algebraic equation. Useful general analytical expressions for the ground-state prop- n analytical expressions can then be obtained for the con- o erties of spherical, cigar-shaped,anddisk-shapedconden- densate ground-state properties [4]. In the simple case c satesthatreduceto thecorrectanalyticalformulasinboth of a spherical trap characterized by an oscillator length : theTFandthe perturbativeregimes,andremainvalidand v a = ~/mω,Eq. (1)leadsintheTFlimitto i r accurateinbetweenthesetwolimitingcases. X We begin by considering a BEC in a spherical trap. In p ar 1mω2r2+gN ψ(r) 2 = µ, 0 r R (2) principle, we start from the TF relation of Eq. (2). How- 2 | | ≤ ≤ ever, since we intend to apply this equation to arbitrarily small condensates, we introduce a lower cutoff radius r , where the condensate radius R = 2µ/~ωar is deter- defined through 1mω2r2 = 3~ω, in order to be consis0- mined from the condition |ψ(r)|2 ≥p0, and the chemical tent with the fact2that th0e cont2ribution from the harmonic potential µ = 1 (15Na/a )2/5~ω follows from the nor- oscillatorenergycannotbesmallerthanthezero-pointen- 2 r malizationofψ(r). ergy. As for the small volume V a3 correspondingto 0 ∼ r In the opposite limit, when N is small enough that the r r ,wedonotaspiretogetapreciseknowledgeofthe 0 ≤ interaction energy can be treated as a weak perturbation, wavefunction therein. Instead, we contentourselveswith oneentersthe(idealgas)perturbativeregime. Inthiscase, an effective condensate density n¯ in that region. As we 0 tothelowestorder,ψ(r)isgivenbytheharmonicoscillator shall see, this is all that is needed to obtain very accurate 2 approximateformulas for mostof the condensateground- 5 20 (a) (b) stateproperties. Thuswestartfromtheansatz R 4 µ 15 1 2 mω2r2+gN ψ(r) 2 = µ, r0 < r R (4a) 3 µ 2 | | ≤ µ 1,5 3 3 TF 0 1 10 R 2~ω+g 6/πn¯0 = µ, 0 ≤ r ≤ r0 (4b) 2 εpot ε 2 RTF with ψ(r) = 0 fopr r > R. A renormalization constant int 5 1 1 κ−1 6/π hasbeenintroducedin Eq. (4b)to guaran- ε 0 tee th≡e cporrect perturbative limit. In this limit µ → 32~ω 00 5 10 1k5in 20 00 0 55 10 5 and R r0 = √3ar. Under these circumstances, only χ 1xχ10 2x10 Eq. (4b→) contributes significantly to the chemical poten- 0 0 tial, and in this case n¯ = N/V . This corresponds to a 0 0 FIG.1:(Coloronline)Theoreticalpredictionfortheground-state uniformsphericalcondensate,definedinthefinite volume propertiesofsphericalcondensates(solidlines).Theopencircles V . In order for this uniform density to producethe same 0 aretheexactnumericalresults.Forcomparisonpurposeswehave chemicalpotentialasthegroundstateoftheharmonicos- alsoincludedtheTFprediction(dashedlines). cillatoroverthevolumeoftheentirespaceitisonlyneces- sary to renormalizethe correspondinginteractionstrength bymultiplyingby 6/π. Equations(4)alsoyieldthecor- satisfiesEq. (5)witharesidualerror[8]smallerthan0.7% rectresultintheTFregime. Thisismainlyaconsequence p foranyχ [0, ). Figures 1(a) and1(b)show, respec- ofthedirectrelationexistingbetweenthenumberofparti- 0 ∈ ∞ clesandthesizeofatrappedBEC.Forlargecondensates, tively,thepredictedchemicalpotential,µ = 1R2,andcon- such that µ ~ω, one has R r and, as a result, the densate radius, obtained from Eq. (6) (solid2lines), along 0 ≫ ≫ relativecontributionfromEq. (4b)tothenormalizationin- with the exact results obtained from the numerical solu- tegralthatdeterminesµbecomesnegligible.Sincewehave tion of the 3D GPE (opencircles). For the numericalcal- renouncedanexplicitexpressionforψ(r)inV ,inthisre- culation we have defined the radius through the condition 0 spect,ourapproachcannotprovidemoreinformationthan ψ(R) 2 = 0.05 ψ(0) 2. With this definition, Eq. (6)re- theTFapproach. OnlywhenR r0 canwehaveasuffi- |produc|esthenum|erical|Rwitharelativeerrorsmallerthan ≫ cientlypreciseknowledgeofthewavefunctionand,inthis 3% for any χ . Most of the error, however, comes from 0 case,itcoincideswiththeTFwavefunction. theregionwhereχ 1(TFlimit)becauseinthatregion 0 The chemical potential follows from the normalization ≫ 2 R R andit rathersatisfies ψ(R) = 0. Theaccu- TF ofψ(r). Afterastraightforwardcalculationoneobtains → | | racywithrespecttothenumericalµisbetterthan0.5%. Astraightforwardcalculationyieldsthemean-fieldinter- 1 5 √3 2 3√3 3 a actionenergyperparticle,ǫint ǫint/~ω Eint/N~ω, R + (κ 1)R κ = N , (5) ≡ ≡ 15 2 − − 2 − 5 a where R R/a , and µ = 1R2(cid:18)is the c(cid:19)hemical proten- ǫint = 8χ1 1085R7+√3(κ−1)R4 tial in uni≡ts of ~ωr. As Eq. (25) shows, the ground-state 0 (cid:20) properties depend on the sole parameter χ Na/a . 6√3 κ 3 R2+9√3 κ 3 . (7) 0 ≡ r − − 5 − 7 When χ 1 (TF limit) the above equation leads to (cid:18) (cid:19) (cid:18) (cid:19)(cid:21) 0 ≫ µ = 1(15χ )2/5, as expected. The χ 1 limit cor- For χ 1, one recovers the TF result, ǫ = (2/7)µ. 2 0 0 ≪ 0 ≫ int responds to the perturbative regime and, in this case, one 2 In the χ 1 limit, using that R = 3+2 2/πχ 0 0 obtainsµ = 3/2+ 2/πχ ,whichisnothingbuttheper- ≪ − 0 (1/π)(1/9+ 2/3π)χ2+O(χ3)isaperturbativesolution turbative result (3). For arbitrary χ , in principle one has 0 0 p p 0 ofEq. (5),oneobtainsǫ = χ ~ω/√2π =gn¯/2,which tosolvenumericallytheabovequinticpolynomialequation p int 0 againisthecorrectresult. Finally,thekineticandpotential (which has only one physically meaningfulreal solution). energies can be readily obtained in terms of the previous This is a simple task thatcan becarried outwith standard resultsbyusingtheexactrelations[5] mathematicalsoftwarepackages.Wehavefound,however, aratheraccurateapproximatesolution.Itcanbeshownthat ǫ E /N = µ/2 (7/4)E /N, (8a) kin kin int theexpression ≡ − ǫ E /N = µ/2 (1/4)E /N. (8b) pot pot int −1 ≡ − 2 1 1 π/2 R = 3+ (15χ0)25 + 25 + 72χ101/15+10 + p2χ0 ! aInndFiǫgp.o1t,(ao)bwtaeinsheodwfrtohmethEeqosr.et(i6ca)–l(p8r)ed(iscotliiodnlfionreǫs)in,ta,lǫokning, (6) withtheexactnumericalresults(opencircles). 3 NextweconsideraBECconfinedinacigar-shapedmag- 5 2 netic trap with oscillator lengths a⊥ = ~/mω⊥ and χ1 = 1.0 az = ~/mωz and an aspect ratio λ p= ωz/ω⊥ ≪ 4 µ 2 1 2. We shall restrict ourselves to the mean-field regime, which rpequires Nλa2/a2 1 [9, 10, 11]. As be- 3 µ µ fore, we startfrom the⊥usual≫TF expression,whichwe as- TF 10 1 0√ λ sume to be valid up to a minimum radial distance r0 = 2 ε z ⊥ pot ε √2a⊥, determined from the condition that the contri- int -1 1 bution from the radial harmonic oscillator energy should ε not be smaller than ~ω . This defines an outer region kin ⊥ V (r ,z): r2/R2+z2/Z2 1 r > r0 , 0 -2 w+hich≡isn{oth⊥ingbutt⊥heusualTFellTipFso≤idal∧dens⊥ityclo⊥u}d, 0 5 1χ0 15 20 0 0a,3 n 0,6 1 1 truncated at r = r0. Note that unlike what happens ⊥ ⊥ with the condensate radius R = 2µ/~ω a , which FIG.2:(Coloronline)Theoreticalpredictionfortheground-state ⊥ ⊥ remains the same, now the axial condensate half-length properties of arbitrary cigar-shaped condensates with λ ≪ 2 Z = 2(µ/~ω 1)a /√λ cpoincides with the TF (solid lines). The open circles are exact numerical results ob- value Z only in⊥ −the limzit µ/~ω 1. For large tainedwithλ=0.2. ThedashedlineistheTFprediction. TF ⊥ condenspates, when µ ~ω (TF re≫gime), this is the ⊥ ≫ only region that contributes significantly. On the con- trary, in the perturbative regime, as µ ~ω most of µ = 1+1(3χ )2/3andZ = λ−1/2(3χ )1/3,inagreement → ⊥ 2 1 1 the contribution comes from the inner cylinder V withpreviousresults[10,11]. Ingeneral,forarbitraryχ , − 1 (r ,z): r r0 z Z . In this case, if a ≡ anapproximatesolutionsatisfyingEq. (10)witharesidual { ⊥ ⊥ ≤ ⊥ ∧ | | ≤ } ≪ a , the transverse dynamics becomes frozen in the radial errorlessthan0.75%foranyχ [0, )isgiven by ⊥ 1 ∈ ∞ groundstateoftheharmonictrapandthecondensatewave functioncanbefactorizedasψ(r ,z) = ϕ(r )φ(z),with −1 ⊥ ⊥ 1 1 1 4 tϕo(ar⊥m)ea=n-fi(πelad2⊥q)u−a1s/i2-1eDxpc(o−ndr⊥e2n/s2aate2⊥.)S.uTbshtiistuctoinrgretshpeonndins √λZ = (15χ1)54 + 13 + 57χ1+345 + (3χ1)43! Eq. (1)andintegratingouttheradialdynamics,onefinds (11) Themean-fieldinteractionenergyǫ ǫ /~ω is 1 int int ⊥ ~ω + mω2z2+g N φ(z)2 = µ, (9) ≡ ⊥ 2 z 1D | | 1 1 where g = g/2πa2 [12], and we have used that µ ǫint = (√λZ)5+ (√λZ)7 . (12) ~ω 1D1~ω to neg⊥lect the axial kinetic energy. No∼te 15χ1 (cid:18) 7 (cid:19) ⊥ ≫ 2 z Forχ 1,Eq. (12)reducestoǫ = (2/7)µ, whilefor that g can be conveniently rewritten as gn¯ with n¯ = 1 int 1/π(r10D)2,indicatingthatonecanaccountfor2thecontr2ibu- χ1 1≫, itleadstoǫint = (2/5)(µ ~ω⊥), whichagain ⊥ are≪thecorrectanalyticallimits. Asfo−rthecondensateden- tionfromtheradialgroundstatebyusingauniformmean 2 densityperunitareanormalizedtounityinV−. Guidedby sity per unit length, n1(z) ≡ N 2πr⊥dr⊥|ψ(r⊥,z)| , thesesimpleideas,wethenproposethefollowingansatz: afterastraightforwardcalculationonefinds R (√λZ)2 z2 (√λZ)4 z2 2 n (z) = 1 + 1 12mω⊥2r⊥2 + 12mωz2z2+gN |ψ(r⊥,z)|2 =µ, r∈ V+ 1 4a (cid:18) − Z2(cid:19) 16a (cid:18) − Z2(1(cid:19)3) 1 ThefirsttermisthecontributionfromV andthusitisthe ~ω + mω2z2+gNn¯ φ(z)2 =µ, r V − ⊥ 2 z 2| | ∈ − onlyonethatcontributessignificantlyintheχ1 1limit. ≪ Onthecontrary,thesecondterm,whichisthecontribution withψ = 0elsewhere. Thenormalizationofψleadsto from V , gives the dominant contribution in the χ 1 + 1 ≫ limit,ingoodagreementwithpreviousresults[11]. 1 1 a Figure 2 shows the theoretical predictions for the (√λZ)5+ (√λZ)3 = Nλ , (10) 15 3 a ground-state properties of arbitrary cigar-shaped conden- ⊥ sates with λ 2, obtained from Eqs. (8) and (11)–(13) where Z Z/a and R R/a . The chemical poten- ≪ z ⊥ (solid lines), alongwith exactnumericalresults (opencir- ≡ ≡ tial µ ≡ µ/~ω⊥ is given by µ = 1 + 21(√λZ)2. Now cles). the relevantparameterdeterminingthe ground-stateprop- Finally, we consider a BEC in a disk-shaped trap with (e1rt0i)esleiasdχs1to≡µN=λa1/(1a5⊥χ. W)2h/e5nanχd1Z≫=1λ(T−F1/r2e(g1i5mχe),)1E/q5.. bλat≫ive2reagnidmaez, w≫hicah. oIncctuhrisscwahseen, iχnthemNeaan/-λfi2ealdpetu1r-, 2 1 1 2 ≡ z ≪ Whenχ 1(mean-fieldquasi-1Dregime), oneobtains thesystemreducestoaquasi-2Dcondensatesatisfying 1 ≪ 4 5 2 χ = 1.0 1 1 2 2~ωz + 2mω⊥2r⊥2 +g2DN|ϕ(r⊥)|2 = µ, (14) 4 1,5 µ 1,5 µ wpmgeκiart−2lhiu1znn¯agit21tiD,olenwn=fghatechrtegoa/rnn¯,√da1nκ2d=π−2p1ar1≡zo/p[2o1as3ez2].t/ihsπeWaifsoeultnlhtoihefweoainrpnmpgrreaomwnpersriaaitanetzte:dgre2enDnsoiatrys- 23 µεpToFt 0 0 ε µ TF 1 1r/ λ√⊥ p int 0,5 1 1 1 ε mω2z2+ mω2r2 +gN ψ(r ,z) 2 =µ, r V kin 2 z 2 ⊥ ⊥ | ⊥ | ∈ + 0 0 0 5 10 15 20 0 0,1 0,2 0,3 1 1 χ a a n ~ω + mω2r2 +gκ−1Nn¯ ϕ(r )2 =µ, r V 2 z 2 2 z 2 ⊥ ⊥ 2 1| ⊥ | ∈ − (15) FIG.3:(Coloronline)Theoreticalpredictionfortheground-state propertiesofarbitrarydisk-shapedcondensateswithλ≫2(solid with ψ = 0 elsewhere. In the above equations, V (r ,z): r2/R2 +z2/Z2 1 z > z andV+ ≡ lines).Theopencirclesareexactnumericalresultsobtainedwith { ⊥ ⊥ TF ≤ ∧ | | 0} − ≡ λ=20. ThedashedlineistheTFprediction. (r ,z): r R z z ,wherez = a ,R = ⊥ ⊥ 0 0 z TF { ≤ ∧ | | ≤ } 2µ/~ω √λa , R = 2(µ/~ω 1/2)√λa , and z ⊥ z ⊥ − Z = 2µ/~ω a . More precisely, one expects κ−1 p z z p 2 → 2/π intheperturbativeregime(χ 1),whileκ−1 1inthpeTFregime(χ 1). Thefin2a≪lresultsareno2tve→ry ξ[2µ (r ) 1] [2µ (r )]3/2 1 spensitive to the specifi2c≫functional form of κ−21. We thus n2(r⊥) = 4zπa⊥az− + z 6π⊥aaz − , (20) proposeoneofthesimplestpossibilities: 2 whereξ (κ 1)and2µ (r ) 1+R (1 r2/R2). ≡ 2− z ⊥ ≡ λ − ⊥ InFig. 3weshowtheground-statepropertiesofarbitrary κ−1(χ ) 2/π +Θ(χ 0.1) 2 2 ≡ 2− disk-shaped condensates with λ 2, obtained from our ≫ p RTF(χ2 = 0.1) analytical formulas [Eqs. (8) and (18)–(20)] (solid lines), 1 2/π 1 , (16) × − − R (χ ) alongwithexactnumericalresults(opencircles). (cid:18) TF 2 (cid:19) (cid:16) p (cid:17) In conclusion, modifying the usual TF approximation where Θ(x) is the Heaviside function and R (χ ) = TF 2 conveniently to account for the zero-point energy contri- (15χ )1/5a is the TF radius. The normalization of ψ 2 ⊥ bution,we havederivedgeneralanalyticalexpressionsfor yields theground-statepropertiesofspherical,cigar-shaped,and disk-shapedcondensatesthatreducetothecorrectanalyti- 4 2 calformulasinboththeTFandthemean-fieldperturbative 1 5 1 R R 1 Na Z + (κ 1) = , (17) regimes and remain valid and accurate in between these 15 8 2− λ2 − 6λ − 15 λ2a z two limiting cases. Mean-fieldquasi-1Dand-2Dconden- whereZ Z/a ,R R/a ,andZ2 R2/λ = 1. The satesappearassimpleparticularcasesofourformulation. z ⊥ chemical≡potentialisµ≡ µ/~ω = 1(1−+R2/λ). This work has been supported by MEC (Spain) and ≡ z 2 FEDERfund(EU)(ContractNo. Fis2005-02886). Forχ 1Eq. (17)leadstotheusualTFresults,while 2 ≫ forχ 1(mean-fieldquasi-2Dregime),oneobtainsµ = 2 ≪ 1/2+(2 2/πχ )1/2 andR = λ1/2(8 2/πχ )1/4. An 2 2 approximatesolutionthatsatisfiesEq. (17)witharesidual errorlesspthan0.95%foranyχ [0, p)isgiven by ∗ Electronicaddress: [email protected] 2 ∈ ∞ † Electronicaddress: [email protected] −1/8 R R/√λ = (1/15χ )8/5+(κ /8χ )2 (18) [1] M.H.Andersonetal.,Science269,198(1995). λ 2 2 2 ≡ [2] K.B.Davisetal.,Phys.Rev.Lett.75,3969(1995). After some calcuhlation one finds the followiing expres- [3] C.C.Bradleyetal.,Phys.Rev.Lett.78,985(1997). sionsforthemean-fieldinteractionenergyǫ ǫ /~ω [4] G.BaymandC.J.Pethick,Phys.Rev.Lett.76,6(1996). int ≡ int z [5] Forareviewsee,forexample,F.Dalfovo,S.Giorgini,L.P. andthecondensatedensityperunitarean (r ): 2 ⊥ Pitaevskii,andS.Stringari,Rev.Mod.Phys.71,463(1999). [6] A.L.Fetter,J.LowTemp.Phys.106,643(1997). 1 8Z7 R6 R4 4R2 8 [7] P.SchuckandX.Vin˜as,Phys.Rev.A61,43603(2000). ǫint = +ξ λ λ λ [8] Given P(R) = χ, we define the residual error associated 8χ2 105 6 − 3 − 15 − 105! withtheapproximatesolutionRεas[P(Rε)−χ]/χ. (19) [9] D.S.Petrovetal.,Phys.Rev.Lett.85,3745(2000). 5 [10] V.Dunjko,V.Lorent,andM.Olshanii,Phys.Rev.Lett.86, [12] M.Olshanii,Phys.Rev.Lett.81,938(1998). 5413(2001). [13] D.S.Petrovetal.,Phys.Rev.Lett.84,2551(2000). [11] C. Menotti and S. Stringari, Phys. Rev. A 66, 043610 (2002).