# Dilute quantum liquid in a K-Rb Bose mixture

V. Cikojević  University of Split, Faculty of Science, Ruđera Boškovića 33,
HR-21000 Split, Croatia  Departament de Física, Universitat Politècnica de
Catalunya, Campus Nord B4-B5, E-08034 Barcelona, Spain  E. Poli  Institut für
Experimentalphysik, Universität Innsbruck, Austria  Dipartimento di Fisica e
Astronomia "Galileo Galilei" and CNISM, Universita di Padova, via Marzolo 8,
35122 Padova, Italy  F. Ancilotto  Dipartimento di Fisica e Astronomia
"Galileo Galilei" and CNISM, Universita di Padova, via Marzolo 8, 35122
Padova, Italy  L. Vranješ-Markić  University of Split, Faculty of Science,
Ruđera Boškovića 33, HR-21000 Split, Croatia  J. Boronat  Departament de
Física, Universitat Politècnica de Catalunya, Campus Nord B4-B5, E-08034
Barcelona, Spain

###### Abstract

A quantum liquid in a heterogeneous mixture of 41K and 87Rb atoms is studied
using the diffusion Monte Carlo method and Density Functional Theory. The
perturbative Lee-Huang-Yang term for a heterogeneous mixture is verified and
it is proved to be valid only near the gas-liquid transition. Based on the
equations of state of the bulk mixture, calculated with diffusion Monte Carlo,
extensions to Lee-Huang-Yang corrected mean-field energy functionals (MF+LHY)
are presented. Using Density Functional Theory, a systematic comparison
between different functionals is performed, focusing on the critical atom
number, surface tension, surface width, Tolman length, and compressibility.
These results are given as a function of the inter-species interaction
strength, within the stability domain of the liquid mixture.

##  I Introduction

A new quantum state of matter has been predicted [1], and experimentally
realized [2, 3, 4] recently in ultracold atomic gases, where a subtle
interplay between inter-species attractive interactions and quantum
fluctuations may result in the formation of self-bound, ultradilute droplets.
Such liquid droplets are fundamentally different from those in classical or
Helium fluids, where they arise instead from the interplay between the short-
range repulsive and long-range attractive components of the interatomic
potential [5]. The existence of self-bound ultra dilute quantum droplets, made
of atoms of a binary mixture of Bose-Einstein condensates, was predicted by
Petrov [1] and has been experimentally confirmed shortly later [2, 3]. These
systems, whose peculiar properties are shared by other systems like dipolar
Bose gases [6, 7, 8], are characterized by ultralow densities, orders of
magnitude lower than that of the prototypical quantum liquid, i.e., liquid
helium. Noticeably, quantum droplets exist at temperatures that are several
orders of magnitude lower than the freezing points of classical liquids.

Mean-field analysis predicts that binary mixtures of Bose-Einstein condensates
become unstable against collapse when the attractive inter-species interaction
overcomes the repulsive contact potential between atoms [9]. However, in the
utradilute liquid phase the mean-field collapse is avoided if beyond-mean-
field first-order perturbative corrections, in the form of the Lee Huang Yang
energy functional [10, 11], are included. This correction is repulsive in
nature, and thus stabilizes the system.

The formation of heteronuclear quantum droplets in an attractive bosonic
mixture of 41K and 87Rb has been observed recently [4]. At variance with the
largely studied mixture of two hyperfine states of 39K, longer lived self-
bound states are observed in the K-Rb mixture, both in free space and in
optical waveguides. The K-Rb mixture has proven to be robust even when the two
components are exposed to different confining potentials. Such long-lived,
self-bound droplets remain localized on a timescale of several tens of
milliseconds, more than a factor of 10 larger than in the 39K mixture [3].

The increased lifetime of this new liquid mixture allows not only for a more
detailed experimental observation and characterization of isolated droplets,
but will also permit the observation of more complex scenarios arising from
the interactions between self-bound droplets. Moreover, the observed reduction
of three-body losses allows for the realization of larger droplets. Recently,
collisions between two droplets have been proposed as a useful experimental
tool to investigate the dynamical properties of self-bound systems [12]. When
two such droplets approach each other with a given relative velocity, they can
either merge in a single droplet (coalescence) or separate into two or more
droplets after the collision (bouncing/fragmentation) [13]. These different
outcomes depend on whether or not the surface tension is large enough to
counterbalance the kinetic energy of the colliding drops. A different
phenomenology is expected for the coalescence dynamics of two droplets
colliding at very small velocities, which, in analogy to previous studies on
helium clusters, could be a probe of their superfluid properties. One
interesting outcome of collision would be the formation of vortices or other
topological structures during the merging process, as a consequence of their
condensate nature. These effects are known to arise during the merging of
superfluid liquid helium nanodroplets [14].

A fruitful comparison with experiments needs accurate and reliable theoretical
schemes, which have not been applied yet on K-Rb mixture. The use of Density
Functional Theory (DFT), in its time-dependent version, is known to allow a
quite accurate description of the dynamics of inhomogeneous superfluid
systems, even at the level of the Local Density approximation. Finite-range
effects permit to widen the applicability of numerical simulations by going
beyond the usual mean-field theory corrected with the LHY term (MF+LHY). A
most natural way to include such effects within the DFT scheme is to use
density functionals built from results of first-principles quantum Monte Carlo
calculations for the homogeneous phase [15, 16]. Importantly, first-principles
finite-$N$ QMC calculations in the low-density regime show good agreement with
the DFT approach [17] only in the vicinity of the gas-liquid transition while,
at larger densities, modifications to MF+LHY approach are necessary [18, 19,
20, 21].

This work is organized as follows. In Sec. II, we introduce the MF+LHY
functional for a K-Rb mixture and predict a range of scattering lengths which
allow self-bound state. In Sec. III, we describe the diffusion Monte Carlo
(DMC) methodology for obtaining the ground-state energies of the bulk mixture.
We report the DMC energies obtained using the short-range and finite-range
potentials in Sec. IV and Sec. V, respectively. In Sec. VI, we apply the
density-functional formalism to the characterization of the droplets,
comparing the MF+LHY and QMC-based functionals for each quantity. Namely, we
provide results for the surface tension, density profile, critical atom
number, surface width, compressibility, and Tolman length for an
experimentally relevant range of scattering parameters. Finally, a summary of
results and conclusions is given in Sec. VII.

##  II MF + LHY equation of state

We consider a homogeneous, heteronuclear Bose-Bose mixture with two components
(with masses $m_{1}$ and $m_{2}$) in a volume $V$ and a total number of bosons
$N=N_{1}+N_{2}$. By neglecting finite-range effects (which will be considered
in Sec. III), only the $s$-wave scattering lengths are used to characterize
the inter-particle interactions. Within the LHY-extended mean-field framework
(abbreviated MF+LHY hereafter), the energy of the system per unit volume is
given by the functional [1]

$\mathcal{E}=\mathcal{E}_{\rm MF}+\mathcal{E}_{\rm LHY}\ ,$ (1)

where the mean-field (MF) and Lee-Huang-Yang (LHY) terms read [9, 22, 23]

$\mathcal{E}_{\rm
MF}=\dfrac{1}{2}g_{11}\rho_{1}^{2}+\dfrac{1}{2}g_{22}\rho_{2}^{2}+g_{12}\rho_{1}\rho_{2}\
,$ (2) $\mathcal{E}_{\rm
LHY}=\dfrac{8m_{1}^{3/2}\left(g_{11}\rho_{1}\right)^{5/2}}{15\pi^{2}\hbar^{3}}\left[1+\left(\dfrac{m_{2}}{m_{1}}\right)^{3/5}\dfrac{g_{22}\rho_{2}}{g_{11}\rho_{1}}\right]^{5/2}\
,$ (3)

Here, $\rho_{i}$, ($i=1,2$) are the number densities of each component of the
mixture, normalized such that $\int_{V}\rho_{i}\,{\rm d}{\bf r}=N_{i}$, and
$g_{ij}=2\pi\hbar^{2}a_{ij}/\mu_{ij}$ are the $ij$-interaction strengths, with
$\mu_{ij}^{-1}=m_{i}^{-1}+m_{j}^{-1}$ the reduced mass.

The mixture of the two species is stable against fluctuations in the
concentration $N_{1}/N_{2}$ if [1]

$\dfrac{\rho_{2}}{\rho_{1}}=\sqrt{\dfrac{g_{11}}{g_{22}}}\ .$ (4)

As pointed out in Refs. [1, 22, 24], it is safe to assume that this optimal
composition is realized everywhere in the system. Thus, the energy functional
$\mathcal{E}_{\rm MF+LHY}[\rho]=\mathcal{E}_{\rm MF}[\rho]+\mathcal{E}_{\rm
LHY}[\rho]$ becomes effectively single-component, and can be written in terms
of the total density $\rho=\rho_{1}+\rho_{2}$. Under this assumption, the MF
and LHY terms read

$\displaystyle\mathcal{E}_{\rm MF}$ $\displaystyle=$
$\displaystyle\frac{2\pi\hbar^{2}\left(2a_{11}m_{2}+a_{12}\sqrt{\frac{a_{11}m_{2}}{a_{22}m_{1}}}\left(m_{1}+m_{2}\right)\right)}{m_{1}m_{2}\left(\sqrt{\frac{a_{11}m_{2}}{a_{22}m_{1}}}+1\right)^{2}}\rho^{2}\hskip
5.69046pt$ (5) $\displaystyle\mathcal{E}_{\mathrm{LHY}}$ $\displaystyle=$
$\displaystyle\frac{256\sqrt{\pi}\hbar^{2}a_{11}^{5/2}}{15m_{1}}$ (6)
$\displaystyle\times$
$\displaystyle\left[\frac{1+\left(\dfrac{m_{2}}{m_{1}}\right)^{{1}/{10}}\sqrt{a_{22}/a_{11}}}{1+\sqrt{\dfrac{m_{2}a_{11}}{m_{1}a_{22}}}}\right]^{{5}/{2}}\rho^{5/2}$

The MF+LHY energy per particle can be compactly written as

$\dfrac{E/N}{|E_{0}^{\rm MF+LHY}|/N}=-3\dfrac{\rho}{\rho_{0}^{\rm
MF+LHY}}+2\left(\dfrac{\rho}{\rho_{0}^{\rm MF+LHY}}\right)^{3/2}\ ,$ (7)

$\rho_{0}^{\rm MF+LHY}$ being the equilibrium density within the MF+LHY
theory, i.e., the density that minimizes the functional $\mathcal{E}_{\rm
MF+LHY}[\rho]$. Explicitly,

$\displaystyle\rho_{0}^{\rm MF+LHY}=$ $\displaystyle\displaystyle\frac{25\pi
m_{2}^{3}\left(2a_{11}m_{2}+a_{12}\sqrt{\frac{a_{11}m_{2}}{a_{22}m_{1}}}\left(m_{1}+m_{2}\right)\right)^{2}}{4096\left(a_{11}m_{2}+a_{22}m_{1}\left(\frac{m_{2}}{m_{1}}\right)^{\frac{3}{5}}\sqrt{\frac{a_{11}m_{2}}{a_{22}m_{1}}}\right)^{5}}$
(8)
$\displaystyle\times\left(\sqrt{\frac{a_{11}m_{2}}{a_{22}m_{1}}}+1\right).$

The energy per particle at equilibrium is $E_{0}^{\rm
MF+LHY}/N=\mathcal{E}_{\rm MF+LHY}[\rho_{0}^{\rm MF+LHY}]/\rho_{0}^{\rm
MF+LHY}$. The fact that the MF+LHY functional can be written in a universal
form (Eq. 7) proves that all the results from Ref. [1] can be applied here
with a proper change of units.

In what follows, we define the hyperfine state $\ket{F=1,m_{F}=1}$ of 41K as
component 1, and the hyperfine state $\ket{F=1,m_{F}=1}$ of 87Rb as component
2. The scattering parameters describing the intra-species repulsion are fixed
and their values are equal to $a_{11}=65a_{0}$ [25] and $a_{22}=100.4a_{0}$
[26]. With those two parameters, the MF+LHY theory predicts a self-bound state
(hereafter called "liquid") for $a_{12}<a_{12}^{c}$, with $a_{12}^{c}$ given
by

$a_{12}^{c}=\dfrac{-2\sqrt{a_{22}/a_{11}}}{\sqrt{m_{2}/m_{1}}\left(1+m_{1}/m_{2}\right)}\approx-75.4a_{0}.$
(9)

In the experiment, accessible values of $a_{12}$ are in the range between
$a_{12}=-80a_{0}$ and $a_{12}=-95a_{0}$ [4].

##  III  Diffusion Monte Carlo

We use the diffusion Monte Carlo (DMC) method to determine the energy per
particle in the homogeneous phase. This method was previously applied in
related problems regarding the study of a Bose-Bose liquid by some of the
authors [27, 15]. DMC is nowadays a well-know method that is able to solve
exactly the imaginary-time Schrödinger equation of the many-particle system,
within some statistical noise. The starting point of DMC is the decomposition
of the imaginary-time evolution operator. In this paper, we use a propagator
that is accurate up to second order in the timestep [28], following the
implementation outlined in Ref. [29].

To reduce the variance in the estimation of energy, we use standard importance
sampling through a trial wavefunction, written as a Jastrow product over pairs
[30]

$\Psi({\bm{\mathrm{R}}})=\prod_{j>i=0}^{N_{1}}f^{(11)}(r_{ij})\prod_{j>i=N_{1}}^{N_{1}+N_{2}}f^{(22)}(r_{ij})\prod_{i,j}^{N_{1},N_{2}}f^{(12)}(r_{ij})\
,$ (10)

where the two-particle correlation functions $f^{(\alpha,\beta)}(r)$
($\alpha,\beta=1,2$) are chosen as

$f^{(\alpha,\beta)}(r)=\begin{cases}f_{\rm 2b}(r)&r<\tilde{R}\\\
B(1-\frac{a_{\alpha,\beta}}{r})&r<R_{v},\\\
C\exp(-\frac{D}{r}+\frac{E}{r^{2}}),&R_{v}<r<L/2\\\ 1,&r>L/2\ .\\\
\end{cases}$ (11)

The two-particle correlation function at short distances $f_{\rm 2b}(r)$ is
the solution to the two-body problem for a given interaction potential.
Throughout the paper we use short-ranged potentials, such that the potential
is zero at distances greater than $\tilde{R}$. The function $f_{\rm 2b}$ is
connected to the asymptotic form $1-a_{\alpha,\beta}/r$, where
$a_{\alpha,\beta}$ is a corresponding scattering length. At $r=R_{v}$, this is
connected to the long-range phononic form
$C\exp(-\frac{D}{r}+\frac{E}{r^{2}})$ [31]. Finally, we impose that the
function is constant (one) at the boundary of the simulation box ($r=L/2$).
This trial wavefunction has only one variational parameter, namely $R_{v}$,
which we optimize variationally finding that, in all cases, $R_{v}=0.45L$. Our
DMC results are unbiased for time-steps $\Delta\tau\lesssim
0.5m_{41}a_{11}^{2}/\hbar^{2}$ and walker number $n_{\rm w}\approx 200$.
Simulations are performed in a cubic box of size $L=(N/\rho)^{1/3}$, with
periodic boundary conditions applied to particle coordinates. For each
density, a set of several calculations of increasing number of particles,
namely $N=130$, 160, 200, 250, and 500 is performed in order to study finite-
size effects. The final energies, corresponding to the thermodynamic limit,
are obtained by an extrapolation to $N\rightarrow\infty$ assuming a correction
which decreases as $N^{-1}$.

##  IV Short-range potentials

In a first approach to the problem, we use in DMC interatomic potentials that
have a very short range, which we call POT-SR in the following. Since we
cannot use a contact interaction in DMC, we model the short-range interaction
by a set of potentials with a range $r_{p}$ satisfying $\rho r_{p}^{3}\ll 1$,
where $\rho$ is the typical number density. Under this criterium, we choose a
hard-core potential for repulsive interactions (between equal species), with a
diameter corresponding to the $s$-wave scattering length [9]

$V_{ii}(r)=\begin{cases}\infty&r<a_{ii},\\\ 0,&\mathrm{otherwise}\\\
\end{cases},\hskip 14.22636pti=1,2.$ (12)

The attraction between different species is modeled by a short-range square-
well potential,

$V_{12}(r)=\begin{cases}-V_{\rm sr},&r<R_{\rm sr},\\\ 0,&\mathrm{otherwise}\\\
\end{cases},$ (13)

where we choose $R_{\rm sr}=a_{11}$. The particular choice of potentials given
in Eqs. (12) and (13) resembles a zero-range case, since the probability of
finding two particles within the diameter $R_{\rm sr}$ is $\rho_{0}R_{\rm
sr}^{3}\approx 2\cdot 10^{-4}$, evaluated at the equilibrium density for the
densest liquid analyzed, corresponding to $a_{12}=-95a_{0}$. By properly
setting $V_{\rm sr}$, we obtain a target scattering length $a_{12}$, which for
a square-well potential is given by [32]

$a_{12}=R_{\rm sr}\left\\{1+\dfrac{\tan\left(K_{0}R_{\rm
sr}\right)}{K_{0}R_{\rm sr}}\right\\}\ ,$ (14)

where $K_{0}=\sqrt{2V_{\rm sr}(m_{1}+m_{2})/(\hbar^{2}m_{1}m_{2})}$. DMC
results of the energy per particle, obtained with POT-SR potentials (Eqs. 12
and 13), are shown in Fig. 1, for $a_{12}=-77a_{0}$, $-85a_{0}$, $-90a_{0}$,
and $-95a_{0}$, which include experimentally accessible values. They are
compared with the predictions of MF +LHY theory. The DMC energy per particle
is well fitted to the form

$\dfrac{E}{N}=\alpha\rho+\beta\rho^{\gamma}\ ,$ (15)

where $\alpha$, $\beta$, and $\gamma$ are the fitting parameters. The
equilibrium density and the coefficients of the fit to the DMC energy per
particle are reported in Table 1. As we can see in Fig. 1, the DMC results of
the energy per particle reproduce well the MF +LHY theory for
$a_{12}=-77a_{0}$, which is close to the critical value,
$a_{12}^{c}=-75.4a_{0}$. However, the DMC results for different $a_{12}$ do
not fit a single line defined by Eq. (7), implying that the universality of
the MF +LHY theory is broken, with the deviation growing as $a_{12}$ becomes
larger. Interestingly, we observe repulsive beyond-LHY contributions when the
fluid enters a more correlated regime. Repulsive beyond-LHY contributions to
the energy have already been observed in symmetric Bose-Bose fluids [27] and
in the liquid 39K mixture [15] for the potentials with a small effective
range. This effect could be due to bosonic pairing between atoms of different
species [20, 19].

Figure 1: DMC energy per particle as a function of density, obtained with the POT-SR models (Eqs. 12 and 13), in a mixture having the optimal composition $\rho_{2}/\rho_{1}=\sqrt{g_{11}/g_{22}}$. Energy per particle and the total density are normalized with respect to the equilibrium values given in Eq. (8). Table 1: Coefficients of fit $E/N=\alpha\rho a_{11}^{3}+\beta(\rho a_{11}^{3})^{\gamma}$ of the DMC energies per particle, obtained with short-range set of potentials POT-SR (Eqs. (12) and (13), reported in Fig. (1). Parameters $\alpha$ and $\beta$ are given in units $\hbar^{2}/(m_{1}a_{11}^{2})=10^{-3}{\rm Kelvin}$, and $\gamma$ is adimensional. $\rho_{0}$ stands for the equilibrium density of a QMC functional. Scattering lengths in the repulsive channel are $a_{11}=65a_{0}$ and $a_{22}=100.4a_{0}$. $~{}a_{12}~{}$ | Method | $\alpha$ | $\beta$ | $\gamma$ | $\rho_{0}a_{11}^{3}$
---|---|---|---|---|---
$-77a_{0}$ | QMC | $-0.056$ | $35.709$ | $1.500$ | $1.08\times 10^{-6}$
| MF+LHY | $-0.057$ | $36.437$ | $1.5$ | $1.09\times 10^{-6}$
$-85a_{0}$ | QMC | $-0.385$ | $15.981$ | $1.395$ | $3.48\times 10^{-5}$
| MF+LHY | $-0.340$ | $36.437$ | $1.5$ | $3.86\times 10^{-5}$
$-90a_{0}$ | QMC | $-0.592$ | $16.344$ | $1.384$ | $7.53\times 10^{-5}$
| MF+LHY | $-0.516$ | $36.437$ | $1.5$ | $8.92\times 10^{-4}$
$-95a_{0}$ | QMC | $-0.772$ | $16.768$ | $1.379$ | $1.26\times 10^{-4}$
| MF+LHY | $-0.693$ | $36.437$ | $1.5$ | $1.61\times 10^{-4}$

##  V Finite-range potentials

The effective range in the 41K-87Rb mixture is not known, but can be estimated
from the combined knowledge of the $C_{6}$ coefficient of the leading term in
the long-range tail of the Van Der Waals interaction potential and the
scattering length obtained from measurements of the Feshbach resonances [33,
34, 35]. Knowing the van der Waals coefficient $C_{6}$, one can estimate the
effective range using a semiclassical approximation [35],

$r_{\rm
eff}=\dfrac{\Gamma(1/4)^{2}}{6\pi}a_{m}\left(1-2\dfrac{a_{m}}{a}+2\left(\dfrac{a_{m}}{a}\right)^{2}\right)\
,$ (16)

with $a_{m}={4\pi R_{\rm vdw}}/{\Gamma(1/4)^{2}}$ the mean scattering length
and $R_{\rm vdw}=\left({2\mu C_{6}}/{\hbar^{2}}\right)^{1/4}/2$ the van der
Waals length, where $\mu$ is the reduced mass. The $C_{6}$ coefficients for
41K and 87Rb are given, along with the effective ranges derived from (16), in
Table 2.

Since only two scattering parameters, namely the scattering length $a$ and the
effective range $r_{\rm eff}$, cannot uniquely define the interaction
potential, we resort to model potentials satisfying the two scattering
criteria. To investigate the role of the "shape" of the interaction dictated
by the higher-order scattering parameters, we performed two independent sets
of calculations, each having different models of the interaction potential. We
call the two models POT-I and POT-II, and both can be written as follows

$V(r)=\begin{cases}V_{0}&0<r<R_{0},\\\ -V_{1},&R_{0}<r<R_{1},\\\
0,&\mathrm{otherwise}.\\\ \end{cases}\hskip 14.22636pti=1,2.$ (17)

This particular form of interaction is convenient because the analytic
expressions for both the $s$-wave scattering length and effective range are
analytically known [36]. The specific values of the interaction parameters in
all three channels for POT-I and POT-II potentials are summarized in Table 3.

In Fig. (2), we report the comparison between the DMC and MF+LHY equations of
state, for three values of $a_{12}$ and using POT-I and POT-II as model
potentials (see Eq. (17) and Table (3)). For all three values of $a_{12}$,
there is a slight increase in the equilibrium density, relative to the MF+LHY
theory. This phenomenon was previously observed in a mixture with symmetric
interactions [27, 18, 17], whereas in a 39K liquid mixture [15] the
significant increase of equilibrium density occurs because in that mixture the
effective range is much larger than in the K-Rb one. In other words, in the K
mixture the finite-range, beyond-LHY, negative energy contributions dominate.
On the other hand, we observe that energies obtained with two different
potentials but with the same $s$-wave scattering length and effective range
collapses to a single equation of state, at least at densities not much larger
than the equilibrium one. Therefore, the range of universality is extended,
but now in terms of two scattering parameters [15].

Table 2: $C_{6}$ coefficients of each channel and the corresponding scattering lengths and effective ranges. Channel |  $C_{6}$ (a. u.) | $a/a_{0}$ | $r_{\rm eff}/a_{0}$
---|---|---|---
K-K (1-1) [37, 25] | 3897 | $65$ | $168$
Rb-Rb (2-2) [26] | 4707 | $100.4$ | $153$
K-Rb (1-2) [38] | 4285 | $-85$ | $795$
|  | $-90$ | $748$
|  | $-95$ | $707$
Table 3:  Parameters of POT-I and POT-II potentials which reproduce both scattering parameters (see Table 2). $R_{0}$ and $R_{1}$ are given in units $a_{11}=65a_{0}$, and $V_{0}$ and $V_{1}$ are given in units $\hbar^{2}/(m_{1}a_{11}^{2})=10^{-3}{\rm Kelvin}$. $a$ | Legend | $R_{0}$ | $R_{1}$ | $V_{0}$ | $V_{1}$
---|---|---|---|---|---
$65a_{0}$ | POT-I | $1.8797$ | $0$ | $1.2123$ | $0$
| POT-II | $2.8386$ | $3.5581$ | $0.9140$ | $0.3687$
$100.4a_{0}$ | POT-I | $2.7084$ | $0$ | $0.3337$ | $0$
| POT-II | $3.8374$ | $4.8201$ | $0.2366$ | $0.0824$
$-85a_{0}$ | POT-I | $0$ | $4.8897$ | $0$ | $0.01865$
| POT-II | $2.7513$ | $3.6684$ | $0.0632$ | $0.1189$
$-90a_{0}$ | POT-I | $0$ | $4.8069$ | $0$ | $0.0205$
| POT-II | $2.7419$ | $3.6558$ | $0.0586$ | $0.1197$
$-95a_{0}$ | POT-I | $0$ | $4.7315$ | $0$ | $0.0222$
| POT-II | $2.7329$ | $3.6438$ | $0.0541$ | $0.1205$

Figure 2: Energy per particle as a function of density, for values of
$a_{12}=-85a_{0}$ (a), $a_{12}=-90a_{0}$ (b) and $a_{12}=-95a_{0}$ (c). DMC
calculations are performed with POT-I and POT-II potentials, which satisfy
both scattering parameters (see Eq. (17) and table (4). Ratio of
concentrations $\rho_{2}/\rho_{1}=\sqrt{g_{11}/g_{22}}$ is kept fixed for each
density, a criteria coming from the mean-field, also verified by QMC to
predict the ground state. Table 4: Coefficients of fit $E/N=\alpha\rho
a_{11}^{3}+\beta(\rho a_{11}^{3})^{\gamma}$ of the DMC energies per particle,
obtained with set of potentials POT-I and POT-II (se Eq. (17 and Table (3)).
Comparison between different functionals are presented in Fig. (2). Parameters
$\alpha$ and $\beta$ are given in units
$\hbar^{2}/(m_{1}a_{11}^{2})=10^{-3}{\rm Kelvin}$, and $\gamma$ is
adimensional. Scattering lengths in the repulsive channel are $a_{11}=65a_{0}$
and $a_{22}=100.4a_{0}$. For comparison, MF+LHY parameters are written in
Table (1).

$~{}a_{12}~{}$ | Potentials | $\alpha$ | $\beta$ | $\gamma$
---|---|---|---|---
-85 | POT-I | -0.37 | 14.087 | 1.393
| POT-II | -0.371 | 14.547 | 1.396
-90 | POT-I | -0.578 | 11.778 | 1.359
| POT-II | -0.573 | 12.894 | 1.369
-95 | POT-I | -0.824 | 9.752 | 1.316
| POT-II | -0.774 | 12.22 | 1.351

##  VI Density-functional results

The stability of the self-bound mixture of 41K-87Rb in free space implies the
presence of a surface and a positive surface tension associated to it. We
studied the surface properties of spherical 41K-87Rb droplets within the DFT
approach in the MF+LHY framework, for a fixed ratio of densities
$\rho_{2}/\rho_{1}=\sqrt{g_{11}/g_{22}}$ corresponding to the equilibrium one
for the homogeneous mixture. We provide additional results using a more
accurate density functional obtained from the ab-initio QMC results discussed
in the previous Section.

We have first studied a planar surface, with the aim of determining the
surface tension $\sigma$ in a range of values of the inter-particle attractive
interaction $a_{12}$ accessible to the experiments. When studying a planar
surface, it is useful to use a slab geometry, i.e., we assume an extended
homogeneous system in the $xy$ plane (with periodic boundary conditions), and
with a finite extension along the $z$-direction. In this direction, two
liquid-vacuum interfaces are formed, with the slab width thick enough to have
a constant density region between the two confining surfaces ("bulk" phase).
This amounts to neglect curvature effects, which will be explicitly considered
later on. The surface width of the density profile along $z$ can be quantified
by the parameter $\Delta$, which measures the width between surface points at
$90\%$ and $10\%$ of the bulk total density.

By defining the following coefficients,

$C_{\rm
K}=\frac{1}{4}\left(\frac{\hbar^{2}}{2m_{1}}+\frac{\hbar^{2}}{2m_{2}}\sqrt{\frac{g_{11}}{g_{22}}}\right)\
,$ (18) $C_{\delta}=g_{11}+g_{12}\sqrt{\frac{g_{11}}{g_{22}}}\ ,$ (19)
$C_{\rho}=\frac{8}{15\pi^{2}}\left(\frac{m_{1}}{\hbar^{2}}\right)^{3/2}g_{11}^{5/2}\left[1+\left(\frac{m_{1}}{m_{2}}\right)^{3/5}\sqrt{\frac{g_{22}}{g_{11}}}\right]\
,$ (20)

the effective single-component energy density of the mixture within the MF+LHY
theory, expressed for simplicity in terms of the density $\rho_{1}$, reads

$\mathcal{E}=C_{\rm
K}\frac{(\nabla\rho_{1})^{2}}{\rho_{1}}+C_{\delta}\rho_{1}^{2}+C_{\rho}\rho_{1}^{5/2}.$
(21)

We recall that knowledge of one density is enough to characterize the whole
droplet because of the underlying assumption
$\rho_{2}/\rho_{1}=\sqrt{g_{11}/g_{22}}$ [22, 24].

Remarkably, the surface tension of the planar interface described by local
energy functionals of the form (21) can be estimated, without any prior
knowledge of the density profile, by calculating the following integral [39],

$\sigma=2\int_{0}^{\rho_{0}}d\rho_{1}\,\sqrt{C_{\rm
K}\left(C_{\delta}\rho_{1}+C_{\rho}\rho_{1}^{3/2}-\mu_{0}\right)}\,,$ (22)

where $\mu_{0}=C_{\delta}\rho_{1}+C_{\rho}\rho_{1}^{3/2}$ is the chemical
potential of a liquid system in equilibrium with the vacuum, evaluated at the
equilibrium density $\rho_{0}$. The density profile can also be obtained by
simple quadrature, solving the implicit equation

$z(\rho)=z_{0}+\int_{\rho_{0}/2}^{\rho}\frac{1}{h(\rho^{\prime})}d\rho^{\prime}$
(23)

where $\rho(z_{0})=\rho_{0}/2$ and

$h(\rho)=-\sqrt{\left(\frac{\rho}{C_{\rm
K}}\right)\left[C_{\delta}\rho^{2}+C_{\rho}\rho^{5/2}-\mu_{0}\rho\right]}$
(24)

Here $\rho$ can have any value in the interval $[0,\rho_{0}]$.

The calculated surface tension for different values of $a_{12}$ is shown in
Fig. 3, and compared with the results obtained from the QMC-based functional.
Notice that relatively small changes in the inter-species interaction strength
cause order-of-magnitude changes in the surface tension, which is highlighted
by the logarithmic scale introduced in the $\sigma$-axis. Surface tension
obtained with the QMC-based functionals POT-SR and POT-I are both below the
predictions of the MF+LHY functional, with POT-SR having larger deviations
from MF+LHY.

Figure 3: Surface tension as a function of the (attractive) inter-species
scattering length $a_{12}$. The solid red line is obtained using the MF+LHY
functional, green squares and blue crosses are the predictions assuming POT-SR
and POT-I functionals, respectively (see Table 1 and 4).

The width of the liquid-vacuum interface profile is strongly dependent on the
inter-species scattering length $a_{12}$. In a droplet this will affect the
overall shape of the droplet itself, depending upon the value of the total
number of particles $N$. For larger values of $N$ the droplet will be
characterized by a central region of fairly uniform density (“bulk") and an
external surface region where the density drops to zero with the distance from
its center, whereas it will be an “all surface", gaussian-like droplet, where
the central bulk region is almost absent, for low values of $N$. By defining
the droplet radius $R$ as

$R=\sqrt[3]{\frac{3N}{4\pi\rho_{\text{bulk}}}},\,$ (25)

in the first case the ratio between the surface width $\Delta$ and the droplet
radius $R$ is $\Delta/R\ll 1$, while in the second case $\Delta/R\gg 1$.

We must notice that not all values of $N$ are allowed in a droplet for a given
$a_{12}$ because small droplets, i.e., those with a number of particle below
some critical value $N_{c}$ become unstable when the kinetic energy dominates
over the interaction energy, eventually causing the evaporation of the
droplet. In order to estimate the critical size $N^{c}$ we make a simple
variational ansatz for the radial density profile, which is a good
approximation for spherical, small droplets:

$\rho_{1}(r)=\frac{N_{1}}{\pi^{3/2}\sigma^{3}}e^{-r^{2}/\sigma^{2}}.$ (26)

We use the above ansatz in the energy functional (21) and impose the condition
for a minimum, $\partial E/\partial\sigma=0$, together with the additional
requirement $E=0$ which marks the line separating stable droplets with
negative total energies from unstable ones with positive energies. Solving for
the total number of atoms of species 1, $N_{1}$, we find the critical droplet
size $N_{c}=N_{c1}(1+\sqrt{g_{11}/g_{22}})$, where

$N_{c1}=-\frac{c}{\lambda^{9/2}(a/\lambda^{2}+b/\lambda^{3})}$ (27)

Here $\lambda=-3b/5a$, with $a,b,c$ being given by the expressions $a=6C_{\rm
K}$, $b=C_{\delta}/(2\pi^{3/2})$ and $c=4C_{\rho}/(5\sqrt{10}\pi^{9/4})$. The
table 5 and the figure 4 show the critical atom number for some values of the
scattering length computed with DFT and QMC methods.

Figure 4: Critical atom number $N_{c}$ as a function of the inter-species
scattering length. The solid red line is obtained using the MF+LHY functional,
green squares and blue crosses are the predictions assuming POT-SR and POT-I
functionals, respectively (see Table 1 and 4).

Figure 5 summarizes our results: the color scales shows the ratio between the
surface width $\Delta$ and the radius $R$ of the droplet, which evaporates for
a number of particles below the critical value (solid black line). The red
region identifies gaussian-like, "all-surface" droplets, while the green and
blue region identifies droplets with a well defined, central "bulk" density.
Notice that the black line belongs to the red region for each $|a_{12}|$ so
the assumption of a gaussian density profile for the variational study of the
critical atom number is indeed justified.

Figure 5:  The ratio $\Delta/R$ is represented as a function of $|a_{12}|$ and
$N$. The red region identifies gaussian-like radial density profiles while the
green and blue region identifies droplets with a central bulk region. The
black solid line marks the total critical number of particles $N_{c}$ below
which the droplet evaporates.

The surface tension for the 41K - 87Rb self-bound mixture has been obtained
for a planar interface (we will call it $\sigma_{0}$ to distinguish it from
the size-dependent surface tension of a droplet), though in real droplets an
interfacial curvature of the surface in contact with the vacuum is present.
The curvature-dependent surface tension can be expressed in term of the so-
called Tolman length $\delta$ [40]. To a first approximation the Tolman length
$\delta$ is independent of the droplet size, and it gives the size-dependent
surface tension in term of the one for a planar surface [41, 42]

$\sigma(R)=\sigma_{0}\left(1-\frac{2\delta}{R}\right).$ (28)

A thermodynamic argument relates the Tolman length to the isothermal
compressibility $\kappa^{-1}$ and surface tension $\sigma_{0}$ [43]

$\delta\approx-\kappa^{-1}\sigma_{0}.$ (29)

The isothermal compressibility $\kappa^{-1}$ of a self-bound quantum mixture
with equilibrium bulk densities $\rho_{1}$ and $\rho_{2}$ in the MF+LHY
approach is given by [22]

$\kappa^{-1}=\left(g_{11}\rho_{1}^{2}+g_{22}\rho_{2}^{2}+2g_{12}\rho_{1}\rho_{2}+\frac{15}{4}\mathcal{E}_{\text{LHY}}\right)^{-1}.$
(30)

Previously, some studies have noticed that the product $\kappa^{-1}\sigma_{0}$
is a fundamental characteristic length in liquid droplets [44, 45, 46, 47, 48]
though the connection with the Tolman length $\delta$ was never explicitly
made. We checked the validity of Eq. (29) by independently computing $\delta$
using the Liquid Drop Model (LDM), i.e., writing the calculated total energy
of a droplet made of $N$ atoms as

$E=aN+bN^{2/3}+cN^{1/3}$ (31)

where the separate bulk, surface and curvature contributions to the total
energy of the droplet are highlighted. By using the relation (28) one can see
that $c=8\pi(3/4\pi\rho_{0})^{1/3}\sigma_{0}\delta$. The coefficient $c$ is in
turn obtained by fitting the calculated energies using the LDM expression
quoted above, allowing to determine the Tolman length $\delta$. Predictions
for the Tolman length are summarized in Table 5.

In Fig. 6, the compressibility computed with the MF+LHY approach is shown for
a set of values of the scattering length and compared to QMC results. In
contrast to 39K droplets [16], the compressibility of QMC-based functional is
higher compared to the MF+LHY, due to relatively smaller values of effective
ranges in a K-Rb mixture.

Table 5: Summary of all the quantities reported in the paper. Mass $m$ is the mass of 41Rb atom, and $a_{0}$ is Bohr radius. $\epsilon_{r}$(POT-SR) and $\epsilon_{r}$(POT-I) are the relative differences $(O^{\rm POT-SR}-O^{\rm MF+LHY})/|O^{\rm MF+LHY}|$ and $(O^{\rm POT-I}-O^{\rm MF+LHY})/|O^{\rm MF+LHY}|$ for observable $O$, given in percentages. Positive (negative) values of $\epsilon_{r}$ mean that QMC functionals predict higher (lower) values of the observable with respect to MF+LHY one. Observable | $a_{12}/a_{0}$ | MF+LHY | POT-SR | POT-I |  $\epsilon_{r}$(POT-SR) |  $\epsilon_{r}$(POT-I)
---|---|---|---|---|---|---
| $-85$ | $2.07\cdot 10^{4}$ | $2.31\cdot 10^{4}$ | $2.26\cdot 10^{4}$ | $12$ |  $9$
$N_{c}$ | $-90$ | $7.28\cdot 10^{3}$ | $8.47\cdot 10^{3}$ | $8.31\cdot 10^{3}$ | $16$ | $14$
| $-95$ | $3.49\cdot 10^{3}$ | $4.44\cdot 10^{3}$ | $4.10\cdot 10^{3}$ | $27$ | $17$
| $-85$ | $1.97\cdot 10^{-20}$ | $1.40\cdot 10^{-20}$ | $1.77\cdot 10^{-20}$ | $-29$ |  $-10$
$\dfrac{\sigma}{\hbar^{2}/(ma_{0}^{4})}$ | $-90$ | $8.55\cdot 10^{-20}$ | $5.45\cdot 10^{-20}$ | $7.42\cdot 10^{-20}$ | $-36$ | $-13$
| $-95$ | $2.40\cdot 10^{-19}$ | $1.34\cdot 10^{-19}$ | $1.95\cdot 10^{-19}$ | $-44$ | $-19$
| $-85$ | $4.17\cdot 10^{4}$ | $4.78\cdot 10^{4}$ | $4.48\cdot 10^{4}$ | $14$ |  $7$
$\Delta/a_{0}$ | $-90$ | $2.23\cdot 10^{4}$ | $2.65\cdot 10^{4}$ | $2.46\cdot 10^{4}$ | $19$ | $11$
| $-95$ | $1.43\cdot 10^{4}$ | $1.80\cdot 10^{4}$ | $1.65\cdot 10^{4}$ | $26$ | $15$
| $-85$ | $-5.53\cdot 10^{3}$ | $-6.57\cdot 10^{3}$ | $-6.16\cdot 10^{3}$ | $-19$ |  $-11$
${\delta}/{a_{0}}$ | $-90$ | $-2.95\cdot 10^{3}$ | $-3.67\cdot 10^{3}$ | $-3.43\cdot 10^{3}$ | $-24$ | $-16$
| $-95$ | $-1.90\cdot 10^{3}$ | $-2.50\cdot 10^{3}$ | $-2.35\cdot 10^{3}$ | $-32$ | $-24$
| $-85$ | $2.81\cdot 10^{23}$ | $4.69\cdot 10^{23}$ | $3.47\cdot 10^{23}$ | $67$ |  $24$
$\dfrac{\kappa^{-1}}{ma_{0}^{5}/\hbar^{2}}$ | $-90$ | $3.45\cdot 10^{22}$ | $6.73\cdot 10^{22}$ | $4.62\cdot 10^{22}$ | $95$ | $34$
| $-95$ | $7.90\cdot 10^{21}$ | $1.86\cdot 10^{22}$ | $1.20\cdot 10^{22}$ | $135$ | $52$
Figure 6: Compressibility $\kappa^{-1}$ as a function of the inter-species
scattering length $a_{12}$ in the self-bound droplet regime. The solid red
line is obtained using the MF+LHY functional, green squares and blue crosses
are the predictions assuming POT-SR and POT-I functionals, respectively (see
Table 1 and 4). Figure 7: Radial density profile for a droplet with
$a_{12}=-90\,a_{0}$: the prediction of MF+LHY (red solid line) are compared
with the results obtained from a density functional reproducing the QMC
results, either at zero-range (POT-SR, green dashed line) and with finite
range effects (POT-I, blue dashed line).  Figure 8: Corrected and normalized
surface tension $\sigma/\sigma_{0}$ as a function of the total number of
particles $N$ in the droplet for a set of values of the inter-species
scattering length $a_{12}$, computed within the MF+LHY framework.

Fig. 7 shows the calculated radial density profile for a droplet with
$a_{12}=-90\,a_{0}$, where a comparison is made between the prediction of the
MF+LHY approach and those from the QMC-based energy functionals. The QMC-based
functional, which relies on short-range model potentials (POT-SR, see Fig. 1),
predicts less dense droplets due to repulsive beyond-LHY energy contributions.
Profiles obtained with the POT-I and MF+LHY functionals are more similar since
the effect of increasing the effective range, included in POT-I through model
potentials, is to increase binding energy and the peak density.

In Fig. 8, we show the dependence of the corrected surface tension $\sigma$ as
a function of the droplet’s total number of particles $N$, computed within the
MF+LHY framework by combining Eq. (28) with Eq. (25).

In Table 5, we summarize the predictions for all observables analyzed in our
work with the MF+LHY and QMC-based functionals for $a_{12}=-85a_{0}$,
$-90a_{0}$ and $-95a_{0}$. It can be seen that QMC functionals show increasing
deviation from the predictions of the MF+LHY functional as $|a_{12}|$
increases. QMC-based density functionals constructed assuming short-range
model potentials (POT-SR) predict larger deviations from the MF+LHY theory.
When the correct effective range is included in the model potentials (POT-I),
the predictions become more similar to those of MF+LHY. Variable most
sensitive to beyond-LHY energy corrections appears to be the compressibility
$\kappa^{-1}$, with a relative difference in the range $\approx 20\%$ to
$130\%$. This could have an impact on the collective excitation modes of a
droplet [49], which will be investigated in a future study.

##  VII Conclusion

We have performed QMC calculations for the ground state of a K-Rb liquid
mixture. We found that the zero-range model potentials used in a QMC
calculation predict significant beyond LHY energy contributions. Otherwise,
when the $s$-wave effective range is included in the model potentials, the
energies fall close to the MF+LHY energies. Using both the MF+LHY and the QMC-
based functionals, we have investigated fundamental surface properties for an
experimentally relevant range of scattering parameters. These properties are
relevant to the ongoing experiments because the observed droplets have a large
surface-to-volume ratio. Upon entering a more correlated (denser) regime, the
differences between the predictions for all quantities with the MF+LHY and the
QMC-based functional grow.

###### Acknowledgements.

This work has been supported by the Ministerio de Economia, Industria y
Competitividad (MINECO, Spain) under grants Nos. FIS2017-84114-C2-1-P and
FIS2017-87801-P (AEI/FEDER, UE), and by the EC Research Innovation Action
under the H2020 Programme, Project HPC-EUROPA3 (INFRAIA-2016-1-730897). V. C.
gratefully acknowledges the computer resources and technical support provided
by Barcelona Supercomputing Center. We also acknowledge financial support from
Secretaria d’Universitats i Recerca del Departament d’Empresa i Coneixement de
la Generalitat de Catalunya, co-funded by the European Union Regional
Development Fund within the ERDF Operational Program of Catalunya (project
QuantumCat, ref. 001-P-001644).

## References

  * Petrov [2015] D. S. Petrov, Quantum mechanical stabilization of a collapsing bose-bose mixture, Phys. Rev. Lett. 115, 155302 (2015).
  * Cabrera  _et al._ [2018] C. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell, Quantum liquid droplets in a mixture of bose-einstein condensates, Science 359, 301 (2018).
  * Semeghini  _et al._ [2018] G. Semeghini, G. Ferioli, L. Masi, C. Mazzinghi, L. Wolswijk, F. Minardi, M. Modugno, G. Modugno, M. Inguscio, and M. Fattori, Self-bound quantum droplets of atomic mixtures in free space, Phys. Rev. Lett. 120, 235301 (2018).
  * D’Errico  _et al._ [2019] C. D’Errico, A. Burchianti, M. Prevedelli, L. Salasnich, F. Ancilotto, M. Modugno, F. Minardi, and C. Fort, Observation of quantum droplets in a heteronuclear bosonic mixture, Phys. Rev. Research 1, 033155 (2019).
  * Barranco  _et al._ [2006] M. Barranco, R. Guardiola, S. Hernández, R. Mayol, J. Navarro, and M. Pi, Helium nanodroplets: An overview, Journal of low temperature physics 142, 1 (2006).
  * Böttcher  _et al._ [2019] F. Böttcher, M. Wenzel, J.-N. Schmidt, M. Guo, T. Langen, I. Ferrier-Barbut, T. Pfau, R. Bombín, J. Sánchez-Baena, J. Boronat, and F. Mazzanti, Dilute dipolar quantum droplets beyond the extended gross-pitaevskii equation, Phys. Rev. Research 1, 033088 (2019).
  * Chomaz  _et al._ [2016] L. Chomaz, S. Baier, D. Petter, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino, Quantum-fluctuation-driven crossover from a dilute bose-einstein condensate to a macrodroplet in a dipolar quantum fluid, Phys. Rev. X 6, 041039 (2016).
  * Bisset  _et al._ [2021] R. N. Bisset, L. A. P. Ardila, and L. Santos, Quantum droplets of dipolar mixtures, Phys. Rev. Lett. 126, 025301 (2021).
  * Pethick and Smith [2008] C. J. Pethick and H. Smith,  _Bose–Einstein condensation in dilute gases_ (Cambridge university press, 2008).
  * Lee  _et al._ [1957] T. D. Lee, K. Huang, and C. N. Yang, Eigenvalues and eigenfunctions of a bose system of hard spheres and its low-temperature properties, Phys. Rev. 106, 1135 (1957).
  * Larsen [1963] D. M. Larsen, Binary mixtures of dilute bose gases with repulsive interactions at low temperature, Annals of Physics 24, 89 (1963).
  * Ferioli  _et al._ [2019] G. Ferioli, S. G., L. Masi, G. Giusti, G. Modugno, M. Inguscio, A. Gallemí, A. Recati, and M. Fattori, Phys. Rev. Lett. 122, 090401 (2019).
  * Cikojević  _et al._ [2021] V. Cikojević, L. V. Markić, M. Pi, M. Barranco, F. Ancilotto, and J. Boronat, Dynamics of equilibration and collisions in ultradilute quantum droplets (2021), arXiv:2104.09102 [cond-mat.quant-gas] .
  * Escartín  _et al._ [2019] J. M. Escartín, F. Ancilotto, M. Barranco, and M. Pi, Vorticity and quantum turbulence in the merging of superfluid helium nanodroplets, Physical Review B 99, 10.1103/physrevb.99.140505 (2019).
  * Cikojević  _et al._ [2020] V. Cikojević, L. V. Markić, and J. Boronat, Finite-range effects in ultradilute quantum drops, New Journal of Physics 22, 053045 (2020).
  * Cikojević  _et al._ [2020] V. Cikojević, L. V. Markić, M. Pi, M. Barranco, and J. Boronat, Towards a quantum monte carlo–based density functional including finite-range effects: Excitation modes of a ${}^{39}\mathrm{K}$ quantum droplet, Phys. Rev. A 102, 033335 (2020).
  * Parisi and Giorgini [2020] L. Parisi and S. Giorgini, Quantum droplets in one-dimensional bose mixtures: A quantum monte carlo study, Phys. Rev. A 102, 023318 (2020).
  * Petrov and Astrakharchik [2016] D. S. Petrov and G. E. Astrakharchik, Ultradilute low-dimensional liquids, Phys. Rev. Lett. 117, 100401 (2016).
  * Ota and Astrakharchik [2020] M. Ota and G. E. Astrakharchik, Beyond Lee-Huang-Yang description of self-bound Bose mixtures, SciPost Phys. 9, 20 (2020).
  * Hu and Liu [2020a] H. Hu and X.-J. Liu, Consistent theory of self-bound quantum droplets with bosonic pairing, Phys. Rev. Lett. 125, 195302 (2020a).
  * Cikojević  _et al._ [2018] V. Cikojević, K. Dželalija, P. Stipanović, L. Vranješ Markić, and J. Boronat, Ultradilute quantum liquid drops, Phys. Rev. B 97, 140502 (2018).
  * Ancilotto  _et al._ [2018] F. Ancilotto, M. Barranco, M. Guilleumas, and M. Pi, Self-bound ultradilute bose mixtures within local density approximation, Phys. Rev. A 98, 053623 (2018).
  * Minardi  _et al._ [2019] F. Minardi, F. Ancilotto, A. Burchianti, C. D’Errico, C. Fort, and M. Modugno, Effective expression of the lee-huang-yang energy functional for heteronuclear mixtures, Physical Review A 100, 063636 (2019).
  * Staudinger  _et al._ [2018] C. Staudinger, F. Mazzanti, and R. E. Zillich, Self-bound bose mixtures, Phys. Rev. A 98, 023633 (2018).
  * D’Errico  _et al._ [2007] C. D’Errico, M. Zaccanti, M. Fattori, G. Roati, M. Inguscio, G. Modugno, and A. Simoni, Feshbach resonances in ultracold 39k, New Journal of physics 9, 223 (2007).
  * Marte  _et al._ [2002] A. Marte, T. Volz, J. Schuster, S. Dürr, G. Rempe, E. G. M. van Kempen, and B. J. Verhaar, Feshbach resonances in rubidium 87: Precision measurement and analysis, Phys. Rev. Lett. 89, 283202 (2002).
  * Cikojević  _et al._ [2019] V. Cikojević, L. V. Markić, G. E. Astrakharchik, and J. Boronat, Universality in ultradilute liquid bose-bose mixtures, Phys. Rev. A 99, 023618 (2019).
  * Sarsa  _et al._ [2002] A. Sarsa, J. Boronat, and J. Casulleras, Quadratic diffusion monte carlo and pure estimators for atoms, The Journal of chemical physics 116, 5956 (2002).
  * Boronat [2002] J. Boronat, Microscopic approaches to quantum liquids in confined geometries, World Scientific , 21 (2002).
  * Jastrow [1955] R. Jastrow, Many-body problem with strong forces, Physical Review 98, 1479 (1955).
  * Reatto and Chester [1967] L. Reatto and G. Chester, Phonons and the properties of a bose system, Physical Review 155, 88 (1967).
  * Stoof  _et al._ [2009] H. T. Stoof, K. B. Gubbels, and D. Dickerscheid,  _Ultracold quantum fields_ (Springer, 2009).
  * Chin  _et al._ [2010] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Feshbach resonances in ultracold gases, Rev. Mod. Phys. 82, 1225 (2010).
  * Gribakin and Flambaum [1993] G. F. Gribakin and V. V. Flambaum, Calculation of the scattering length in atomic collisions using the semiclassical approximation, Phys. Rev. A 48, 546 (1993).
  * Flambaum  _et al._ [1999] V. Flambaum, G. Gribakin, and C. Harabati, Analytical calculation of cold-atom scattering, Phys. Rev. A 59, 1998 (1999).
  * Jensen  _et al._ [2006] L. M. Jensen, H. Nilsen, and G. Watanabe, Bcs-bec crossover in atomic fermi gases with a narrow resonance, Physical Review A 74, 043608 (2006).
  * Bohn  _et al._ [1999] J. L. Bohn, J. P. Burke, C. H. Greene, H. Wang, P. L. Gould, and W. C. Stwalley, Collisional properties of ultracold potassium: Consequences for degenerate bose and fermi gases, Phys. Rev. A 59, 3660 (1999).
  * Thalhammer  _et al._ [2009] G. Thalhammer, G. Barontini, J. Catani, F. Rabatti, C. Weber, A. Simoni, F. Minardi, and M. Inguscio, Collisional and molecular spectroscopy in an ultracold bose-bose mixture, New Journal of Physics 11, 055044 (2009).
  * Stringari and Treiner [1987] S. Stringari and J. Treiner, Surface properties of liquid 3he and 4he: A density-functional approach, Phys. Rev. B 36, 8369 (1987).
  * Tolman [1949] R. C. Tolman, The effect of droplet size on surface tension, The journal of chemical physics 17, 333 (1949).
  * Blokhuis and Kuipers [2006] E. Blokhuis and J. Kuipers, Thermodynamic expressions for the tolman length, The Journal of chemical physics 124, 74701 (2006).
  * Rehner and Gross [2018] P. Rehner and J. Gross, Surface tension of droplets and tolman lengths of real substances and mixtures from density functional theory, The Journal of Chemical Physics 148, 164703 (2018), https://doi.org/10.1063/1.5020421 .
  * Bartell [2001] L. S. Bartell, Tolman’s $\delta$, surface curvature, compressibility effects, and the free energy of drops, The Journal of Physical Chemistry B 105, 11615 (2001).
  * Mayer [1963] S. Mayer, A molecular parameter relationship between surface tension and liquid compressibility, The Journal of Physical Chemistry 67, 2160 (1963).
  * Helfrich [1973] W. Helfrich, Elastic properties of lipid bilayers: theory and possible experiments, Zeitschrift für Naturforschung C 28, 693 (1973).
  * Fisk and Widom [1969] S. Fisk and B. Widom, Structure and free energy of the interface between fluid phases in equilibrium near the critical point, The Journal of Chemical Physics 50, 3219 (1969).
  * Egelstaff and Widom [1970] P. Egelstaff and B. Widom, Liquid surface tension near the triple point, The Journal of Chemical Physics 53, 2667 (1970).
  * Mon and Stroud [1981] K. Mon and D. Stroud, Surface widths of simple liquids and an empirical law of freezing, The Journal of Chemical Physics 74, 2078 (1981).
  * Hu and Liu [2020b] H. Hu and X.-J. Liu, Collective excitations of a spherical ultradilute quantum droplet, Phys. Rev. A 102, 053303 (2020b).