Quantum Phases of Self-Bound Droplets of Bose-Bose Mixtures
QQuantum Phases of Self-Bound Droplets of Bose-Bose Mixtures
Junqiao Pan,
1, 2
Su Yi,
1, 2, 3, 4, ∗ and Tao Shi
1, 3, † CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics,Chinese Academy of Sciences, Beijing 100190, China School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China CAS Center for Excellence in Topological Quantum Computation,University of Chinese Academy of Sciences, Beijing 100049, China Shenzhen Institute for Quantum Science and Engineering,Southern University of Science and Technology, Shenzhen 518055, China (Dated: March 2, 2021)We systematically investigate the ground-state properties of self-bound droplets of quasi-two-dimensional binary Bose gases by using the Gaussian state theory. We find that quantum dropletsconsists two macroscopic squeezed phases and a macroscopic coherent phase. We map out the phasediagram and determine all phase boundaries via both numerical and nearly analytical methods. Inparticular, we find three easily accessible signatures for the quantum phases and the stablizationmechanism of the self-bound droplets by precisely measuring their radial size. Our studies indicatethat binary droplets represent an ideal platform for in-depth investigations of the quantum natureof the droplet state.
Introduction .—Quantum droplets in both dipolar [1–4] and binary condensates [5–9] represent a new stateof matter emerging in the mean-field unstable regime.They have attracted great interests in studying theirfascinating properties, such as self-bound feature [2, 3,5, 7], collective excitations [10, 11], soliton to droplettransition [6], droplet-droplet collision [12], supersolidstates [13–15], Goldstone mode [16], and so forth (seerecent reviews [17–19] and references therein). A widelyadopted treatment for droplet states is the extend Gross-Pitaevskii equation (EGPE) which perturbatively incor-porates quantum fluctuation into the Gross-Pitaevskiiequation in terms of the Lee-Huang-Yang correction [20–24]. Although EGPE has provided satisfactory explana-tions to experimental observations, there are still discrep-ancies with experimental measurements [4, 5, 18].Theoretical methods beyond EGPE have also been em-ployed to study different aspects of the droplet states [25–34]. Among them, we studied the dipolar droplets usingthe Gaussian-state theory (GST) in which quantum fluc-tuation is included in a self-consistent manner [35]. Com-pared to other theoretical approaches beyond the EGPEtheory, GST is computationally efficient and can be ap-plied to realistic systems with large number of atoms. In-terestingly, we found that, besides coherent state, dipolardroplets also contained two new macroscopic squeezedphases which were characterized by large second-ordercorrelation and asymmetric atom-number distribution(AND), in striking contrast with those of the coherent-state-based droplets [35, 36]. Whereas the asymmetricAND was experimentally observed [2, 4], verifying thepredication on the second-order correlation function re-quires further experiments beyond simple density mea- ∗ Electronic address: [email protected] † Electronic address: [email protected] surements.In this Letter, we study the quantum phases of self-bound droplets in quasi-two-dimensional (quasi-2D) bi-nary condensates. Although this system has same quan-tum phases as those in dipolar droplets, here we discoverthree readily accessible signatures for the determinationof the quantum states and the stablization mechanisms.Specifically, due to the multiple quantum phases contain-ing in the droplet states, the radial size ( σ ) versus atomnumber ( N ) curve is of W shape, as opposed to the V-shaped curve expected for single quantum phase. There-fore, the observation of double dips on the σ - N curvemay rule out the pure coherent explanation of the dropletstate. Moreover, the critical atom number (CAN), i.e.,the minimal number of atoms required for forming a self-bound state, is determined by the quantum states of thecondensate and its precise measurement of the CAN al-lows us to distinguish squeezed state from coherent one.Finally, the dip atom number (DAN), i.e., the atom num-ber at the dip of the σ - N curve, depends on both quan-tum phase and stabilization mechanism, which makes ita quantitative criterion for stability mechanism. Formulation .—We consider a ultracold gas of N Katoms with N ↑ atoms being in the hyperfine state |↑(cid:105) ≡| F, m F (cid:105) = | , − (cid:105) and N ↓ in |↓(cid:105) ≡ | , (cid:105) . The totalHamiltonian of the system, H = H + H B + H B , consists of the single-, two-, and three-particle terms. Insecond-quantized form, the single-particle part reads H = (cid:88) α (cid:90) d r ˆ ψ † α ( r ) h α ˆ ψ α ( r ) , (1)where ˆ ψ α ( r ) ( α = ↑ , ↓ ) are the field operators and h α = − (cid:126) ∇ / (2 M ) − µ α with M being the mass of the atomsand µ α the chemical potential introduced to fixed the a r X i v : . [ c ond - m a t . qu a n t - g a s ] F e b atom number in the α th component. The two-body (2B)interaction Hamiltonian takes the form H B = (cid:88) αβ g αβ (cid:90) d r ˆ ψ † α ( r ) ˆ ψ † β ( r ) ˆ ψ β ( r ) ˆ ψ α ( r ) , (2)where g αβ = 4 π (cid:126) a αβ /M characterize the 2B interac-tion strengths with a αβ being scattering length betweencomponents α and β . The scattering lengths are tunablethrough Feshbach resonance and for scenario of interestto quantum droplets, the intra- and inter-species scatter-ing lengths satisfy a ↑↑ , a ↓↓ > a ↑↓ <
0, respectively.Finally, the three-body (3B) interaction Hamiltonian canbe expressed as H B = g (cid:88) αβγ (cid:90) d r ˆ ψ † α ( r ) ˆ ψ † β ( r ) ˆ ψ † γ ( r ) ˆ ψ γ ( r ) ˆ ψ β ( r ) ˆ ψ α ( r ) , (3)where g is the 3B coupling constant which, for sim-plicity, is assumed to be independent of the spin com-ponents. Moreover, since we are only interested in theground states of the system, we assume that g is realand positive. The value of g shall be determined byfitting the experimental data.We study self-bound droplet states using the GST.Specifically, we assume that the many-body wave func-tion of a binary gas adopts the variational ansatz [36–38] | Ψ GS (cid:105) = e ˆΨ † ( r )Σ z ( r , r (cid:48) )Φ( r (cid:48) ) e i ˆΨ † ( r ) ξ ( r , r (cid:48) ) ˆΨ( r (cid:48) ) / | (cid:105) , (4)where ˆΨ( r ) ≡ (cid:18) ˆ ψ ( r )ˆ ψ † ( r ) (cid:19) with ˆ ψ ( r ) = (cid:18) ˆ ψ ↑ ( r )ˆ ψ ↓ ( r ) (cid:19) is the fieldoperator expressed in the Nambu basis,Σ z ( r , r (cid:48) ) = (cid:18) I ⊗ δ ( r − r (cid:48) ) 00 − I ⊗ δ ( r − r (cid:48) ) (cid:19) (5)with I being a 2 × r ) = (cid:104) ˆΨ( r ) (cid:105) = (cid:18) φ ( c ) φ ( c ) ∗ (cid:19) with φ ( c ) = (cid:32) φ ( c ) ↑ φ ( c ) ↓ (cid:33) is thevariational parameter describing the coherent part of thecondensates, and ξ ( r , r (cid:48) ) is the variational parameter thatcharacterizes the squeezed part of the condensate. Theground state wave function can be found by minimiz-ing energy E = (cid:104) Ψ GS | H | Ψ GS (cid:105) through imaginary-timeevolution [36, 37, 39].Physically, it is more convenient to factorize the Gaus-sian state wave function into a multimode squeezed co-herent state [35, 39], i.e., | Ψ GS (cid:105) = e √ N ( c ) ( ˆ c † − ˆ c ) ∞ (cid:89) i =1 e ξ i ( ˆ s † i − ˆ s i ) | (cid:105) , (6)where N ( c ) = (cid:82) d r (cid:2) φ ( c ) ( r ) (cid:3) † φ ( c ) ( r ) is the atom num-ber in the coherent mode and ˆ c = (cid:82) d r (cid:2) ¯ φ ( c ) ( r ) (cid:3) † ˆ ψ ( r )with ¯ φ ( c ) ( r ) = φ ( c ) ( r ) / √ N ( c ) being the normalized mode function for the coherent component. In the squeez-ing operators, ˆ s i = (cid:82) d r (cid:104) ¯ φ ( s ) i ( r ) (cid:105) † ˆ ψ ( r ) is the annihila-tion operator for the i th squeezed mode ¯ φ ( s ) i ≡ (cid:32) φ ( s ) i, ↑ φ ( s ) i, ↓ (cid:33) that satisfy (cid:82) d r (cid:104) ¯ φ ( s ) i ( r ) (cid:105) † ¯ φ ( s ) j ( r ) = δ ij . Furthermore,sinh ξ i = (cid:113) N ( s ) i with N ( s ) i being the occupation numberin the i th squeezed mode. The total number of atomsin the squeezed modes is N ( s ) = (cid:80) i N ( s ) i and the totalsqueezed density is n ( s ) ( r ) = (cid:80) i N ( s ) i | ¯ φ ( s ) i ( r ) | . For con-venience, it is always assumed that N ( s ) i are sorted indescending order with respect to the index i . Mode j isnotably populated if N ( s ) j /N ( s ) ≥ .
1% [35].In the main text, we focus on the self-bound droplets inquasi-2D geometry achieved by imposing a harmonic con-finement, V z ( z ) = M ω z z /
2, along the z axis [5]. When ω z is sufficiently large, the motion of atoms along the z direction is frozen to the ground state of V z ( z ). In addi-tioin, we always fix the atom number ratio between twocomponents at N ↑ /N ↓ = (cid:112) a ↓↓ /a ↑↑ by following the ex-periments [5, 7]. Numerical results can be convenientlychecked using the virial relation, E k + E B + 2 E B = 0,where E k , E B , and E B are kinetic, 2B, and 3B ener-gies, respectively [39]. We point out that results for 3Ddroplets are presented in the Supplemental Material [39]. Three-body coupling strength .—Let us first fix the valueof g for K atom. Previously, g of Dy atom was deter-mined by fitting the AND of the droplet states [35]. Herewe instead determine g by fitting the radial size σ of thequasi-2D droplets. Following the experiment [5], we ex-tract σ by fitting the total density with a two-dimensionalGaussian. In Fig. 1(a), we plot the radial size σ as a func-tion of the effective scattering length δa ≡ a ↑↓ + √ a ↑↑ a ↓↓ (or, equivalently, the magnetic field B ) computed with g = 6 . × − (cid:126) m / s. Good agreement between thenumerical and the experimental results is achieved. Wepoint out that varying the value of g does not changethe shape of the red line; it only shifts the red line ver-tically. Therefore, to achieve better agreement for thefitting, one may have to make g spin dependent.To further demonstrate the validity of the fitted g , weshow, in Fig. 1(b) and (c), the N dependence of σ GST under different magnetic fields. As can be seen, althoughdiscrepancies exist at small N , good agreements are stillachieved for large N . Particularly, the main feature of theexperimental results is captured by the numerical simu-lations (see below for details). Hereafter, we shall alwaysuse the fitted g for all numerical results presented below. Phase diagram .—Given that N ↑ /N ↓ = (cid:112) a ↓↓ /a ↑↑ , it isfound numerically that the coherent and squeezed modessatisfy φ ( c ) ↑ ( r ) /φ ( c ) ↓ ( r ) = (cid:112) N ↑ /N ↓ and φ ( s ) i, ↑ ( r ) /φ ( s ) i, ↓ ( r ) = (cid:112) N ↑ /N ↓ for any i , respectively. As a result, N ( c ) α /N α isindependent of the spin component α , which allows us todefine the coherent fraction as f c ≡ N ( c ) /N to character-ize the property of the condensate, instead of considering -5.69 -4.95 -4.18 -3.37 -2.530510 (a) (b) (c) FIG. 1: (color online). (a) Radial size versus reduced scat-tering length for N = 1 . × . Solid line is computed usingthe GST with g = 6 . × − (cid:126) m / s. Empty circles witherror bars and dashed line (both are extracted from Ref. [5])represent the experimental data and the EGPE result, respec-tively. (b) and (c) show radial size as functions of atom num-ber for various magnetic fields. Empty circles (extracted fromRef. [5]) are experimental data. More comparisons betweenthe numerically computed radial size and the experimentaldata can be found in the Supplemental Material [39]. the coherent fraction of individual spin component.Figure 2(a) shows the distribution of f c on the ( δa, N )parameter plane. Similar to dipolar droplets [35], thereare four distinct regions: expanding gas (EG), squeezed-vacuum state (SVS), squeezed-coherent state (SCS), andcoherent state (CS). The EG region lies below the CAN(solid line or (cid:52) ), in which the attractive two-body inter-action is insufficient for the gas to form a self-bound state.Both SVS and SCS phases contain a single squeezedmode except SCS phase also consists of a coherent com-ponent. The SVS-to-SCS transition (dash-dotted line or ♦ ) breaks the Z symmetry of the condensate and is athird-order phase transition. The CS phase is dominatedby the coherent component with f c > .
75 and the SCS-to-CS transition is of first order. Interestingly, as shownby the solid and the dash-dotted lines, the boundaries forthe self-bound state and the SVS-to-SCS transition canbe found nearly analytically with high accuracy [39]. Wemay also determine the boundary of SCS-to-CS transi-tion by analyzing the stability of Bogliubov excitation inthe CS phase [39].The phase diagram is most conveniently reproduced byplotting the number of the notably-populated squeezedmodes S . As shown in Fig. 2(b), both SVS and SCSphases are featured by single squeezed mode; while CSphase may contains a large number of squeezed modeseven if the squeezed component becomes negligibly small. FIG. 2: (color online). Phase diagram. (a) Distribution of thecoherent fraction f c on the ( δa, N ) parameter plane. Markers (cid:52) , ♦ , and (cid:53) denote numerical results for CAN, SVS-to-SCStransition boundary, and SCS-to-CS transition boundary, re-spectively. The solid and dash-dotted lines show the corre-sponding analytic results. Empty circles ( (cid:35) ) with error barare experimental data for CAN [5]. (b)-(d) demonstrate thenumber of the notably-populated squeezed modes S , the dis-tribution of the peak density n p (units of m − ), and the radialsize σ GST (units of µ m) on the ( δa, N ) parameter plane, re-spectively. Physically, squeezing in SVS and SCS phases are macro-scopic quantum state originating from two-body attrac-tion [35, 36]; while squeezing in CS phase representsquantum depletion induced by 3B repulsion [35]. Similarto dipolar droplets, the presence of squeezing also signif-icantly modifies the AND of these quantum states [39].The properties of these quantum phases can be furtherexplored by examining the peak density n p and the ra-dial size σ GST as plotted in Fig. 2(c) and (d), respec-tively. Among them, σ GST exhibits visible difference ascompared to the distributions of other quantities.
Radial size .—To explore the properties of σ GST , weplot the N dependence of σ GST for δa = − . a inFig. 3. As a comparison, the radial size numerically com-puted via EGPE, σ EGPE , is also plotted. Surprisingly, σ GST is W-shaped function of N , in striking contrast tothe V-shaped σ EGPE . To understand its origin, we notethat the total energy per atom can be expressed as [39] ε = E ( σ ) N ∝ σ − ˜ g Nσ + ˜ g ν N ν − σ ν − (for ν > , (7)where E = E k + E B + E B are the total energy. On theright hand side, the first term is the kinetic energy, thesecond term is contributed by the 2B interaction, and thelast term is the energy associated with the stabilizationforce, i.e., ν = 3 for 3B repulsion and 5 / g ( >
0) and ˜ g ν ( >
0) are re-
FIG. 3: (color online). Radial size versus mean atom numberfor δa = − . a . Vertical dotted lines mark the locationsof two phase transitions. duced interaction parameters (RIPs). From Eq. (7), theequilibrium size can be immediately found to be σ = √ N (cid:20) ( ν − g ν N ˜ g N − (cid:21) / [2( ν − (for N > / ˜ g ) , (8)which is a V-shaped function of N . More specifically, σ diverges as N approaches the CAN N cri ≡ g , (9)grows asymptotically as σ ≈ √ N [( ν − g ν / ˜ g ] / [2( ν − for large N , and reaches its minimal value at the DAN N dip ≡ ν − ν − g . (10)In contrast, we note that the peak density, n p ∝ N/σ = { (˜ g N − / [( ν − g ν N ] } / ( ν − , is always a monotoni-cally increasing function of N and converges to a constantproportional to { ˜ g / [( ν − g ν ] } / ( ν − for large N . It isworthwhile to point out that σ ( N ) depends not only onthe stabilization mechanism through ν but also on themany-body wave function via the RIPs (˜ g , ˜ g ν ) [39].To carry out quantitative comparisons, we assume thatall density profiles are Gaussian functions. Then for puresqueezed and coherent states with 3B interaction, we ob-tain the RIPs (˜ g ( s )2 , ˜ g ( s )3 ) and (˜ g ( c )2 , ˜ g ( c )3 ), respectively, sat-isfying ˜ g ( s )2 = 3˜ g ( c )2 and ˜ g ( s )3 = 15˜ g ( c )3 [35, 39]. In addition,for the coherent-state-based EGPE theory, the RIPs are(˜ g ( c )2 , ˜ g ( c )5 / ) [39]. Correspondingly, we plot σ ( c,s ) ν ( N ) us-ing Eq. (8) in Fig. 3, where the subscript denotes thestabilization mechanism and the superscript denotes thequantum states. As can be seen, σ ( s )3 agrees very wellwith σ GST around the left dip where the condensate is aSVS. This agreement also indicates that the density pro-file of a SVS is well approximated by a Gaussian func-tion. Around the right dip, the quantum state of the con-densate is more complicated than a pure coherent state. Therefore, σ ( c )3 only exhibits rough agreement with σ GST for large N where the quantum state is dominated by co-herent component. Still, the remaining discrepancy canbe explained as the density profile of the coherent stateis significantly deviated from a Gaussian function [39].As to σ ( c )5 / , it agree well with σ EGPE , which confirms thevariational calculation. One should note that, although σ ( c )5 / and σ ( c )3 have the same CANs, they define differentDANs.It is now clear that the W-shaped σ GST stems fromthe multiple quantum phases in the self-bound droplets.Remarkably, this feature can be vaguely seen in the ex-perimental data [5]. In fact, as shown in Fig. 1(b) and(c), after the completion of the left dip, the experimentaldata of the radial size starts to decrease again as one fur-ther increases N . Since radial size will eventually grow as √ N in the liquid phase, the drop after the first dip thensignals the existence of the second dip and, consequently,the W shape. As the coherent-state-based EGPE theoryonly predicts a V-shaped radial size curve, a clear obser-vation of the double dips may qualitatively exclude theEGPE explanation of the droplet states.More interestingly, variational analyses can even leadto quantitative criteria for quantum phases and stabiliza-tion mechanisms. To see this, let us examine N cri , theCAN extracted analytically from Eq. (8). At first sight,it is surprised that N cri is independent of ˜ g ν . This can beunderstood as follows. Because ε k and ε B have the samepower dependence on σ , the 2B attraction can be can-celed out exactly by the kinetic energy when N = N cri .Then, in the absence of stabilization force, a self-boundstate can be of arbitrary radial size. Once the stabiliza-tion force is turned on, in order to have a self-bound stateof the same atom number, the radial size of the droplethas to be infinite to make the repulsive interaction energyvanish. This universality of N cri is highly nontrivial as itallows our theory to be examined by experimental dataregardless of the accuracy of the fitted g . We point outthat the above analysis is not applicable in 3D geometrywhere a self-bound droplet always has a finite size [39].Still, the CAN depends on the quantum state of thecondensate through the RIP ˜ g . Then corresponding tothe pure squeezed and coherent states, we have distinctCANs N ( s )cri and N ( c )cri which satisfy N ( c )cri = 3 N ( s )cri . Theprecise measurements of CAN can therefore be used todistinguish the quantum state of the condensates. InFig. 2(a), we show the experimentally measured CANversus the reduced scattering length. From the right-most three data points, it seems that the quantum stateat the self-bound boundary is coherent state; while thequantum state for the other three data point is unclear.However, after carefully examining these data, we findthat high precision measurements are still needed to ac-curately determine the quantum phase [39]. In fact, toprecisely determine CAN, the experimental data shouldbe fitted using proper curves corresponding to appropri-ate stabilization mechanisms. In particular, the experi-mental data should be of V shape close to the self-boundboundary.Finally, we note that even the stabilization mechanismcan be identified using DAN N dip , a quantity dependingon ν . From Eq. (10), it is clear that, corresponding tothree σ ( c,s ) ν ( N ), there are three distinct DANs that sat-isfy N ( c )dip , = 3 N ( s )dip , and N ( c )dip , / = (9 / N ( s )dip , , where N ( s )dip , = 2 / ˜ g ( s )2 . Therefore, precise measurement of DANallows us to distinguish not only the quantum phase butalso the stabilization mechanism. Moreover, N dip is alsoindependent of ˜ g ν , which makes the determination of thestabilization mechanism independent of the fitting of g . Conclusion .—We have provided a detailed phase dia-gram for self-bound droplets of quasi-2D binary conden-sates beyond the coherent-state based mean field theory.Among all quantum phases, SVS and SCS are two macro-scopic quantum squeezed states that demand experimen-tal verification. More importantly, other than measuringAND or second-order correlation function, we have foundthree easily accessible experimental signatures: i) the ob- servation of the double dips on the σ ( N ) curve rule outthe explanations of the liquid droplets based on singlequantum state; ii) the precise measurements of N cri or N dip further allow us to distinguish the squeezed and thecoherent states; iii) we may even determine the stabiliza-tion mechanism from N dip . More remarkably, both N cri and N dip are independent of the strength of the stabi-lization force, which makes the comparison between thetheoretical and the experimental results immune to anyunknown parameters. Acknowledgement. —We thank enlightening discus-sions with H. Hu, X.-J. Liu, C. R. Cabrera, C. Navarrete-Benlloch, T. Pfau, and H. P. Buechler. This work wassupported by National Key Research and DevelopmentProgram of China (Grant No. 2017YFA0304501), by theNSFC (Grants No. 11974363, and No. 12047503), by theStrategic Priority Research Program of CAS (Grant No.XDB28000000), and by the Open Project of Shenzhen In-stitute of Quantum Science and Engineering (Grant No.SIQSE202006). [1] H. Kadau, M. Schmitt, M. Wenzel, C. Wink, T. Maier,I. Ferrier-Barbut, and T. Pfau, Nature , 194(2016), ISSN 1476-4687, URL https://doi.org/10.1038/nature16485 .[2] M. Schmitt, M. Wenzel, F. B¨ottcher, I. Ferrier-Barbut,and T. Pfau, Nature , 259 (2016), ISSN 1476-4687,URL https://doi.org/10.1038/nature20126 .[3] L. Chomaz, S. Baier, D. Petter, M. J. Mark, F. W¨achtler,L. Santos, and F. Ferlaino, Physical Review X ,041039 (2016), URL https://link.aps.org/doi/10.1103/PhysRevX.6.041039 .[4] F. B¨ottcher, M. Wenzel, J.-N. Schmidt, M. Guo, T. Lan-gen, I. Ferrier-Barbut, T. Pfau, R. Bomb´ın, J. S´anchez-Baena, J. Boronat, et al., Physical Review Research , 033088 (2019), URL https://link.aps.org/doi/10.1103/PhysRevResearch.1.033088 .[5] C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor,P. Thomas, P. Cheiney, and L. Tarruell, Science , 301 (2018), URL http://science.sciencemag.org/content/359/6373/301.abstract .[6] P. Cheiney, C. R. Cabrera, J. Sanz, B. Naylor,L. Tanzi, and L. Tarruell, Physical Review Letters , 135301 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.135301 .[7] G. Semeghini, G. Ferioli, L. Masi, C. Mazzinghi, L. Wol-swijk, F. Minardi, M. Modugno, G. Modugno, M. In-guscio, and M. Fattori, Physical Review Letters ,235301 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.235301 .[8] C. D’Errico, A. Burchianti, M. Prevedelli, L. Salasnich,F. Ancilotto, M. Modugno, F. Minardi, and C. Fort,Physical Review Research , 033155 (2019).[9] A. Burchianti, C. D’Errico, M. Prevedelli, L. Salasnich,F. Ancilotto, M. Modugno, F. Minardi, and C. Fort, Con-densed Matter , 21 (2020).[10] G. Natale, R. van Bijnen, A. Patscheider, D. Petter,M. Mark, L. Chomaz, and F. Ferlaino, Physical review letters , 050402 (2019).[11] L. Tanzi, S. Roccuzzo, E. Lucioni, F. Fam`a, A. Fioretti,C. Gabbanini, G. Modugno, A. Recati, and S. Stringari,Nature , 382 (2019).[12] G. Ferioli, G. Semeghini, L. Masi, G. Giusti,G. Modugno, M. Inguscio, A. Gallem´ı, A. Re-cati, and M. Fattori, Physical Review Letters ,090401 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.122.090401 .[13] F. B¨ottcher, J.-N. Schmidt, M. Wenzel, J. Hertkorn,M. Guo, T. Langen, and T. Pfau, Physical Review X , 011051 (2019), URL https://link.aps.org/doi/10.1103/PhysRevX.9.011051 .[14] L. Tanzi, E. Lucioni, F. Fam`a, J. Catani, A. Fioretti,C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno,Physical Review Letters , 130405 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.122.130405 .[15] L. Chomaz, D. Petter, P. Ilzh¨ofer, G. Natale, A. Traut-mann, C. Politi, G. Durastante, R. M. W. V. Bijnen,A. Patscheider, M. Sohmen, et al., Physical Review X , 021012 (2019), URL https://link.aps.org/doi/10.1103/PhysRevX.9.021012 .[16] M. Guo, F. B¨ottcher, J. Hertkorn, J.-N. Schmidt,M. Wenzel, H. P. B¨uchler, T. Langen, and T. Pfau, Na-ture , 386 (2019).[17] Z.-H. Luo, W. Pang, B. Liu, Y.-Y. Li, and B. A.Malomed, Frontiers of Physics , 32201 (2021),URL http://journal.hep.com.cn/fop/EN/abstract/article_28427.shtml .[18] F. B¨oettcher, J.-N. Schmidt, J. Hertkorn, K. Ng, S. Gra-ham, M. Guo, T. Langen, and T. Pfau, Reports onProgress in Physics , 012403 (2021), URL https://doi.org/10.1088/1361-6633/abc9ab .[19] M. Guo and T. Pfau, Frontiers of Physics , 32202(2021), URL https://link.springer.com/article/10.1007/s11467-020-1035-8 . [20] T. D. Lee, K. Huang, and C. N. Yang, Physical Review , 1135 (1957), URL https://link.aps.org/doi/10.1103/PhysRev.106.1135 .[21] R. Sch¨utzhold, M. Uhlmann, Y. Xu, and U. R. Fischer,Int. J. Mod. Phys. B , 3555 (2006), ISSN 0217-9792,URL https://doi.org/10.1142/S0217979206035631 .[22] A. R. P. Lima and A. Pelster, Physical Review A , 041604 (2011), URL https://link.aps.org/doi/10.1103/PhysRevA.84.041604 .[23] A. R. P. Lima and A. Pelster, Physical Review A , 063609 (2012), URL https://link.aps.org/doi/10.1103/PhysRevA.86.063609 .[24] D. S. Petrov, Physical Review Letters , 155302(2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.115.155302 .[25] H. Saito, Journal of the Physical Society of Japan ,053001 (2016), https://doi.org/10.7566/JPSJ.85.053001,URL https://doi.org/10.7566/JPSJ.85.053001 .[26] A. Macia, J. S´anchez-Baena, J. Boronat, and F. Maz-zanti, Phys. Rev. Lett. , 205301 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.205301 .[27] F. Cinti, A. Cappellaro, L. Salasnich, and T. Macr`ı, Phys.Rev. Lett. , 215302 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.215302 .[28] R. Bombin, J. Boronat, and F. Mazzanti, Phys. Rev.Lett. , 250402 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.250402 .[29] V. Cikojevi´c, L. V. c. v. Marki´c, G. E. Astrakharchik, andJ. Boronat, Phys. Rev. A , 023618 (2019), URL https://link.aps.org/doi/10.1103/PhysRevA.99.023618 .[30] M. Ota and G. E. Astrakharchik, SciPost Phys. , 20 (2020), URL https://scipost.org/10.21468/SciPostPhys.9.2.020 . [31] A. Boudjemˆaa and N. Guebli, Phys. Rev. A ,023302 (2020), URL https://link.aps.org/doi/10.1103/PhysRevA.102.023302 .[32] H. Hu, J. Wang, and X.-J. Liu, Phys. Rev. A ,043301 (2020), URL https://link.aps.org/doi/10.1103/PhysRevA.102.043301 .[33] H. Hu and X.-J. Liu, Phys. Rev. Lett. ,195302 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.125.195302 .[34] Q. Gu and L. Yin, Phys. Rev. B , 220503 (2020),URL https://link.aps.org/doi/10.1103/PhysRevB.102.220503 .[35] Y. Wang, L. Guo, S. Yi, and T. Shi, Physical ReviewResearch , 043074 (2020), URL https://doi.org/10.1103/PhysRevResearch.2.043074 .[36] T. Shi, J. Pan, and S. Yi, arXiv preprintarXiv:1909.02432 (2019), URL https://arxiv.org/abs/1909.02432 .[37] T. Shi, E. Demler, and J. I. Cirac, Annals ofPhysics , 245 (2018), ISSN 0003-4916, URL .[38] T. Guaita, L. Hackl, T. Shi, C. Hubig, E. Dem-ler, and J. I. Cirac, Physical Review B ,094529 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.094529 .[39] See Supplemental Material for the imaginary-time evo-lution in Gaussian-state theory, the atom-number distri-bution of a Gaussian state, the derivation of the virialrelation, the density profiles of the squeezed and coher-ent atoms, the boundaries of the SVS-to-SCS and theSCS-to-CS transitions, the variational analysis for radialsize, and the properties of 3D droplets. Supplemental Material
In the Supplemental Material, we present the derivation for the imaginary-time evolution in Gaussian-state theory(Sec. SM-I), the atom number distribution of a Gaussian state (Sec. SM-II), the virial relation (Sec. SM-III), thedensity profile (Sec. SM-IV), the boundary of the SVS-to-SCS transition (Sec. SM-V), the boundary of the SCS-to-CStransition (Sec. SM-VI), the variational analysis for radial size (Sec. SM-VII), and the main properties of the 3Ddroplets (Sec. SM-VIII).
SM-I. IMAGINARY-TIME EVOLUTION
Here we show the details of how to use the imaginary-time evolution to get the Gaussian mean-field ground state.To this end, we assume that the many-body wave function of a binary gas adopts the variational ansatz [1–4] | Ψ GS (cid:105) = e ˆΨ † ( r )Σ z ( r , r (cid:48) )Φ( r (cid:48) ) e i ˆΨ † ( r ) ξ ( r , r (cid:48) ) ˆΨ( r (cid:48) ) / | (cid:105) , (SM1)where ˆΨ( r ) ≡ (cid:18) ˆ ψ ( r )ˆ ψ † ( r ) (cid:19) with ˆ ψ ( r ) = (cid:18) ˆ ψ ↑ ( r )ˆ ψ ↓ ( r ) (cid:19) is the field operator expressed in the Nambu basis,Σ z ( r , r (cid:48) ) = (cid:18) I ⊗ δ ( r − r (cid:48) ) 00 − I ⊗ δ ( r − r (cid:48) ) (cid:19) (SM2)with I being a 2 × r ) = (cid:104) ˆΨ( r ) (cid:105) = (cid:18) φ ( c ) φ ( c ) ∗ (cid:19) with φ ( c ) = (cid:32) φ ( c ) ↑ φ ( c ) ↓ (cid:33) is the variationalparameter describing the coherent part of the condensates, and ξ ( r , r (cid:48) ) is the variational parameter that characterizesthe squeezed part of the condensate. We point out that the multiplication in the exponents in Eq. (SM1) representshort-hand notations which can be computed by first performing the matrix multiplications in the Nambu space andthen integrating over the repeated spatial coordinates. As an example, we haveˆΨ † ( r )Σ z ( r , r (cid:48) )Φ( r (cid:48) ) = (cid:88) α (cid:90) d r d r (cid:48) (cid:104) ˆ ψ † α ( r ) δ ( r − r (cid:48) ) φ ( c ) α ( r (cid:48) ) − ˆ ψ α ( r ) δ ( r − r (cid:48) ) φ ( c ) ∗ α ( r (cid:48) ) (cid:105) . (SM3)It should be noted that there exists a gauge redundancy in ξ ( r , r (cid:48) ) [1, 2]. Namely, different ξ ’s may lead to thesame Gaussian state. To remove this redundancy, we introduce the covariance matrixΓ( r , r (cid:48) ) ≡ (cid:18) Γ ( r , r (cid:48) ) Γ ( r , r (cid:48) )Γ ( r , r (cid:48) ) Γ ( r , r (cid:48) ) (cid:19) = (cid:104){ δ ˆΨ( r ) , δ ˆΨ † ( r (cid:48) ) }(cid:105) , where δ ˆΨ = ˆΨ − Φ is the fluctuation and { , } here denotes the anticommutator of the vectorial operator where indetail Γ ( r , r (cid:48) ) = (cid:32) (cid:104){ δ ˆ ψ †↑ ( r (cid:48) ) , δ ˆ ψ ↑ ( r ) }(cid:105) (cid:104){ δ ˆ ψ †↑ ( r (cid:48) ) , δ ˆ ψ ↓ ( r ) }(cid:105)(cid:104){ δ ˆ ψ †↓ ( r (cid:48) ) , δ ˆ ψ ↑ ( r ) }(cid:105) (cid:104){ δ ˆ ψ †↓ ( r (cid:48) ) , δ ˆ ψ ↓ ( r ) }(cid:105) (cid:33) (SM4)and Γ ( r , r (cid:48) ) = (cid:18) (cid:104){ δ ˆ ψ ↑ ( r (cid:48) ) , δ ˆ ψ ↑ ( r ) }(cid:105) (cid:104){ δ ˆ ψ ↑ ( r (cid:48) ) , δ ˆ ψ ↓ ( r ) }(cid:105)(cid:104){ δ ˆ ψ ↓ ( r (cid:48) ) , δ ˆ ψ ↑ ( r ) }(cid:105) (cid:104){ δ ˆ ψ ↓ ( r (cid:48) ) , δ ˆ ψ ↓ ( r ) }(cid:105) (cid:19) (SM5)and Γ = Γ † , Γ = Γ T . It can be shown that Γ and ξ are connected through Γ = SS † , where S = e i Σ z ξ is asymplectic matrix satisfying S Σ z S † = Σ z . So it is convenient to take the elements of φ ( c ) ( r ) and Γ( r , r (cid:48) ) as variationalparameters.To proceed further, we substitute the expansion ψ ( r ) = φ ( c ) ( r ) + δψ ( r ) into the Hamiltonian H defined in maintext, which leads to single-particle Hamiltonian H = (cid:88) α = ↑ , ↓ (cid:90) d r (cid:16) φ ( c ) α ( r ) † h α φ ( c ) α ( r ) + δ ˆ ψ † α ( r ) h α φ ( c ) ( r ) + φ ( c ) α ( r ) † h α δ ˆ ψ α ( r ) + δ ˆ ψ † α ( r ) h α δ ˆ ψ α ( r ) (cid:17) (SM6)two-body interaction Hamiltonian H B = (cid:88) α,β = ↑ , ↓ g αβ (cid:90) d r (cid:104) | φ ( c ) α ( r ) | | φ ( c ) β ( r ) | + (cid:16) δ ˆ ψ † α ( r ) φ ( c ) α ( r ) | φ ( c ) β ( r ) | + δ ˆ ψ † α ( r ) δ ˆ ψ † β ( r ) φ ( c ) β ( r ) φ ( c ) α ( r )+ 2 δ ˆ ψ † α ( r ) δ ˆ ψ † β ( r ) δ ˆ ψ β ( r ) φ ( c ) α ( r ) + H.c. (cid:17) + 2 δ ˆ ψ † α ( r ) δ ˆ ψ α ( r ) | φ ( c ) β ( r ) | + 2 δ ˆ ψ † α ( r ) δ ˆ ψ β ( r ) φ ( c ) ∗ β ( r ) φ ( c ) α ( r )+ δ ˆ ψ † α ( r ) δ ˆ ψ † β ( r ) δ ˆ ψ β ( r ) δ ˆ ψ α ( r ) (cid:105) (SM7)and three-body interaction Hamiltonian H B H B = g (cid:88) α,β,γ = ↑ , ↓ (cid:90) d r (cid:104) | φ ( c ) α ( r ) | | φ ( c ) β ( r ) | | φ ( c ) γ ( r ) | + (cid:16) δ ˆ ψ † α ( r ) φ ( c ) α ( r ) | φ ( c ) β ( r ) | | φ ( c ) γ ( r ) | + 3 δ ˆ ψ † α ( r ) δ ˆ ψ † β ( r ) φ ( c ) β ( r ) φ ( c ) α ( r ) | φ ( c ) γ ( r ) | + 3 δ ˆ ψ † α ( r ) δ ˆ ψ β ( r ) φ ( c ) ∗ β ( r ) φ ( c ) α ( r ) | φ ( c ) γ ( r ) | + δ ˆ ψ † α ( r ) δ ˆ ψ † β ( r ) δ ˆ ψ † γ ( r ) φ ( c ) γ ( r ) φ ( c ) β ( r ) φ ( c ) α ( r ) + 6 δ ˆ ψ † α ( r ) δ ˆ ψ † β ( r ) δ ˆ ψ α ( r ) φ ( c ) β ( r ) | φ ( c ) γ ( r ) | + 3 δ ˆ ψ † α ( r ) δ ˆ ψ † β ( r ) δ ˆ ψ γ ( r ) φ ( c ) ∗ γ ( r ) φ ( c ) β ( r ) φ ( c ) α ( r ) + 3 δ ˆ ψ † α ( r ) δ ˆ ψ † β ( r ) δ ˆ ψ † γ ( r ) δ ˆ ψ γ ( r ) φ ( c ) β ( r ) φ ( c ) α ( r )+ 3 δ ˆ ψ † α ( r ) δ ˆ ψ † β ( r ) δ ˆ ψ † γ ( r ) δ ˆ ψ γ ( r ) δ ˆ ψ β ( r ) φ ( c ) α ( r ) + H.c. (cid:17) + 3 δ ˆ ψ † α ( r ) δ ˆ ψ α ( r ) | φ ( c ) β ( r ) | | φ ( c ) γ ( r ) | + 3 δ ˆ ψ † α ( r ) δ ˆ ψ † β ( r ) δ ˆ ψ β ( r ) δ ˆ ψ α ( r ) | φ ( c ) γ ( r ) | + 6 δ ˆ ψ † α ( r ) δ ˆ ψ † β ( r ) δ ˆ ψ γ ( r ) δ ˆ ψ β ( r ) φ ( c ) ∗ γ ( r ) φ ( c ) α ( r )+ δ ˆ ψ † α ( r ) δ ˆ ψ † β ( r ) δ ˆ ψ † γ ( r ) δ ˆ ψ γ ( r ) δ ˆ ψ β ( r ) δ ˆ ψ α ( r ) (cid:105) (SM8)After applying the Wick’s theorem with respect to the Gaussian state | Ψ GS (cid:105) , we obtain the mean-field Hamiltonian H MF = E + ( δ ˆ ψ † η + η ∗ δ ˆ ψ ) + 12 : δ ˆΨ † H δ ˆΨ : , (SM9)where : ˆ O : denotes the normal-ordered operator defined with respect to the Gaussian state, E ≡ (cid:104) Ψ GS | H | Ψ GS (cid:105) isthe expectation value of the Hamiltonian and is composed of the kinetic energy E k = (cid:88) α = ↑ , ↓ (cid:90) d r (cid:16) φ ( c ) ∗ α ( r ) h α φ ( c ) α ( r ) + (cid:104) δ ˆ ψ † α ( r ) h α δ ˆ ψ α ( r ) (cid:105) (cid:17) , (SM10)the two-body interaction energy E B = (cid:88) α,β = ↑ , ↓ g αβ (cid:90) d r (cid:104) | φ ( c ) α ( r ) | | φ ( c ) β ( r ) | + 2Re (cid:16) (cid:104) δ ˆ ψ α ( r ) δ ˆ ψ β ( r ) (cid:105) φ (c) ∗ β ( r ) φ (c) ∗ α ( r ) (cid:17) + |(cid:104) δ ˆ ψ α ( r ) δ ˆ ψ β ( r ) (cid:105)| + (cid:16) | φ ( c ) α ( r ) | + (cid:104) δ ˆ ψ † α ( r ) δ ˆ ψ α ( r ) (cid:105) (cid:17) (cid:104) δ ˆ ψ † β ( r ) δ ˆ ψ β ( r ) (cid:105) + (cid:16) φ ( c ) ∗ α ( r ) φ ( c ) β ( r ) + (cid:104) δ ˆ ψ † α ( r ) δ ˆ ψ β ( r ) (cid:105) (cid:17) (cid:104) δ ˆ ψ † β ( r ) δ ˆ ψ α ( r ) (cid:105) (cid:105) , (SM11)and the tree-body interaction energy E B = g (cid:88) α,β,γ = ↑ , ↓ (cid:90) d r (cid:104) | φ ( c ) α ( r ) | | φ ( c ) β ( r ) | | φ ( c ) γ ( r ) | + 6Re (cid:16) (cid:104) δ ˆ ψ † α ( r ) δ ˆ ψ † β ( r ) (cid:105) φ (c) β ( r ) φ (c) α ( r ) | φ (c) γ ( r ) | (cid:17) + 6 (cid:104) δ ˆ ψ † α ( r ) δ ˆ ψ β ( r ) (cid:105) φ ( c ) ∗ β ( r ) φ ( c ) α ( r ) | φ ( c ) γ ( r ) | + 3 (cid:104) δ ˆ ψ † α ( r ) δ ˆ ψ α ( r ) (cid:105)| φ ( c ) β ( r ) | | φ ( c ) γ ( r ) | + 2 (cid:104) δ ˆ ψ α ( r ) δ ˆ ψ β ( r ) (cid:105) (cid:16) (cid:104) δ ˆ ψ † α ( r ) δ ˆ ψ † β ( r ) (cid:105)(cid:104) δ ˆ ψ † γ ( r ) δ ˆ ψ γ ( r ) (cid:105) + 2 (cid:104) δ ˆ ψ † α ( r ) δ ˆ ψ † γ ( r ) (cid:105)(cid:104) δ ˆ ψ † β ( r ) δ ˆ ψ γ ( r ) (cid:105) (cid:17) + 6Re (cid:16) φ (c) β ( r ) φ (c) α ( r ) (cid:16) (cid:104) δ ˆ ψ † α ( r ) δ ˆ ψ † β ( r ) (cid:105)(cid:104) δ ˆ ψ † γ ( r ) δ ˆ ψ γ ( r ) (cid:105) + 2 (cid:104) δ ˆ ψ † α ( r ) δ ˆ ψ † γ ( r ) (cid:105)(cid:104) δ ˆ ψ † β ( r ) δ ˆ ψ γ ( r ) (cid:105) (cid:17)(cid:17) + (cid:16) | φ ( c ) α ( r ) | + (cid:104) δ ˆ ψ † α ( r ) δ ˆ ψ α ( r ) (cid:105) (cid:17) (cid:16) |(cid:104) δ ˆ ψ β ( r ) δ ˆ ψ γ ( r ) (cid:105)| + (cid:104) δ ˆ ψ † γ ( r ) δ ˆ ψ γ ( r ) (cid:105)(cid:104) δ ˆ ψ † β ( r ) δ ˆ ψ β ( r ) (cid:105) + |(cid:104) δ ˆ ψ † γ ( r ) δ ˆ ψ β ( r ) (cid:105)| (cid:17) + 2 (cid:16) φ ( c ) ∗ β ( r ) φ ( c ) α ( r ) + (cid:104) δ ˆ ψ † β ( r ) δ ˆ ψ α ( r ) (cid:105) (cid:17) (cid:16) (cid:104) δ ˆ ψ † α ( r ) δ ˆ ψ † γ ( r ) (cid:105)(cid:104) δ ˆ ψ γ ( r ) δ ˆ ψ β ( r ) (cid:105) + (cid:104) δ ˆ ψ † α ( r ) δ ˆ ψ γ ( r ) (cid:105)(cid:104) δ ˆ ψ † γ ( r ) δ ˆ ψ β ( r ) (cid:105) + (cid:104) δ ˆ ψ † α ( r ) δ ˆ ψ β ( r ) (cid:105)(cid:104) δ ˆ ψ † γ ( r ) δ ˆ ψ γ ( r ) (cid:105) (cid:17)(cid:105) . (SM12)The linear driving term η = (cid:18) η ↑ η ↓ (cid:19) is composed of η α = h α φ ( c ) α + (cid:88) β = ↑↓ g αβ (cid:16)(cid:16) | φ ( c ) β | + (cid:104) δ ˆ ψ † β δ ˆ ψ β (cid:105) (cid:17) φ ( c ) α + (cid:104) δ ˆ ψ α δ ˆ ψ β (cid:105) φ ( c ) ∗ β + (cid:104) δ ˆ ψ † β δ ˆ ψ α (cid:105) φ ( c ) β (cid:17) + g (cid:88) β,γ = ↑ , ↓ (cid:104)(cid:16) | φ ( c ) β | | φ ( c ) γ | + 2Re (cid:16) (cid:104) δ ˆ ψ † β δ ˆ ψ † γ (cid:105) φ (c) γ φ (c) β (cid:17) + 2 (cid:104) δ ˆ ψ † γ δ ˆ ψ β (cid:105) φ (c) ∗ β φ (c) γ + 2 (cid:104) δ ˆ ψ † β δ ˆ ψ β (cid:105)| φ (c) γ | + |(cid:104) δ ˆ ψ β δ ˆ ψ γ (cid:105)| + (cid:104) δ ˆ ψ † γ δ ˆ ψ γ (cid:105)(cid:104) δ ˆ ψ † β δ ˆ ψ β (cid:105) + |(cid:104) δ ˆ ψ † γ δ ˆ ψ β (cid:105)| (cid:17) φ ( c ) α + 2 (cid:16) (cid:104) δ ˆ ψ α δ ˆ ψ β (cid:105)| φ ( c ) γ | + (cid:104) δ ˆ ψ α δ ˆ ψ β (cid:105)(cid:104) δ ˆ ψ † γ δ ˆ ψ γ (cid:105) + (cid:104) δ ˆ ψ † γ δ ˆ ψ β (cid:105)(cid:104) δ ˆ ψ γ δ ˆ ψ α (cid:105) + (cid:104) δ ˆ ψ † γ δ ˆ ψ α (cid:105)(cid:104) δ ˆ ψ γ δ ˆ ψ β (cid:105) (cid:17) φ ( c ) ∗ β + 2 (cid:16) (cid:104) δ ˆ ψ † β δ ˆ ψ α (cid:105)| φ ( c ) γ | + (cid:104) δ ˆ ψ † β δ ˆ ψ † γ (cid:105)(cid:104) δ ˆ ψ γ δ ˆ ψ α (cid:105) + (cid:104) δ ˆ ψ † β δ ˆ ψ γ (cid:105)(cid:104) δ ˆ ψ † γ δ ˆ ψ α (cid:105) + (cid:104) δ ˆ ψ † β δ ˆ ψ α (cid:105)(cid:104) δ ˆ ψ † γ δ ˆ ψ γ (cid:105) (cid:17) φ ( c ) β (cid:105) . (SM13)The matrix elements of H = (cid:18) E ∆∆ † E ∗ (cid:19) are defined as E = (cid:18) E ↑↑ E ↑↓ E ↓↑ E ↓↓ (cid:19) and ∆ = (cid:18) ∆ ↑↑ ∆ ↑↓ ∆ ↓↑ ∆ ↓↓ (cid:19) , where E αβ = (cid:88) γ = ↑ , ↓ (cid:104) h α + g αγ (cid:16) | φ ( c ) γ | + (cid:104) δ ˆ ψ † γ δ ˆ ψ γ (cid:105) (cid:17)(cid:105) δ αβ + g αβ (cid:16) φ ( c ) α φ ( c ) ∗ β + (cid:104) δ ˆ ψ † β δ ˆ ψ α (cid:105) (cid:17) + g (cid:88) γ = ↑ , ↓ (cid:88) γ (cid:48) = ↑ , ↓ (cid:18) | φ ( c ) γ | | φ ( c ) γ (cid:48) | + | φ ( c ) γ | (cid:104) δ ˆ ψ † γ (cid:48) δ ˆ ψ γ (cid:48) (cid:105) + φ ( c ) ∗ γ (cid:48) φ ( c ) γ (cid:104) δ ˆ ψ † γ δ ˆ ψ γ (cid:48) (cid:105) + Re (cid:16) φ (c) ∗ γ φ (c) ∗ γ (cid:48) (cid:104) δ ˆ ψ γ (cid:48) δ ˆ ψ γ (cid:105) (cid:17) + 12 (cid:16) (cid:104) δ ˆ ψ † γ δ ˆ ψ γ (cid:105)(cid:104) δ ˆ ψ † γ (cid:48) δ ˆ ψ γ (cid:48) (cid:105) + (cid:104) δ ˆ ψ † γ δ ˆ ψ γ (cid:48) (cid:105)(cid:104) δ ˆ ψ † γ (cid:48) δ ˆ ψ γ (cid:105) + (cid:104) δ ˆ ψ † γ δ ˆ ψ † γ (cid:48) (cid:105)(cid:104) δ ˆ ψ γ (cid:48) δ ˆ ψ γ (cid:105) (cid:17)(cid:17) δ αβ + (cid:16) | φ ( c ) γ | + (cid:104) δ ˆ ψ † γ δ ˆ ψ γ (cid:105) (cid:17) (cid:16) φ ( c ) α φ ( c ) ∗ β + (cid:104) δ ˆ ψ † β δ ˆ ψ α (cid:105) (cid:17) + φ ( c ) γ φ ( c ) ∗ β (cid:104) δ ˆ ψ † γ δ ˆ ψ α (cid:105) + φ ( c ) α φ ( c ) ∗ γ (cid:104) δ ˆ ψ † β δ ˆ ψ γ (cid:105) + φ ( c ) α φ ( c ) γ (cid:104) δ ˆ ψ † γ δ ˆ ψ † β (cid:105) + φ ( c ) ∗ γ φ ( c ) ∗ β (cid:104) δ ˆ ψ γ δ ˆ ψ α (cid:105) + (cid:104) δ ˆ ψ † γ δ ˆ ψ α (cid:105)(cid:104) δ ˆ ψ † β δ ˆ ψ γ (cid:105) + (cid:104) δ ˆ ψ α δ ˆ ψ γ (cid:105)(cid:104) δ ˆ ψ † γ δ ˆ ψ † β (cid:105) (cid:35) (SM14)and ∆ αβ = g αβ (cid:16) φ ( c ) α φ ( c ) β + (cid:104) δ ˆ ψ α δ ˆ ψ β (cid:105) (cid:17) + g (cid:88) γ = ↑ , ↓ (cid:104)(cid:16) | φ ( c ) γ | + (cid:104) δ ˆ ψ † γ δ ˆ ψ γ (cid:105) (cid:17) (cid:16) φ ( c ) α φ ( c ) β + (cid:104) δ ˆ ψ α δ ˆ ψ β (cid:105) (cid:17) + φ ( c ) β φ ( c ) γ (cid:104) δ ˆ ψ † γ δ ˆ ψ α (cid:105) + φ ( c ) α φ ( c ) γ (cid:104) δ ˆ ψ † γ δ ˆ ψ β (cid:105) + (cid:104) δ ˆ ψ α δ ˆ ψ γ (cid:105) φ ( c ) ∗ γ φ ( c ) β + (cid:104) δ ˆ ψ β δ ˆ ψ γ (cid:105) φ ( c ) ∗ γ φ ( c ) α + (cid:104) δ ˆ ψ α δ ˆ ψ γ (cid:105)(cid:104) δ ˆ ψ † γ δ ˆ ψ β (cid:105) + (cid:104) δ ˆ ψ β δ ˆ ψ γ (cid:105)(cid:104) δ ˆ ψ † γ δ ˆ ψ α (cid:105) (cid:105) . (SM15)Based on the GST, the ground-state solutions, φ ( c ) ( r ) and Γ( r , r (cid:48) ), can be obtained by numerically evolving theimaginary-time equations of motion [1, 2], ∂ τ Φ = − Γ (cid:18) ηη ∗ (cid:19) , (SM16a) ∂ τ Γ = Σ z H Σ z − Γ H Γ , (SM16b)which converge when the imaginary time τ is sufficiently large. SM-II. ATOM-NUMBER DISTRIBUTION
Following the approach used in dipolar droplets [4], here we derive the atom-number distribution of a Gaussianstate describing a binary condensate. To this end, we expand the coherent mode ˆ c in terms of the squeezed modes asˆ c = ∞ (cid:88) i =1 α i ˆ s i + α ⊥ ˆ c ⊥ (SM17)0 FIG. SM1: (color online). Atom number distributions P ( n ) for the states with mean atom numbers (from left to right) N = 6 . × , 3 × , and 10 . The y axis of the blue line (red lines) is on the left (right). The inset shows the coherentfraction as a function of mean atom number where the markers denote the mean atom numbers used to plot P ( n )’s. Thereduced scattering length is δa = − . a . where α i = (cid:82) d r (cid:104) ¯ φ ( s ) i ( r ) (cid:105) † ¯ φ ( c ) ( r ), α ⊥ = 1 − (cid:80) i | α i | , and ˆ c ⊥ represents the mode that is perpendicular to all ˆ s i . TheGaussian state wave function can now be expressed as | Ψ GS (cid:105) = e √ N ( c ) α ⊥ ( ˆ c †⊥ − ˆ c ⊥ ) ∞ (cid:89) i =1 e √ N ( c ) ( α i ˆ s † i − α ∗ i ˆ s i ) e ξ i ( ˆ s † i − ˆ s i ) | (cid:105) . (SM18)The AND for mode ˆ c ⊥ is p ⊥ ( (cid:96) ) = e − N ( c ) | α ⊥ | (cid:16) N ( c ) | α ⊥ | (cid:17) (cid:96) /(cid:96) !and that for mode ˆ s i is [5] p i ( (cid:96) ) = (cid:0) tanh ξ i (cid:1) (cid:96) (cid:96) ! cosh ξ i (cid:12)(cid:12)(cid:12) H (cid:96) (cid:104) γ i (sinh 2 ξ i ) − / (cid:105)(cid:12)(cid:12)(cid:12) exp (cid:20) − N ( c ) | α i | − N ( c ) (cid:0) α ∗ i + α i (cid:1) tanh ξ i (cid:21) , (SM19)where γ i = √ N ( c ) α i cosh ξ i + √ N ( c ) α ∗ i sinh ξ i and H (cid:96) ( x ) are Hermite polynomials. The AND of | Ψ GS (cid:105) can then beexpressed into a recursive form as follows. Let P i ( n ) denote the AND of the state containing first i squeezed modes,the AND containing first i + 1 squeezed modes can then be expressed as P i +1 ( n ) = n (cid:88) (cid:96) =0 P i ( n − (cid:96) ) p i +1 ( (cid:96) ) (SM20)where P ( n ) ≡ p ⊥ ( n ). Applying this equation successively will eventually leads to the total AND P ( n ).Figure SM1 shows three typical ANDs corresponding to states in SVS, SCS, and CS phases. Due to the squeezing, P ( n ) is asymmetric with a long tail at the large n side, which was observed in dipolar droplet experiments [6, 7].However, as f c increases, P ( n ) becomes more symmetric. Other features of the AND include that P ( n ) has a peakat n ≈ N ( c ) and the mean value of atom number is (cid:88) n nP ( n ) = N ( c ) + N ( s ) = N. SM-III. VIRIAL RELATION
Here we derive a virial relation among different contributions to the total energy [8]. To this end, we shall make useof the fact that different spin component share same normalized density profile. Hence we can use ¯ φ ( c ) ( r ) to represent1 FIG. SM2: (color online). ε k , ε B , and ε B as functions of N for δa = − . a . Energy is in units of (cid:126) / ( ML ) with L = 6 µ m.Vertical dotted lines mark the locations of two phase transitions. the normalized mode function for the coherent component and for i th squeezed mode the normalized mode functionsare ¯ φ ( s ) i ( r ). If total energy E = E k + E B + E B is stationary for any variation of ¯ φ ( c ) ( r ) and ¯ φ ( s ) i ( r ) around theground state solution, we can choose scaling transformations of the form φ ( x, y ) → (1 + λ ) / φ [(1 + λ ) x, y ] and insertkinetic energy, 2B and 3B interaction energies, we have δE k =[(1 + λ ) − E xk δE B = λE B δE B =[(1 + λ ) − E B , (SM21)where E xk = − M (cid:90) d r (cid:32) N ( c ) ¯ φ ( c ) ∗ ( r ) ∂ ∂x ¯ φ ( c ) ( r ) + (cid:88) i N ( s ) i ¯ φ ( s ) ∗ i ( r ) ∂ ∂x ¯ φ ( s ) i ( r ) (cid:33) . (SM22)By imposing the energy variation δE = δE k + δE B + δE B to vanish at first order in λ we can find2 E xk + E B + 2 E B = 0 . (SM23)Analogous expressions can be obtained by impose same scaling transform in y directions. By summing over the x andthe y directions, we can find the virial relation E k + E B + 2 E B = 0 (SM24)for our two-dimensional system.In Fig. SM2, we plot the N dependence of the energies, where ε k = E k /N , ε B = E B /N , and ε B = E B /N . Ascan be seen, unlike ε B and ε B which are monotonic functions of N , the variation of ε k is negatively correlated tothat of the radial size σ GST because ε k ∝ /σ . Furthermore, the phase transitions can also be visualized from thebehavior of the energies. In particular, it appears that the SVS-to-SCS transition occurs when ε k = ε B , which willbe proven analytically in Sec. SM-V. SM-IV. DENSITY PROFILE
Here we present the normalized density profiles of the coherent and the squeezed atoms, i.e., ¯ n ( c ) ( ρ ) = | ¯ φ ( c ) ( ρ ) | and¯ n ( s ) ( ρ ) = n ( s ) ( ρ ) /N ( s ) , here ρ = (cid:112) x + y and we have used the fact that both ¯ n ( c ) and ¯ n ( s ) are axially symmetric.For convenience, we introduce the width of the coherent/squeezed atoms w ( c,s ) = (cid:20)(cid:90) dxdyρ ¯ n ( c.s ) ( x, y ) (cid:21) / (a) 0 5 10 15 20036912 10 (b) FIG. SM3: (color online). Normalized density profiles of the coherent (a) and the squeezed atoms (b) with δa = − . a and N = 9 . × ( (cid:35) ), 2 . × ( ♦ ), 3 . × ( (cid:52) ), 4 . × ( (cid:53) ), and 9 × ( (cid:3) ). Insets in (a) and (b) show the N dependenceof w ( c ) and w ( s ) , respectively. Vertical dotted lines mark the locations of two phase transitions. Symbols on w ( c,s ) mark the N values used to plot the corresponding density profiles. to characterize the geometry of the states. Clearly, w ( c,s ) are better defined than the radial size when ¯ n ( c,s ) ( ρ ) deviatesignificantly from Gaussian function.In Fig. SM3, we plot ¯ n ( c,s ) ( ρ ) for various N ’s in SCS and CS phases. Immediately to the right of the SVS-to-SCStransition, e.g. N = 9 . × , it is found that ¯ n ( c ) is identical to ¯ n ( s ) . This observation is understandable as thecoherent atoms originates from the three-body repulsion which is simply proportional to ¯ n ( s ) when N ( c ) is essentiallyzero. When N is further increased, ¯ n ( c ) ( ρ ) changes in according with the size of the coherent part as shown inFig. SM3(a). Particularly, in the large N limit, the top of ¯ n ( c ) ( ρ ) becomes flatted and significantly deviates from aGaussian function, which is in consistent to the liquid nature of the condensate. As to the profiles of the squeezedatoms, ¯ n ( s ) ( ρ = 0) continuously decreases with N such that ¯ n ( s ) ( ρ ) is no longer a monotonic decreasing function of ρ in the CS phase, which further demonstrates the squeezing here is quantum depletion due to the three-body repulsion. SM-V. SVS-TO-SCS TRANSITION
Here we derive the equation describing the boundary between the SVS and the SCS phases. To this end, we notethat the 2B and 3B interaction energies for a SCS are, respectively, E B = 4 π (cid:126) δaM √ a ↑↑ a ↓↓ ( √ a ↑↑ + √ a ↓↓ ) (cid:90) d r (cid:104) N ( c )2 | ¯ φ ( c ) ( r ) | + 6 N ( c ) N ( s ) | ¯ φ ( c ) ( r ) | | ¯ φ ( s ) ( r ) | + 3 N ( s )2 | ¯ φ ( s ) ( r ) | (cid:105) ,E B = g (cid:90) d r (cid:104) N ( c )3 | ¯ φ ( c ) ( r ) | + 15 N ( c )2 N ( s ) | ¯ φ ( c ) ( r ) | | ¯ φ ( s ) ( r ) | + 45 N ( c ) N ( s )2 | ¯ φ ( c ) ( r ) | | ¯ φ ( s ) ( r ) | +15 N ( s )3 | ¯ φ ( s ) ( r ) | (cid:105) . (SM25)In the vicinity of the SVS-to-SCS transition, the squeezed and the coherent modes have the same density profile (seethe main text for details), i.e., ¯ n ( r ) ≡ ¯ n ( s ) ( r ) = ¯ n ( c ) ( r ). Then, the interaction energies reduce to E B = ω B (cid:16) N ( c )2 + 6 N ( c ) N ( s ) + 3 N ( s )2 (cid:17) ,E B = ω B (cid:16) N ( c )3 + 15 N ( c )2 N ( s ) + 45 N ( c ) N ( s )2 + 15 N ( s )3 (cid:17) , (SM26)where ω = 4 π (cid:126) δaM √ a ↑↑ a ↓↓ ( √ a ↑↑ + √ a ↓↓ ) (cid:90) d r ¯ n ( r ) and ω = g (cid:90) d r ¯ n ( r ) . (SM27)To find the phase boundary, we consider a pure squeezed state with all N atoms occupying the squeezed mode.The total energy is then E ( N ) = ε k N + 3 ω B N + 15 ω B N . (SM28)3Now, we add δN atom to the condensate. If all newly added atoms go to the squeezed mode, we have E (SVS)0 ( N + δN ) = ε k ( N + δN ) + 3 ω B ( N + δN ) + 15 ω B ( N + δN ) ; (SM29)while if all newly added atoms go to the coherent mode, we have E (SCS)0 ( N + δN ) = ε k ( N + δN ) + ω B (3 N + 6 N δN + δN ) + ω B (15 N + 45 N δN + 15 N δN + δN ) . (SM30)If the transition happens at N , we must have E (SCS)0 ( N + δN ) < E (SVS)0 ( N + δN ), which, to the second order in δN ,leads to − ω B < ω B N. This result says that the atom number at the third order phase transition, N ∗ , satisfies − ω B = 15 ω B N ∗ . Makinguse of Eq. (SM28), we prove − E B = 3 E B and, subsequently, E k = E B on the boundary of the SVS-to-SCS transition. Finally, combined with the virial relation and Eq. ( ?? ), we obtain N ∗ = 3˜ g ( s )2 . (SM31) SM-VI. SCS-TO-CS TRANSITION
In this section, we derived a criterion for the first-order SCS-to-CS transition based on the stability of Bogoliubovexcitation in the the presence of 3B interactions. Specifically, in the CS phase, we may neglect the contributions fromterms (cid:104) δ ˆ ψ † ( r (cid:48) ) δ ˆ ψ ( r ) (cid:105) and (cid:104) δ ˆ ψ ( r (cid:48) ) δ ˆ ψ ( r ) (cid:105) such that Eqs. (SM13)-(SM15) reduce to η α = h (cid:48) α + (cid:88) β g (cid:48) αβ | φ ( c ) β ( x, y ) | + g (cid:48) (cid:88) β,γ | φ ( c ) β ( x, y ) | | φ ( c ) γ ( x, y ) | φ ( c ) α ( x, y ) E αβ = (cid:32) h (cid:48) α + (cid:88) γ g (cid:48) αγ | φ ( c ) γ ( x, y ) | (cid:33) δ αβ + g (cid:48) αβ φ ( c ) α ( x, y ) φ ( c ) ∗ β ( x, y )+ g (cid:48) (cid:88) γ (cid:88) γ (cid:48) | φ ( c ) γ ( x, y ) | | φ ( c ) γ (cid:48) ( x, y ) | δ αβ + | φ ( c ) γ ( x, y ) | φ ( c ) α ( x, y ) φ ( c ) ∗ β ( x, y ) ∆ αβ = g (cid:48) αβ φ ( c ) α ( x, y ) φ ( c ) β ( x, y ) + g (cid:48) (cid:88) γ | φ ( c ) γ ( x, y ) | φ ( c ) α ( x, y ) φ ( c ) β ( x, y ) , (SM32)where h (cid:48) α = − (cid:126) ( ∂ x + ∂ y ) / (2 M ) − µ α is the single particle term and g (cid:48) αβ = g αβ / ( √ πa z ) and g (cid:48) = g / ( √ πa z ) arethe interaction parameters of quasi-2D system. Meanwhile, the covariance matrix Γ reduces to a unit matrix. As aresult, the steady-state condition ∂ τ Φ = 0 leads to the Gross-Pitaevskii equations h (cid:48) α + (cid:88) β g (cid:48) αβ | φ ( c ) β ( x, y ) | + g (cid:48) (cid:88) β,γ | φ ( c ) β ( x, y ) | | φ ( c ) γ ( x, y ) | φ ( c ) α ( x, y ) = 0 . (SM33)For simplicity, we consider a homogeneous gas. The chemical potentials can then be solved to be µ α = (cid:88) β g (cid:48) αβ n β + g (cid:48) (cid:88) β,γ n β n γ , (SM34)where n α = | φ ( c ) α | is the density of α th component.4 FIG. SM4: (color online). Peak condensate density n p as a function of atom number N for δa = − . a . The dashed linemarks the critical condensate density n ∗ . The vertical dotted line denotes the boundary between the SCS and the CS phases. In the momentum space, the mean-field Hamiltonian becomes H MF = E + (cid:88) k
12 : δ ˆΨ † k H k δ ˆΨ k : (SM35)where δ ˆΨ k is the Fourier transform of the fluctuation operator δ ˆΨ and H k = (cid:18) E k ∆ k ∆ † k E Tk (cid:19) with E k = (cid:18) (cid:15) k + ( g (cid:48)↑↑ + g (cid:48) (cid:80) β n β ) n ↑ ( g (cid:48)↑↓ + g (cid:48) (cid:80) β n β ) √ n ↑ n ↓ ( g (cid:48)↑↓ + g (cid:48) (cid:80) β n β ) √ n ↑ n ↓ (cid:15) k + ( g (cid:48)↓↓ + g (cid:48) (cid:80) β n β ) n ↓ (cid:19) (SM36)and ∆ k = (cid:18) ( g (cid:48)↑↑ + g (cid:48) (cid:80) β n β ) n ↑ ( g (cid:48)↑↓ + g (cid:48) (cid:80) β n β ) √ n ↑ n ↓ ( g (cid:48)↑↓ + g (cid:48) (cid:80) β n β ) √ n ↑ n ↓ ( g (cid:48)↓ , ↓ + g (cid:48) (cid:80) β n β ) n ↓ (cid:19) . (SM37)Here (cid:15) k = (cid:126) k / (2 M ). After diagonalizing H k , we obtain the Bogoliubov excitation spectra ξ ± ( k ) = 12 (cid:104) ξ ↑ + ξ ↓ ± (cid:113) ( ξ ↑ − ξ ↓ ) + 16( g (cid:48)↑↓ + g (cid:48) n t ) n ↑ n ↓ ε k (cid:105) (SM38)where n t = n ↑ + n ↓ is the total condensate density, and ξ α = (cid:15) k [ (cid:15) k + 2( g (cid:48) αα + g (cid:48) n t ) n α ] (SM39)are the familiar quasiparticle energies in a single-component condensate of spin- α atoms. Apparently, the CS phasebecomes unstable when quasiparticle energy in the lower branch is imaginary, i.e.,( ξ ↑ + ξ ↓ ) < ( ξ ↑ − ξ ↓ ) + 16( g (cid:48)↑↓ + g (cid:48) n t ) n ↑ n ↓ (cid:15) k . (SM40)After some straightforward calculations, we find n t < n ∗ ≡ (cid:114) π a z g π (cid:126) M a ↑↓ − a ↑↑ a ↓↓ ( a ↑↑ + a ↓↓ − a ↑↓ ) , (SM41)which can be interpreted as the criterion for the CS-to-SCS transition. In Fig. SM4, we plot the numerically obtainedpeak condensate density n p as a function of N . As can be seen, the boundary for the SCS-to-CS transition determinedwith the criterion n p = n ∗ is in good agreement with numerical result.5 (a) (b) (c) (d) FIG. SM5: (color online). Radial size σ versus atom number N for various magnetic fields. Solid line is computed using theGST with g = 6 . × − (cid:126) m / s. Empty circles are experimental data extracted from Ref. [10]. SM-VII. VARIATIONAL ANALYSIS FOR RADIAL SIZE
Let ¯ φ ( r ) be a normalized mode function, the kinetic energy and the two-body interaction energies can be generallyexpressed as E k = − (cid:126) N M (cid:90) d r ¯ φ ∗ ( r ) ∇ ¯ φ ( r ) , (SM42) E B = g (cid:48) N (cid:90) d r | ¯ φ ( r ) | . (SM43)We further assume that the energy associated with the stabilization force is E νB = g (cid:48) ν N ν (cid:90) d r | ¯ φ ( r ) | ν (for ν > . (SM44)To proceed, we assume that the mode function adopts a Gaussian form, i.e.,¯ φ ( r ) = 1 √ πσ e − ( x + y ) / (2 σ ) √ πa z ) / e − z / (2 a z ) (SM45)where a z = (cid:112) (cid:126) / ( M ω z ) is the width of the axial harmonic oscillator. The total energy of the self-bound droplets canthen be expressed as ε = (cid:126) M (cid:18) σ − ˜ g Nσ + ˜ g ν N ν − σ ν − (cid:19) . (SM46)To derive the reduced interaction parameters ˜ g ( c,s )2 and ˜ g ( c,s ) ν , we first consider a pure coherent state for which thetwo-body and the thee-body interactions are E ( c )2 B = 4 π (cid:126) δaM √ a ↑↑ a ↓↓ ( √ a ↑↑ + √ a ↓↓ ) N (cid:90) d r | ¯ φ ( r ) | ,E ( c )3 B = g N (cid:90) d r | ¯ φ ( r ) | . The reduced interaction parameters can be evaluated to be˜ g ( c )2 = 4 √ π √ a ↑↑ a ↓↓ ( √ a ↑↑ + √ a ↓↓ ) | δa | a z , ˜ g ( c )3 = M g √ π (cid:126) a z . (SM47)Next, for single-mode squeezed vacuum state, we have E ( s )2 B = 3 4 π (cid:126) δaM √ a ↑↑ a ↓↓ ( √ a ↑↑ + √ a ↓↓ ) N (cid:90) d r | ¯ φ ( r ) | ,E ( s )3 B = 15 g N (cid:90) d r | ¯ φ ( r ) | . (a) (b) FIG. SM6: (color online). (a) σ r , σ ( s ) r , and σ ( c ) r versus N . (b) ε k , ε B , and ε B as functions of N . Here the s -wave scatteringlength is fixed at δa = − . a and energy is in units of (cid:126) / ( ML ) with L = 6 µ m. From left to right, two vertical dottedlines mark the critical atom numbers for the third- and the first-order phase transitions, respectively. which lead to ˜ g ( s )2 = 3˜ g ( c )2 and ˜ g ( s )3 = 15˜ g ( c )3 . Finally, for the EGPE theory [9], ˜ g ( c )5 / can be similarly evaluated as˜ g ( c )5 / = 102475 π (cid:32) √ a ↑↑ a ↓↓ (cid:112) a ↓↓ /a ↑↑ (cid:33) / (cid:112) / a / z π / f (cid:32) a ↑↓ a ↑↑ a ↓↓ , (cid:114) a ↓↓ a ↑↑ (cid:33) (SM48)where f ( x, y ) = (cid:88) η = ± √ (cid:16) y + η (cid:112) (1 − y ) + 4 xy (cid:17) / . (SM49)In Fig. SM5, we compare, the numerically computed σ GST curves with experimental data [10] not shown in the maintext. Apparently, some data do not have the signature V shape close to the critical atom number for the self-bounddroplet. Therefore, high precision data for radial size are still needed.
SM-VIII. THREE-DIMENSIONAL SELF-BOUND DROPLETS
For three-dimensional (3D) gases, the properties of the self-bound droplet states [11] can be studied similarly usingthe Gaussian-state theory. Here we also obtain SVS, SCS, and CS phases as those found in two-dimensional droplets.Figure SM6(a) shows the typical behavior of the radial size σ r as a function atom number N . Apparently, σ r ( N ) alsopossesses a double-dip structure. In addition, the SVS-to-SCS and the SCS-to-CS transitions are of third and firstorder, respectively.Analytically, we may also carry out similar analyses to those in two-dimensional gases. Here, we only present themain results for 3D binary droplets. The virial relation in 3D takes the form2 E k + 3 E B + 6 E B = 0 . (SM50)With the assumption of a Gaussian mode function¯ φ ( r ) = 1( √ πσ ) / e − r / (2 σ ) (SM51)with σ r being the radial size, the energy per atom can then be expressed as ε = 3 (cid:126) M (cid:18) σ − ˜ g Nσ + ˜ g N σ (cid:19) , (SM52)where the reduced interaction parameters depend on the quantum phases of the condensate. Specifically, for SVS, ˜ g and ˜ g becomes ˜ g ( s )2 = 8 √ π | δa |√ a ↑↑ a ↓↓ ( √ a ↑↑ + √ a ↓↓ ) and ˜ g ( s )3 = 10 M √ π (cid:126) g , (SM53)7 FIG. SM7: (color online). Peak condensate density n p as a function of atom number N for δa = − . a . The dashed linemarks the critical condensate density n ∗ . The vertical dotted line denotes the boundary between the SCS and the CS phases. respectively. And similar to the 2D case, ˜ g ( c )2 = ˜ g ( s )2 / g ( c )3 = ˜ g ( s )3 /
15 for CS.The equilibrium radial size is the solution of equation ∂ε /∂σ r = 0 that satisfies ∂ ε /∂σ r >
0. In Fig. SM6, weplot the equilibrium radial sizes σ ( s ) r ( N ) and σ ( c ) r ( N ) obtained with the reduced interactions parameters (˜ g ( s )2 , ˜ g ( s )3 )and (˜ g ( c )2 , ˜ g ( c )3 ), respectively. As can be seen, they are in rough agreement with the numerical results at SVS and CSregimes. In particular, the critical atom number of the self-bound droplet is determined by the equations: ∂ε /∂σ r = 0and ∂ ε /∂σ r = 0, which leads to N cri = 64 √ ˜ g / (27˜ g ). We then find two critical atom numbers N ( s )cri = 2 π ( √ a ↑↑ + √ a ↓↓ ) a ↑↑ a ↓↓ δa (cid:115) M g √ π (cid:126) and N ( c )cri = 2 π ( √ a ↑↑ + √ a ↓↓ ) a ↑↑ a ↓↓ δa (cid:115) M g √ π (cid:126) (SM54)for SVS and CS phases, respectively.In Fig. SM6, we plot the energy contributions as a function of N . It can be analytically shown that the SVS-to-SCStransition occurs when 2 E k = | E B | or at N ∗ = 4˜ g (cid:114) g π ( √ a ↑↑ + √ a ↓↓ ) a ↑↑ a ↓↓ δa (cid:115) M g √ π (cid:126) . (SM55)Finally, for the SCS-to-CS transition, the phase boundary is located at roughly n ∗ = 1 g π (cid:126) M a ↑↓ − a ↑↑ a ↓↓ ( a ↑↑ + a ↓↓ − a ↑ , ↓ ) . (SM56)Again, as shown in Fig. SM7, it can be verified that the phase boundary given by the criterion n p = n ∗ is in goodagreement with the numerical result. [1] T. Shi, E. Demler, and J. I. Cirac, Ann. Phys. (N Y) , 245 (2018).[2] T. Shi, J. Pan, and S. Yi, arXiv:1909.02432 (2019).[3] T. Guaita, L. Hackl, T. Shi, C. Hubig, E. Demler, and J. I. Cirac, Phys. Rev. B , 094529 (2019).[4] Y. Wang, L. Guo, S. Yi, and T. Shi, Phys. Rev. Research , 043074 (2020).[5] See, e.g., C. Gerry and P. Knight, Introductory Quantum Optics (Cambridge University Press, 2005).[6] M. Schmitt, M. Wenzel, F. B¨ottcher, I. Ferrier-Barbut, and T. Pfau, Nature , 259 (2016).[7] F. B¨ottcher, M. Wenzel, J.-N. Schmidt, M. Guo, T. Langen, I. Ferrier-Barbut, T. Pfau, R. Bomb´ın, J. S´anchezBaena, J.Boronat, and F. Mazzanti, Phys. Rev. Research , 033088 (2019).[8] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. , 463 (1999).[9] D. S. Petrov, Phys. Rev. Lett. , 155302 (2015).[10] Cabrera, CR and Tanzi, L and Sanz, J and Naylor, B and Thomas, P and Cheiney, P and Tarruell, L, Science , 301(2018)[11] G. Semeghini, G. Ferioli, L. Masi, C. Mazzinghi, L. Wolswijk, F. Minardi, M. Modugno, G. Modugno, M. Inguscio, andM. Fattori, Phys. Rev. Lett.120