Borromean droplet in three-component ultracold Bose gases
BBorromean droplet in three-component ultracold Bose gases
Yinfeng Ma,
1, 2
Cheng Peng,
1, 2 and Xiaoling Cui
1, 3, ∗ Beijing National Laboratory for Condensed Matter Physics,Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China Songshan Lake Materials Laboratory, Dongguan, Guangdong 523808, China (Dated: February 2, 2021)Borromean ring refers to a peculiar structure where three rings are linked together while any twoof them are unlinked. Here we propose the realization of its quantum mechanical analog in a many-body system of three-component ultracold bosons. Namely, we identify the
Borromean droplet ,where only the ternary bosons can form a self-bound droplet while any binary subsystems cannot.Its formation is facilitated by an additional attractive force induced by the density fluctuation of athird component, which enlarges the mean-field collapse region in comparison to the binary case andrenders the formation of Borromean droplet after incorporating the repulsive force from quantumfluctuations. Outside the Borromean regime, the competition between ternary and binary dropletsleads to an interesting phenomenon of droplet phase separation, manifested by double plateaus inthe density profile. We further show that the transition between different droplets and gas phase canbe conveniently tuned by boson numbers and interaction strengths. The study reveals the possibilityof Borromean binding in the many-body world and sheds light on more intriguing many-body boundstate formed in multi-component systems.
Introduction.
As an exotic topological structure, Bor-romean ring has gained much attention in literature fromdifferent branches of science. It has been successfullyconstructed in both biology[1] and chemistry[2], and itsquantum mechanical analog, the Borromean binding, hasalso been reported in nuclear physics in terms of halonuclei[3, 4] and in ultracold gases as the Efimov effect[5–7]. In these occasions, the Borromean binding refers tothe trimer formation in few-body clusters where no dimeris present, such as the Efimov trimer observed in thenegative scattering length side[8–15]. Theoretical studieshave found that the Borromean trimer can be supportedby fine-tuning the pairwise short-range potential[16–20]or by modifying the single-particle dispersion[21]. In gen-eral, to afford such a peculiar bound state calls for strin-gent requirements that are usually difficult to meet inreality. As a result, the study of Borromean binding hasonly been limited in the few-body level, and its occur-rence in a many-body system has not been explored upto date.In the many-body world, droplet represents a typicalclass of bound state that has been well studied in heliumliquid[22, 23] and also proposed in Bose-Einstein systemsin early times by Huang[24]. Recently, there has been arevived study of droplet in dilute atomic gases, knownas the quantum droplet, following a pioneer theory byPetrov[25]. Such quantum droplet is stabilized by themean-field attraction and the Lee-Huang-Yang repulsionin three dimension from quantum fluctuations, and so farit has been successfully observed in dipolar gases[26–32]and in binary Bose gases of alkali atoms[33–36]. Theoret-ical study of quantum droplet has recently been extendedto low dimensions[37–42], Bose-Fermi mixtures[43–48],dipolar mixtures[49, 50] etc.In this work, we report the discovery of Borromean binding in a many-body system, namely, the
Borromeandroplet in three-component boson mixtures. Here theterm “Borromean” means that only ternary bosons canform the droplet while any binary subsystems cannot.The physical origin of such peculiar situation is thatthe third component can induce an additional attrac-tive force to the system via density fluctuations. Thisattraction further intensifies the mean-field collapse ascompared to the binary systems and renders the forma-tion of Borromean droplet after incorporating the force ofLee-Huang-Yang repulsion. Moreover, outside the Bor-romean regime, the competition between ternary andbinary droplets leads to an interesting phenomenon ofdroplet phase separation, characterized by two distinctplateaus in the droplet density profile. We further showthat the emergence of these different droplets can be con-veniently tuned by the boson number of each componentand the coupling strengths in-between. Experimentalrelevance of our results as well as extension to multi-component systems will also be discussed.
Model.
We start with the Hamiltonian for three-component bosons H = (cid:82) d r H ( r ), with ( (cid:126) = 1) H ( r ) = (cid:88) i =1 , , Ψ † i ( r )( − ∇ m i )Ψ i ( r )+ (cid:88) ij g ij † i Ψ † j Ψ j Ψ i ( r ) . (1)Here r is the coordinate; m i and Ψ i are respectively themass and field operator of boson species i ; g ij is the s-wave coupling strength between species i and j .For a homogeneous system with uniform densities { n i } ( i = 1 , , (cid:15) mf = 12 (cid:88) i,j =1 g ij n i n j (2) a r X i v : . [ c ond - m a t . qu a n t - g a s ] F e b Following the standard Bogoliubov theory to treatquantum fluctuations[51], we obtain the Lee-Huang-Yang(LHY) energy per volume as: (cid:15)
LHY = (cid:90) d k π ) (cid:88) i =1 ( E i k − (cid:15) i k − g ii n i ) + (cid:88) ij m ij g ij n i n j k . (3)Here E i k ( i = 1 , ,
3) are the Bogoliubov spectra[52].To describe a droplet with inhomogeneous densities, weadopt the local density approximation(LDA) and writethe total LHY energy as E LHY = (cid:82) d r (cid:15) LHY ( n i ( r )), with n i ( r ) = | Ψ i ( r ) | . This leads to three coupled Gross-Pitaevskii(GP) equations for { Ψ i ( r ) } ( i = 1 , , i∂ t Ψ i = − ∇ m i + (cid:88) j g ij | Ψ j | + ∂(cid:15) LHY ∂n i Ψ i , (4)The ground state can be approached by the imaginarytime evolution of above equations.In this work, to facilitate discussions while keeping theessence of physics, we consider the equal mass case m i ≡ m and the coupling strengths with following symmetry g = g ≡ g, g = g ≡ g (cid:48) . (5) Mean-field stability.
We first analyze the mean-field stability against density fluctuations of three-component(ternary) bosons, and compare it with thetwo-component(binary) cases. The stability is deter-mined by the second-order variation of (cid:15) mf with re-spect to small change of local densities δn i : δ (cid:15) mf = (cid:80) ij g ij δn i δn j , which gives: δ (cid:15) mf = g − g δn − + g + g δn + g δn + √ g (cid:48) δn δn + , (6)with δn ± ≡ ( δn ± δn ) / √ diagonalized fluctuationmodes for components 1 and 2. Eq.6 clearly shows thatthe fluctuation of component 3 will interfere with δn + and result in two new eigen-modes, while δn − is left un-changed. The mean-field stability requires δ (cid:15) mf > δn i , which leads to the following condition for a sta-ble ternary system: g > | g | ; g > g (cid:48) g + g ; (7)Note that the first condition ensures the stability of(1&2) system, while the second one is due to the inter-ference between (3) and (1&2) and ensures the stabilityof (1&2&3). Importantly, compared to the stability con-dition g (cid:48) < gg for (2&3) or (1&3), the requirementin (7) is more stringent. Therefore, there exists a finiteparameter window g (cid:48) g ∈ ( g + g , g ) , (8) (cid:1)(cid:3) (cid:2) (cid:3)(cid:1)(cid:6)(cid:1)(cid:5)(cid:1)(cid:4)(cid:2)(cid:4) (cid:5)(cid:5)(cid:1)(cid:3)(cid:2)(cid:5)(cid:6)(cid:5)(cid:5)(cid:5) (cid:5)(cid:5)(cid:1)(cid:4)(cid:2) (cid:2) (cid:1)(cid:2) (cid:1)(cid:2) (cid:5) (cid:1) / gg ¢ FIG. 1. (Color online) Mean-field phase diagram for three-component bosons under couplings (5) and g = g . The grayarea marks the mean-field stable region(“S”), smaller thanthat for binary subsystems (bounded by red square). Afterincorporating the LHY repulsion, Borromean droplet can takeplace in region I and II(A). The dashed lines separating II(A)and II(B) is determined by C min = 0 at the droplet-gas tran-sition (see text). such that all binary subsystems are stable against densityfluctuations while (1&2&3) is not.The intensified mean-field instability of ternary bosons,as compared to all binary subsystems, can be attributedto the additional attractive force brought by the thirdcomponent. To see this efficiently, let us consider a spe-cial case with g = 0 and start from the subsystem (1&3)whose stability condition is g (cid:48) < gg . This conditioncan be re-formulated as g + g ind >
0, with g ind = − g (cid:48) /g the induced interaction (to component 3) by the densityfluctuation of component 1[51]. Now if add the com-ponent 2 to (1&3), the fluctuation of component 2 willinduce an additional attraction to component 3 and now g ind = − g (cid:48) /g is doubled. Thus g needs to be more re-pulsive than in (1&3) case in order to stabilize (1&2&3).For a finite g , the fluctuations of 1 and 2 will coupletogether and give g ind = − g (cid:48) / ( g + g ), again moreattractive than the binary cases. We have checked thatsuch enhanced g ind in ternary system robustly applies formore general coupling strengths beyond (5).In Fig.1, we plot out the mean-field phase diagram of(1&2&3) system in ( g (cid:48) , g ) plane taking a fixed g = g .The mean-field stable region, as required by (7) and la-beled as “S”, is shown to be smaller than the stable re-gion of binary subsystems (bounded by red square). Forother regions in the diagram, a homogeneous (1&2&3)system will undergo a collapse or phase separation dueto density fluctuations, as determined by the eigen-modesof (6)[52]. Among them, there are four regions, labeledas I,II,III,IV in Fig.1, that all the three components un- (cid:2) (cid:3) (cid:4) (cid:5) (cid:6)(cid:3)(cid:2)(cid:3)(cid:2)(cid:2) (cid:1)(cid:4)(cid:2) (cid:3)(cid:3) (cid:4) (cid:1) (cid:2)(cid:1)(cid:3) (cid:2) (cid:3) (cid:2) (cid:3) (cid:4) (cid:5) (cid:6)(cid:2)(cid:1)(cid:2)(cid:3)(cid:2)(cid:1)(cid:3)(cid:3) (cid:1)(cid:6)(cid:2) (cid:2) (cid:3) (cid:4) (cid:5) (cid:6)(cid:2)(cid:1)(cid:2)(cid:3)(cid:2)(cid:1)(cid:3)(cid:3)(cid:3)(cid:2) (cid:1)(cid:5)(cid:2) FIG. 2. (Color online.) Droplet region(colored) in the N - C plane for three typical points (marked by “*”) on the verticalline in Fig.1, which have a fixed g (cid:48) /g = −√ / g /g = − . a ); − . b ); − c ). Red dashed lines show C min when the system has minimal total energy at given N . In (b),the yellow area marks the regime for droplet phase separation. dergo collapse simultaneously. This offers a possibilityfor droplet formation when further incorporating the re-pulsive force from quantum fluctuations. Of particularlyinteresting is regions I and II, as analyzed below. Borromean droplet.
It is obvious that for region I,which satisfies (8), the ternary system (1&2&3) can forma self-bound droplet while any binary subsystem cannot.By definition, this is the
Borromean droplet . Meanwhile,since the actual droplet formation also depends on par-ticle numbers { N i } , we will show that by tuning { N i } within certain range, such intriguing droplet can also ex-ist in other regions such as II(A). Given the couplingsymmetry in Eq.(5), we have N = N for the groundstate and there left two tunable parameters for { N i } : to-tal number N = 2 N + N and number ratio C = N /N .To explore essential properties of Borromean droplet,we carry out full simulations of the GP equations (4)to search for ground state with a fixed g (cid:48) /g = −√ / g /g in regions I and II, i.e., following thevertical line in Fig.1. In Fig.2(a-c), we show the areaof droplet formation in the N - C plane for three typicalvalues of g /g , where we also show the value of C whenthe energy reaches minimum for each N , denoted as C min (dashed lines in Fig.2). One can see that for g /g inregion I (Fig.2(a)), a ternary droplet can be supportedwhen N is beyond a critical number, N t,c , where a gas todroplet transition occurs. For all N > N t,c , the dropletcan only survive for C within a narrow window around C min . This characterizes the Borromean nature of thedroplet, i.e., its formation cannot extend to C = 0 (whenthe third component is absent).Interestingly, the Borromean droplet can extend topart of the region II. As an example, for the parame-ters considered in Fig.2(b), we can see that as increasing N , a ternary droplet first emerges at N t,c with a finite C . As N is further increased to N b,c , the binary (1&2)droplet appears at C = 0. It is then followed that theBorromean droplet is stabilized within the number win-dow N ∈ ( N t,c , N b,c ), where N is large enough to sup- FIG. 3. (Color online.) (a) C min (solid line) as a functionof g /g at the critical number N t,c , in comparison with C (0)min (dashed line) from Eq.(9). (b) Phase diagram in the N - g plane. The Borromean droplet (“BD”) occurs in region Ifor N > N t,c and in II(A) for N ∈ ( N t,c , N b,c ). For N < N t,c the system is in gas phase; for
N > N b,c the binary dropletcan also exist. Dashed lines show the function fits of N (0) t,c (Eq.11) and N (0) b,c (from Ref.[25]). Here g (cid:48) /g = −√ / port a ternary droplet but still small for the binary one.The Borromean droplet will vanish when go deep intoregion II. As shown in Fig.2(c) for large attractive g ,as increasing N the droplet solution first emerges at N b,c with C = 0. In this case, a binary droplet is more favoredthan a ternary one.Given the different behaviors of droplet formation inregion II, we have separated this region into II(A) andII(B) in Fig.1 — the former can support the Borromeandroplet (within certain number window) while the lattercannot. Their boundary (dashed line in Fig.1) is deter-mined by the zero crossing of C min at critical N t,c . InFig.3(a), we show C min at N t,c as a function of g /g witha given g (cid:48) /g = −√ /
2. We can see that C min continu-ously decreases as g becomes more attractive, and re-duces to zero at g /g ≈ −
4, which separates region II(A)from II(B) on the vertical line in Fig.1. In fact, C min canbe estimated from the minimization of (cid:15) mf , which gives C (0)min = g + g − g (cid:48) g − g (cid:48) . (9)In Fig.3(a), we can see that Eq.(9) can well fit C min inregion I, but deviate visibly as entering region II(A) whenthe system is far away from the mean-field collapse line.In Fig.3(b), we further map out the phase dia-gram highlighting the Borromean droplet(“BD”) in the( g , N ) plane, taking a fixed g (cid:48) /g = −√ /
2. To summa-rize, the Borromean droplet occurs at
N > N t,c in regionI and N ∈ ( N t,c , N b,c ) in region II(A). We can see thatboth N t,c and N b,c decrease as g gets more attractive.Now we analytically estimate N t,c . First, we inves-tigate the equilibrium density of Borromean droplet byenforcing the zero pressure P = (cid:80) i n i ∂(cid:15)/∂n i − (cid:15) = 0,where (cid:15) = (cid:15) mf + (cid:15) LHY . Utilizing (cid:15) mf = gn f with f = C / Cg (cid:48) /g + 1 + g /g , and the LHY energy at themean-field collapse line[52] (cid:15) LHY = 8 / (15 π )( gn ) / f ,with f = (1 + g /g + C ) / + (1 − g /g ) / , we obtainthe density of component- i in the ternary droplet n (0) t,i = η i π a (cid:18) f f (cid:19) . (10)Here a = mg/ (4 π ), η = η = 1 and η = C (0)min . Further,based on (9,10) and the single-mode assumption Ψ i ( r ) = (cid:113) n (0) t,i ψ ( r ), the coupled GP equations (4) can be reducedto a single one similar to that in the binary case[25]. Thisresults in the following critical number at the transitionbetween ternary droplet and gas phase: N (0) t,c = (2 + C ) / (cid:18) (cid:19) / N c π f f , (11)with ˜ N c = 18 .
65 (at the vanishing of droplet solution).In Fig.2(c), we show that (11) fits numerical N t,c qualita-tively well over the parameter regime considered. When C = 0, (10,11) recover the results for binary droplets[25]. Droplet phase separation.
Outside the Borromeanregime, both the ternary and binary droplets can sur-vive and will directly compete with each other. Here wewill demonstrate an interesting phenomenon of dropletphase separation as below.We consider the region II(A) and with a large N ( >N b,c ), i.e., above the “BD” region in Fig.3(b). In thiscase, as shown in Fig.1(b), the droplet solution can ap-pear in a reasonably broad window of C staring fromzero and its energy minimum occurs at finite C min .Among these C -values, C = 0 and C = C min representtwo typical solutions corresponding to, respectively, thebinary(1&2) and ternary(1&2&3) droplets. We find thatfor certain intermediate C , two types of droplet can co-exist in the form of phase separation, as shown in Fig.4,manifested by two different plateaus in the density pro-file that well fit the equilibrium densities of ternary andbinary droplets. Specifically, the (1&2&3) droplet occu-pies the central part with C ∼ C min and a higher density,while the (1&2) droplet stays at edge with C = 0 and alower density. Such distribution is believed to lower thesurface energy the most. (cid:1) (cid:2)(cid:1)(cid:1) (cid:3)(cid:1)(cid:1) (cid:4)(cid:1)(cid:1) (cid:5)(cid:1)(cid:1)(cid:2)(cid:1) (cid:1)(cid:3) (cid:2)(cid:1) (cid:1)(cid:2) ( 0 )b , 2( 0 )b , 1( 0 )t , 3 ( 0 )t , 2( 0 )t , 13 21 nnn nnn nn === (cid:1) (cid:1) FIG. 4. (Color online.) Density profile displaying thephase separation between ternary (1&2&3) and binary (1&2)droplets. Here g (cid:48) /g = −√ / , g /g = − . , N = 10 , C =0 .
2, corresponding to triangular point in Fig.2(b). Horizontallines show function fits to the equilibrium densities of ternary(Eq.10) and binary (from Ref.[25]) droplets. The length anddensity units are, respectively, a and 1 /a ( a = mg/ (4 π )). Here we estimate the parameter regime in (
N, C ) planeto support such droplet phase separation. For given N and C , we have N = N C/ ( C + 2) and N = N = N/ ( C + 2). Since the full number of component 3, to-gether with part of 1&2 components, occupy at the cen-ter to form (1&2&3) droplet with number ratio C min ,we can obtain the total number of ternary droplet as N t = N ( C min + 2) /C min . The rest 1&2 componentswill be left to form the binary droplet with total num-ber N b = 2 N − N /C min . The appearance of twotypes of droplets thus require N t > N t,c and N b > N b,c ,which set the constraint for allowed N - C values. For anygiven N ( > N b,c ), the allowed C is within certain window( C L , C H ), as shown by yellow region in Fig.2(b). For verylarge N , we have C L ∼ C H ∼ C min , and thus thedroplet phase separation can occur for any C ∈ (0 , C min ). Summary and discussion.
In summary, we have re-vealed intriguing physics of droplet formation in three-component boson mixtures, including the enhanced den-sity fluctuations towards mean-field collapse, the occur-rence of Borromean droplet, as well as the competitionand phase separation between different types of droplets.Though we have focused on the equal mass case withcertain coupling symmetries (5), the underlying physicsrevealed here can be extended to more general case ofmass ratios and coupling strengths.For the experimental detection of our results, espe-cially regions I and II(A) in Fig.1, one would need thethree-component bosons to hold (i) all repulsive intra-species couplings and (ii) at least two inter-species cou-plings to be attractive. For the three | F = 1 (cid:105) hy-perfine states of K atom, (ii) can be satisfied easilywhile (i) can be difficult given the generic negative back-ground couplings[53, 54]. However, as more and moreFeshbach resonances are explored between different al-kali bosons, such as K- Rb[55], Na- Rb[56], K- Cs[57], etc, it would be promising in future to probethe three-component droplets in hetero-nuclear mixtures.Finally, we discuss the possibility of droplet formationwith high-order Borromean structure, i.e., the
Brunnian ring[58–60]. By definition, the n -th order Brunnian ringdescribes n rings linked together while any subsystem of m ( < n ) rings is unlinked. Here we remark that its quan-tum mechanical analog in a many-body system, the n -thorder Brunnian droplet, can also exist. It means thatonly the n -component bosons form a self-bound statewhile any m ( < n )-component cannot. The underlyingmechanism resembles that of Borromean droplet revealedin this work, i.e., an extra component will bring addi-tional attractive force to the system via density fluctua-tions. Take the example of a simple case where g ii = g ( i = 2 , ...n ) and all inter-species couplings are zero ex-cept g i = g (cid:48) , the induced interaction to component 1is g ind = − ( n − g (cid:48) /g , which gets more attractive if n is larger. As a result, for gg /g (cid:48) ∈ ( n − , n − n -component mixture undergoes the mean-field col-lapse while any subsystem of m ( < n )-component is sta-ble. Further joining the LHY repulsive force, this willsupport the formation of n -th order Brunnian droplet,and the Borromean droplet discussed in this work can beseen as the n = 3 special case within this generalization. Acknowledgment.
The work is supported by the Na-tional Key Research and Development Program of China(2016YFA0300603, 2018YFA0307600), the National Nat-ural Science Foundation of China (No.12074419), and theStrategic Priority Research Program of Chinese Academyof Sciences (No. XDB33000000). ∗ [email protected][1] C. Mao, W. Sun and N. C. Seeman, Nature , 137(1997).[2] K. S. Chichak, S. J. Cantrill, A. R. Pease, S.-H. Chiu, G.W. V. Cave, J. L. Atwood, and J. F. Stoddart, Science , 1308 (2004).[3] M. V. Zhukov, B. V. Danilin, D. V. Fedorov, J. M. Bang,I. S. Thompson, and J. S. Vaagen, Phys. Rep. , 151(1993).[4] D. V. Fedorov, A. S. Jensen, and K. Riisager, Phys. Rev.C , 201 (1994).[5] E. Braaten and H.-W. Hammer, Phys. Rep. , 259(2006).[6] C.H. Greene, P. Giannakeas, and J. P´erez-R´ıos, Rev.Mod. Phys. , 035006 (2017).[7] P. Naidon, S. Endo, Rep. Prog. Phys. , 056001(2017).[8] T. Kraemer, M. Mark, P. Waldburger, J. G. Danzl, C.Chin, B. Engeser, A. D. Lange, K. Pilch, A. Jaakkola,H.-C. N¨agerl and R. Grimm, Nature , 315 (2006).[9] T. B. Ottenstein, T. Lompe, M. Kohnen, A. N. Wenz,and S. Jochim, Phys. Rev. Lett. , 203202 (2008). [10] J. R. Williams, E. L. Hazlett, J. H. Huckans, R. W.Stites, Y. Zhang, and K. M. O’Hara, Phys. Rev. Lett. , 130404 (2009).[11] M. Zaccanti, B. Deissler, C. D’Errico, M. Fattori, M.Jona-Lasinio, S. M¨uller, G. Roati, M. Inguscio and G.Modugno, Nat. Phys. , 586 (2009).[12] N. Gross, Z. Shotan, S. Kokkelmans and L. Khaykovich,Phys. Rev. Lett. , 163202 (2009).[13] S. E. Plooack, D. Dries and R. G. Hulet, Science ,1683 (2009).[14] M. Berninger, A. Zenesini, B. Huang, W. Harm, H.-C.N¨agerl, F. Ferlaino, R. Grimm, P. S. Julienne and J. M.Hutson, Phys. Rev. Lett. , 120401 (2011).[15] R. J. Wild, P. Makotyn, J. M. Pino, E. A. Cornell andD. S. Jin, Phys. Rev. Lett. , 145305 (2012).[16] J.-M. Richard and S. Fleck, Phys. Rev. Lett. , 1464(1994).[17] S. Moszkowski, S. Fleck, A. Krikeb, L. Theussl, J.M.Richard, and K. Varga, Phys. Rev. A , 032504 (2000)[18] E. Nielsen, D. V. Fedorov, and A. S. Jensen, Few-BodySystems , 15 (1999);[19] A. G. Volosniev, D. V. Fedorov, A. S. Jensen, and N. T.Zinner, Eur. Phys. J. D , 95 (2013).[20] A. G. Volosniev, D. V. Fedorov, A. S. Jensen, and N. T.Zinner, arxiv: 1312.6535.[21] X. Cui, W. Yi, Phys. Rev. X , 031026 (2014).[22] S. Grebenev, J. P. Toennies, A. F. Vilesov, Science ,2083 (1998).[23] M. Barranco, R. Guardiola, S. Hernandez, R. Mayol, J.Navarro, J. Low Temp. Phys. , 1(2006)[24] K. Huang, Phys. Rev. , 765 (1959); Phys. Rev. ,1129 (1960).[25] D.S. Petrov, Phys. Rev. Lett. , 155302 (2015).[26] I. Ferrier-Barbut, H. Kadau, M. Schmitt, M. Wenzel, andT. Pfau, Phys. Rev. Lett. , 215301 (2016).[27] M. Schmitt, M. Wenzel, F. B¨ottcher, I. Ferrier-Barbut,and T. Pfau, Nature , 259 (2016).[28] I. Ferrier-Barbut, M. Schmitt, M. Wenzel, H. Kadau, andT. Pfau, J. Phys. B , 214004 (2016).[29] L. Chomaz, S. Baier, D. Petter, M.J. Mark, F. W¨achtler,L. Santos, and F. Ferlaino, Phys. Rev. X , 041039(2016).[30] L. Tanzi, E. Lucioni, F. Fama, J. Catani, A. Fioretti, C.Gabbanini, R. N. Bisset, L. Santos, and G. Modugno,Phys. Rev. Lett. , 130405 (2019).[31] F. B¨ottcher, J.-N. Schmidt, M. Wenzel, J. Hertkorn, M.Guo, T. Langen, and T. Pfau, Phys. Rev. X , 011051(2019).[32] L. Chomaz, D. Petter, P. Ilzh¨ofer, G. Natale, A. Traut-mann, C. Politi, G. Durastante, R.M.W. van Bijnen, A.Patscheider, M. Sohmen, M.J. Mark, and F. Ferlaino,Phys. Rev. X , 021012 (2019).[33] C.R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas,P. Cheiney, and L. Tarruell, Science , 301 (2018).[34] P. Cheiney, C. R. Cabrera, J. Sanz, B. Naylor, L. Tanzi,L. Tarruell, Phys. Rev. Lett. , 135301 (2018).[35] G. Semeghini, G. Ferioli, L. Masi, C. Mazzinghi, L. Wol-swijk, F. Minardi, M. Modugno, G. Modugno, M. Ingus-cio, M. Fattori, Phys. Rev. Lett. , 235301 (2018).[36] C. D’Errico, A. Burchianti, M. Prevedelli, L. Salasnich,F. Ancilotto, M. Modugno, F. Minardi, and C. Fort,Phys. Rev. Research , 033155 (2019).[37] D. S. Petrov and G. E. Astrakharchik, Phys. Rev. Lett. , 100401 (2016). [38] D. Edler, C. Mishra, F. W¨achtler, R. Nath, S. Sinha, andL. Santos, Phys. Rev. Lett. , 050403 (2017).[39] K. Jachymski and R. Oldziejewski, Phys. Rev. A ,043601 (2018).[40] P. Zin, M. Pylak, T. Wasak, M. Gajda, and Z. Idziaszek,Phys. Rev. A , 051603(R) (2018).[41] T. Ilg, J. Kumlin, L. Santos, D. S. Petrov, and H. P.B¨uchler, Phys. Rev. A , 051604(R) (2018).[42] X. Cui, Y. Ma, arxiv:2010.10723.[43] X. Cui, Phys. Rev. A , 023630 (2018).[44] S. Adhikari, Laser Phys. Lett , 095501 (2018).[45] D. Rakshit, T. Karpiuk, M. Brewczyk, and M. Gajda,SciPost Phys. , 079 (2019).[46] D. Rakshit, T. Karpiuk, P. Zin, M. Brewczyk, M. Lewen-stein, and M. Gajda, New J. Phys. , 073027 (2019).[47] M. Wenzel, T. Pfau and I. Ferrier-Barbut, PhysicaScripta , 10 (2018).[48] J.-B. Wang, J.-S. Pan, X. Cui, W. Yi, Chin. Phys. Lett. , 076701 (2020).[49] Joseph C. Smith, D. Baillie, and P. B. Blakie,arxiv:2007.00366.[50] R. N. Bisset, L. A. Pe˜ n a Ardila, and L. Santos,arxiv:2007.00404.[51] C. J. Pethick and H. Smith, Bose-Einstein Condensationin Dilute Gases , Cambridge University Press, 2002.[52] See supplementary material for the more details on thederivation of mean-field instability and the LHY energy.[53] M. Lysebo and L. Veseth, Phys. Rev A , 032702 (2010).[54] S. Roy et al, Phys. Rev. Lett. , 053202 (2013).[55] L. Wacker et al , Phys. Rev. A , 053602 (2015).[56] F. Wang et al , J. Phys. B , 015302 (2015).[57] M. Gr¨obner et al , Phys. Rev. A , 022715 (2017).[58] H. Brunn, Math. Phys. Klasse, , 77(1892).[59] H. Debrunner, Duke Math. J., , 17 (1961); D.E. Pen-ney, Duke Math. J., ,31(1969); T. Yanagawa, Osaka J.Math., ,127 (1964).[60] Nils A. Baas, D. V. Fedorov, A. S. Jensen, K. Riisager,A. G. Volosniev, N. T. Zinner, Physics of Atomic Nuclei, , 361 (2014). Supplemental Materials
In this Supplemental Material, we provide more details on the derivations of mean-field instability and the LHYenergy for three-component bosons.
I. Mean-field instability against density fluctuations
As discussed in the main text, the mean-field stability in the presence of density fluctuations is determined by thesecond-order variation of mean-field energy δ (cid:15) mf , as shown by Eq.(6) in the main text. Further diagonalizing δ (cid:15) mf leads to: δ (cid:15) mf = ˜ g δ ˜ n + ˜ g δ ˜ n + ˜ g δ ˜ n , (12)where the eigen-fluctuation-energy reads ˜ g = g − g , (13)˜ g = g + g + g − ∆4 , (14)˜ g = g + g + g + ∆4 , (15)with ∆ = (cid:112) g (cid:48) + ( g + g − g ) . The according eigen-fluctuation density reads δ ˜ n = 1 √ δn − δn ) , (16) δ ˜ n = (cid:115) g (cid:48) ∆(∆ + ( g + g − g )) (cid:18) δn + δn − ∆ + ( g + g − g )2 g (cid:48) δn (cid:19) , (17) δ ˜ n = (cid:115) g (cid:48) ∆(∆ − ( g + g − g )) (cid:18) δn + δn + ∆ − ( g + g − g )2 g (cid:48) δn (cid:19) , (18)Here δ ˜ n is just δn − in the main text, and δ ˜ n , δ ˜ n are the two eigen-modes due to coupling between δn + and δn . (cid:1)(cid:4) (cid:1)(cid:3) (cid:2) (cid:3) (cid:4)(cid:1)(cid:5)(cid:2)(cid:5)(cid:6) (cid:7)(cid:6)(cid:6)(cid:6) (cid:1) (cid:1) (cid:4)(cid:3) (cid:2) (cid:1)(cid:2) (cid:1)(cid:2) (cid:9)(cid:1)(cid:2)(cid:9) (cid:5)(cid:6)(cid:6)(cid:8) (cid:6)(cid:6) FIG. 5. Mean-field phase diagram for three-component bosons with couplings g ii = g ( i = 1 , , , g = g = g (cid:48) . Themean-field stable region is labeled as “S”. The rest ones, labeled as “I,II,III,IV” and “A,B,C,D”, are all mean-field unstableregions. To ensure the mean-field stability, all the three eigen-modes should have positive energies, i.e., all ˜ g i >
0. Thisleads to the stability condition presented as Eq.7 in the main text. In Fig.5, we show the mean-field phase diagram
Region Label g /g g (cid:48) /g − (1,2);1 −
2A (1 , ∞ ) ( −∞ , − g (cid:48) c ) C;PSB (1 , ∞ ) ( − g (cid:48) c , g (cid:48) c ) S;PSC (1 , ∞ ) ( g (cid:48) c , ∞ ) PS;PSD ( −∞ ,
1) ( g (cid:48) c , ∞ ) PS;CS (-1,1) ( − g (cid:48) c , g (cid:48) c ) S;SI ( − ,
1) ( − , − g (cid:48) c ) C;CII ( −∞ , −
1) ( − ,
0) C;CIII ( − ,
1) ( −∞ , −
1) C;CIV ( −∞ , −
1) ( −∞ , −
1) C;CTABLE I. Mean-field instability for each region in Fig.5. The last column lists the tendency of the system due to densityfluctuations between 3 and (1,2) (“3 − (1,2)”), and between 1 and 2 (“1 − g (cid:48) c = (cid:113) g /g . for three-component bosons with couplings g ii = g ( i = 1 , , , g = g = g (cid:48) . The mean-field stable region is labeledas “S”, and the rest unstable regions are labeled as “I,II,III,IV” and “A,B,C,D”. Next we will analyze the fate of ahomogeneous ternary mixture under density fluctuations in these unstable regions. • Region A: ˜ g < , ˜ g < , ˜ g >
0. To lower the energy, the system tends to have finite δ ˜ n , δ ˜ n and zero δ ˜ n .Thus one has δn (cid:54) = δn and δn + δn ∝ δn . This implies that the density fluctuations lead to a collapsebetween 3 and (1,2) mixture, while (1,2) itself tends to phase separate. • Region B: ˜ g < , ˜ g > , ˜ g >
0, which leads to a finite δ ˜ n while δ ˜ n = δ ˜ n = 0. Hence we have δn = − δn , δn = 0, which means that (1,2) phase separate, and 3 remains stable. • Region C: ˜ g < , ˜ g < , ˜ g >
0, which leads to finite δ ˜ n , δ ˜ n while δ ˜ n = 0. This gives δn (cid:54) = δn , δn + δn ∝ − δn . Therefore there is a phase separation between 3 and (1,2), and 1, 2 themselves also phaseseparate. • Region D: ˜ g > , ˜ g < , ˜ g >
0, which leads to finite δ ˜ n and δ ˜ n = δ ˜ n = 0. This gives δn = δn , δn + δn ∝ − δn , and therefore there is a phase separation between 3 and (1,2), while 1, 2 themselves collapse. • Region I, II, III, IV: ˜ g > , ˜ g < , ˜ g >
0, which, similar to region D, leads to finite δ ˜ n and δ ˜ n = δ ˜ n = 0.However, g (cid:48) in this region has an opposite sign with region D. Therefore we have δn = δn and δn + δn ∝ δn .This means that 3 and (1,2) tend to collapse, and 1,2 themselves also collapse. These are the regions that thethree components all undergo collapse under density fluctuations.Above results are summarized in Table.I, where we have listed the parameter regime for each region and thetendency of three components with respect to density fluctuations, including collapse(C), phase separation(PS) andstable(S). II. Lee-Huang-Yang energy for three-component bosons
Based on the standard Bogoliubov theory, we expand the field operator as:Ψ i = √ n i + (cid:88) k (cid:54) =0 √ V exp( i kr ) θ i k (19)Here θ i k is the fluctuation operator for component- i boson at momentum k . Then the Hamiltonian canbe transformed to the bilinear form H/V = (1 /V ) (cid:80) k φ † h k φ + (cid:15) mf − (1 /V ) (cid:80) i k ( (cid:15) i k + g ii n i ), where φ = (cid:16) θ k , θ k , θ k , θ † − k , θ † − k , θ † − k (cid:17) T , and h k = (cid:15) k + g n g √ n n g √ n n g n g √ n n g √ n n g √ n n (cid:15) k + g n g √ n n g √ n n g n g √ n n g √ n n g √ n n (cid:15) k + g n g √ n n g √ n n g n g n g √ n n g √ n n (cid:15) k + g n g √ n n g √ n n g √ n n g n g √ n n g √ n n (cid:15) k + g n g √ n n g √ n n g √ n n g n g √ n n g √ n n (cid:15) k + g n (20)The LHY energy can then be derived as Eq.3 in the main text, with E i k ( i = 1 , ,
3) the three Bogoliubov excitationenergies. After straightforward algebra, we find that E i k are the three roots of following equation x + bx + cx + d = 0 , (21)where b = − (cid:88) i ω i , (22) c = (cid:88) i 3) are the Bogoliubov spectra for the individual components. Under the equalmass case m = m = m = m, (cid:15) k = (cid:15) i k and coupling symmetry g = g = g and g = g = g (cid:48) , we can analyticallywrite down the three Bogoliubov energies at the mean-file collapse line g (cid:48) = − (cid:112) g ( g + g ) / E k = (cid:15) k E k = (cid:104) (cid:15) k + (cid:15) k (cid:16) g ( n + n ) + g n − (cid:112) g ( n − n ) + ( g n + 2 g n )( g n + 2 g n ) (cid:17)(cid:105) / E k = (cid:104) (cid:15) k + (cid:15) k (cid:16) g ( n + n ) + g n + (cid:112) g ( n − n ) + ( g n + 2 g n )( g n + 2 g n ) (cid:17)(cid:105) / (25)One can see that at the mean-file collapse line, one mode E k becomes quadratic while the other modes are still linearat low energy. In this case, LHY energy can be simplified as ε LHY = √ π (cid:104) ( α + β ) / + ( α − β ) / (cid:105) (26)where α = g ( n + n ) + g n , β = (cid:112) g ( n − n ) + ( g n + 2 g n )( g n + 2 g n ). When g = g, n = n and n /n = C , the LHY energy is given by (cid:15) LHY = π f , with f = (1 + g /g + C ) / + (1 − g /g ) / . This is usedto estimate n (0) ii