Theory for Self-Bound States of Dipolar Bose-Einstein Condensates

We investigate the self-bound states of dipolar Dy condensates with the Gaussian-state ansatz which improves the conventional coherent-state ansatz with multimode squeezed coherent states. We show that the self-bound states consist of the experimentally observed self-bound liquid phase and the unobserved self-bound gas phase. The numerically obtained gas-liquid boundary is in good agreement with experimental data. Our theory also allows one to extract the real part of the three-body coupling constant of the Dy atoms from the particle number distribution of the condensates. In particular, we results show that the self-bound states are stabilized by the short-range three-body repulsion. Our study shed a different light to understand the self-bound droplets of Bose-Einstein condensates.

Introduction.-Since the observations of the stable liquid droplets in dipolar condensates in the mean-field collapse regime [1,2], this novel quantum phase has received increasing attention in the community. To date, there have rapid experimental progresses in exploring the fundamental properties of these quantum droplets, such as the stabilization mechanism [2,3], the collective excitations [2,4], the self-bound states [5,6], the striped states [7], and the metastable states with supersolid characteristics [8][9][10]. Aside from dipolar systems, liquid droplets were also observed in attractive bosonic mixtures [11][12][13]. In addition, the droplet-droplet collisions in Bose-Bose mixtures were also experimentally studied [14].
From the theoretical side, the stabilization mechanism of the liquid droplets is of fundamental importance. It is now widely accepted that, unlike normal liquids, the twobody attraction is balanced at large density by quantum fluctuation [2,3]. To the first order, the quantum fluctuation can be incorporated into the Gross-Pitavskii equation (GPE) through the Lee-Huang-Yang (LHY) correction to the condensate energy [15] and leads to the extended Gross-Pitaevskii equation (EGPE) [16][17][18]. Although the EGPE explains many experimental observations, the discrepancy between the theoroetical and experimental results still exists [6]. In fact, there exists an intrinsic flaw to apply the EGPE to the droplet state as the LHY correction is obtained under the assumption that the collective excitation upon the coherent ground state is stable. Moreover, the EGPE cannot explain the asymmetric particle number distribution (PND) obtained in experiments [5,6].
In this work, we revisit the stabilization-mechanism problem of the self-bound dipolar droplets by using the * Electronic address: syi@itp.ac.cn † Electronic address: tshi@itp.ac.cn Gaussian-state theory (GST) [19][20][21] in which the manybody ground state wave functions of the condensates take the form of a multimode squeezed coherent state. As a result, the quantum fluctuation is treated in a selfconsistent manner within the mean-field framework. A critical ingredient in our study is the conservative threebody interaction. The inclusion of the squeezing components allows us to obtain the unknown real part of the three-body coupling strength by fitting the experimentally measured PND [5], which is subsequently used to map out the phase diagram on the parameter plane consisting of the relative dipolar interaction strength and the particle number. Originating from the competition between the long-range two-body attraction and the shortrange three-body repulsion, there exist two self-bound phases, i.e., the self-bound liquid (SBL) and the selfbound gas (SBG) phases. Remarkably, the numerical results for the critical particle number of the SBL phase are in good agreement with the experimental values. Although the three-body repulsion was previously considered the force that balances the two-body attraction at higher densities in Refs. [22][23][24], these GPE-based studies implicitly assume that the many-body ground state wave function of the condensate is a coherent state. Consequently, the estimated real part of the three-body coupling constant is much larger than the expected value [3]. While in our work, the appearance of the squeezed component in the self-bound states gives rise to unusual enhancement to the three-body coupling constant such that the two-body attraction is balanced at a rather smaller three-body coupling strength.
Formulation.-We consider a ultracold gas of N 164 Dy atoms in the free space. The total Hamiltonian of the system, H = H 0 + H 2B + H 3B , consists of the single-, two-, and three-particle terms. In the secondquantized form, the single-particle Hamiltonian reads H 0 = drψ † (r)h 0ψ (r), where ψ(r) is the field operator and h 0 = −∇ 2 /(2M ) − µ with M being the mass of the atoms and µ the chemical potential. The two-particle Hamiltonian takes the form Here, we assume that two atoms interact via the contact and the dipolar interactions, i.e., where a s is the s-wave scattering length, µ 0 is the vacuum permeability, µ = 9.93 µ B is the magnetic dipole moment of Dy atom with µ B being the Bohr magneton, and θ r is the polar angle of r. For convenience, we introduce the relative dipolar strength ε dd = µ 0 µ 2 M/(12π 2 a s ). Finally, the three-body interaction term can be expressed as where g 3 is the three-body coupling constant. Although the real part of g 3 is difficult to measure experimentally, as shall be shown, it can be estimated through our theory. Without loss of generality, we assume that g 3 is real since we shall only focus on the ground state of the system. To proceed, let us briefly outline the GST which assumes that the many-body wave function of the condensate takes the following variational ansatz [19][20][21] where the coherent part of the condensate is described by the mean value Φ(r) = Ψ (r) = φ c (r), φ * c (r) T ofΨ(r) ≡ ψ (r),ψ † (r) T defined in the Nambu basis, Σ z (r, r ) = σ z ⊗ δ(r − r ) with σ z being the Pauli matrix, and the matrix ξ(r, r ) characterizes the squeezing of the condensate. It should be noted that, for shorthand notation, the products in the exponents of Eq. (4) should be understood as the matrix multiplication in the Nambu space and the integration in the coordinate space [25]. To remove the redundancy [19,20], instead of ξ we introduce the covariance matrix Γ(r, r ) = {δΨ(r), δΨ † (r )} as variational parameters, where {·, ·} denotes the anticommutator and δΨ =Ψ − Φ is the fluctuation field. It can be verified that Γ and ξ are connected through Γ = SS † with S = e iΣ z ξ satisfying SΣ z S † = Σ z . Furthermore, the matrix elements of Γ in the 2 × 2 Nambu representation are related to the normal and the anomalous correlation functions, G(r, r ) = δψ † (r )δψ(r) and F (r, r ) = δψ(r)δψ(r ) , through the relations Γ 11 (r, r ) = Γ T 22 (r, r ) = 2G(r, r ) + δ(r − r ) and Γ 12 (r, r ) = Γ † 21 (r, r ) = 2F (r, r ). The Wick's theorem gives rise to the mean-field Hamiltonian [25] H MF = E 0 + (δψ † η + H.c.) + 1 2 : δΨ † HδΨ : which converge for a sufficiently large imaginary time τ . The numerical procedures for evolving Eq. (6) are detailed in the Supplementary Materials. General properties of the ground states.-Before presenting the results, let us first analyze the general structure of the many-body wave function for the condensate.
To this aim, we diagonalize the normal and the anomalous correlation functions into G(r, r ) ≈ S j=1 N s,jφs,j (r)φ * s,j (r ) and F (r, r ) ≈ S j=1 N s,j (N s,j + 1)φ s,j (r)φ s,j (r ) (via the Takagi diagonalization), where S is the total number of squeezed modes, N s,j is the atom number in the jth squeezed modē φ s,j satisfying the orthonormal condition drφ * s,iφ s,j = δ ij . Without loss of generality, we always assume that N s,j are sorted in descending order with respect to the index j. The total number of atoms in the squeezed modes is N s = j N s,j . Now, the ground state wave function can be expressed as where N c = dr|φ c (r)| 2 is the atom number in the coherent state,b c = drψ(r)φ * c (r) withφ c (r) = φ c (r)/ √ N c being the normalized mode function for the coherent component, and ξ j = arcsinh N s,j is the squeezing parameter in the modeb s,j = drψ(r)φ * s,j (r). Of particular interest, in the presence of the squeezed component, the PND and equivalently the quantum statistics of |Ψ GS may differ from that of a coherent state significantly. To see this, we expand the coherent modeb c using the squeezed modes aŝ where α j = drφ * s,j (r)φ c (r), α 2 ⊥ = 1 − j |α j | 2 , andb ⊥ represents the mode that is perpendicular to allb s,j . The variational ground state (7) can now be expressed as

The PND for modeb
where γ j = α j cosh ξ j + α * j sinh ξ j and H (x) are Hermite polynomials. The PND of |Ψ GS can then be expressed into a recursive form as follows. Let P j (n) denote the PND of the state containing first j squeezed modes, the PND containing first j + 1 squeezed modes can then be expressed as where P 0 (n) ≡ p ⊥ (n). Applying the relation (10) successively will eventually leads to the total PND P S (n). In general, P S (n) is peaked at roughly n ≈ N c . To the left of this peak, the PND is mainly determined by Poission distribution of the coherent state such that P S (n) increase abruptly as n approaches N c . While to the right of the peak, P S (n) has a long tail at large n [26]. As a result, P S (n) becomes asymmetric around N c , which is in striking contrast to the PND of a coherent state.
Determination of g 3 .-In Ref. [5], the PNDs of the condensates were experimentally measured under various magnetic fields, or equivalently, ε dd 's. Because the measured PNDs were asymmetric around the mean particle number, they were fitted with a function that is the convolution of a Gaussian and a Maxwell-Boltzmann distribution in Ref. [5]. Here we attribute the measured PNDs to that of the Gaussian state |Ψ GS and fit the measured PNDs with P S (n). Specifically, for every PND obtained under a given magnetic field in Ref. [5], we have carried out extensive fittings by systematically varying the values of N and g 3 . It is found that the experimentally measured PNDs are broader than that obtained with P S (n), which might originate from the detection noise and the thermal fluctuations of the condensates in the evaporation. The best fit, as shown in Fig. 1, is achieved under the magnetic field B = 6.624 G with fitting parameters N * = 790 and g * 3 = 6.73×10 −40 m 6 /s, where N * roughly agrees with the experimentally measured mean number of atoms. Interestingly, the value of g * 3 is about one order of magnitude smaller that used in previous work [22] and is therefore more favorable [3]. Unless otherwise stated, the value of the three-body coupling strength will always be fixed to g * 3 in below. Phase diagram.-In Fig. 2(a), we map out the coherent-component fraction f c ≡ N c /N of condensate in the ε dd N parameter plane. To reveal more details, we plot the typical ε dd and N dependence of f c in Fig. 2(b) and the ε dd dependence of the kinetic energy E kin , the two-body interaction energy E 2B , and the three-body interaction energy E 3B in Fig. 2(c). The detailed expressions for these energy components can be found in the Supplementary Materials [25]. We note that for the parameters covered by our numerical calculations, we always have E 2B < 0. Moreover, since the peak density n peak and the radial 1/e 2 width σ ρ (or the axial 1/e 2 width σ z ) are of interest to experimentalists, we also map out the values of log n peak and σ ρ on the ε dd N parameter plane in Fig. 3. An immediate observation is that all three figures share the same structure. Namely, the parameter plan is divided into three regions representing the expanding-gas (EG) phase, the SBG phase, and the SBL phase. The condensate of the EG phase is not a selfbound state, it is therefore not of interest to the present work.
In the SBG phase, the attractive interaction or the atom number is large enough for the condensate to form a self-bound state. Because the density is still relatively low [see Fig. 3(a)], the condensate remains in the gas phase. A important feature of the SBG phase is that the squeezed component is always dominated by a single squeezed mode (see Fig. 5). More interestingly, the SBG phase can be further divided into two subregions by the dotted line in Fig. 2(a): i) To the left of this boundary, the coherent-component fraction vanishes completely such that the state of the condensate is a single-mode squeezed vacuum, in analog to that of a trapped weakly attractive condensate [21]. ii) To the right of the dotted line, the ground state of the condensate is a squeezed coherent state in which f c can be as high as around 0.64. Analyzing the ε dd dependence of the ground-state energy E 0 (= E kin + E 2B + E 3B ) further reveals that the transition between the squeezed vacuum and the squeezed coherent states is of third order, since the third derivative of E 0 with respect to ε dd is discontinuous. This transition reflects the fact that the appearance of the coherent component breaks the Z 2 symmetry possessed by the pure squeezed states. Physically, increasing ε dd or N enhances both two-body and three-body interactions, which shrinks the condensate and leads to a higher density (see Fig. 3). It turns out that the coherent condensate is more energetically favorable in the high density regime, where the squeezed atoms are turned into the coherent ones. This mechanism can be further confirmed by deliberately lowering the value of g 3 in numerical calculations. It is found that for the same set of ε dd and N , the coherent-component fraction monotonically increases when g 3 is lowered.
As one further increases ε dd or N , the three-body repulsion becomes determinative such that f c experiences an abrupt increase and is now dominant in the condensate (f c 0.8). The system then enters the SBL phase through a first-order phase transition where, as shown in Fig. 3, the n peak (σ ρ ) of the condensate is significantly increased (reduced) compared that in the SBG phase. Remarkably, the SBG-SBL phase boundary is in good agreement with the experimentally measured critical atom number [5]. As shown in Fig. 4, we have also numerically computed the critical atom number for the recent 162 Dy experiment [6], which is in very good agreement with experimental results compared to EGPE with LHY corrections.
To gain more insight into the gas-liquid transition, we examine the number of the squeezed modes by introducing the so-called notably-populated squeezed modes (NPSMs). Specifically, the jth squeezed mode is an NPSM if N s,j /N s ≥ 0.5%. In Fig. 5, we present the distribution of the number of the NPSMs in the ε dd N plane, where the SBL-SBG boundary is sharply defined. The SBG phase is dominated by a single squeezed mode; while across the gas-liquid boundary, the coherent component becomes dominant and the number of NPSMs increases abruptly. Physically, the squeezing in the SBG phase originates from the overall attractive two-body interactions, same as that in the condensate with pure swave interaction [21]. On the other hand, the squeezing in the SBL phase represents the depletion of the condensate that is induced by the three-body repulsion. Indeed, it is verified that the fraction of the squeezed component increases if we increase g 3 or, equivalently, decrease ε dd .
Density profiles.-We now examine the density profiles which consists of the coherent part n c (r) = |φ c (r)| 2 and the squeezed part n s (r) = G(r, r). Due to axial symmetry of the condensate, both n c and n s are cylindrically symmetric. In Fig. 6, we plot the density profiles along the radial ρ = x 2 + y 2 and the z directions for both coherent and squeezed components with N = 700 and for various ε dd 's. In addition, because the density is symmetric about the z = 0 plane, only the density profiles for z ≥ 0 are plotted.
For the coherent component, n c always reaches its peak value at the center of the condensate. It then monotonically decreases along both the radial and the axial directions. In addition, as ε dd increases, the peak value of n c also grows along with the increase of f c . As to the squeezed component, similar to the behavior of n c , n s also reaches its peak at the condensate center when ε dd is small. However, for large ε dd , the value of n s at the condensate center drops such that both n s (ρ, 0) and n s (0, z) are not monotonic. Physically, this can be understood by noting the total density, n c + n s , is always maximized at the center of the condensate. As a result, the three-body repulsion is strongest at the center of the condensate which locally turns the squeezed component into the coherent one. Now, the structure of the a liquid droplet can be roughly described as follows: the highdensity coherent component forms the liquid core of the condensate; while the low-density squeezed component constitutes a gas shell surrounding the liquid core.
Effective equations for the self-boudn states.-In principle, we may derive effective equations for the groundstate spatial modes φ c and {φ s,j } S j=1 by minimizing the ground-state energy E 0 under the constraints that the unnormalized squeezed modes φ s,j = N s,jφs,j are mutually orthogonal. This however becomes very cumbersome deep inside the SBL phase since the number of the squeezed modes could be very large. We notice that deeply inside the SBL phase most of atoms are in the coherent state, while in the SBG phase the squeezed component is dominated by a single mode, i.e., S = 1. As a result, in the deep SBL phase, the ground state is described by a single macroscopic wavefunction φ c (r) that obeys the conventional GPE. In the SBG phase, both coherent and squeezed condensates contain the macroscopic number of particles, which are characterized by two macroscopic wavefunctions χ(r) ≡ (φ c (r), φ s (r)) T , where φ s (r) ≡ φ s,1 (r) describes the mode with largest squeezing. The conditions δE 0 /δχ † = 0 then give rise to the squeezed-coherent-state GPEs where I is a 2 × 2 identity matrix, M (2) (r) = dr U (r − r ) n 1 (r ) n sc (r ) n sc (r ) n 3 (r ) (12) and are the mean fields induced by the two-and three-body interactions, respectively. Here we have introduced the notations n q = |φ c | 2 + q|φ s | 2 and n sc = 2Re(φ * s φ c ). We note that the diagonal elements of M (2) are contributed by the intra-state interaction and by inter-state Hartree interaction; while the off-diagonal terms originate from the inter-state Fock and pairing interactions.
To appreciate the physical significance of Eqs. (11), let us consider the simplest cases by assuming that the condensate is either in a pure coherent state or in a singlemode squeezed vacuum state. For either case, Eqs. (11) reduce to a unified form: h 0 + c 2 dr U (r − r )|ϕ(r )| + c 3 g3 2 |ϕ(r)| 4 ϕ(r) = 0, where c 2 = c 3 = 1 if ϕ ≡ φ c and c 2 = 3 and c 3 = 15 if ϕ ≡ φ s . As can be seen, although both two-and three-body interaction strengths of the squeezed state are enhanced compared to those of the coherent state, the three-body interaction strength gains more increment. As shown in the expression of E 3B in the Supplemental Material [25], this unusually large enhancement to the three-body interaction strength originates from the fact that the combination of squeezing and three-body interaction yields more terms when applying the Wick's theorem.
We have performed extensive numerical calculations to check the validity of these effective equations. It is found that away from the gas-liquid boundary the results obtained with the effective equations are in perfect agreement with those computed using Eqs. (6). While close to the gas-liquid boundary, because number of the squeezed modes is rapidly increased, the effective equations becomes less accurate.
Discussion and conclusion.-As analyzed in our previous work [21], if we sequentially solve Eqs. (6a) and (6b) iteratively starting with a pure coherent state, the EGPE is merely the resulting equation after the first iteration. Therefore, EGPE is applicable only when the condensate is dominated by the coherent component, which excludes its application to the SBG phase. Even for the coherentcomponent-dominant SBL phase one should be very careful about the validity of EGPE. In fact, at its current form, the EGPE only takes into account the quantum fluctuation induced by the two-body interactions; while the quantum fluctuation in the SBL phase is induced by the three-body repulsion. To further confirm this, we switch off the three-body repulsion by setting g 3 = 0 in numerical calculations, it is found that condensates always collapse for all parameters covered by Fig. 2(a). This result clearly suggests that the quantum fluctuation associated with two-body interactions is insufficient to balance the two-body attractions.
We note that three-body repulsion was previously considered as the source that stabilizes the two-body attraction in Refs. [22][23][24]. However, to match the experimental measurements, the estimated three-body coupling constant using the GPE is much higher than the expected value [5]. Within the GST, the appearance of the squeezed component leads to the unusual enhancement to the three-body coupling constant such that the twobody attraction is balanced at a rather smaller g 3 . For parameters deep inside the SBL phase where the squeezing becomes negligibly small, it has been verified that the GPE with three-body repulsion can lead to the same results as those obtained with Eqs. (6).
In summary, we have studied the self-bound states of dipolar condensates by using the GST which allows us to treat the quantum fluctuation in a self-consistent manner. Aside from the self-bound liquid droplet phase, we predict the existence of the SBG phase. We have also shown that the gas-liquid transition is induced by competition between the long-range two-body attraction and the short-range three-body repulsion. Unlike the conventional theory based on coherent states, squeezed component emerges in the ground states of both phases with distinct physical origins: it is induced by the two-body attraction (three-body repulsion) in the SBG (SBL) phase. The appearance of the squeezed component in the selfbound states is crucial as (particularly in the SBL phase) it allows the three-body repulsion to balance the twobody attraction even with a relatively small three-body coupling constant. For experimental detections, we note that the squeezed-component-dominant SBG phase can be detected by measuring the second-order correlation function of the gas as it exhibits distinct quantum statistics from that of the coherent-component-dominant SBL phase [21]. In the future, we shall further study the dynamical formation of the self-bound states. We believe our study opens up a broad avenue to explore the brand new macroscopic quantum states in cold atom systems.

Supplemental Material
This supplemental material is divided in three sections. In Sec. SM1, we derive the mean-field Hamiltonian, i.e., Eq. (5) in the main text. In Sec. SM2, the EOM in the imaginary time evolution is projected in the angular momentum basis e imϕ / √ 2π, which can be solved efficiently by the Hankel transformation. In Sec. SM3, we derive the coupled GP-like equation for coherent and squeezed condensates by minimizing the ground state energy.