## Effective expression of the Lee-Huang-Yang energy functional for heteronuclear mixtures

- F. Minardi, 1, 2, 3 F. Ancilotto, 4, 5 A. Burchianti, 1, 2 C. D'Errico, 1, 2 C. Fort, 1, 2 and M. Modugno 6, 7 1 Istituto Nazionale di Ottica, CNR-INO, 50019 Sesto Fiorentino, Italy 2 LENS and Dipartimento di Fisica e Astronomia, Universit` di Firenze, 50019 Sesto Fiorentino, Italy a 3 Dipartimento di Fisica e Astronomia, Universit` a di Bologna, 40127 Bologna, Italy 4 Dipartimento di Fisica e Astronomia 'Galileo Galilei' and CNISM, Universit` di Padova, 35131 Padova, Italy a 5 CNR-IOM Democritos, 26534136 Trieste, Italy 6

Depto. de Fis´ ıca Te´rica e Hist. o de la Ciencia, Universidad del Pais Vasco UPV/EHU, 48080 Bilbao, Spain 7 IKERBASQUE, Basque Foundation for Science, 48013 Bilbao, Spain

(Dated: December 6, 2019)

We consider a homogeneous heteronuclear Bose mixture with contact interactions at the meanfield collapse, i.e. with interspecies attraction equal to the mean geometrical intraspecies repulsion. We show that the Lee-Huang-Yang (LHY) energy functional is accurately approximated by an expression that has the same functional form as in the homonuclear case. The approximated energy functional is characterized by two exponents, which can be treated as fitting parameters. We demonstrate that the values of these parameters which preserve the invariance under permutation of the two atomic species are exactly those of the homonuclear case. Deviations from the exact expression of LHY energy functional are discussed quantitatively and a specific application is described.

## I. INTRODUCTION

Bose-Einstein condensates (BEC) with ultracold gases are so dilute that interatomic interactions are generally weak, making the consolidated framework of mean-field (MF) theory a reliable and accurate approach in many cases. However, even with ultracold gases, there are instances where quantum fluctuations become relevant, thus producing sizeable deviations from the MF theory. Going beyond the MF analysis, N. Bogoliubov derived the many-body ground state of a weakly interacting gas of bosons [1], assuming that most particles occupy the same quantum state, i.e. belong to the condensed fraction, except for a small fraction ('quantum depletion'). The leading correction to the MF energy for identical bosons with hard-sphere interactions was later obtained by Lee-Huang-Yang (LHY) [2] in the homogeneous case. These results were then extended to the case of twocomponent bosons by D. Larsen [3].

quantum fluctuations could have striking effects, i.e. a BEC binary mixture [14]: when both components have repulsive intraspecies interactions but the interspecies interactions are attractive and sufficiently large, the MF energy may vanish or even become negative. In this case, the LHY energy correction dictates the behavior of the system, allowing for the existence of 'quantum droplets', i.e. stable states self-bound by interatomic attractive interactions. Such states exhibit the properties of a liquid even if their density is several orders of magnitude lower than in ordinary liquids [14-20]. Similarly, LHY corrections also stabilize droplets with attractive dipolar forces [21, 22].

A well-known example where quantum fluctuations are important in ultracold atomic gases is provided by optical lattices that suppress the single-particle kinetic energy and increase the interatomic interactions [4-6]. Another case is when the interatomic interactions are enhanced by increasing the scattering length in proximity of a Feshbach resonance. This route was firstly followed in fermionic condensates, where accurate measurements of the frequency shifts in collective oscillations [7] and of density profiles [8, 9] showed deviations from the MF predictions. Similar effects were later observed in strongly interacting bosonic gases in the excitation spectrum of 85 Rb [10] and in the equation of state of 7 Li [11]. The momentum distribution of the quantum (and thermal) depletion has been measured in a BEC of metastable He [12] and the total depletion density in a BEC of 39 K [13].

Recently, D. Petrov proposed a new setting where

To describe this novel system, the MF Gross-Pitaevskii equation (GPE) has been modified to incorporate the extra energy due to quantum fluctuations. For mixtures, the two-component LHY energy calculated in the homogeneous case is added to the MF energy functional to obtain a generalized GPE (g-GPE), that can be also used in the inhomogeneous case within the local density approximation [16-18, 20, 23]. While for an homonuclear mixture of equal mass components the LHY energy functional is analytically known in a closed and simple form [14], to date no analytic expression is known for a heteronuclear mixture and the LHY term must be calculated numerically [23]. This implies that the LHY energy must be tabulated for a large number of values of the ratio of two components densities n /n 2 1 , prior to the numerical integration of the g-GPE.

In this paper, guided by symmetry considerations and asymptotic behavior, we derive a simple and appealing functional form for the heteronuclear mixtures at the MF collapse, as detailed below. We verify that this function accurately approximates the numerical integral and, as such, can be efficiently used in the numerical solution of the g-GPE and to extend analytical results existing for homonuclear mixtures to the heteronuclear case.

In Section II, after summarizing known results to introduce the notation, we review the known analytical LHY energy in the special case of a largely imbalanced mixture, where the lighter species is the majority density component. In other words, this corresponds to the case of heavy impurities [3]. In Section III we propose an approximated formula for the LHY energy, which is the main result of this work, and we show its limits of applicability by comparing it with the exact numerical results of the LHY integral. In Secton IV we justify our proposed approximation by means of a Taylor expansion of the LHY integral around the homonuclear case. In Section V we show the usefulness of an analytic formula with a concrete example: we apply our analytic formula to the case of heteronuclear mixtures in the so-called LHY 'quantum fluids' regime, i.e. when the MF interaction is null and the mixture is governed by the LHY term alone.

## II. LHY ENERGY FUNCTIONAL

The LHY energy of a Bose mixture can be written as

$$E _ { \text{LHY} } = \int d r \mathcal { E } _ { \text{LHY} } ( n _ { 1 } ( r ), n _ { 2 } ( r ) ) \,, \quad ( 1 ) \quad T _ { \text{ is ol} } \\ \omega \quad \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot$$

where E LHY is the energy density given by [14, 23]

$$\mathcal { E } _ { \text{LHY} } & = \frac { 8 } { 1 5 \pi ^ { 2 } } \left ( \frac { m _ { 1 } } { h ^ { 2 } } \right ) ^ { 3 / 2 } ( g _ { 1 1 } n _ { 1 } ) ^ { 5 / 2 } f \left ( \frac { m _ { 2 } } { m _ { 1 } }, \frac { g _ { 1 2 } ^ { 2 } } { g _ { 1 1 } g _ { 2 2 } }, \frac { g _ { 2 2 } n _ { 2 } } { g _ { 1 1 } n _ { 1 } } \right ) \\ & = \kappa m _ { 1 } ^ { 3 / 2 } ( g _ { 1 1 } n _ { 1 } ) ^ { 5 / 2 } f ( z, u, x ),$$

with n i ( r )( i = 1 2) being the density of each component, , and κ ≡ 8 / (15 π 2 /planckover2pi1 3 ), z ≡ m /m 2 1 , u ≡ g 2 12 / g ( 11 g 22 ), x ≡ g 22 n / g 2 ( 11 n 1 ). The interaction strengths are proportional to the scattering lengths a 11 , a 22 , a 12 : g ii = 4 π /planckover2pi1 2 a ii /m , g i 12 = 2 π /planckover2pi1 2 a 12 /m 12 , with m 12 the reduced mass. We note that the LHY energy density is independent of the sign of g 12 , but the formation of quantum droplets occurs only for negative g 12 . Therefore, as done in Ref. [14], here we consider the case of a mixture with interspecies attractive interaction, namely g 12 &lt; 0, with g 11 &gt; 0 and g 22 &gt; 0.

For the homonuclear case, z = 1, the function is given in Ref. [14]:

$$\overset { \dots } { \text{on is} } \\ \text{one.} \quad \ f ( 1, u, x ) = \sum _ { \pm } \left ( 1 + x \pm \sqrt { ( 1 - x ) ^ { 2 } + 4 u x } \right ) ^ { 5 / 2 } / ( 4 \sqrt { 2 } ) \\ \text{so that } f ( 1 \ 1 \ r ) \ = \ ( 1 \ + \ r ) ^ { 5 / 2 } \quad \text{For the heteronymeter}$$

/negationslash so that f (1 , 1 , x ) = (1 + x ) 5 / 2 . For the heteronuclear case, z = 1, a simple expression similar to the one above is not available and approximate or numerical methods are needed.

The expression of f ( z, u, x ) for heteronuclear mixtures is obtained from Ref. [14, 23]:

$$f ( z, u, x ) = \frac { 1 5 } { 3 2 } \int _ { 0 } ^ { \infty } k ^ { 2 } \mathcal { F } ( k, z, u, x ) \, d k \,, \quad ( 4 )$$

$$\overline { \mathcal { F } ( k, z, u, x ) } & = \left \{ \frac { 1 } { 2 } \left [ k ^ { 2 } \left ( 1 + \frac { x } { z } \right ) + \frac { 1 } { 4 } k ^ { 4 } \left ( 1 + \frac { 1 } { z ^ { 2 } } \right ) \right ] + \left [ \frac { 1 } { 4 } \left [ \left ( k ^ { 2 } + \frac { 1 } { 4 } k ^ { 4 } \right ) - \left ( \frac { x } { z } k ^ { 2 } + \frac { 1 } { 4 z ^ { 2 } } k ^ { 4 } \right ) \right ] ^ { 2 } + u _ { z } ^ { x } k ^ { 4 } \right ] ^ { 1 / 2 } \right \} ^ { 1 / 2 } ( 5 ) \\ & \quad + \left \{ \frac { 1 } { 2 } \left [ k ^ { 2 } \left ( 1 + \frac { x } { z } \right ) + \frac { 1 } { 4 } k ^ { 4 } \left ( 1 + \frac { 1 } { z ^ { 2 } } \right ) \right ] - \left [ \frac { 1 } { 4 } \left [ \left ( k ^ { 2 } + \frac { 1 } { 4 } k ^ { 4 } \right ) - \left ( \frac { x } { z } k ^ { 2 } + \frac { 1 } { 4 z ^ { 2 } } k ^ { 4 } \right ) \right ] ^ { 2 } + u _ { z } ^ { x } k ^ { 4 } \right ] ^ { 1 / 2 } \right \} ^ { 1 / 2 } \\ & \quad - \frac { 1 + z } { 2 z } k ^ { 2 } - ( 1 + x ) + \frac { 1 } { k ^ { 2 } } \left [ 1 + x ^ { 2 } z + 4 u x \frac { z } { 1 + z } \right ].$$

Since the dimensionless function f ( z, u, x ) is weakly dependent on u , as shown in Ref. [14], the latter is usually replaced by its value at MF collapse, where g 12 + √ g 11 g 22 = 0, i.e. u = 1:

## A. Limit of heavy impurities

$$\mathcal { E } _ { \text{LHY} } \simeq \kappa m _ { 1 } ^ { 3 / 2 } ( g _ { 1 1 } n _ { 1 } ) ^ { 5 / 2 } f ( z, 1, x ). \quad \ ( 6 )$$

The expression of f ( z, 1 , x ) coincides with that given in Ref. [23], where it is remarked that, since the integral converges due to the cancellation of divergent terms, its evaluations requires some care.

Lacking a closed form like Eq. (3), an analytic expression of the LHY energy is given in Ref. [3] for the special case of a repulsive mixture with a 11 = a 12 = a 22 &gt; 0, where the density of the heavy component is much lower

FIG. 1. Values of f ( z, 1 , x ) as a function of x , for z = 41 87, / as obtained from numerical integration of Eq. (4)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00804v2-figure-d66f8e72-0103-415f-bf46-34f878e03e91.jpg)

than that of the light component, i.e. z &gt; 1 and x /lessmuch 1:

$$\mathcal { E } _ { \text{LHY} } & = \frac { 8 } { 1 5 \pi ^ { 2 } } \left ( \frac { m _ { 1 } } { \hslash ^ { 2 } } \right ) ^ { 3 / 2 } ( g _ { 1 1 } n _ { 1 } ) ^ { 5 / 2 } f _ { L } ( z, x ) \quad ( 7 ) \quad 1. \ In n \\ f _ { r } ( z \ r ) & = 1 + r \frac { 1 5 \frac { z + 1 } { 2 } \left [ \frac { z ^ { 2 } } { 2 } \, \text{acc} \, \sqrt { z ^ { 2 } - 1 } - 1 \right ] \quad \ \text{conve} \, } \quad \ \kappa m _ { 2 } ^ { 3 / \colon }$$

$$f _ { L } ( z, x ) = 1 + x \frac { 1 5 } { 1 6 } \frac { \lambda ^ { n } } { z - 1 } \left [ \frac { z ^ { 2 } } { \sqrt { z ^ { 2 } - 1 } } \arctan \sqrt { z ^ { 2 } - 1 } - 1 \right ] \quad \begin{array} { c } \text{conve} \\ \kappa m _ { 2 } ^ { 3 / 2 } \\ \text{Hav} \\ \text{also } r \end{array}$$

Even if the above equation was derived for all positive scattering lengths, we can also apply it to the case of -a 12 = a 11 = a 22 &gt; 0 since the LHY energy is independent of the sign of a 12 , as shown by Eq. (2). Note, however, that the condition of all equal scattering lengths corresponds to u = (1 + z ) 2 / (4 z ), thus u = 1 only for z = 1.

In the general case, however, the numerical integration of the g-GPE requires the knowledge of the LHY term for generic values of x , and the above Eq. (7) is obviously insufficient. However, we will use it to confirm the validity of our approximated formula proposed below. For the purpose, it is useful to consider f L ( z, x ) around z = 1, where it can be compared to f ( z, 1 , x ):

$$f _ { L } ( z, x ) \simeq 1 + \frac { 5 } { 2 } x + \frac { 3 } { 2 } x ( z - 1 ) + O [ z - 1 ] ^ { 2 }. \quad ( 9 ) \quad \text{while}.$$

## III. EFFECTIVE ANSATZ

We introduce here the main result of our work, namely an approximate analytic expression for f ( z, 1 , x ). The exact function, defined by Eqs. (4, 5), can only be evaluated by points, through numerically integration of the above equations. As an example, in Fig. 1 we plot f (41 / 87 1 , , x ), corresponding to the case of a 41 K87 Rb mixture.

When seeking for an approximate expression for f ( z, 1 , x ), it is useful to recall that in the case of a single component the LHY correction takes the form [2, 24]

$$\mathcal { E } _ { \text{LHY} } = \kappa m ^ { 3 / 2 } \left ( g n \right ) ^ { 5 / 2 }. \quad \quad ( 1 0 ) \quad \text{As}$$

Thus, since for x → 0 the LHY energy has to converge to the value of the single species 1, in this limit f ( z, 1 , x ) →

FIG. 2. Numerical values of integrals (points) compared to the approximate function (1 + z 3 / 5 x ) 5 / 2 (red dotted-dashed lines), for z = 0 9 0 5 . , . , 0 1. . We restrict to z ≤ 1 without loss of generality .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00804v2-figure-fe5ef6ce-406d-439c-8e23-e856772f89c0.jpg)

1. In the opposite limit, i.e. for x /greatermuch 1, the LHY energy converges to the value of the single species 2, E LHY → κm 3 / 2 2 ( g 22 n 2 ) 5 / 2 , implying f ( z, 1 , x ) → z 3 / 2 x 5 / 2 .

Having established the above asymptotic behavior, we also recall that E LHY must be invariant under permutation of the two species. It is therefore convenient to rewrite Eq. (2) in a more symmetric form with respect to the species index:

$$\text{how-} \\ \text{ngths} \\ \text{lv for} \quad \text{s} \quad. \. \. \. \. \. \. \. \. \. \.$$

Now it is easy to see that invariance under permutation of indices i = 1 2 implies ,

$$f ( z, u, x ) = z ^ { 3 / 2 } x ^ { 5 / 2 } f \left ( \frac { 1 } { z }, u, \frac { 1 } { x } \right ). \quad \ \ ( 1 2 ) \\.. \quad. \quad. \quad. \quad. \quad.$$

The above properties suggest that one could approximate f ( z, 1 , x ) with the following function

$$F _ { q } ( x, z ) = \left [ 1 + \left ( z ^ { 3 / 2 } x ^ { 5 / 2 } \right ) ^ { 1 / q } \right ] ^ { q }, \quad \ ( 1 3 ) \\ \tilde { \tilde { h } } \, \text{hac a functional form that reach]} \, \text{that of the}$$

which has a functional form that resembles that of the homonuclear case Eq. (3), except for the presence of the factor z 3 / 2 , accounting for the mass ratio, and of the exponent q (1 /q ), that is here considered as a fitting parameter. In the specific case shown of Fig. 1, a best fit to f ( z, 1 , x ) returns q = 2 4612(2) (with the relative residu-. als being below 2% for all values of 0 01 . &lt; x &lt; 100). In general the q parameter returned by the fit depends on z , but for a broad range of z values it is very close to the value q = 5 2, which characterizes the function / f in the homonuclear mixture, Eq. (3) [25]. Thus, we find convenient to approximate f ( z, 1 , x ) with the following, simple generalization of the expression for the homonuclear case:

$$f \left ( z, 1, x \right ) \simeq F _ { 5 / 2 } ( z, x ) \equiv ( 1 + z ^ { 3 / 5 } x ) ^ { 5 / 2 }. \quad ( 1 4 )$$

As shown in Fig. 2, this form, in spite of its simplicity, provides a very accurate representation of the exact (numerical) values of f ( z, 1 , x ), in a wide range of mass

FIG. 3. Top: plot of F 5 / 2 ( z, x ) (solid line) and f (3) ( z, x ) (dashed line), as defined in text. Bottom: relative deviation | f (3) ( z, x ) -F 5 / 2 ( z, x ) | /f (3) ( z, x ) (percentage). Here z = 41 87. /

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00804v2-figure-5c58c24c-b42a-41b7-8e9a-5b92a6041d17.jpg)

ratios. In addition, as requested, Eq. (14) is consistent with above Eq. (9) when we consider x, | z - | /lessmuch 1 1.

Finally, by combining Eq. (6) and Eq. (14), we arrive at the following effective expression for the LHY energy functional

$$\mathcal { E } _ { \text{LHY} } \simeq \frac { 8 m _ { 1 } ^ { 3 / 2 } ( g _ { 1 1 } n _ { 1 } ) ^ { 5 / 2 } } { 1 5 \pi ^ { 2 } \hslash ^ { 3 } } \left [ 1 + \left ( \frac { m _ { 2 } } { m _ { 1 } } \right ) ^ { 3 / 5 } \frac { g _ { 2 2 } n _ { 2 } } { g _ { 1 1 } n _ { 1 } } \right ] ^ { 5 / 2 }, \ ( 1 5 ) \quad \text{(x, t} }$$

which constitutes one of the main results of this work.

## IV. PERTURBATIVE EXPANSION

In order to quantify the deviation of the exact f ( z, 1 , x ) from F 5 / 2 ( z, x ), we note that the latter is only function of ξ ≡ z 3 / 5 x (and not of x and z separately). Then we define the auxiliary function

$$g ( z, \xi ) \equiv f ( z, 1, z ^ { - 3 / 5 } \xi ) ; \quad \quad ( 1 6 ) \quad \text{for} \\ \text{plic}$$

if f ( z, 1 , x ) was a function of ξ alone, then g z, ξ ( ) would be independent on z . Although this is not the case in general, close to z = 1, g z, ξ ( ) is almost independent of z. To show this, we write the Taylor expansion:

$$g ( z, \xi ) = \sum _ { n = 0 } ^ { + \infty } G _ { n } ( \xi ) ( z - 1 ) ^ { n }, \quad \ \ ( 1 7 ) \quad \ e r o i m u _ { \substack { \text{$\text{$mul$} \\ \text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$
$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$                                                                                                                                                                                                        
    g ( z, \xi ) = \sum _ { n = 0 } ^ { + \infty } G _ { n } ( \xi ) ( z - 1 ) ^ { n }, \text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{	$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text	$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\test$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text
$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text{$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text'
$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text'                                                                                                                                                                                                        }                                                                                                                                                                                                        g ( z, \xi ) = \sum _ { n = 0 } ^ { + } G _ { n } ( \xi ) ( z - 1 ) ^ { n }, \text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{#\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{#\text{$\text{$\text{$\text{$\text{$$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$$

and compute analytically the coefficients G n from the integral Eq. (5)

$$G _ { 0 } ( \xi ) = ( 1 + \xi ) ^ { 5 / 2 }$$

$$G _ { 1 } ( \xi ) = 0$$

$$\Phi _ { 2 } ( \xi ) & = - \frac { 1 2 } { 3 5 } \xi \sqrt { 1 + \xi } \quad \quad \ ( 2 0 ) \\ G _ { 3 } ( \xi ) & = \frac { 4 \xi ( 4 7 \xi + 4 3 ) } { \lambda \swarrow \sqrt { \sum _ { \dots } } }.$$

$$G _ { 3 } ( \xi ) = \frac { \stackrel { \cdot \cdot \cdot } { 4 \xi ( 4 7 \xi + 4 3 ) } } { 5 2 5 \sqrt { \xi + 1 } }. \quad \quad ( 2 1 )$$

The vanishing first derivative, G 1 ( ξ ) = 0, confirms that g z, ξ ( ) depends weakly on z around z = 1. Therefore,

$$\begin{array} { l l } \frac { y ( x, y ) \, \exp e n u s \, \text{warn} \, y \, \text{on} \, \text{or} u u \, z - 1. \, \text{inereive}, \\ \infty \\ 1 0 0 \\ f ( z, 1, x ) = ( 1 + z ^ { 3 / 5 } x ) ^ { 5 / 2 } - \frac { 1 2 } { 3 5 } z ^ { 3 / 5 } x \sqrt { 1 + z ^ { 3 / 5 } x } ( z - 1 ) ^ { 2 } \\ \omega _ { i } \\ \omega _ { i } \, \text{dev} _ { i } \, & + \frac { 4 z ^ { 3 / 5 } x ( 4 7 z ^ { 3 / 5 } x + 4 3 ) } { 5 2 5 \sqrt { z ^ { 3 / 5 } x + 1 } } ( z - 1 ) ^ { 3 } + O \left ( ( 1 - z ) ^ { 4 } \right ) \\ \text{Here} \quad & \equiv f ^ { ( 3 ) } ( z, x ) + O \left ( ( z - 1 ) ^ { 4 } \right ) \\ \quad & \text{The above equation shows that our ansatz function} \\ F _ { 5 / 2 } ( z, x ) \, \text{is the zeroth-order term of the Taylor expan} \end{array}$$

The above equation shows that our ansatz function F 5 / 2 ( z, x ) is the zeroth-order term of the Taylor expansion around z = 1. A comparison between F 5 / 2 ( z, x ) and the third-order truncated Taylor series, f (3) ( z, x ), is shown in Fig. 3 for z = 41 87. / The lower panel shows that the first term, F 5 / 2 ( z, x ), approximates f ( z, 1 , x ) better for small and large values of x , as expected since it is determined to match the asymptotic behaviours. Anyway, the maximum relative deviation is below 4%, proving that F 5 / 2 ( z, x ) suffices for a reliable estimation of the LHY integral.

## V. A LEE-HUANG-YANG FLUID WITH HETERONUCLEAR MIXTURES

Recently the so called 'Lee-Huang-Yang fluid' has been proposed, consisting in a dilute quantum mixture where mutual interactions and atom numbers are tuned such that the MF interactions cancel out, and the resulting system is governed only by quantum fluctuations [26, 27]. This system has been studied in the case of a homonuclear mixture ( m 1 = m 2 ) in Ref. [26], where explicit expressions for important quantities such as, e.g., the healing length and the frequency of monopole oscillation frequency, are obtained.

To prove the usefulness of our analytical, albeit approximate, form for the LHY energy discussed in Section III, we provide here a generalization of results derived in Ref. [26] to the experimentally relevant case of heteronuclear mixtures, by explicitly using our ansatz formula Eq. (14).

The total energy of a Bose-Bose binary mixture, with densities n , n 1 2 and interacting via contact potentials, is given by

$$E = \int d { \mathbf r } \left \{ \sum _ { i = 1, 2 } \left [ \frac { \hslash ^ { 2 } } { 2 m _ { i } } | \nabla \psi _ { i } | ^ { 2 } + V _ { i } | \psi _ { i } | ^ { 2 } + \frac { 1 } { 2 } g _ { i i } | \psi _ { i } | ^ { 4 } \right ] + g _ { 1 2 } | \psi _ { 1 } | ^ { 2 } | \psi _ { 2 } | ^ { 2 } + \mathcal { E } _ { L H Y } \right \} \\ & \quad \text{[}$$

where ψ i and V i are the wavefunction and the external potential of the i -th component ( i = 1 2), respectively, , and E LHY is the LHY energy density contribution Eq. (2). We assume in the following a symmetric harmonic trap, V 1 ( r ) = V 2 ( r ) = m ω r / 1 2 1 2 2 ≡ V ho . The associated harmonic length is a ho = /planckover2pi1 /m ω 1 1 .

(namely, the MF collapse condition), and atom numbers such that N /N 2 1 = √ g 11 /g 22 , only the LHY energy term remains in Eq. (23), where the two component wavefunctions satisfy the condition ψ 2 = ψ 1 ( g 11 /g 22 ) 1 / 4 .

√ The atom numbers and interactions can be tuned to cancel completely the mean-field terms in the above energy functional. In fact, by choosing g 12 = - √ g 11 g 22

By defining | ψ | 2 = | ψ 1 | 2 + | ψ 2 | 2 and α = √ g 11 /g 22 = z 1 / 2 √ a 11 /a 22 , one has | ψ 1 | 2 = (1+ α ) -1 | ψ | 2 and | ψ 2 | 2 = α (1+ α ) -1 | ψ | 2 . Thus the total energy functional Eq. (23) reduces to:

$$E = \int d \mathbf r \left [ \frac { \hslash ^ { 2 } } { 2 } | \nabla \psi | ^ { 2 } \left ( \frac { 1 } { m _ { 1 } ( 1 + \alpha ) } + \frac { \alpha } { m _ { 2 } ( 1 + \alpha ) } \right ) + V _ { h o } | \psi | ^ { 2 } + \mathcal { E } _ { L H Y } \right ] = \int d \mathbf r \left [ \frac { \hslash ^ { 2 } } { 2 m ^ { * } } | \nabla \psi | ^ { 2 } + V _ { h o } | \psi | ^ { 2 } + \mathcal { E } _ { L H Y } \right ] \quad ( 2 4 )$$

where

$$m ^ { * } = m _ { 1 } \left [ \frac { 1 + \sqrt { z a _ { 1 1 } / a _ { 2 2 } } } { 1 + \sqrt { a _ { 1 1 } / z a _ { 2 2 } } } \right ]. \quad ( 2 5 ) \quad \begin{array} { c } \text{enc} \\ C _ { L f } \\ \sqrt { n } \end{array} \\ \text{by approximating} \, \mathcal { E } _ { \text{LHY} } \text{ by the effective expression} \quad \text{v}$$

Then, by approximating E LHY by the effective expression in Eq. (15), we can write homonuclear LHY fluid of Ref. [26], the only difference being that the coefficient 256 √ π /planckover2pi1 2 | a 12 | 5 / 2 / (15 m ) ≡ C LHY is replaced by C ∗ LHY , m by m ∗ , and ω 1 by ω ∗ ≡ √ m /m ω 1 ∗ 1 .

$$\mathcal { E } _ { \text{LHY} } & = \frac { 8 } { 1 5 \pi ^ { 2 } } \left ( \frac { m _ { 1 } } { h ^ { 2 } } \right ) ^ { 3 / 2 } \left ( \frac { g _ { 1 1 } n } { 1 + \alpha } \right ) ^ { 5 / 2 } \left ( 1 + z ^ { 3 / 5 } \alpha ^ { - 1 / 2 } \right ) ^ { 5 / 2 } \\ & = C _ { \text{LHY} } ^ { * } | \psi | ^ { 5 } & ( 2 6 ) \\ \text{with} =$$

with

$$C _ { \text{LHY} } ^ { * } \equiv \frac { 2 5 6 \sqrt { \pi } } { 1 5 } \frac { \hslash ^ { 2 } } { m _ { 1 } } \left [ a _ { 1 1 } \frac { 1 + z ^ { 1 / 1 0 } \sqrt { a _ { 2 2 } / a _ { 1 1 } } } { 1 + \sqrt { z a _ { 1 1 } / a _ { 2 2 } } } \right ] ^ { 5 / 2 } \cdot \ ( 2 7 ) & & \text{and t} \\ \text{The expression in Eq. (26) makes the above Eq. (24) } \quad \text{and t}$$

With the above replacements the results of Ref. [26] can be easily generalized to the case of different masses, z = 1. For instance, we readily obtain the healing length in the homogeneous case, V ho ( r ) = 0,

/negationslash

$$\xi _ { L H Y } ^ { - 2 } = \frac { 2 5 6 \sqrt { \pi } n ^ { 3 / 2 } a _ { 1 1 } ^ { 5 / 2 } } { 3 } \left [ \frac { 1 + \sqrt { z a _ { 1 1 } / a _ { 2 2 } } } { 1 + \sqrt { a _ { 1 1 } / z a _ { 2 2 } } } \right ] \times \\ \left [ \frac { 1 + z ^ { 1 / 1 0 } \sqrt { a _ { 2 2 } / a _ { 1 1 } } } { 1 + \sqrt { z a _ { 1 1 } / a _ { 2 2 } } } \right ] ^ { 5 / 2 }, \\ \text{and the frequency $\omega$ of monopole collective oscillation in}$$

The expression in Eq. (26) makes the above Eq. (24) formally identical to the total energy functional for the

√ and the frequency ω of monopole collective oscillation in the weakly interacting limit, i.e. N 3 / 2 | a 12 /a ho | 5 / 2 /lessmuch 1 (where N = N 1 + N 2 is the total atom number):

$$\omega / \omega ^ { * } = 2 + \frac { 6 4 \sqrt { 2 } N ^ { 3 / 2 } } { 5 \sqrt { 5 } \, \pi ^ { 7 / 4 } } \left ( \frac { a _ { 1 1 } } { a _ { h o } } \right ) ^ { 5 / 2 } \left [ \frac { 1 + \sqrt { z a _ { 1 1 } / a _ { 2 2 } } } { 1 + \sqrt { a _ { 1 1 } / z a _ { 2 2 } } } \right ] ^ { 1 3 / 8 } \left [ \frac { 1 + z ^ { 1 / 1 0 } \sqrt { a _ { 2 2 } / a _ { 1 1 } } } { 1 + \sqrt { z a _ { 1 1 } / a _ { 2 2 } } } \right ] ^ { 5 / 2 }.$$

As expected, in the case of equal masses z = 1 the above expressions reduce to the corresponding expressions (10) and (15) of Ref. [26].

In order to estimate the effect of mass imbalances on the properties of the LHY fluid, we evaluate the deviation of the monopole oscillation frequency from the case of ideal gases, defined as ( ω/ω ∗ -2) /N 3 / 2 ≡ B , for various heteronuclear Bose mixtures which are currently realized in experiments of tunable dual-species condensates, as summarized in Table I. Experimentally, the frequencies of collective modes oscillations are measured with high accuracy and have been used throughout as a powerful

| Species 1,2   | Species 1,2   |   a 11 /a 0 |   a 22 /a 0 |   z |   10 7 B | n 0 /N 4 / 13 [cm - 3   | Ref.      |
|---------------|---------------|-------------|-------------|-----|----------|-------------------------|-----------|
| 39 K          | 39 K          |        86   |        34   | 1   |     2.4  | 4 . 0 · 10 13           | [16],[18] |
| 23 Na         | 39 K          |        52   |        13   | 1.7 |     0.32 | 8 . 9 · 10 13           | [28]      |
| 41 K          | 87 Rb         |        65   |       100.4 | 2.1 |     8.2  | 3 . 1 · 10 13           | [20]      |
| 39 K          | 87 Rb         |        10   |       100.4 | 2.2 |     0.76 | 7 . 9 · 10 13           | [29]      |
| 23 Na         | 87 Rb         |        54.5 |       100.4 | 3.8 |     5    | 3 . 0 10 13             | [30]      |

·

TABLE I. Bose mixtures experimentally realized and amenable to collapse at g 12 = - √ g 11 g 22 . For the purpose of comparison, we assume the same harmonic potential for all species, corresponding to an harmonic frequency of 100 Hz for 87 Rb, i.e. ω / j (2 π ) = 100 √ m Rb /m j Hz for species j .

tool to study many-body dynamics. With typical atom numbers of N ∼ 10 , the deviations 4 ω/ω ∗ -2 should be readily measurable.

In the Thomas-Fermi regime, i.e. the opposite limit N 3 / 2 | a 12 /a ho | 5 / 2 /greatermuch 1, we obtain the total peak density from Eq. (9) of Ref. [26]:

$$n _ { 0 } & = \frac { A ^ { 4 / 3 } 3 ^ { 2 / 3 } } { 3 2 ( 2 \pi ) ^ { 1 / 3 } } \left ( \frac { N } { a _ { h o } ^ { * 6 } } \right ) ^ { 4 / 1 3 } \left ( \frac { 1 5 m ^ { * } C _ { L H Y } ^ { * } } { 2 5 6 \sqrt { \pi } \hslash ^ { 2 } } \right ) ^ { - 6 / 1 3 } \\ & \simeq 0. 0 7 8 \left ( \frac { N } { a _ { h o } ^ { 6 } } \right ) ^ { 4 / 1 3 } \left [ a _ { 1 1 } \frac { 1 + z ^ { 1 / 1 0 } \sqrt { a _ { 2 2 } / a _ { 1 1 } } } { 1 + \sqrt { z a _ { 1 1 } / a _ { 2 2 } } } \right ] ^ { - 1 5 / 1 3 } \quad \text{We} \\ & \quad \text{for co} \\ & \quad \text{carefu} \\ & \quad \text{Unive} \\ & \quad \text{Fund}$$

where A /similarequal 1 815. . This quantity is relevant for the lifetime of the atomic LHY quantum fluid that is limited by 3-body inelastic collisions, whose rate is proportional to n 2 0 . In Table I we report also the peak densities evaluated

- [1] N. Bogolubov, J. Phys. (USSR) 11 , 23 (1947).
- [2] T. D. Lee, K. Huang, and C. N. Yang, Phys. Rev. 106 , 1135 (1957).
- [3] D. M. Larsen, Annals of Physics 24 , 89 (1963).
- [4] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81 , 3108 (1998).
- [5] M. Greiner, O. Mandel, T. Esslinger, T. W. H¨ ansch, I. Bloch, and T. W. Ha, Nature 415 , 39 (2002).
- [6] K. Xu, Y. Liu, D. E. Miller, J. K. Chin, W. Setiawan, and W. Ketterle, Phys. Rev. Lett. 96 , 180405 (2006).
- [7] A. Altmeyer, S. Riedl, C. Kohstall, M. J. Wright, R. Geursen, M. Bartenstein, C. Chin, J. H. Denschlag, and R. Grimm, Phys. Rev. Lett. 98 , 040401 (2007).
- [8] Y.-i. Shin, A. Schirotzek, C. H. Schunck, and W. Ketterle, Phys. Rev. Lett. 101 , 070404 (2008).
- [9] N. Navon, S. Nascimb` ene, F. Chevy, and C. Salomon, Science 328 , 729 (2010).
- [10] S. B. Papp, J. M. Pino, R. J. Wild, S. Ronen, C. E. Wieman, D. S. Jin, and E. A. Cornell, Phys. Rev. Lett. 101 , 135301 (2008).

for the different Bose mixtures.

## VI. CONCLUSIONS

We have shown that the Lee-Huang-Yang energy functional for a heteronuclear Bose mixture at the MF collapse can be accurately approximated by an expression that has a functional form similar to the simple expression of the homonuclear case, see Eq. (15). Two different ansatz functions, with one or two exponents as fitting parameters, have been considered [25]. Remarkably, we found that the values of these fitting parameters which preserve the invariance under permutation of the two species are exactly those of the homonuclear case. A quantitative analysis of the deviations from the exact expression of LHY energy functional indicates that the simple expression we propose is effective in a wide range of mass ratios, and may be useful for describing current and future experiments. Our formula greatly simplifies the numerical integration of the generalized GPE's and allows further analytic study, as demonstrated in the case of a Lee-Huang-Yang fluid with heteronuclear mixtures.

## ACKNOWLEDGMENTS

We are grateful to M. Prevedelli and L. Salasnich for comments and suggestions, and to L. Marmugi for carefully reading the manuscript. We acknowledge support by the Spanish Ministry of Science, Innovation and Universities and the European Regional Development Fund FEDER through Grant No. PGC2018-101355-BI00 (MCIU/AEI/FEDER,UE), the Basque Government through Grant No. IT986-16, and the Fondazione CR Firenze through project 'SUPERACI-Superfluid Atomic Circuits'.

- [11] N. Navon, S. Piatecki, K. G¨ unter, B. Rem, T. C. Nguyen, F. Chevy, W. Krauth, and C. Salomon, Phys. Rev. Lett. 107 , 135301 (2011).
- [12] R. Chang, Q. Bouton, H. Cayla, C. Qu, A. Aspect, C. I. Westbrook, and D. Cl´ ement, Phys. Rev. Lett. 117 , 235303 (2016).
- [13] R. Lopes, C. Eigen, N. Navon, D. Cl´ment, R. P. Smith, e and Z. Hadzibabic, Phys. Rev. Lett. 119 , 190404 (2017).
- [14] D. S. Petrov, Phys. Rev. Lett. 115 , 155302 (2015).
- [15] A. Cappellaro, T. Macr` ı, G. F. Bertacco, and L. Salasnich, Scientific Reports 7 , 13358 (2017).
- [16] C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell, Science 359 , 301 (2018).
- [17] P. Cheiney, C. R. Cabrera, J. Sanz, B. Naylor, L. Tanzi, and L. Tarruell, Phys. Rev. Lett. 120 , 135301 (2018).
- [18] G. Semeghini, G. Ferioli, L. Masi, C. Mazzinghi, L. Wolswijk, F. Minardi, M. Modugno, G. Modugno, M. Inguscio, and M. Fattori, Phys. Rev. Lett. 120 , 235301 (2018).
- [19] G. Ferioli, G. Semeghini, L. Masi, G. Giusti, G. Modugno, M. Inguscio, A. Gallem´ ı, A. Recati, and M. Fat-

tori, Phys. Rev. Lett. 122 , 090401 (2019).

- [20] C. D'Errico, A. Burchianti, M. Prevedelli, L. Salasnich, F. Ancilotto, M. Modugno, F. Minardi, and C. Fort, arXiv:1908.00761 [cond-mat.quant-gas].
- [21] I. Ferrier-Barbut, H. Kadau, M. Schmitt, M. Wenzel, and T. Pfau, Phys. Rev. Lett. 116 , 215301 (2016).
- [22] L. Chomaz, S. Baier, D. Petter, M. J. Mark, F. W¨chtler, a L. Santos, and F. Ferlaino, Phys. Rev. X 6 , 041039 (2016).
- [23] F. Ancilotto, M. Barranco, M. Guilleumas, and M. Pi, Phys. Rev. A 98 , 053623 (2018).
- [24] F. Dalfovo, S. Giorgini, and S. Stringari, Rev. Mod. Phys. 71 , 463 (1999).
- [25] We arrive to the same results considering the slightly different trial function (1+ z α x ) β and, (i) either determining α and β exponents as two independent fitting parameters, in this case the fit values of α, β depend slightly on the mass ratio z ; (ii) or requiring invariance under permutation of the two species, that fixed α = 3 5 and / β = 5 2, / independent of z .
- [26] N. B. Jørgensen, G. M. Bruun, and J. J. Arlt, Phys. Rev. Lett. 121 , 173403 (2018).
- [27] The vanishing of the mean-field term could even make possible the investigation of beyond-LHY corrections in such a system, as suggested in Ref. [26]. Their quantitative characterization would require to extend the single component results [31] (for example using the approach of [32]) to the case of a two-component mixture.
- [28] T. A. Schulze, T. Hartmann, K. K. Voges, M. W. Gempel, E. Tiemann, A. Zenesini, and S. Ospelkaus, Phys. Rev. A 97 , 023623 (2018).
- [29] L. Wacker, N. B. Jørgensen, D. Birkmose, R. Horchani, W. Ertmer, C. Klempt, N. Winter, J. Sherson, and J. J. Arlt, Phys. Rev. A 92 , 053602 (2015).
- [30] F. Wang, X. Li, D. Xiong, and D. Wang, J. Phys. B: At. Mol. Opt. Phys. 49 , 015302 (2016).
- [31] A. Fetter and J. Walecka, Quantum Theory of ManyParticle Systems (Dover Publications, New York, 1971).
- [32] R. S. Christensen, J. Levinsen, and G. M. Bruun, Phys. Rev. Lett. 115 , 160401 (2015).
F. Minardi F. Ancilotto A. Burchianti C. D'Errico C. Fort M. Modugno
Effective expression of the Lee-Huang-Yang energy functional for heteronuclear mixtures
We consider a homogeneous heteronuclear Bose mixture with contact
interactions at the mean-field collapse, i.e. with interspecies attraction
equal to the mean geometrical intraspecies repulsion. We show that the
Lee-Huang-Yang (LHY) energy functional is accurately approximated by an
expression that has the same functional form as in the homonuclear case. The
approximated energy functional is characterized by two exponents, which can be
treated as fitting parameters. We demonstrate that the values of these
parameters which preserve the invariance under permutation of the two atomic
species are exactly those of the homonuclear case. Deviations from the exact
expression of LHY energy functional are discussed quantitatively and a specific
application is described.