Table Of Content

Stability of Inhomogeneous Multi-Component Fermi Gases D. Blume,1,2 Seth T. Rittenhouse,3 J. von Stecher,3 and Chris H. Greene3 1Department of Physics and Astronomy, Washington State University, Pullman, Washington 99164-2814 2 JILA, University of Colorado, Boulder, CO 80309-0440 3Department of Physics and JILA, University of Colorado, Boulder, CO 80309-0440 (Dated: February 4, 2008) Two-component equal-mass Fermi gases, in which unlike atoms interact through a short-range two-body potential and like atoms do not interact, are stable even when the interspecies s-wave scattering length becomes infinitely large. Solving the many-body Schrödinger equation within a 8 hyperspherical framework and by Monte Carlo techniques, this paper investigates how the prop- 0 erties of trapped two-component gases change if a third or fourth component are added. If all 0 interspeciesscatteringlengthsareequalandnegative,ourcalculationssuggest thatboththree-and 2 four-component Fermi gases become unstable for a certain critical set of parameters. The relevant n lengthscaleassociatedwiththecollapseissetbytheinterspeciesscatteringlengthandwearguethat a thecollapseis,similartothecollapseofanattractivetrappedBosegas,amany-bodyphenomenon. J Furthermore,weconsiderathree-componentFermigasinwhichtwointerspeciesscatteringlengths 5 are negative while the other interspecies scattering length is zero. In this case, the stability of the 2 Fermi system is predicted to depend appreciably on the range of the underlying two-body potential. Wefindparameter combinations for which thesystem appears to becomeunstable for a finite ] negative scattering length and parameter combinations for which the system appears to be made r e up of weakly-bound trimers that consist of one fermion of each species. h t PACSnumbers: o . t a I. INTRODUCTION scenario, the atomic masses of the different components m differ. In either of these realizations of multi-component - Fermi gases, all or some of the interspecies scattering d Overthe past decade orso, the field ofultracoldgases lengths may be tunable thanks to the possible existence n has seen tremendous breakthroughs. After reaching de- o of magnetic or optical Fano-Feshbach resonances. This generacyinBose[1,2]andFermi[3]gases,therealization c may open the possibility to experimentally investigate [ of an atom laser [4, 5, 6], of the Mott-insulator transi- the stability and to study, provided an extended stable tion[7],andoftheconversionfromanatomictoamolec- regime exists, the behaviors of multi-component Fermi 1 ular gas [8] followed. Many of the present-day studies v gases as a function of the s-wave scatterig length. take advantage of the tunability of the atom-atom scat- 1 Using two different theoretical frameworks,this paper 1 teringlengthinthe vicinity ofaso-calledFano-Feshbach considers three- and four-component Fermi gases, and 0 resonance [9, 10]. As an external magnetic field is tuned compares their behaviors with those of two-component 4 through its resonance value, the sign of the scattering Fermi gases. In particular, we ask how the stability . length changes [11, 12]. Exactly on resonance, the scat- 1 of two-component Fermi gases changes when a third or 0 tering length is infinitely large, allowing for the study of fourth component are added. It is now well established 8 strongly-correlated systems. Experimentally, the speed thattrappedtwo-componentFermigasesarestable even 0 ofthe magnetic field rampcanbe changed,allowingadi- whentheinterspeciess-wavescatteringlengthisnegative : abatic ramps, for example, and the ramp itself can be v and its magnitude is infinitely large [19, 20, 21, 22, 23]. Xi reversed. It is this versatility that made possible the Thestabilityofinhomogeneousaswellasofhomogeneous experimental study of the BCS-BEC crossover,using ul- two-component Fermi gases with attractive short-range r tracold atomic Fermi gases trapped in two different hy- a interactions and arbitrary interspecies scattering length perfine states [13, 14]. can be attributed to the Pauli exclusion principle (also Using present-day technology, the realization of de- referredto asFermipressure),whichintroduceseffective generate multi-component atomic Fermi gases appears repulsive intraspecies interactions that more than com- possible in principle. The occupation of more than two pensate the attractive interspecies interactions. In con- different hyperfine states of the same species requires, trast, homogeneous Bose gases with negative scattering neglecting for the moment possible losses, only moder- lengths are unstable; they can, however,be stabilizedby ate changes of current set-ups. A particularly promising an external confining potential as long as the product of candidate appears to be 6Li [15], and the coexistence of thenumberofbosonsandthes-wavescatteringlengthis threehyperfinestateshasalreadybeendemonstratedfor less negative than a certain critical value [24, 25, 26]. 40K [16]. Alternatively, a number of goups are presently Section II introduces the Hamiltonian of trapped pursuing the simultaneous trapping of three different multi-componentFermigasesaswellas asimple “count- atomicspecies[17,18]. Intheformerscenario,theatomic ingargument”thatturnsouttobequiteusefulinunder- massesofallcomponentsareequal,whereasinthe latter standing the stability of multi-component Fermi gases. 2 Section III investigates the stability of three- and four- component equal-massFermi gases within a hyperspher- TABLE I: Number Natt of attractive interactions, number ical framework, focussing on the large and small parti- Nrep of effectively repulsive interactions and ratio Nrep/Natt for finite and infinite N for a χ-component Fermi gas (χ=2 cle number limits. The physical picture emerging from through4) inwhichallinterspeciesinteractionsareequal(or the hyperspherical treatment is further investigated in resonant). Sec. IV, which solves the many-body Schrödinger equa- tionforshort-rangeinteractionsusingMonteCarlotech- χ=2 χ=3 χ=4 niques. Finally, Sec. V summarizes and connects our Natt 14N2 13N2 38N2 results with those available in the literature [20, 27, 28, Nrep N2 `N2 −1´ N2 `N3 −1´ N2 `N4 −1´ 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40]. Nrep/Natt (N finite) NN−2 N2−N3 N3−N4 Nrep/Natt (N →∞) 1 21 13 II. HAMILTONIAN AND GENERAL CONSIDERATIONS analysis. ThenumberN ofattractivepairinteractions att V , where again α=β, is given by The Hamiltonian H for an atomic Fermi gas with χ αβ 6 components under external spherically symmetric har- N2χ 1 monic confinement is given by N = − , (2) att 2 χ χ Nα ~2 1 andthenumberN ofeffectivelyrepulsiveinteractions, H = 2 + m ω2~r2 + rep −2m ∇~rαi 2 α α αi i.e., the number of like fermion pairs, by α=1i=1(cid:18) α (cid:19) XX χ Nα Nβ N = N N −χ. (3) + Vαβ(~rαi ~rβj ). (1) rep 2 χ | − | α<βi=1j=1 XXX Table I summarizes the values of N , N and att rep Here, the number of atoms of the αth component is de- N /N for χ = 2, 3 and 4. The ratio N /N (re- rep att rep att noted by N , and the total number of atoms is given by ported in the fourth and fifth row of Table I for finite α N, N = χ N . In Eq. (1), ~r denotes the position and infinite N, respectively) decreases with increasing α=1 α αi vector of the ith atom of the αth component, measured χ and approaches 1/(χ 1) in the large N limit, indi- P − withrespecttothecenterofthetrap,andm andω re- cating that the Fermi pressure becomes less important α α spectively the atomic mass and the angular frequency of compared to the interspecies interactions as χ increases. the αth component. The potential V describes the in- Anotherinterestingscenarioariseswheneachcomponent αβ teractionbetweenanatomoftheαthandanatomofthe is occupied by exactly one fermion, i.e., when χ = N. βthcomponent. Thisworkconsidersazero-rangepoten- In this case, no effectively repulsive interactions exist, tial(see towardsthe endofthis sectionandSec.III)and i.e., N /N =0, and the system’s ground state is the rep att apurelyattractiveshort-rangepotential(seeSec.IV). In same as that of the corresponding N-boson system. As both cases, we characterize the strengths of the V by pointed out already in the introduction, Bose gases with αβ the s-wave scattering lengths a . We assume that the negative scattering length become unstable in the limit αβ two-body interactions are independent of spin, implying that the absolute value of the scattering length becomes thatthenumberofatomsineachspinstateisconserved. large. This, together with the fact that two-component Throughout this work, like atoms are taken to be non- Fermigasesarestable for allscatteringlengths, suggests interacting, implying a =0 for α=1, ,χ. that there exists a critical χ-value beyond which multi- αα ··· In the most general case, the Hamiltonian given in component Fermi gases with large negative interspecies Eq. (1) has χ different m , ω and N (α = 1, ,χ), scattering length are unstable. Sections III and IV show α α α ··· andχ(χ 1)/2differentV (α,β =1, ,χandα=β), that χ =3. αβ cr − ··· 6 resulting in a tremendously large parameter space. To In addition to Fermi gases in which all interspecies in- reduce the parameter space, we first consider the case teractions V (α,β = 1, ,χ and α = β) are equal, αβ ··· 6 where all m , all ω , all N , and all a (α = β) are we consider the scenario in which only a subset of inter- α α α αβ 6 equal. A four-component gas of this type could, e.g., species interactions are “turned on”. In particular, we be realized by equally populating and trapping the four consider χ-component Fermi gases with χ 1 equal and − spin states of a fermionic atom whose ground state has non-zero (or resonant) a . For the three-component αβ vanishing total electronic angular momentum J but a system, we take a = a and a = 0. The number 12 23 31 non-vanishing nuclear spin I of 3/2. In this case, the of attractive interactions is in this case by a factor of scattering lengths between the different spin substates 2/3 smaller than in the case with three resonant inter- are equal and s-wave scattering between two atoms in actions, thus increasing the ratio of N /N from 1/2 rep att the same spin substate are forbidden by symmetry. to 3/4in the large N limit. For the four-componentsys- BeforesolvingtheSchrödingerequationfortheHamil- tem, two different “non-trivial” possibilities for turning tonian given in Eq. (1), we present a simple counting on only χ 1 interactions exist: (i) a = a = a 12 13 14 − 3 and a = a = a = 0, and (ii) a = a = a and which can be written as (k (0)a )3 1 for fermions, 23 34 24 12 23 34 F s | | ≪ a = a = a = 0 [the configuration a = a = a is violated and the instability prediction, Eq. (6), can- 13 14 24 12 23 34 and a = a = a = 0, e.g., is equivalent to (i); the not be trusted. Indeed, Eq. (6) disagrees with the 13 14 24 configuration a = a = a and a = a = a = 0, experimental finding that two-component Fermi gases 12 23 31 14 24 34 e.g., is trivial in the sense that it can be broken up into are stable [13, 14, 47, 48, 49, 50]. To resolve this a three-component system with all resonant interactions disagreement, two independent studies recently intro- and a single non-interacting component]. For the non- duced density-dependent “renormalization schemes” of trivial configurations, the number of attractive interac- the scattering length a entering into Eq. (4) [51, 52]. s tions is half as large for the case of three resonant inter- Using these “renormalization schemes”, the (modified) actions as in the case of six resonant interactions, thus density-dependentpseudopotentialcanbe appliedto de- increasing the ratio of N /N from 1/3 to 2/3 in the scribe strongly-interacting Fermi gases and predicts, in rep att large N limit. This counting analysis indicates that the agreement with experimental [13, 14, 47, 48, 49, 50] and values of the ratio N /N for large N for three- and other theoretical [19, 20, 21, 22, 53] works, that two- rep att four-component gases with χ 1 resonant a are be- component Fermi gases are stable even in the strongly- αβ − tween those for two- and three-component Fermi gases interactingunitaryregime. ThefactthattheFermipseu- with χ(χ 1)/2 resonant a . Thus, the question arises dopotentialhas to be modified when applied to fermions αβ − whether multi-component Fermi gases with χ 1 reso- was already suggested earlier (see, e.g., Refs. [44, 54]), − nantinteractionsarestablefor allscatteringlengths and andiswellknowninthenuclearphysicscommunity(see, N values, or whether they become unstable for a certain e.g., Refs. [55, 56] and references therein). critical parameter combination. The instability prediction given in Eq. (6) for two- To analyze the stability of Fermi gases quantitatively, component Fermi gases can be readily generalized to one may attempt to describe the interspecies atom- multi-component Fermi gases with χ(χ 1)/2 and with − atominteractionsby a zero-rangeFermipseudopotential χ 1resonantinteractions(in the caseofχ 1resonant − − V (~r) [41], interactions, we consider the non-trivial scenario intro- δ duced above, in which none of the components can be 2π~2as separated off). The right hand side of Eq. (6) has to V (~r)= δ(~r), (4) δ µ be multiplied by 1/(χ 1) in the former case, and by − χ/[2(χ 1)] in the latter case. Since the diluteness cri- which is directly proportional to the s-wave scattering − terium is fullfilled approximatelyin the all resonantcase length as. Here, µ denotes the reduced mass of the two at the predicted collapse point [e.g., (k (0)a )3 0.23 F s interacting atoms. This pseudopotential has been em- ≈− and 0.067for χ=3and4,respectively],the prediction ployed successfully to predict many properties of dilute − that multi-component Fermi gases with all resonant in- Bose gases. For example, the interaction potential given teractionsbecomeunstableforafinite a ,derivedusing s in Eq. (4) together with a Hartree wave function cor- | | the “bare” interaction given in Eq. (4), may in fact be rectly predicts that trapped Bose gases with negative a s correct. Sections III and IV confirm this. In the case become unstable if [24, 42] with χ 1 resonant interactions, the bare Fermi pseu- − a dopotential, Eq. (4), predicts that the collapse occurs at (N 1) s . 0.575, (5) (k (0)a )3 0.76 and 0.54 for χ = 3 and 4, respec- − a − F s ho ≈ − − tively. This suggests that the bare pseudopotential can- where aho denotes the oscillator length, aho = not be used and that the instability prediction derived ~/(2µω). In general, the true atom-atom interaction usingitmaynotbe correct(see Sec.IVforamany-body canbereplacedinthelong-wavelengthlimitbythepseu- analysis). p dopotential given in Eq. (4) provided the system is dilute, i.e., if n(0)a 3 1, where n(0) denotes the peak s | | ≪ density. Assuming N is not too small, one finds that III. HYPERSPHERICAL FRAMEWORK n(0)a 3 ismuchsmallerthanonewhen(N 1)a equals s s | | − 0.575a ; consequently, the pseudopotential predicts − ho A. N-particle system the instability of dilute Bose gases correctly. Appliedtotwo-componentFermigases,the barepseu- This section investigates the stability of three- and dopotentialemployedwithinahypersphericalframework four-component Fermi gases within a hyperspherical predicts that the system becomes unstable if [43] framework [37, 43, 57, 58]. Throughout this section, we k (0)a . 1.22, (6) assume that all angular trapping frequencies are equal, F s − i.e.,ω =ωforα=1, ,χ. Togaininsightintothesys- α ··· where k (0) denotes the noninteracting Fermi wave vec- tem’s behavior, we employ hyperspherical coordinates: F tor at the trap center. (A mean-field analysis predicts the3N coordinatesaredividedintoahyperradiusRand a slightly more negative critical value of k (0)a . 3N 1 hyperangles, collectively denoted by Ω [59]. In F s − π/2 [44, 45, 46].) However, at the critical parameter the following, only the definition of the hyperradius R, − combination k (0)a = 1.22, the diluteness criterium, which can be thought of as a measure of the size of the F s − 4 system, is needed, strongly-interacting limit [37, 64]. For infinitely large interspecies scattering lengths a (α = β), the adia- 1 αβ 6 1 χ Nα 2 batic approximation introduced above is— for a class of R= m ~r2 . (7) states that arise assuming boundary conditions consis- M α αi α=1i=1 ! tentwithcontactinteractionswhenthedistancebetween XX eachpairofparticlesgoestozero[64],referredto asuni- As before,the positionvectors~r aremeasuredwithre- αi versal states in the following— exact and the universal spect to the center of the trap. In Eq. (7), M denotes effective potentials V (R) are given by the total mass of the system, i.e., M = χ Nα m . ν α=1 i=1 α Using these coordinates, the many-body wave function ~2(1+C ) ν P P V (R)= . (11) Ψ(R,Ω)canbeexpandedintermsofasetofΩ-dependent ν 2MR2 channel functions Φ (Ω;R), which depend parametri- ν cally on R, and R-dependent weight functions Fνn(R), InEq.(11),theCν denoteR-independentconstantsthat arise from the integration over the 3N 1 hyperangles. − Ψ(R,Ω)= R(1−3N)/2F (R)Φ (Ω;R). (8) The functional form given in Eq. (11) has also been de- νn ν rived by explicitly solving the many-body Schrödinger νn X equation for a class of variational wave functions using Inthe adiabatichypersphericalapproximation[57,60, hypersphericalcoordinates[37]. Usingdimensionalargu- 61,62,63],whichneglectsoff-diagonalcouplingelements ments, Eq. (11) can be understood as follows: The term and additionally restricts the sum in Eq. (8) to a single of V (R) that accounts for the particle-particle interac- ν term, the solution of the many-body Schrödinger equa- tions has— because of the absence of any other length tionreducestodetermininganR-dependenteffectivepo- scale in the problem— to scale in the same way with R tentialVν(R),whichincludespartofthekineticenergyas as the 1/R2 term that accounts for the kinetic energy wellasthe effects ofthe two-bodyinteractions,andthen contribution. The total universal effective hyperradial solving an effective one-dimensional radial Schrödinger potential curves at unitarity thus have a simple func- equation in the hyperradius R, tionalform: V (R)dominatesatsmallRandapproaches ν plus or minus infinity in the R 0 limit, depending on ~2 ∂2 → −2M ∂R2 +Vν(R)+Vtrap(R) Fνn(R)= Vwheth(Rer)tdhoemqiunaantetistyat1l+argCeνRis.positive or negative, and (cid:20) (cid:21) trap E F (R), (9) For a two-component unitary Fermi gas, a number of νn νn C valuesare known. ForN 6,selectedC with ν 2 ν ν where have been obtained using a≤correlated Gaussian (C≤G) basis set expansion-type approach [65]. For N 30, 1 Vtrap(R)= 2Mω2R2. (10) upper bounds for the C0 have been obtained by th≤e FN- DMC method [65]. In general, the calculation of the effective potentials In the large N limit, the value of C for equal-mass 0 V (R) is, at least for many-body systems, just as hard two-component Fermi gases at unitarity can be related ν assolvingthemany-bodySchrödingerequationitself. In to the universal parameter β of the homogeneous sys- certain circumstances, however, the effective potentials temusingthe hypersphericalframeworkofRef.[37]that or their functional form can be obtained more easily. employsadensity-dependentscatteringlength[52],lead- For weakly-interacting dilute equal-mass two- ing to C = β [37]. Alternatively, this relationship can 0 component Fermi gases, e.g., the effective hyperradial be derived by applying the local density approximation potential V (R) consists of two parts [37, 43]: The first to the trapped system (see, e.g., Ref. [65]). It is gener- 0 part, which is proportional to c /R2, arises from the allybelievedthatthe mostaccuratevalueforβ hasbeen kin kinetic energy operator, and the second part, which obtained by the FN-DMC method [22, 66], β = 0.58. − is proportional to c /R3, accounts for the particle- Since1+β >0,theuniversalhyperradialpotentialcurve int particle interactions. c is positive, and c is directly V (R)+V (R) for large N (shown by a solid line in kin int 0 trap proportional to the s-wave scattering length a [43]. If Fig. 1 using β = 0.58)is repulsive at small R, prevent- s − a is sufficiently small (a < 0), the small R region, ing the collapse of the two-componentunitary Fermi gas s s | | where the attractive c /R3 term dominates over the towards cluster-like bound states. The hyperspherical int repulsive c /R2 term, is separated by a “barrier” potential curve picture reveals an intuitive way of un- kin from the R region where the repulsive c /R2 term derstandingthestabilityoftheenergeticallylowest-lying kin dominatesoverthe attractivec /R3 term. Thisbarrier gas-like state of two-component unitary Fermi gases. int prevents the Fermi gas from collapsing to a cluster-like In Fig. 1, the hyperradial potential V(R) and the hy- many-body bound state, thus explaining the stability perradius R are scaled by E and R , respectively, NI NI of weakly-interacting dilute equal-mass two-component where E denotes the energy of the non-interacting NI Fermi gases (see also Ref. [51]). χ-component Fermi system and R the square root NI Analytical expressions for V (R) for equal-mass two- of the expectation value of R2 calculated for the non- ν component Fermi gases have also been obtained in the interacting χ-component Fermi system, i.e., R = NI 5 4, respectively [70]. The resulting hyperradial potential curves V (R)+V (R) are shown in Fig. 1 in the large 0 trap 4 N limit by dotted and dashed lines for χ = 3 and 4, respectively. The attractive small-R behavior is due to thenegativevaluesof(1+C ),andsuggeststhatunitary 0 NI 2 three-andfour-componentFermigaseswithallresonant E / interationsareunstable. NotethatthecoefficientC0 and ) N /N both scale as 1 χ with χ, reflecting the fact NI 0 att rep − R thatV0(R)isinpartdeterminedbythepairinteractions. R/ Interestingly, both the bare pseudo-potential (see the V( -2 discussion towards the end of Sec. II) and the density- dependent pseudo-potential (used in this section within thehypersphericalframework)predictthatχ-component -40 1 2 3 Fermigases,χ>2,withallresonantinteractionsbecome R/R unstable for a finite value of the interspecies scattering NI length. The critical value predicted by the bare interaction is presumably not negative enough while that for FIG. 1: (Color online) Hyperradial potential curve V (R)+ 0 the density-dependent interaction would presumably be Vtrap(R) as a function of the hyperradius R for χ = 2 (solid toonegative(weexpectthattherenormalizationscheme line),χ=3(dottedline)andχ=4(dashedline)inthelarge N limit. Allinteractionsbetweenunlikeatomsarecharacter- originally developed for the two-component Fermi gas ized byan infinitescattering length,and thecoefficient C is “over-renormalizes”the scattering length). 0 takentobe(χ−1)βwithβ =−0.58. Bothlengthandenergy The analysis above can be extended to the case where are scaled bythecorresponding valuesof thenon-interacting only χ 1 interspecies scattering lengths are infinite; as system (see text). before,−we focuss on what we defined in Sec. II as “non- trivial” scenarios. In the case of χ 1 resonant inter- − actions, the bare pseudopotential is expected to fail (see R2 [37,43]. Usingthevirialtheorem[64,67],R NI NI Sec. II). Using, just as above, density-dependent zero- h i can be related to E , NI range two-body interactions to describe the V with p αβ non-zero a , the coefficient C becomes in the large N ~ ENI limit 2(χ α1β)β;justasinthea0ll-resonantcase,C scales RNI = Mω ~ω . (12) χ − 0 r r with χ in the same way as Natt/Nrep. Using β = 0.58, − (1+C )ispositiveforχ=3and4(thecorrespondingpo- In Fig. 1, E and R have been evaluated by writing 0 NI NI E as(λ+3N/2)~ωandusingtheleadingorderofλfor tentialcurvesareshowninFig.2); the repulsivesmall-R NI large closed-shell systems, λ=(3N)4/3/4. We note that behavior suggests that three- and four-component unitary Fermi gases with χ 1 resonant interactions are thescaledquantityR/R remainsfiniteinthelimitthat NI − stable. We notethattheanalysisinhypersphericalcoor- N goes to infinity, allowing the thermodynamic limit to dinates, which employs density-dependent interactions, be taken. has two short-comings: i) The derivation of the effec- For all-resonant three- and four-component gases at tive potential curves assumes the existence of a single unitarity, the functional form given in Eq. (11) remains length scale, the hyperradius R; since some of the inter- valid for universal states [37, 64]. For these systems, species scattering length are zero while others are res- the values of C are not known (the only exception be- ν onant, a more accurate description would presumably ing the N = 3 system [68, 69], see below). It has been introduce an additional length scale and be based on argued [37], parametrizing V by a density-dependent αβ a “hyperradial vector” as opposed to a single hyperra- zero-rangetwo-bodypotential[52]andneglectingtheex- dius. ii)The derivationneglectsthe possibleexistenceof istence of weakly-bound trimer or tetramer states, that weakly-bound trimer (see the next section) or tetramer the coefficient C is in the large N limit determined by 0 states, or said differently, the existence of non-universal theparameterβoftheunitarytwo-componentFermigas, states. Section IVB shows that the behaviors of multi- i.e.,C =(χ 1)β. Duetothe lackofbenchmarkresults 0 − componentFermigaseswithχ>2 maydiffer depending for all-resonant three- and four-component Fermi gases, on whether or not such states exist. there is some ambiguity in how the simplest adaptation ofthe renormalizationscheme,originallydesignedto describe the physics of two-component Fermi gases [52], is applied to three- and four-component Fermi gases, B. Three-particle system thusintroducingsomeuncertaintyabouttheexactvalues of C that describe unitary three- and four-component This section discusses the behaviors of three distin- 0 Fermigases. UsingC =(χ 1)β andβ = 0.58[22,66], guishablefermionsinteractingthroughzero-rangepoten- 0 − − (1+C )isnegativeforthree-andfour-componentgases, tials with equal masses and equal trapping frequencies. 0 (1+C ) = 0.16 and (1+C ) = 0.74 for χ = 3 and For these systems, the effective hyperradial potential 0 0 − − 6 radius R is defined without the center-of-mass motion andscaledbythehyperradiusR ofthenon-interacting NI 4 system without the center-of-mass motion. For the discussion that follows, the difference between the hyperradius defined without the center-of-mass and that de- NI 2 fined in Eq. (7) is irrelevant. For the scattering lengths E shown, the hyperradial potentials show a barrier around / ) R 0.2R (i.e., R 0.3a ), whichseparatesthe large NI 0 ≈ NI ≈ ho R R region where gas-like states live from the small R re- R/ gion. As as becomes more negative, the barrier height V( -2 decreases and moves to slightly larger R values. The critical scattering length a at which the barrier disap- c pears (a 0.39a ) provides a first estimate for when c ho -4 the three-≈bo−dy system becomes “unstable” against the 0 1 2 3 R/R formation of tightly-bound trimer states with negative NI energy,assumingsuchstatesexist;theconditionsforthe existence of negative energy states are discussed below. FIG. 2: (Color online) Hyperradial potential curve V (R)+ 0 A more accurate estimate would account for the zero- Vtrap(R)asafunctionofthehyperradiusRforχ=3(dotted pointenergyofthequantumsystem,resultinginasome- line)andχ=4(dashedline)inthelargeN limitforχ−1res- onant interactions (the scattering lengths a , β = 2,··· ,χ, whatlessnegativecriticalscatteringlengththatisinfair 1β are infinite and all other scattering lengths are zero). The agreement with the GP prediction given in Eq. (5) (see coefficientC istakentobe 2(χ−1)β withβ=−0.58. Both Ref. [57] for a discussion of the N-boson system). The 0 χ length and energy are scaled by the corresponding values of lifetime ofthegas-likethree-bodystatecanbeestimated the non-interacting system (see text). The plotting range of by calculating the tunneling rate through the potential theaxis is thesame as in Fig. 1 for ease of comparison. barrierasafunctionofa [57]. Whenthebarrierheightis s large compared to the energy of the non-interacting system (small a ), the lifetime of the gas-like trimer state s | | is much larger than 1/ω; in this case, the system can be considered stable (it is, in fact, metastable). The tunneling rate is appreciable only when the barrier height 4 becomescomparabletotheenergyofthenon-interacting system, i.e., when a approaches a . NI s c E / 2 Toobtaintheenergyspectrumintheadiabaticapprox- ) NI imation,wesolvetheone-dimensionalradialSchrödinger R equation [Eq. (9)]. To this end, we add a term V (R) R/ 0 SR thatisrepulsiveforR.R andnegligiblysmallforR& ( c V R to V (R), V (R)=2(1 √3tanh2(R/R ))/(2mR2). c ν SR c − -2 VSR(R) cures the divergence of the Vν(R) at small R, which is related to the Thomas collapse [73]. When the 0 0.5 1 1.5 2 2.5 barrierheightislargecomparedtotheenergyofthenon- R/R interacting system and the short-range length scale R NI c is chosen sufficiently small, this “ad-hoc” regularization FIG. 3: (Color online) Solid, dash-dotted and dashed lines should affect the eigenenergies of the gas-like states at show the hyperradial potential curve V0(R)+Vtrap(R) as a most weakly. The eigenenergies of molecular-like states function of the hyperradius R, scaled by RNI, for three dis- (states livinginthe smallR region),however,shouldde- tinguishable fermions interacting through zero-range interac- pend comparatively strongly on the particular value of tionswithas=−0.1aho,−0.15aho and−0.2aho,respectively. Rc. Thus, some properties of the system are expected Note that R is defined without the center-of-mass motion, to be non-universal. Indeed, our calculations for differ- implying ENI =3~ω. ent R show that the first state with negative energy c appears when the scattering length a is approximately s equal to 8.4R . For R = 0.001a , e.g., the first neg- c c ho curves Vν(R) can be obtained analytically for any com- ative ener−gy state appears for as 0.0084aho, a value bination of scattering lengths a12, a23 and a31 [71, 72]. muchlessnegativethanac. Atthis≈sc−atteringlength,the Figure 3 shows V (R)+V (R) for the three-particle barrierheightislargecomparedtotheenergyofthenon- 0 trap system as a function of R for three different negative interacting system and the gas-like system is— because scattering lengths a , i.e., a = 0.1a , 0.15a and of the existence of a negative trimer state— metastable. s s ho ho − − 0.2a (here, a = a with α = β) as a function of IfR ischosensothattheabsolutevalueofthescattering ho s αβ c − 6 the hyperradius R. Throughout this section, the hyper- length at which the first bound state appears is compa- 7 rable to orlargerthan a , then the gas-likestate would IV. MONTE CARLO TREATMENT c | | bethetruegroundstateofthesystemfor a < a . For s c | | | | realisticdilute alkaligaseswithall-resonantinteractions, This section discusses the solutions of the many-body we expect that three-body bound states appear for scat- SchrödingerequationfortheHamiltoniangiveninEq.(1) teringlengthslessnegativethanac; consequently,trimer obtained by Monte Carlo techniques as a function of the as well as larger three-component systems with as & ac scatteringlength. SectionIVAintroducesthe VMC and are expected to be metastable. The picture developed FN-DMC methods, and Sec. IVB presents our results. herewithintheadiabaticapproximationremainsqualita- tively valid if the coupling between channels is included. A. Variational and fixed-node diffusion Monte In addition to the all-resonant three-particle system, Carlo method we consider the trapped three-particle system with two resonant interactions (a = a = a and a = 0). In s 12 23 31 this case,the hyperradialpotential barrierdisappears at Usingthe VMC andFN-DMC methods[77], wedeter- mine solutions of the many-body Schrödinger equation amorenegativescatteringlengththanintheall-resonant case,i.e., ata 1.00a . Ananalysisofthe tunneling for two different purely attractive short-rangemodel po- c ho ≈− tentials V : a square well interaction potential VSW, probabilitysuggeststhatthelifetime ofthetwo-resonant αβ αβ system is, considering the same a , enhanced by a factor s ofabouttencomparedtotheall-resonantsystem. Alife- V for r <R VSW(r)= − αβ,0 αβ,0 (13) time enhancement is expected since the ratio Natt/Nrep αβ (0 for r >Rαβ,0, (see Sec. II) is smaller for the three-component system with two resonant interactions than for that with three and a Gaussian interaction potential VG, αβ resonant interactions. Calculating the trimer energies in the adiabatic ap- VG(r)= V exp r2 . (14) proximation,we find thatthe firstthree-body state with αβ − αβ,0 −2Rα2β,0! negative energy appears at a 2500R . For exam- s c ≈ − ple, a negative energy state appears for as ≈ −2.5aho if Both potentials depend on a depth Vαβ,0, Vαβ,0 ≥ 0, Rc =0.001aho and for as 0.25aho if Rc =0.0001aho. and a range Rαβ,0. In our calculations, we set all Rαβ,0 ≈− In the latter case, the disappearing of the hyperradial to the same value, Rαβ,0 aho, and adjust the depths ≪ potential barrier is accompanied by a “collapse” of the Vαβ,0 until the s-wavescattering lengths aαβ assume the metastable gas-like state to a cluster-like state. In the desired values. In Sec. IVB, Vαβ,0 and Rαβ,0 are chosen formercase,incontrast,thedisappearingofthehyperra- so that the potential Vαβ supports no two-body s-wave dial potential barrier is accompaniedby the lowest-lying boundstateandsothataαβ 0. Theuseoftwodifferent ≤ gas-like state evolving smoothly to a cluster-like state functionalformsfortheinteractionpotentialsVαβ allows with negative energy. Although the short-range param- forexploringthe dependence oftheresultsonthedetails eter R cannot be straightforwardly related to the pa- of the short-range physics. c rameters of typical atom-atom potentials, it seems plau- We now discuss how the solutions of the many-body sible that realistic three-component Fermi systems with Schrödinger equation are obtained by the variational tworesonantinteractionscouldfallintoeitherofthetwo Monte Carlo method. The functional form of the vari- categoriesdiscussedhere. SectionIVinvestigatesthebe- ational wave function ψT is chosen on physical grounds haviors of three-component many-body systems and at- (see below), and ψT is parametrized in terms of a set of tempts to determine how the system’s behaviors depend parametersp~, which are optimized so as to minimize the onthe range ofthe underlying finite-rangetwo-body po- energy expectation value E(p~) [77], tential (this range can, roughly speaking, be connected ψ (p~)H ψ (p~) with the short-range parameter Rc). E(p~)= h T | | T i. (15) ψ (p~)ψ (p~) T T We note that the negative scattering length a at h | i s which the first weakly-bound trimer state appears (as- Sinceψ dependsingeneralonthe3N coordinatesofthe T suming no deeply-bound states exist), can be estimated system,theintegrationontherighthandsideofEq.(15) using the following “rule of thumb” [74, 75, 76]: a = is high-dimensional; this implies that standard tech- s | | R exp(π/s ), where s is the coefficient that determines niquessuchastheSimpsonrule[78],whichscale,roughly c 0 0 the energy spectrum of the lowest hyperradial potential speaking, exponentially with the number of degrees of curve at unitarity (s = 1.006 for a = a = a = a , freedom cannot be used for its evaluation. Instead, we 0 s 12 23 31 and s = 0.413 for a = a = a and a = 0 [68]). evaluate the high-dimensional integral by Monte Carlo 0 s 12 23 31 This rule of thumb predicts a negative scattering length techniquesusingstandardMetropolissampling[79]. The at which the first weakly-bound trimer state appears of stochastic evaluation of the integral introduces a statis- 23R and 2000R for three andtwo resonantinterac- ticaluncertainty,whichcanbe reducedbyincreasingthe c c − − tions,respectively. Thesevaluesaretobecomparedwith number of samples K included in the evaluation of the 8.4R and 2500R found numerically (see above). expectation value (the errorbar decreases as 1/√K with c c − − 8 increasing K). We denote the energy expectation value particle harmonic oscillator wave function. We choose obtained by the VMC method by E . to write ψ in terms of the product φ ϕ instead of the VMC T α i Thevariationalwavefunctionψ formulti-component harmonicoscillatorfunctionsthemselves,becausethisal- T Fermi gases is written as (see Refs. [21, 22, 53, 65, 80, lowsthe widths b ofthe Gaussianstobe variedwithout α 81,82]forMonteCarlostudiesoftwo-componentsystems changing the nodal surface of ψ . T andRef.[35]foraMonteCarlostudyofthree-component In addition to the VMC method, we apply the FN- systems) DMC method [77, 83]. The FN-DMC method uses the variationalwavefunction ψ asa so-calledguiding func- T χ Nα Nβ χ Nα tion and determines the energy of a state whose nodal ψT = fαβ(|~rαi−~rβj|)× ϕα(~rαi) surface is identical to that of ψT. It can be shown that αY<βiY=1jY=1 αY=1iY=1 the FN-DMC energy provides an upper bound to the χ lowest eigenenergy of the eigenstate that has the same × A(φ1(~rα1),··· ,φNα(~rαNα)),(16) symmetry as ψT [83]. If the nodal surface of ψT coin- α=1 cides with that of the true eigenfunction, then the FN- Y DMC method results, within the statistical uncertainty, where denotes the anti-symmetrizer and f , ϕ and αβ α A in the exact eigenenergy. In general, the quality of the φ denote two- or one-body functions whose functional i FN-DMC energy depends crucially on the construction form is discussed in detail in the following. of the nodal surface of ψ . As in the VMC method, The f denote two-body correlation factors. For T αβ expectation values calculated by the FN-DMC method small r, f coincides with the zero-energy scatter- αβ have a statistical uncertainty, which can be reduced by ing wave function for two particles interacting through increasing the computational efforts. V . Beyond some matching value R , which is αβ αβ,m treated as a variational parameter, f is taken to be αβ c +c exp( c r),wherec denotesavariational αβ,1 αβ,2 αβ αβ parameter and c− and c are determined by the B. Monte Carlo results αβ,1 αβ,2 condition that f and its derivative are continuous at αβ r =Rαβ,m. The Rαβ,m andcαβ are optimized under the We first consider three- and four-component Fermi constraint that f 0 for all r. gases with one atom per component and equal masses αβ ≥ The ϕα denote Gaussian orbitals with widths bα, and trapping frequencies. If all interspecies scattering ϕ (~r)=exp( r2/(2b2)). The parameters b controlthe lengths areequal,these Fermigasesaredescribedby the α − α α size of the Fermi cloud, and are optimized variationally. sameHamiltonianasthecorrespondingBosegasandthe In the non-interacting case, b = a . For a < 0, systemshouldbecomeunstabletowardstheformationof α ho αβ the system becomes more compact due to the attractive negative energy states, characterized by small interpar- interactions, resulting in a smaller optimal value of b . ticle distances, when the inequality given in Eq. (5) is α If all mα and all bα are the same, i.e., mα = m and fulfilled. ForN =3and4,this impliescriticalscattering b = b, then the product χ Nα ϕ (~r ) reduces to lengths of a 0.29a and 0.19a , respectively. α α=1 i=1 α αi s ≈− ho − ho exp( NR2/(2b2))andthevariationalparameterbispro- To investigate how this instability or collapse arises − Q Q portionaltothemeanhyperradiusofthevariationalwave within the many-body Monte Carlo framework, we per- function. form DMC calculations for the three- and four-particle The φ (~r) are given by H (x)H (y)H (z), where systems interacting through a square well potential with i nx ny nz (n ,n ,n ) = (0,0,0) for i = 1, (n ,n ,n ) = (1,0,0) R =0.01a (thesubscriptsαandβ ofR havebeen x y z x y z 0 ho αβ,0 for i=2, (n ,n ,n )=(0,1,0) for i=3 and so on, and dropped for notational convenience). The depth is cho- x y z the H denote Hermite polynomials of degree n. sen so that a takes the desired value and so that the n s The anti-symmetry and nodal surface of ψ are de- potential supports no two-body bound states. For small T termined by the χ Slater determinants [second line of a ,theDMC methodresultsinenergiesthatagreevery s | | Eq. (16)], which contain the one-body functions φ . The well with the mean-field Gross-Pitaevskii (GP) energy i nodal surface of ψ coincides with that of the non- for both N = 3 and 4. As a becomes more negative T s interacting multi-component Fermi gas. For small a , (ofthe orderof 0.06a [84]), the DMC walkerssample αβ ho | | − a < 0, this is expected to be a very good approxima- atfirstexclusivelypositive energyconfigurations;after a αβ tion. It has been shownthat the nodalsurfaceused here certain propagation time some walkers sample configu- also provides reasonably tight bounds for the energies of rations with negative energy that correspond to tightly- trapped two-component gases [65, 81, 82] all the way to boundstates. Theexistenceoftightly-boundthree-body unitarity. Whether this holds true for multi-component states,whichdependonthedetailsoftheunderlyingtwo- Fermi gases with χ>2 should be addressed in more de- bodypotential,forscatteringlengthsa thatarelessneg- s tailinfollow-upwork. Notethattheguidingfunctionψ ative than the critical scattering lengths determined at T used throughout this work does not build in any “pair- the mean-field level or within the hyperspherical frame- ing physics” from the outset (see Sec. V for further dis- work was already pointed out in Sec. IIIB. Our calcula- cussion). The product φ (~r)ϕ (~r), i = 1, ,N , with tions here show that the DMC algorithm “sees” these α i α ··· b =a coincides with the ith (non-normalized) single- three-body bound states (as well as four-body bound α ho 9 gases [86] (the second part of the Appendix in Ref. [85] provides a detailed discussion). As in those earlier studies, we interpret the existence of the energy barrier as an indication that the Hilbert space is divided into two 4 nearly orthogonal spaces if a is quite a bit less nega- s tive than the critical scattering length predicted at the NI 2 mean-field level: Gas-like states live in one region of the E Hilbert space,and cluster-likebound states in the other. / C M While this reasoning leads to an intuitive understand- V 0 ing of the stability and decay of Bose gases with neg- E ative scattering length, the question arises why the en- ergybarrierdisappearsfor scatteringlengths a that are -2 s notably less negative than the critical scattering length predicted by the GP equation. We believe that this is 0 0.5 1 1.5 2 2.5 b/a due to the existence of tightly-bound states with nega- ho tive energy. If a is sufficiently large,the overlapof the s | | variationalwavefunctionwitheigenfunctions ofcluster- FIG. 4: (Color online) Circles, squares and triangles show like bound states increases with decreasing b. For non- the variational energy EVMC for N = 3 atoms interacting vanishing overlap, E provides a rather poor upper through a square well potential (range R = 0.01a ) with VMC 0 ho as = −0.05aho, −0.1aho and −0.15aho, respectively, as a boundtothe truegroundstateenergyofthe systemand function of the variational parameter b [see the discussion notanupperboundtotheenergeticallylowest-lyinggas- followingEq.(16)]. EVMC isscaledbytheenergyENI ofthe like state. This implies that our VMC calculations do non-interacting system, ENI = 4.5~ω. Dotted lines connect notallowforareliabledeterminationofthe criticalscat- data points for ease of viewing. tering length. Instead, they indicate that the existence of the energy barrier gives rise to the stability of the gas-like state and that the disappearing of the energy states) if the simulation time is sufficiently long and a barrier qualitatively explains how the instability of Bose s is sufficiently large. If a is not too large, the DMC|al-| gases and multi-component Fermi gases with one atom s gorithm does— if the i|nit|ial walker configurations corre- per component arises when as becomes too negative. spond to gas-likestates— not “know”about the tightly- We note that the picture developed here based on our bound states with negative energy. VMCandDMCcalculationsiscloselyrelatedtotheanal- To gain further insight, we treat the N =3 system by ysisbasedonthehypersphericalformulationpresentedin the VMC method. For a given a , we fix the variational Sec. III and in Refs. [37, 57, 87, 88]. In the VMC cal- s parameters R and c (their exact values are not impor- culations, the role of the hyperradius R is played by the m tant for the discussion that follows) and vary b (we drop Gaussianwidthb[seediscussionfollowingEq.(16)]. Fur- the subscripts α and β of R and c , and the sub- thermore, variational mean-field calculations have been αβ,m αβ script α of b ). Circles, squares and triangles in Fig. 4 interpreted in much the same way [89, 90]. α show the variational energy EVMC for N =3 as a func- We now investigate the behaviors of multi-component tion of b for as = 0.05aho, 0.1aho and 0.15aho, re- Fermigaseswithmorethanoneatompercomponentfor − − − Gspaeucstsiviaenly.wiFdotrhabsis=la−rg0e.r05tahhaon, EabVoMutCaihnocr(etahseessyifsttehme warheicrhesaolnlainntt;erasspebceifeosres,cawtteesreintgalαeβng=thassa.αWβ ewfiitrhstαso6=lvβe expands compared to its optimal size, thereby reducing the many-body Schrödinger equation for a square well the attractionfelt betweenthe atoms)and if b is smaller potential with range R = 0.01a by the FN-DMC αβ,0 ho than about aho (the system shrinks compared to its op- method. Circles in Fig. 5 show the FN-DMC energy for timal size, thereby increasing the kinetic energy). For a four-component Fermi gas with four atoms per com- as = 0.1aho and 0.15aho (squares and triangles in ponent (this corresponds to a closed shell) as a func- − − Fig. 4), the variational energy assumes a local minimum tion of the s-wave scattering length a . The energy de- s at b ≈ aho, shows a local maximum at b ≈ 0.25aho for creases approximately linearly with decreasing as, sug- as = 0.1aho and at b 0.5aho for as = 0.15aho, and gesting that this weakly-interacting system can be de- − ≈ − becomes negative for yet smaller b. We refer to the local scribedto a goodapproximationperturbatively. Assum- maximumasan“energybarrier”thatseparatesthelocal ingzero-rangeinteractionswithscatteringlengthsa and s minimum at b aho from a global minimum that exists treating the system within first order perturbation the- ≈ at smaller b values. Finally, for more negative scattering ory, we find E 36~ω+ca /a , where c=93~ω/√2π, s ho lengths (notshownin Fig.4)no energybarrierexists for for N =16. Th≈is perturbatively calculated energy (solid the variational wave functions considered. line in Fig. 5) describes the FN-DMC energies (circles) The energy barrier discussed here for small s-wave very well. By additionally treating the N = 17 fermion interacting Bose gases also exists in three-dimensional system with N = 5 and N = N = N = 4 (squares 1 2 3 4 dipolar Bose gases [85] and one-dimensional Bose in Fig. 5), we find that the chemical potential decreases, 10 40 4 ν h 38 NI E 0 /C /C M M D 36 V E E -4 34 -0.04 -0.02 0 -8 0 0.5 1 1.5 2 2.5 a / a b/a s ho ho FIG.5: (Coloronline)CirclesandsquaresshowtheFN-DMC FIG. 6: (Color online) Circles and squares show the varia- energy EDMC for a four-component Fermi gas with N = 16 tional energy EVMC for a three-component Fermi gas with and 17 atoms, respectively, interacting througha squarewell N = 12 atoms interacting through a square well potential potential with range R0 =0.01aho as a function of the inter- (range R0 = 0.01aho) with as = −0.05aho and −0.1aho (all speciesscatteringlengthas(allinterspeciesscatteringlengths interspecies scattering lengths are equal), respectively, as a are equal). While the N = 16 system has four atoms per function of the variational parameter b. EVMC is scaled by component, the N = 17 system has one component with the energy ENI of the non-interacting system, ENI = 27~ω. five atoms and three components with four atoms. Dotted Dotted lines connect data points for ease of viewing. lines show a linear fit to theFN-DMC energies: The slope is 37.7~ω/a forN =16and40.7~ω/a forN =17. Forcom- ho ho parison,thesolidlineshowstheenergyforN =16calculated within first order perturbation theory, E = 36~ω+cas/aho with c≈37.1~ω. 4 just as the energy, approximately linearly with increasing as . Forscatteringlengthsmorenegativethanabout NI | | E 0 0.045a [84],theDMCwalkerssample—justasinthe − ho /C case of the small Bose systems— negative energy config- M urations. We find that the three-component Fermi gas V E -4 with N = 12 behaves similarly. At the FN-DMC level, the instability ofmulti-componentFermigasesis accompanied by the existence of cluster-like states with negative energy. To investigate whether these cluster-like -8 0 0.5 1 1.5 2 2.5 states “live” for sufficiently small as , as in the case of b/a Bose gases, in a Hilbert space that| is|approximately or- ho thogonal to the Hilbert space where the gas-like states FIG. 7: (Color online) Circles and squares show the vari- live, we perform a series of VMC calculations. ational energy EVMC for a four-component Fermi gas with Circles and squares in Fig. 6 show the variational en- N = 16 atoms interacting through a square well potential ergy EVMC for the three-component Fermi system with (range R0 = 0.01aho) with as =−0.05aho and −0.07aho, re- N =12asafunctionofbforas = 0.05ahoand 0.1aho, spectively,asafunctionofthevariationalparameterb. EVMC − − respectively (as before, b = bα and as = aαβ where is scaled by the energy ENI of the non-interacting system, α=β). The lowestvariationalenergyforas = 0.05aho ENI = 36~ω. Dotted lines connect data points for ease of is6obtained for b a , E = 26.09(4)~−ω. The viewing. ho VMC ≈ FN-DMC energy for this interspecies scattering length, E = 26.08(4)~ω, agrees with E within error DMC VMC bars, indicating that our variational wave function— iorofthethree-andfour-componentsystemsisverysim- assuming the nodal surface is adequate— captures the ilar, and resembles that discussed above for Bose gases: key physics of the system. For comparison, circles and For small a , E increases if the Gaussian width b s VMC | | squares in Fig. 7 show the variational energy E for is larger than about a and if b is smaller than about VMC ho thefour-componentFermisystemwithN =16asafunc- a . For somewhatmore negative interspecies scattering ho tionofbfora = 0.05a and 0.07a ,respectively. A lengthsa ,E showsalocalminimumatb a and s ho ho s VMC ho − − ≈ comparisonofFigs.6and7showsthattheoverallbehav- alocalmaximumatb 0.3a ,andE becomesneg- ho VMC ≈