Beyond Lee-Huang-Yang description of self-bound Bose mixtures

We investigate the properties of self-bound ultradilute Bose-Bose mixtures, beyond the Lee-Huang-Yang description. Our approach is based on the determination of the beyond mean-field corrections to the phonon modes of the mixture in a self-consistent way and calculation of the associated equation of state. The newly obtained ground state energies show excellent agreement with recent quantum Monte Carlo calculations, providing a simple and accurate description of the self-bound mixtures with contact type interaction. We further show numerical results for the equilibrium properties of the finite size droplet, by adjusting the Gross-Pitaevskii equation. Finally, we extend our analysis to the one-dimensional mixtures where an excellent agreement with quantum Monte Carlo predictions is found for the equilibrium densities.


Introduction
In classical physics the formation of a liquid droplet, i.e. of a self-bound state, typically arises from the interplay between the short-range repulsive and long-range attractive components of the interatomic potential. In quantum fluids, the formation of self-bound droplets has been intensively investigated in liquid Helium, with experimental observations of strongly interacting nanodroplets in both 3He and 4He [1,2]. Recently, it has been pointed out that an ultradilute self-bound state of matter can be formed in mixtures of ultracold atomic gases, though the underlying physics is different [3]. For binary mixtures of Bose-Einstein condensates (BECs), mean-field analysis predicts the system to become unstable against collapse when the attractive inter-species interaction overcomes the repulsive interaction between identical atoms [4]. However, in the utradilute liquid phase, the mean-field collapse is avoided as quantum fluctuations stabilize the system. The liquid droplets formed in ultracold atomic gases are fundamentally different from those in classical or Helium fluids, since it arises from beyond mean-field effects and exhibits extreme diluteness. The observation of such ultradilute liquids has been first achieved in dipolar Bose gases [5,6], where the formation mechanism is the same, arising from an interplay between the attractive dipolar interaction and quantum fluctuations [7]. More recently, the liquid phase has been also observed in attractive Bose-Bose mixtures, both in trap-confined [8,9] and free-space configurations [10]. These experimental works found overall good agreement with the theory developed in the seminal work of Petrov [3].
Although the theory of Petrov [3] reckons success in explaining the stabilization mechanism and providing the energy functional, it is known that the model suffers from a serious conceptual problem. Indeed, in the relevant regime of droplet formation, the theory predicts a purely imaginary phonon velocity for the low-lying excitation spectrum, and thus a complex energy functional. In the original work [3], this problem is contoured by using the velocities calculated at the threshold point, and explicitly putting to zero the value of the suspicious phonon velocity. While such approximation is justified at the mean-field collapse point where the droplet is yet to be formed, its validity for any finite droplet is questionable. As a matter of fact, quantum Monte Carlo (QMC) method from Ref. [11] showed that the accuracy of predictions of Ref. [3] becomes worse as one increases the attractive inter-species interaction. Another issue concerns the disagreement between experiment and theory for the critical number of atoms and the droplet size, as reported in Ref. [8]. In this regard, theoretical works based on beyond mean-field variational [12] and QMC [13] approaches pointed out the crucial role played by finite-range effects. Recently, many theoretical works have been devoted to the investigation of the liquid phase in low-dimensional systems, motivated by the enhanced role of quantum fluctuations [14][15][16]. In particular, one-dimensional (1D) binary mixtures experimentally constitute a perfect playground due to enhanced stability, as the three-body recombination rate is greatly suppressed, and accessibility of a wide regime of interactions, as the coupling constant can even take infinite values without compromising the stability of the system. Also theoretically, 1D geometry is appealing since the energy functional does not suffer from the aforementioned imaginary part and pseudopotential interaction can be used in QMC simulations.
The aim of this paper is to provide a description of the symmetric droplet in binary mixtures of bosons, going beyond the Lee-Huang-Yang (LHY) framework. This is achieved in a phenomenological way, by explicitly including higher order corrections to the Bogoliubov speed of sound in the LHY energy. Although calculated in an approximated way, the resulting beyond LHY correction to the equation of state is found to deeply modify the equilibrium properties of the symmetric mixture. In particular, we find a strong dependence of the equilibrium density on the value of interactions, in excellent agreement with available QMC simulations. We further investigate the equilibrium properties of the finite-size droplet within the local density approximation, and extend our analysis to the 1D mixtures.
This paper is organized as follows: in Sec. 2 we introduce the beyond LHY theory for the droplet, based on the calculation of second order terms in the long wavelength modes of the excitation spectrum. In Sec. 3 we report results obtained in the thermodynamic limit N → ∞ and V → ∞ with N/V = const. These results are compared with available QMC calculations. Section 4 is devoted to the numerical analysis of the finite-size droplet, using a generalized Gross-Pitaevskii equation. Extension of the analysis to the one-dimensional mixture is discussed in Sec. 5.

Theory
We consider a uniform binary mixture of bosons with equal masses (m 1 = m 2 = m). In terms of the single-particle creation and annihilation operators in each component,â † i,k and a i,k (i = 1, 2), the Hamiltonian including all two-body collisions takes the form: where ε k = 2 k 2 /(2m) and we have assumed a contact-type interactions between particles, characterized by coupling constants g ij . The ground state energy of the system is obtained by diagonalizing the Hamiltonian (1). This is achieved by applying the Bogoliubov prescription and replacingâ i,k andâ † i,k by the total number of atoms in each component √ N i , as well as appropriated canonical transformations. The details of the calculation can be found elsewhere [17,18] leading to the following form: where α † k and β † k are, respectively, the creation operators for the quasiparticles in the density and spin channels obeying Bose statistics. The excitation spectrum of the system reads E ±,k = ε 2 k + 2mc ± ε k with c ± the sound velocities defined hereafter. The ground state energy becomes The first term of Eq. (3) describes the mean-field internal energy, whereas the second one corresponds to the contribution from quantum fluctuations and is often referred to as the LHY term [19]. Within the Bogoliubov theory, the long wavelength modes of the excitation spectrum are given by the linear phonons. For a symmetric mixture g 11 = g 22 = g, the speed of sound is [4] c 2 ±,B = 1 2m g(n 1 + n 2 ) ± g 2 (n 1 − n 2 ) 2 + 4g 2 12 n 1 n 2 , with n i = N i /V the atomic density of each component. The main idea of our work is to extend the LHY description in a perturbative way, by evaluating the sound velocities beyond the Bogoliubov formula (4) in a self-consistent manner and to obtain the correction to the ground state energy Eq. (3). The calculation of higher order terms for the excitation spectrum can be achieved either microscopically, by developing the second-order Beliaev theory for the mixtures [20,21], or in a less fastidious macroscopic approach, by means of thermodynamic relations. Indeed, it is known for a single-component weakly interacting Bose gas that at T = 0, the long wavelength phonon mode is related to the compressibility κ as c = (mnκ) −1 [22].
In an analogous way, one can relate the sound mode in the density and spin channel of the Bose mixtures, to the compressibility κ + and spin susceptibility κ − of the system, respectively [23,24]. These quantities are obtained from the energy (3) according the thermodynamic relation yielding the speed of sound c ± = (mnκ ± ) −1 , with n = n 1 + n 2 the total atoms density. In what follows, we evaluate both the speed of sound and the associated ground state energy for the three-dimensional (3D) mixtures. The extension of LHY theory in lower dimension follows essentially the same path and we will discuss as an example the one-dimensional (1D) mixture in the last part of this paper. In 3D, the LHY contribution in Eq. (3) exhibits an ultraviolet divergence, arising from an approximate relation between the coupling constant g ij = 4π 2 a ij /m and s-wave scattering length a ij in the first Born approximation. This is conveniently solved by a proper renormalization of the coupling constant [25]: g ij → g ij (1 + g/V k m/( k) 2 ). Then, the momentum sum in Eq. (3) can be turned into an integral which can be performed analytically resulting in The regime of interest corresponds to repulsive intra-species interaction g > 0 and attractive inter-species interaction g 12 < 0, with a small disbalance |δg|/g 1 where δg = g + g 12 . For a system satisfying the inequality δg < 0, the mean-field field theory would result in energy ∝ n 2 given by the first two terms of Eq. (6) and would predict a collapse of a homogeneous state towards bright soliton formation. The beyond mean-field theory eliminates the mechanical instability as the quantum fluctuations generate a repulsive term ∝ n 5/2 . The interplay between such attractive and repulsive forces is at the heart of droplet formation.
However, energy (6) suffers from the presence of a dynamical instability. This can be easily seen for the unpolarized mixture (n 1 = n 2 = n/2), for which the Bogoliubov speed of sound (4) in the density channel takes the value c +,B = δgn/2m, therefore being purely imaginary in the liquid phase. In the original work of Petrov [3], this issue is circumvented by putting δg = 0 in Eq. (4): The imaginary phonon in the Bogoliubov theory indicates that not only the equation of state (6) needs the presence of a beyond mean-field LHY term to be stabilized, but also the sound velocity requires additional higher-order correction in order to be well defined. In the unpolarized configuration, one immediately finds from Eqs. (4) and (6) the compressibilities of the mixture: It is noteworthy that both susceptibilities are complex in the liquid phase. Although non-zero imaginary part arises naturally in the perturbative approach, one notices that the imaginary component is of order |δg| 5/2 , and can be safely neglected for the experimentally relevant parameter range |δg|/g 1. This is equivalent to neglecting fluctuations in the density channel, while preserving those in the spin channel. Therefore we obtain the following beyond Bogoliubov expressions for the speed of sound: We show in Fig. 1 a comparison between the Bogoliubov sound velocity Eq. (4) and higher order sound velocity (10). One striking feature which our theory predicts is that the velocity of the density mode becomes real above a certain density when beyond mean-field corrections are included, while it is a purely imaginary quantity in the Bogoliubov description. Another  (4), while the black dashed line is obtained by setting g 12 /g = −1 in the LHY term as it was originally done in Ref. [3]. The green dashed and blue dotted-dashed lines correspond, respectively, to the inclusion of higher order term in the density (10a) and spin (10b) sound velocity. The red solid line includes both density and spin corrections. The circle indicates the spinodal density.
feature is that even in the higher order description, the speed of sound becomes imaginary in a small window of density na 3 (δg/g) 2 (see inset of Fig. 1). However, this instability has a physical nature and it defines the spinodal point below which the uniform liquid is unstable towards the formation of multiple droplets. That is, at zero pressure, the liquid is self-bound and it stays at the equilibrium density which corresponds to the position of the minimum in the equation of state. If positive (negative) pressure is applied, the density of the liquid increases (decreases) with respect to the equilibrium density and the energy increases. If the applied pressure is large and positive the energy eventually becomes positive, still the homogeneous system remains stable. On the contrary, for large negative pressures the homogeneous shape can no longer be sustained and the liquid fragments into droplets each having density close to the equilibrium one.
Once we have shown that higher-order corrections remove the unphysical instability associated with the complex values of the speed of density mode, it is useful to investigate if the predictions for the ground-state energy can be also improved.

Energy analysis
We now recalculate the LHY term using the beyond-Bogoliubov sound velocities Eq. (10) and improve the equation of state (3). The resulting energy is shown in Fig. 2 with a red solid line. It has a shape typical for a liquid with the minimum associated with the equilibrium density. It is instructive to compare our results with the ones obtained using the original prescription from Ref. [3], Eq.  (12) and density (11), for different values of δg/g. The lowest black solid line is the universal result from Ref. [3]. The color solid lines are the predictions from the beyond LHY theory, while the color dots are QMC results from Ref. [11].
the Bogoliubov sounds (4) (black dotted line). We remind that in the latter case, the energy is complex, and we only show its real part in Fig. 2. Taking as reference the LHY energy with the Bogoliubov dispersion law, one can see that inclusion of higher order terms in the density sound (10a) (top green dashed line) changes only slightly the behavior of the energy. Instead, the inclusion of higher order terms in the spin sound (10b) (blue dashed-dotted line) strongly suppresses the energy. Thus, we conclude that the main correction to the equation of state arises from quantum fluctuations in the spin channel. It is worth noticing that the inclusion of density sound (10a) leads to a spinodal point below which the uniform liquid is unstable against density fluctuations (filled circle in Fig. 2). Although our theory provides a higher-order correction to the speed of sound, still not all second-order terms are taken into account as it would be in Beliaev theory [21]. Thus it is important to verify the validity of our results through a direct comparison with available Monte-Carlo calculations [11]. For this purpose, it is convenient to normalize both the energy and the density to their equilibrium values obtained within the approach of Petrov [3]: n 0 = 25π(a + a 12 ) 2 16384a 5 .
(12) Figure 3 shows E/|E 0 | calculated for different values of interaction disbalance δg/g, as a function of density n/n 0 . In the rescaled form, the energies evaluated within the original Petrov theory collapse to a single curve (shown with a black solid line in Fig. 3) which is independent of the specific value of δg/g. Predictions of our theory are shown with color lines and should be confronted with QMC results taken from Ref. [11]. The agreement between our beyond LHY theory and QMC results is surprisingly good, especially in the most interesting  7) and Eq. (4), respectively. The red solid line is the calculation within the beyond LHY theory, relying on Eq. (10), while the black dots are QMC results from Ref. [11].
region around the minimum of energy. This hints that that although we do not perform a systematic calculation of the third-order terms of the perturbative theory, in practice the contributions which we miss are small in the considered case of a symmetric mixture. The agreement of our theory with QMC results is further emphasized in Fig. 4 where we show the calculated equilibrium density (corresponding to the minimum of energy in Fig. 3) as a function of δg/g. In particular, for the largest value of |δg| our theory predicts a decrease for the equilibrium density of about 50% in respect to the prediction from Ref. [3]. For comparison, we also show the equilibrium density obtained using the Bogoliubov sounds(4) in the LHY energy, which exhibits a lower value.

Finite size droplet
Following the success of the presented theory in removing the instability in the speed of sound and significantly improving the equation of state of a homogeneous liquid, we aim at an improved description of finite-size droplets. A common path to do so [3,7,11] is to improve the energy functional. Differently from single-component gas, there is a separation of scales in the considered binary mixtures. That is, a finite-size droplet changes its shape at the distances of the order of the "large" healing length defined by the chemical potential, ξ ∝ 2 /(m|µ|), while the main contributions to the LHY terms in Eq. (3) arise from "short" distances ∝ 1/c + [3]. Under these specific conditions it is possible to incorporate higher-order terms locally as non-linear terms in the Gross-Pitaevskii equation (GPE) for the droplet. Following the notation of Ref. [3] we introduce the rescaled coordinater = r/ξ with ξ = 6 2 /(|δg|mn 0 ), where the equilibrium density n 0 is given in Eq. (12). Then, the energy functional associated to the equation of state (6) with the beyond Bogoliubov speed of sound (10) is written as where η = 6ξ 3 /(n 2 0 |δg|) and Φ is normalized as dr|Φ| 2 = N/n 0 . In Eq. (13) we have assumed c + = 0 to hold, so as to avoid the imaginary part of the energy functional. This is motivated from the analysis of the previous section, in which neglecting the density mode was found not to alter greatly the behavior of the equation of state. The energy functional Eq. (13) reduces to the one used in Ref. [3] in the limit δg → 0.
Before discussing the numerical results, it is insightful to notice that the leading order correction introduced by the beyond Bogoliubov sound in the energy functional is of attractive three-body nature. Expanding the quantum fluctuations term in Eq. (13), the energy functional yields term E ∝ K 3 n 3 0 |Φ| 6 /3!. Such cubic dependence can be interpreted as corresponding to three-body interactions with the strength given by An estimate using typical experimental values for 39 K with a = √ a 11 a 22 48a 0 where a 0 is the Bohr radius and a 12 /a = −0.115 provides |K 3 |/3! ∼ 10 −41 m 6 /s, therefore being of the same order as the three-body loss rate measured in real experiment [10]. Thus, the effective three-body interactions (15) might be of the same order as the three-body terms which are not included in the model Hamiltonian (1). It was proposed that inclusion of three-body interactions on its own might lead to stabilization of a droplet [27,28].
We show in Fig Fig. 5 are calculated for (a) N/n 0 ξ 3 = 3000 where a bulk region is formed in the center which is a hallmark of a liquid, and (b) N/n 0 ξ 3 = 30 corresponding to typical experimental conditions [10]. The most crucial effect is that the central density of the droplet is decreased while its size is simultaneously increased, in agreement with diminishing equilibrium density found in a homogeneous liquid, see Fig. 4. It is interesting to note that also in the experiment of Ref. [8], the measured data for the droplet size was found to be larger than the prediction from Ref. [3]. For a sufficiently small number of atoms, the density profile of the liquid phase is well described by a Gaussian function, and one can extract the width of the droplet from a fitting. The obtained result is shown in the inset of Fig. 5(b), where the size of the self-bound mixture is found to be systematically larger than the value σ 0 predicted from Petrov's theory. At the same time, σ 0 increases when δg → 0, so the actual droplet size is larger for a smaller disbalance. Another experimentally relevant quantity is the critical atoms number for the droplet formation. Below a certain number of atoms, the droplet state becomes unstable and it evaporates. The critical number of atoms for the unstable phase can be conveniently investigated by means of variational approach. Close to the critical number, one can indeed safely use the Gaussian ansatz and assume the density profile of the gas to be [28] with σ the waist. Then for a fixed value of N , the equilibrium state corresponds to the value of σ for which the energy functional Eq. (13) is minimized. We have verified that while for sufficiently large number of atoms there is a global minimum in E(σ) at finite value of σ, corresponding to the droplet state, the later becomes a local minimum with E(σ) > 0 as one crosses the metastable point N ≤ N meta . Further decreasing N one reaches the critical number N c below which the energy minimum at finite σ vanishes and the ground state corresponds to a gas, σ → ∞. The metastable number as well as the critical number of atoms as a function of δg/g is reported in Fig. 6. We briefly note that within the variational approach, the theory of Petrov predictsÑ meta = 24 andÑ c = 19.6, thus slightly larger than the values Figure 6: Critical number of atoms as a function of δg/g. The solid and dashed lines are the predicted atom numbers for the unstable (N c ) and metastable (N meta ) droplet solution, respectively. The top red curves are evaluated within our theory, while the bottom black curves are calculated using the universal approach of Ref. [3].
N meta = 22.55 andÑ c = 18.65 reported in Ref. [3], calculated from GPE (14). We find that the inclusion of beyond LHY terms in the energy functional is responsible for shifting the critical number to higher values. Experimentally, the critical number of atoms for the droplet state of Bose mixtures in both trap-confined [8] and free space [10] configuration has been measured. While in the free space measurement N c was found to lie near the prediction of Ref. [3] for the metastable state (= 22.55), the trap confined measurement showed a deviation of N c to lower value. Recently, this deviation was accounted for the effects of finite-range in the interaction potential, which is neglected in the contact type s-wave description [12,13].

1D Mixtures
Although the true Bose-Einstein condensation associated with an off-diagonal long-range order does not exist in one-dimensional geometry, it is known [29] that the Bogoliubov theory is quantitatively correct for predicting the energy in the regime where the coherence is sustained for distances large compared to the mean interparticle distance [30]. Thus, one can still use the Bogoliubov theory Eq. (2) to study the mixture, and the momentum sum in Eq. (3) can be evaluated straightforwardly, yielding the result with the interaction coupling constant related to the s-wave scattering length according to g ij = −2 2 /(ma ij ). It is worth noticing that in the 1D mixture, quantum fluctuations have opposite contribution (notice the negative sign) as compared to the 3D case Eq. (6). Consequently, droplets are formed in the dominantly repulsive regime [14] where δg = g + g 12 > 0.  (7) of Ref. [3] while the dotted line is the beyond mean-field result with the Bogoliubov speed of sound Eq. (4). The red solid line is the prediction from our theory. Black dots are the QMC results from Ref. [15]. The dashed and solid vertical lines in panel (b) indicate the critical values for δg/g above which the liquid becomes metastable in respect to the molecular gas and atomic gas states, respectively.
Therefore the beyond mean-field energy functional (17) does not suffer from any complex sound velocities, in contrast to what happens in the 3D geometry. Nevertheless, it is instructive to obtain higher order corrections to the speed of sound, following a similar procedure as in the 3D case. After evaluating the compressibilities from Eqs. (5) and (17), we obtain the following expressions for the sound velocities: (18b) One can see that in 1D, higher-order corrections to the Bogoliubov speed of sound are responsible for a dynamic instability as n|a| → 0. However, a peculiarity of 1D geometry is that the mean-field regime corresponds to the large density, n|a| 1. In other words, the instability is predicted in the regime where the Bogoliubov theory is not applicable. It has been found in QMC simulations [15] that the liquid evaporates for δg/g > 0.54, in agreement with threshold value where effective interaction between dimers becomes attractive [31] while three-dimer interaction is still repulsive [32].
We show in Fig. 7(a) the equation of state of a one-dimensional Bose mixture for δg/g = 0.3 as predicted from our new theory and compared with QMC calculations from Ref. [15]. The energy is normalized by the binding energy of dimers composed from atoms from different components, B = − 2 /(ma 2 12 ). Finally, Fig. 7(b) shows the equilibrium density in the liquid phase, obtained from the minimum of the ground state energy. Surprisingly, we find that our theory essentially coincides with the QMC results up to extremely strong interactions, δg/g 0.5. By increasing δg/g, we find that the liquid phase becomes first metastable with respect to a molecular gas composed of N/2 free dimers, E(n eq ) > N ε B /2, with the threshold value (δg/g) Meta,1 = 0.356 indicated by the dashed vertical line in Fig. 7(b). Further increase in interactions make the liquid metastable in respect to the atomic gas state, E(n eq ) > 0, with the threshold value (δg/g) Meta,2 = 0.626 shown by a solid vertical line. Eventually the minimum in the energy of the liquid disappears for δg/g > 0.64 making the liquid phase mechanically unstable. The exact threshold value is δg/g > 0.54 as found from QMC [15] and few-body calculations [31]. In other words, our predictions for the equilibrium density turn out to be very precise up to almost the threshold value where the liquid state disappears.

Conclusions
In conclusion, we have developed beyond-LHY theory for the description of self-bound quantum droplet in attractive mixtures of BECs. Our theoretical approach is based on selfconsistent inclusion of higher order terms in the sound velocities, calculated in a perturbative way. The corrections brought to the speed of sound in the density channel is shown to yield a real part, in contrast to the prediction of Bogoliubov theory which is purely imaginary. The new sound velocities are used in turn to improve the energy functional. For the mixtures in three dimensions, our approach is found to describe accurately the equation of state in the liquid phase, predicting an equilibrium density in close agreement with available ab-initio calculations. We further investigate finite-size droplets by means of Gross-Pitaevskii equation, and calculate experimentally measurable quantities such as the size of the droplet and the critical number of atoms. As well, we construct beyond-Bogoliubov theory for 1D geometry. We find an excellent agreement with quantum Monte Carlo values for the equilibrium density while predictions of the Bogoliubov theory are rather imprecise. Finally, our theory predicts disappearance of the liquid phase for strong interactions.
A natural extension of this work consists in investigating the asymmetric mixture. Indeed, on-going experiments use mixtures of 39 K atoms in different hyperfine states, with g 11 = g 12 [8,10]. Recently the realization of self-bound Bose mixtures with different atomic component m 1 = m 2 has been also reported [33]. On the other hand, the theory developed in this work can also be extended for the dipolar gas [5,6], which has attracted much interest these last years, due to the richness of its phase diagram [34][35][36].