Weber number and the outcome of binary collisions between quantum droplets

In most circumstances, degenerate states of dilute atomic gases are realized using samples trapped by different means. The question of whether they can also exist in isolation, without a trap potential, was theoretically answered in the affirmative within two scenarios. The first one corresponds to a homogeneous atomic gas where the two-body scattering length is negative and much larger than the effective two-body interaction radius; then the contribution to the ground state energy due to the three–body correlations could become dominant and provide a self–trapping mechanism [1]. The second scenario corresponds to a condensed Bose-Bose mixture where the interspecies attraction becomes stronger than the geometrical average of the intra–species repulsions [2]. The stabilization mechanism is then provided by quantum fluctuations that are a direct manifestation of beyond mean-field effects. The implementation of the latter mechanism was the basis for the experimental observation of self-bound droplets of ultra cold atoms in homo–nuclear [3, 4] and hetero–nuclear[5] mixtures of alkali atoms, and in dipolar condensates [7, 8, 9, 10]. In the latter case, it was required to drive a Bose-Einstein condensate (BEC) into the strongly interacting regime.

The first correction to the mean-field contribution to the ground state energy of a homogeneous weakly repulsive Bose gas is given by the Lee-Huang-Yang (LHY) term [11]. This term is negligible in many circumstances. However, there are situations where the LHY and the mean field contributions can be of the same order without leaving the weakly interacting regime. For a Bose-Bose mixture where the inter– and intra–species coupling constants are independently controlled, the mean–field term, which is proportional to the square of the density $n^{2}$, could be negative and the LHY contribution, which is proportional to $n^{5/2}$, could be positive. Then, the quantum LHY repulsion could neutralize the mean-field attraction and stabilize the system against collapse [2]. For a large enough number of atoms, the resulting self–trapping quantum state is characterized by a core with a quasi uniform density and a thin free surface similar to that of a liquid droplet. A less intuitive, though more precise description of the self–trapping mechanism can be obtained via Quantum Monte Carlo simulations [12], they predict similar general properties for the quantum droplets, though different specific scenarios of their occurrence. The general properties are: generic features of the collective modes [9], the presence of striped states in confined geometries [13], and the existence of arrays of phase coherent droplets with transient super solid properties [14, 10, 15].

For small atom numbers, the quantum droplet has no collective modes with energy lower than the particle emission threshold so that self-evaporation can occur [2, 16]. In this work we consider droplets with small and large atom numbers. For the latter, self-evaporation is not expected and the resulting quantum state is closer to the liquid droplet expectations, i.e., it exhibits properties nearby a classical liquid, a helium nano droplet, or an atomic nuclei within the liquid drop model. Note, however, that in this regime three–body recombination collisions may reduce considerably the lifetime of the quantum droplets created by a binary mixture of Bose-Bose ultra cold atoms. Thus, we take into consideration recent experimental results [5] which demonstrate that, for an adequate choosing of hetero-nuclear mixtures, quantum droplets with lifetimes in the order of tens of milliseconds can be produced.

The collisions between droplets is an interesting mechanism both to study the quantum droplet excitations, and to realize quantum simulations of analog many body systems. Up to now, this mechanism has been used to observe the crossover between compressible and incompressible regimes of quantum droplets [17]. In the incompressible regime surface waves occur and surface tension should be a relevant parameter to understand the conditions required for different results of quantum droplet collisions.

Our analisis is structured as follows. Following the ideas implemented by the nuclear physics community to study the expectations of the nuclear liquid drop model [20, 21, 22], expressions for the low energy surface excitations are used as a variational ansatz. These expressions include effects directly due to the thickness of the surface of the droplet, and allow to quantify the surface tension of the ground state droplet. We evaluate numerically – within the extended Gross-Pitaevskii equation (EGPE) that incorporates the LHY terms in three dimensions– the ground state for a finite number of atoms, ranging from $10^{4}-10^{7}$ in Bose-Bose mixtures. The analysis is performed for binary mixtures of homo- and hetero–nuclear atomic gases within the conditions on the interaction strengths that give rise to quantum droplets. In the case of homo–nuclear BEC mixtures, we also obtain the ground state within an effective formalism that considers finite-range effects as developed in Refs. [12, 18]. General features of the resulting ground states are then evaluated. They include the surface tension derived from the different surface excitation models described in the previous Section, and the dependence of the droplets radius and saturation energy with the number of atoms in the droplets. Whenever it is possible, the results are compared to previous theoretical predictions. In the next Section, atom losses due to self-evaporation and three–body recombination effects are analyzed. Finally, the outcomes of frontal binary collisions are explored. For classical liquid droplets the collisions are usually investigated in terms of the impact parameter and of the Weber number. The latter is a measure of the relative importance of the inertia of a fluid in terms of the kinetic energy compared to its surface tension [23]. In this work, several regimes for binary front collisions that range from the coalescence of the quantum droplets to their disintegration into smaller droplets are identified in terms of the Weber number.

We are interested in obtaining expressions for the surface tension of a quantum droplet. To that end we shall follow the formalism developed by Berstch in the study of capillary waves in a superfluid [20], and which was later extended to Fermi systems with spherical symmetry [21] having in mind the liquid drop model of nuclei. This formalism is based on the random phase approximation (RPA) and the Thouless variational theorem. In our case, the random phase approximation is equivalent to the linearization of the time-dependent equations that define the dynamics of the system in the vicinity of a variational minimum; Thouless variational theorem guarantees that actual excitation energies are a lower bound to the excitation energies that result from using, as an ansatz, known expressions of excitations with harmonic time dependence and satisfying adequate boundary conditions.

Let us consider a binary mixture of Bose gases composed by $N^{(a)}$ and $N^{(b)}$ atoms of each species. Let us assume that the state of the system is approximately described by the Hartree many-body wavefunction

$$\Psi_{N}(\vec{r}_{1},...,\vec{r}_{N^{(a)}};\vec{r}^{\prime}_{1},...,\vec{r}^{\prime}_{N^{(b)}})=\Big{[}\prod_{i=1}^{N^{(a)}}\chi^{(a)}(\vec{r}_{i})\Big{]}\Big{[}\prod_{i=1}^{N^{(b)}}\chi^{(b)}(\vec{r}^{\prime}_{i})\Big{]}.$$

(1)

As order parameters we consider

$$\psi^{(\alpha)}(\vec{r},t)=\sqrt{N^{(\alpha)}}\chi^{(\alpha)}(\vec{r},t),\quad\quad\alpha=a,b,$$

(2)

$\psi^{(\alpha)}$ evolves according to the equation

	$\displaystyle i\hbar\partial_{t}\psi^{(\alpha)}(\vec{r},t)$	$\displaystyle=$	$\displaystyle\hat{H}\psi^{(\alpha)}(\vec{r},t)$		(3)
		$\displaystyle=$	$\displaystyle\hat{H}_{0}^{\alpha}\psi^{(\alpha)}+\Big{(}\int dr^{3}dr^{\prime 3}U_{\alpha\alpha}(\vec{r},\vec{r}^{\prime};\rho^{(\alpha)}(\vec{r}^{\prime},t))+U_{\alpha\beta}(\vec{r},\vec{r}^{\prime};\rho^{(\beta)}(\vec{r}^{\prime},t))\Big{)}\psi^{(\alpha)}(\vec{r},t),\quad\alpha\neq\beta$		(4)
	$\displaystyle\hat{H}_{0}^{\alpha}$	$\displaystyle=$	$\displaystyle-\frac{\hbar^{2}}{2m_{\alpha}}\nabla^{2}+V^{\alpha}_{ext}(\vec{r})$		(5)

In Eq. (4), $V^{\alpha}_{ext}(\vec{r})$ represents an external potential and $\rho^{(\alpha)}=\psi^{(\alpha)*}\psi^{(\alpha)}$ the density of $\alpha$–atoms. The potentials $U_{\alpha\beta}$ may incorporate, in an effective way, interactions beyond the standard mean field approximation. In the case of the extended Gross-Pitaeviskii formalism, the density dependent interactions are superpositions of contact terms, and have the structure,

$$U_{\alpha\beta}(\vec{r},\vec{r}^{\prime};\rho^{(\gamma)}(\vec{r}^{\prime},t))=\delta(\vec{r}-\vec{r}^{\prime})\sum_{i}g^{i}_{\alpha\beta}(\rho^{(\gamma)}(\vec{r},t))^{s_{i}},\quad\quad\alpha,\beta,\gamma=a,b,$$

(6)

with $s_{i}$ a real number. The stationary solutions of Eq. (4) define the chemical potentials $\mu_{\alpha}$,

$$\psi^{(\alpha)}_{0}(\vec{r},t)=e^{-i\mu_{\alpha}t/\hbar}\phi^{(\alpha)}_{0}(\vec{r}),\quad\quad\alpha=a,b.$$

(7)

Let us consider collective Bogolubov excitations with a harmonic time dependence,

$$\psi^{(\alpha)}_{exc}(\vec{r},t)=e^{-i\mu_{\alpha}t/\hbar}\Big{(}\phi^{(\alpha)}_{0}(\vec{r})+\sqrt{N^{(\alpha)}}\sum_{q}(u_{q}^{(\alpha)}(\vec{r})e^{-i\omega_{q}t}+v_{q}^{(\alpha)*}(\vec{r})e^{-i\omega_{q}t})\Big{)}.$$

(8)

In this equation $q$ is a label that determine the characteristics of the excitation derived, for example, from its geometry, excitation energy etc. The excitation functions belong to the space expanded by a normalized and complete basis set $\{u_{q}^{a},v_{q}^{a},u_{q}^{b},u_{q}^{b}\}$, according to the scalar product

$$\int d^{3}r[u_{q}^{(\alpha)}(\vec{r})u_{p}^{(\alpha)*}(\vec{r})-v_{q}^{(\alpha)}(\vec{r})v_{p}^{(\alpha)*}(\vec{r})]=\frac{\omega_{q}}{|\omega_{q}|}\delta_{qp},\quad\omega_{q}\neq 0.$$

(9)

Excitations with $\omega_{q}=0$ are assumed to be orthogonal to $\phi^{(\alpha)}_{0}$. Within the linear response approximation around the stationary function $\psi^{(\alpha)}_{0}(\vec{r},t)$ , Eq. (4) is equivalent to linear equations with the structure

$${\mathcal{M}}\begin{pmatrix}u_{q}^{(a)}\\ v_{q}^{(a)}\\ u_{q}^{(b)}\\ v_{q}^{(b)}\end{pmatrix}=\hbar\omega_{q}\begin{pmatrix}u_{q}^{(a)}\\ -v_{q}^{(a)}\\ u_{q}^{(b)}\\ -v_{q}^{(b)}\end{pmatrix}.$$

(10)

Thouless variational theorem establishes that the energy $\hbar\omega_{exact}$ of harmonic exact solutions of Eq. (4) are a lower bound of $\hbar\omega_{q}$. This condition is written in our case as

$$\hbar\omega_{exact}\leq\frac{1}{2}\frac{\sum_{\alpha,\beta=a,b}\langle\eta_{q}^{(\alpha)}|\hat{\Xi}^{\alpha\beta}|\eta_{q}^{(\beta)}\rangle+\sum_{\alpha=a,b}\langle\zeta_{q}^{(\alpha)}|\hat{\Delta}^{\alpha}|\zeta_{q}^{(\alpha)}\rangle}{\sum_{\alpha=a,b}(\langle\zeta_{q}^{(\alpha)}|\eta_{q}^{(\alpha)}\rangle)}=\hbar\omega_{\eta\zeta}$$

(11)

where

	$\displaystyle\eta_{q}^{(\alpha)}$	$\displaystyle=$	$\displaystyle\sqrt{N^{(\alpha)}}(u_{q}^{(\alpha)}+v_{q}^{(\alpha)*}),\quad\zeta_{q}^{(\alpha)}=\sqrt{N^{(\alpha)}}(u_{q}^{(\alpha)}-v_{q}^{(\alpha)*}),$		(12)
	$\displaystyle\hat{\Delta}^{\alpha}$	$\displaystyle=$	$\displaystyle\hat{H}_{0}^{\alpha}+\int dr^{3}\int dr^{\prime 3}\big{(}U_{\alpha\alpha}(\vec{r},\vec{r}^{\prime};\rho^{(\alpha)}_{0}(\vec{r}^{\prime}))+U_{\alpha\beta}(\vec{r},\vec{r}^{\prime};\rho^{(\beta)}_{0}(\vec{r}^{\prime})\big{)}-\mu_{\alpha},$		(13)
	$\displaystyle\hat{\Xi}^{\alpha\beta}\eta_{q}^{(\beta)}(\vec{r})$	$\displaystyle=$	$\displaystyle\delta_{\alpha\beta}\hat{\Delta}^{\alpha}\eta_{q}^{(\beta)}(\vec{r})+2\phi^{(\alpha)}_{0}(\vec{r})\int d^{3}r^{\prime}U_{\alpha\beta}(\vec{r},\vec{r}^{\prime};\rho^{(\beta)}_{0})\phi^{(\beta)}_{0}(\vec{r}^{\prime})\eta_{q}^{(\beta)}(\vec{r}^{\prime}).$		(14)

Our ansatz for the excitation functions depend functionally on $\phi_{0}$ and its derivatives

$$\eta_{q}^{(\alpha)}=\eta_{q}^{(\alpha)}[\phi_{0}^{(\alpha)},\partial\phi_{0}^{(\alpha)}],\quad\zeta_{q}^{(\alpha)}=\gamma_{\alpha}\tilde{\zeta}_{q}^{(\alpha)}[\phi_{0}^{(\alpha)},\partial\phi_{0}^{(\alpha)}].$$

(15)

The $\gamma^{\alpha}$ parameters are determined by requiring to yield a minimum for $\omega_{\eta\zeta}$, the result is

	$\displaystyle\gamma_{a}$	$\displaystyle=$	$\displaystyle C_{a}\sqrt{\frac{AB_{b}}{B_{a}(B_{b}C_{a}^{2}+B_{a}C_{b}^{2})}},\quad\gamma_{b}=\frac{B_{a}C_{b}}{B_{b}C_{a}}\gamma_{a}$		(16)
	$\displaystyle A$	$\displaystyle=$	$\displaystyle\sum_{\alpha,\beta=a,b}\langle\eta_{q}^{(\alpha)}|\hat{\Xi}^{\alpha\beta}|\eta_{q}^{(\beta)}\rangle,\quad B_{\alpha}=\langle\tilde{\zeta}^{(\alpha)}_{q}|\hat{\Delta}^{\alpha}|\tilde{\zeta}^{(\alpha)}_{q}\rangle,\quad C_{\alpha}=\langle\tilde{\zeta}_{q}^{\alpha}|\eta_{q}^{\alpha}\rangle.$		(17)

The optimal variational excitation energy satisfies the equation

$$\hbar^{2}\omega_{opt}^{2}=\frac{AB_{a}B_{b}}{B_{a}C_{b}^{2}+B_{b}C_{a}^{2}}.$$

(18)

We shall consider functions $\eta$ and $\tilde{\zeta}$ with the simple structure,

$$\eta^{\alpha}_{q}=\vec{F}_{q}^{\alpha}(\vec{r})\cdot\xi^{(\alpha)}\vec{\nabla}\phi_{0}^{(\alpha)},\quad\tilde{\zeta}^{\alpha}_{q}=\Lambda^{\alpha}_{q}(\vec{r})\phi_{0}^{(\alpha)}.$$

(19)

The functions $\vec{F}^{\alpha}(\vec{r})$ and $\Lambda^{\alpha}(\vec{r})$ are usually chosen to incorporate geometric properties of the droplets surface. The parameter $\xi^{(\alpha)}$ corresponds to a characteristic scale length for the $\alpha$-fluid. A direct calculation shows that, in case the dynamics is governed by an extended Gross Pitaevskii equation characterized by contact interactions, Eq.(6), the contribution of the interaction terms of the excitation functions Eq. (19) to $A$ and $B_{\alpha}$ is null. However, information about such interaction terms is still contained in the stationary solution $\phi_{0}^{(\alpha)}$.

The particular selection

$$\vec{F}_{q}^{\alpha}(\vec{r})=\xi^{(\alpha)}\nabla\Lambda_{q}^{\alpha}(\vec{r}),$$

(20)

allows a hydrodynamic interpretation of $\zeta^{(\alpha)}_{q}/\phi_{0}^{(\alpha)}$ as a velocity field asociated to an infinitesimal change of the density in the direction given by $\vec{\nabla}\zeta^{(\alpha)}_{q}/\phi_{0}^{(\alpha)}$ whenever the interspecies interaction $U_{\alpha\beta}$ is negligible.

If a sphere delimits the boundary of the $\alpha$-fluid there are two alternative configurations. In the case the fluid is contained within the sphere (droplet) or outside it (bubble). Then, singularities at the origin or at infinity are avoided by taking $\Lambda^{\alpha}$ as the adequate solution of Laplace equation in spherical coordinates

	$\displaystyle\Lambda_{\ell m}^{\alpha}(\vec{r})$	$\displaystyle=$	$\displaystyle\mathcal{A}_{\ell m}Y_{\ell m}(\theta,\varphi)\Big{(}\frac{r}{\xi^{(\alpha)}}\Big{)}^{\ell}\quad\quad\quad{\mathrm{droplet}}$		(21)
		$\displaystyle=$	$\displaystyle\mathcal{A}_{\ell m}Y_{\ell m}(\theta,\varphi)\Big{(}\frac{r}{\xi^{(\alpha)}}\Big{)}^{-(\ell+1)}\quad{\mathrm{bubble}},$		(22)

with $\mathcal{A}_{\ell m}$ a normalization constant to be evaluated according to Eq. (9) depending on the selection of the $\vec{F}^{\alpha}_{\ell m}(\vec{r})$ function. Two reasonable options are identified,

Ansatz 1.

$$\vec{F}^{\alpha}_{\ell m}(\vec{r})=\mathcal{A}^{(1)}_{\ell m}Y_{\ell m}(\theta,\varphi)\hat{r}.$$

(23)

Ansatz 2.

$$\vec{F}^{\alpha}_{\ell m}(\vec{r})=\xi^{(\alpha)}\vec{\nabla}\Lambda_{\ell m}^{\alpha}.$$

(24)

We now focus on the particular case where $\phi_{0}^{(\alpha)}$ and the external potential $V_{ext}$ do not depend on $\theta$ and $\varphi$. Then, the corresponding expressions for the $A$, $B_{\alpha}$ and $C_{\alpha}$ functionals that determine the excitation energy, Eq. (18), for a dynamics determined by an extended Gross-Pitaevskii equation are

Ansatz 1.

	$\displaystyle A$	$\displaystyle=$	$\displaystyle|\mathcal{A}^{(1)}_{\ell m}|^{2}\sum_{\alpha=a,b}\xi^{(\alpha)2}\Big{[}-\frac{1}{2}\int dr\partial_{r}\rho^{(\alpha)}\partial_{r}V^{\alpha}_{ext}+\frac{\hbar^{2}}{2m_{\alpha}}(\ell(\ell+1)-2\big{)}\int dr(\partial_{r}\phi_{0}^{(\alpha)})^{2}\Big{]}$		(25)
	$\displaystyle B_{\alpha}$	$\displaystyle=$	$\displaystyle-|\mathcal{A}^{(1)}_{\ell m}|^{2}\frac{\hbar^{2}}{2m_{\alpha}}\ell\int drr\partial_{r}\rho^{(\alpha)}\Big{(}\frac{r}{\xi^{(\alpha)}}\Big{)}^{2\ell}\quad\quad\quad\quad\quad{\mathrm{droplet}}$		(26)
		$\displaystyle=$	$\displaystyle|\mathcal{A}^{(1)}_{\ell m}|^{2}\frac{\hbar^{2}}{2m_{\alpha}}(\ell+1)\int drr\partial_{r}\rho^{(\alpha)}\Big{(}\frac{r}{\xi^{(\alpha)}}\Big{)}^{-(2\ell+2)}\quad{\mathrm{bubble}}$		(27)
	$\displaystyle C_{\alpha}$	$\displaystyle=$	$\displaystyle|\mathcal{A}^{(1)}_{\ell m}|^{2}\frac{\xi^{(\alpha)}}{2}\int drr^{2}\partial_{r}\rho^{(\alpha)}\Big{(}\frac{r}{\xi^{(\alpha)}}\Big{)}^{\ell}\quad\quad\quad\quad\quad\quad{\mathrm{droplet}}$		(28)
		$\displaystyle=$	$\displaystyle|\mathcal{A}^{(1)}_{\ell m}|^{2}\frac{\xi^{(\alpha)}}{2}\int drr^{2}\partial_{r}\rho^{(\alpha)}\Big{(}\frac{r}{\xi^{(\alpha)}}\Big{)}^{-(\ell+1)}\quad\quad\quad\quad{\mathrm{bubble}}$		(29)

Ansatz 2.

	$\displaystyle A$	$\displaystyle=$	$\displaystyle|\mathcal{A}^{(2)}_{\ell m}|^{2}\sum_{\alpha=a,b}\frac{\xi^{(\alpha)4}}{2}\Big{[}-\int drr^{2}\partial_{r}\rho_{0}^{(\alpha)}\partial_{r}V_{ext}(\vec{r})\Big{(}\frac{r}{\xi^{(\alpha)}}\Big{)}^{2\ell}$
		$\displaystyle+$	$\displaystyle\frac{\hbar^{2}}{m_{\alpha}}\ell^{2}(\ell-1)(2\ell+1)\int drr^{-2}(\partial_{r}\phi_{0}^{(\alpha)})^{2}\Big{(}\frac{r}{\xi^{(\alpha)}}\Big{)}^{2\ell}\Big{]}\quad\quad\quad\quad\quad\quad{\mathrm{droplet}}$
		$\displaystyle=$	$\displaystyle|\mathcal{A}^{(2)}_{\ell m}|^{2}\sum_{\alpha=a,b}\frac{\xi^{(\alpha)4}}{2}\Big{[}-\int drr^{2}\partial_{r}\rho_{0}^{(\alpha)}\partial_{r}V_{ext}(\vec{r})\Big{(}\frac{r}{\xi^{(\alpha)}}\Big{)}^{-2(\ell+1)}$
		$\displaystyle+$	$\displaystyle\frac{\hbar^{2}}{m_{\alpha}}(\ell+1)^{2}(\ell+2)(2\ell+1)\int drr^{-2}(\partial_{r}\phi_{0}^{(\alpha)})^{2}\Big{(}\frac{r}{\xi^{(\alpha)}}\Big{)}^{-2(\ell+1)}\Big{]}\quad{\mathrm{bubble}}$		(31)
	$\displaystyle B_{\alpha}$	$\displaystyle=$	$\displaystyle-|\mathcal{A}^{(2)}_{\ell m}|^{2}\frac{\hbar^{2}}{2m_{\alpha}}\ell\int dr\;r\partial_{r}\rho_{0}^{(\alpha)}\Big{(}\frac{r}{\xi^{(\alpha)}}\Big{)}^{2\ell}\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad{\mathrm{droplet}}$		(32)
		$\displaystyle=$	$\displaystyle|\mathcal{A}^{(2)}_{\ell m}|^{2}\frac{\hbar^{2}}{2m_{\alpha}}((\ell+1))\int dr\;r\partial_{r}\rho_{0}^{(\alpha)}\Big{(}\frac{r}{\xi^{(\alpha)}}\Big{)}^{-2(\ell+1)}\quad\quad\quad\quad\quad\quad{\mathrm{bubble}}$		(33)
	$\displaystyle C_{\alpha}$	$\displaystyle=$	$\displaystyle-\frac{m_{\alpha}\xi^{(\alpha)2}}{\hbar^{2}}B_{\alpha}$		(34)

For quantum droplets it can be assumed that the ground state order parameters are proportional to each other,

$$\phi^{(a)}_{0}=\sqrt{\bar{\beta}}\phi^{(b)}_{0}=\sqrt{\frac{N^{(a)}}{N^{(b)}}}\phi^{(b)}_{0}$$

(35)

the last equality being a consequence of Eq. (2). Since, they can be produced without an external potential, and the ground state depends only on the radial coordinate, the excitation energies become

Ansatz 1.

$$\hbar^{2}\omega^{2}\leq-\ell(\ell-1)(\ell+2)\frac{\hbar^{2}}{m_{a}m_{b}}\frac{N^{(a)}m_{b}+N^{(b)}m_{a}}{N^{(a)}m_{a}+N^{(b)}m_{b}}\frac{\int dr(\partial_{r}\phi^{(a)}_{0})^{2}\int dr\partial_{r}\rho^{(a)}_{0}r^{2\ell+1}}{\left(\int dr\partial_{r}\rho^{(a)}_{0}r^{\ell+2}\right)^{2}}$$

(36)

Ansatz 2.

$$\hbar^{2}\omega^{2}\leq-\ell(\ell-1)(2\ell+1)\frac{\hbar^{2}}{m_{a}m_{b}}\frac{N^{(a)}m_{b}+N^{(b)}m_{a}}{N^{(a)}m_{a}+N^{(b)}m_{b}}\frac{\int dr\,(\partial_{r}\phi^{(a)}_{0})^{2}r^{2\ell-2}}{\int dr\,\partial_{r}\rho^{(a)}_{0}r^{2\ell+1}},$$

(37)

whenever an equal length scale

$$\xi^{(a)}=\xi^{(b)}=\xi,$$

(38)

is chosen for both species. In these equations we notice the presence of the total mass of the droplet

$$M_{T}=N^{(a)}m_{a}+N^{(b)}m_{b}.$$

(39)

In his seminal work [25], Lord Rayleigh considered a classical spherical droplet with constant density of particles $\rho_{0}$ with mass $m$ and a radius $R_{0}$. For the frequency of excitation $\omega_{exc}$, he showed that

$$\omega_{exc}^{2}=\ell(\ell-1)(\ell+2)\frac{\sigma}{m\rho_{0}R_{0}^{3}}$$

(40)

and identified $\sigma$ as the surface tension. Lord Rayleigh obtained such an structure for $\omega_{exc}$ analyzing the lowest energy surface excitations of a classical droplet via a variational theorem. The term

$$m\rho_{0}R_{0}^{3}=\frac{3}{4\pi}M_{T}$$

is directly related to the total mass of the droplet $M_{T}$.

We observe that the dependence on the multipole parameter $\ell$ both for ansatz 1 and 2, shows that monopole excitations are not surface excitations, and that the $\ell=1$ contribution is null. This is expected since this term corresponds to a translation of the whole droplet. However for $\ell\geq 2$, the estimation value for $\epsilon=\hbar\omega_{exc}$ of the spherical quantum droplet, and as a consequence of the surface tension, is different for each ansatz:

• (i)

  For ansatz 1, it resembles that of Rayleigh and that of the liquid drop model of nuclei in the limit of $N\rightarrow\infty$ nucleons [26].

• (ii)

  For ansatz 2, it is similar to that predicted for atomic nuclei by Bertsch [21] and by Casas and Strigari [22] using the random phase approximation and the density-density Green’s function formalism; in both cases the second ansatz was taken for the excitation modes.

A variational criteria based on the excitation energy can be applied to select between the two excitation ansatz. Following Lord Rayleigh ideas, the corresponding surface tension estimation would be given by the expressions

Ansatz 1.

$$\sigma^{(1)}_{\ell}=-\frac{\hbar^{2}}{M}\frac{\int dr(\partial_{r}\phi^{(a)}_{0})^{2}\int dr\partial_{r}\rho^{(a)}_{0}r^{2\ell+1}}{\left(\int dr\partial_{r}\rho^{(a)}_{0}r^{\ell+2}\right)^{2}}$$

(41)

Ansatz 2.

$$\sigma^{(2)}_{\ell}=-\frac{\hbar^{2}}{M}\frac{\int dr\,(\partial_{r}\phi^{(a)}_{0})^{2}r^{2\ell-2}}{\int dr\,\partial_{r}\rho^{(a)}_{0}r^{2\ell+1}}$$

(42)

with

$$M=\frac{4\pi}{3}\frac{m_{a}m_{b}}{N^{(a)}m_{b}+N^{(b)}m_{a}}.$$

Consider a homogeneous weakly repulsive Bose gas composed by hard spheres of diameter “$\mathrm{a}$” having a density $n$. The energy correction to the mean field ground state energy introduced by Lee, Huang and Yang [11] is the first term in a power series expansion in the parameter $(n{\mathrm{a}}^{3})^{1/2}$. The analogous term for a mixture of two ultracold Bose gases with atom masses $m_{a}$ and $m_{b}$, densities $n_{a}$ and $n_{b}$, and scattering lengths $\mathrm{a}_{\alpha\beta}$, $\alpha,\beta=a,b$, is [5]:

	$\displaystyle\epsilon_{LHY}=$	$\displaystyle\frac{1}{4\pi^{2}}\left(\frac{m_{a}}{\hbar^{2}}\right)^{3/2}(g_{aa}n_{a})^{5/2}\int^{\infty}_{0}k^{2}\mathcal{F}(k,z,u,x)dk$		(43)
		$\displaystyle\mathcal{F}(k,z,u,x)=$
		$\displaystyle\bigg{(}\frac{1}{2}\left[k^{2}\big{(}1+\frac{x}{z}\big{)}+\frac{k^{4}}{4}\big{(}1+\frac{1}{z^{2}}\big{)}\right]+\bigg{[}\frac{1}{4}\big{[}k^{2}+\frac{k^{4}}{4}-\big{(}\frac{x}{z}k^{2}+\frac{1}{4z^{2}}\big{)}\big{]}^{2}+\frac{ux}{z}k^{4}\bigg{]}^{1/2}\bigg{)}^{1/2}+$
		$\displaystyle\bigg{(}\frac{1}{2}\left[k^{2}\big{(}1+\frac{x}{z}\big{)}+\frac{k^{4}}{4}\big{(}1+\frac{1}{z^{2}}\big{)}\right]-\bigg{[}\frac{1}{4}\big{[}k^{2}+\frac{k^{4}}{4}-\big{(}\frac{x}{z}k^{2}+\frac{1}{4z^{2}}\big{)}\big{]}^{2}+\frac{ux}{z}k^{4}\bigg{]}^{1/2}\bigg{)}^{1/2}$
		$\displaystyle-\frac{1+z}{2z}k^{2}-(1+x)+\frac{1}{k^{2}}\big{[}1+x^{2}z+4ux\frac{z}{1+z}\big{]}.$

In this equation the coupling strength factors $g_{\alpha\beta}=2\pi\hbar^{2}\mathrm{a}_{\alpha\beta}/m_{\alpha\beta}$, with $m_{\alpha\beta}=m_{\alpha}m_{\beta}/(m_{\alpha}+m_{\beta})$, were introduced as well as the notation,

$$\frac{m_{b}}{m_{a}}=z,\quad\frac{g^{2}_{ab}}{g_{aa}g_{bb}}=u,\quad\frac{g_{bb}n_{b}}{g_{aa}n_{a}}=x.$$

Let us define

	$\displaystyle f(z,u,x)$	$\displaystyle=:$	$\displaystyle\frac{15}{32}\int^{\infty}_{0}k^{2}\mathcal{F}(k,z,u,x)dk$		(44)
		$\displaystyle\overset{k=\tan(t)}{=}$	$\displaystyle\frac{15}{32}\int^{\pi/2}_{0}\left(\frac{\sin^{2}(t)}{\cos^{4}(t)}\right)\mathcal{F}(\tan(t),z,u,x)dt.$

In spite of being a combination of terms that are individually divergent, in general, $f(z,u,x)$ converges [24]. Consider the case $g_{aa}>0$, $g_{bb}>0$, and $g_{ab}<0$ and define $\delta g=g_{ab}+\sqrt{g_{aa}g_{bb}}$. Under the weak coupling condition $n\mathrm{a}^{3}<<1$ it results  [2] that: (i)the main contribution to the finite integral comes from terms with $k=1/\xi_{h}=\sqrt{mgn}$; (ii) long-wavelength instabilities expected from a mean field theory approach are cured by quantum fluctuations at shorter wavelengths $k\sim 1/\xi_{h}\gg\sqrt{m|\delta g|n}$; (iii) the equillibrum densities are determined by $g_{\alpha\alpha}$, $n_{b}^{(0)}/n_{a}^{(0)}=\sqrt{g_{aa}/g_{bb}}$. These conditions define the quantum droplet state.

For equal masses[2]

$$f(1,u,x)=\frac{1}{4\sqrt{2}}\Big{[}\big{(}1+x+\sqrt{(1-x)^{2}+4ux}\big{)}^{5/2}+\big{(}1+x-\sqrt{(1-x)^{2}+4ux}\big{)}^{5/2}\Big{]},$$

(45)

which becomes a complex valued expression for $u<1$ and the effective LHY terms yield a non Hermitean EGPE. Within the quantum droplets regime $u\sim 1$, it results that [5]

$$f(z,1,x)\simeq(1+z^{3/5}x)^{5/2}.$$

(46)

This approximate expression for the LHY term in the extended Gross-Pitaevskii equation facilitates the numerical simulations. However, a direct numerical calculation of $f(z,u,x)$ is also feasible. In Fig. 1 the results of a numerical calculation of $f(z,u,x)$ and the result of using the approximate expression Eq. (46) are illustrated for $u\sim 1$ and several $z$ values, notice the logarithmic scale on (a) and (c) axes. The value $u=169/150$ in Figs. 1(a) and (b) approximates that expected for homo–nuclear binary mixtures of ³⁹K under feasible experimental conditions [3, 4] ($\mathrm{a}_{aa}=\mathrm{a}_{bb}=48.57\mathrm{a}_{0}$, $\mathrm{a}_{ab}=-51.86\mathrm{a}_{0}$), or hetero–nuclear binary mixtures of ⁴¹K and ⁸⁷Rb, also under feasible experimental conditions [6] ($\mathrm{a}_{KK}=62.0\mathrm{a}_{0}$, $\mathrm{a}_{RbRb}=100.4\mathrm{a}_{0}$ and $\mathrm{a}_{KRb}=-82.0\mathrm{a}_{0}$). Notice the small relative value of the imaginary part of $f$ with respect to the real one; this fact is frequently used to consider a Hermitean EGPE obtained by just neglecting the imaginary terms derived from $f$.

The parameters

$$\xi=\hbar\sqrt{\frac{3}{2}\frac{\sqrt{g_{bb}}/m_{a}+\sqrt{g_{aa}}/m_{b}}{|\delta g|\sqrt{g_{aa}}n_{a}^{(0)}}},\quad\tau=\frac{3\hbar}{2}\frac{\sqrt{g_{aa}}+\sqrt{g_{bb}}}{|\delta g|\sqrt{g_{aa}}n_{a}^{(0)}}$$

(47)

are identified as the natural units of length and time respectively. In these equations, the saturation density $n_{\alpha}^{(0)}$ of a quantum droplet in the limit of an infinite number $N$ of atoms for the homo–nuclear case is

$$n_{\alpha}^{(0)}=\frac{25\pi}{1024}\frac{(\mathrm{a}_{ab}+\sqrt{\mathrm{a}_{aa}\mathrm{a}_{bb}})^{2}}{\mathrm{a}_{aa}\mathrm{a}_{bb}\sqrt{\mathrm{a}_{\alpha\alpha}}(\sqrt{\mathrm{a}_{aa}}+\sqrt{\mathrm{a}_{bb}})^{5}},$$

(48)

while for the hetero–nuclear case

$$n_{\alpha}^{(0)}=\frac{25\pi}{1024}\frac{1}{f^{2}(m_{b}/m_{a},1,\sqrt{g_{bb}/g_{aa}})}\frac{1}{\mathrm{a}_{\alpha\alpha}^{3}}\frac{\delta g^{2}}{g_{aa}g_{bb}}.$$

(49)

The self-trapping regime requires a minimum number of atoms $N_{c}^{(\alpha)}$ which can be evaluated by calculating the ground state numerically. There is a transient regime characterized by a set of additional critical numbers $N_{ql}^{(\alpha)}$. For $N_{c}^{(\alpha)}<N^{(\alpha)}<N_{ql}^{(\alpha)}$ the atomic density exhibits a surface with a thickness similar to the radius of the droplet. In this regime the concept of surface excitations is questionable and self-evaporation occurs. For higher $N$ values, the ground state corresponds to that expected for a quantum liquid: the internal density is almost uniform and the thickness of its surface is much smaller than the droplet radius. Within the LHY formalism, the values of $N_{c}^{(\alpha)}$ and $N_{ql}^{(\alpha)}$ depend only on the $z,u,x$ variables that determine the $f$ function.

Most of our calculations will be limited to the case $u\sim 1$ using Eq. (46), so that the extended Gross-Pitaevskii equation (EGPE) for the binary Bose-Bose mixtures is taken as

		$\displaystyle i\hbar\partial_{t}\Psi_{\alpha}=$		(50)
		$\displaystyle\Big{(}-\frac{\hbar^{2}}{2m_{\alpha}}\nabla^{2}+g_{\alpha\alpha}|\Psi_{\alpha}|^{2}+g_{\alpha\beta}|\Psi_{\beta}|^{2}+$
		$\displaystyle\frac{4}{3\pi^{2}}\frac{m_{\alpha}^{3/5}g_{\alpha\alpha}}{\hbar^{3}}(m_{\alpha}^{3/5}g_{\alpha\alpha}|\Psi_{\alpha}|^{2}+m_{\beta}^{3/5}g_{\beta\beta}|\Psi_{\beta}|^{2})^{3/2}\Big{)}\Psi_{\alpha}.$

Within the LHY approximation and a symmetric mixture $N^{(a)}$ = $N^{(b)}=N/2$, $m_{a}=m_{b}$ and $\mathrm{a}_{aa}=\mathrm{a}_{bb}$, the equation of state of the droplet reads [2]

$$\frac{E}{E_{0}}=-3\Big{(}\frac{\rho}{\rho_{0}}\Big{)}+2\Big{(}\frac{\rho}{\rho_{0}}\Big{)}^{3/2}$$

(51)

with the equilibrium energy $E_{0}$ and density $\rho_{0}$ given by

$$\rho_{0}=\frac{25\pi(\mathrm{a}_{aa}+\mathrm{a}_{ab})^{2}}{1638\mathrm{a}_{aa}^{5}},\quad\frac{E_{0}}{N}=-\frac{25\pi^{2}\hbar^{2}|\mathrm{a}_{aa}+\mathrm{a}_{ab}|^{3}}{49152m\mathrm{a}_{aa}^{5}}.$$

Alternative formulations to that of LHY corresponds to performing quantum Monte Carlo [18] or density functional methods [12, 19] to study dilute Bose-Bose mixtures of atoms with attractive interspecies and repulsive intraspecies interactions at $T=0$. Such calculations have been reported for homo–nuclear mixtures considering diverse finite range interaction potentials. The results are considered to be valid on the bulk, and can be compared to those resulting from LHY in terms of the equation of state for $u=1$. Second order density Monte Carlo calculations [12] for homo–nuclear mixtures of atoms in two different hyperfine states give rise to a generalization of Eq. (51) with the structure

	$\displaystyle\frac{E}{N}$	$\displaystyle=$	$\displaystyle\frac{|E_{0}|}{N}\Big{[}-3\Big{(}\frac{\rho}{\rho_{0}}\Big{)}+\beta\Big{(}\frac{\rho}{\rho_{0}}\Big{)}^{\gamma}\Big{]}$		(52)
	$\displaystyle\beta$	$\displaystyle=$	$\displaystyle 1.956\frac{\mathrm{a}_{ab}}{\mathrm{a}_{aa}}+\Big{(}0.231+0.236\frac{\mathrm{a}_{ab}}{\mathrm{a}_{aa}}\Big{)}\frac{\mathrm{R}(\mathrm{a}_{ab},r_{\mathrm{eff}})}{\mathrm{a}_{aa}}$
	$\displaystyle\gamma$	$\displaystyle=$	$\displaystyle 1.83+0.32\frac{\mathrm{a}_{ab}}{\mathrm{a}_{aa}}+0.030\Big{(}1+\frac{\mathrm{a}_{ab}}{\mathrm{a}_{aa}}\Big{)}\frac{\mathrm{R}(\mathrm{a}_{ab},r_{\mathrm{eff}})}{\mathrm{a}_{aa}}$		(53)

and $\mathrm{R}$ a function which depends on an effective range parameter $r_{\mathrm{eff}}$. This function can be accessed numerically [12]. The critical number of atoms $N_{c}^{(\alpha)}$ and $N_{ql}^{(\alpha)}$ will now also depend on $\mathrm{R}$. These calculations have a universal character similar to that inferred from EGPE based on LHY approximation whenever the effective range is included as a parameter. Using the local density approximation, the following equation for the corresponding order parameter $\Psi$ is found,

		$\displaystyle i\hbar\partial_{t}\Psi(\vec{x})=\Big{(}-\frac{\hbar^{2}}{2m}\nabla^{2}+\frac{25\pi^{2}\hbar^{2}|\mathrm{a}_{aa}+\mathrm{a}_{bb}|^{3}}{49152m\mathrm{a}_{aa}^{5}}\times$		(54)
		$\displaystyle\Big{[}-6\frac{N^{2}|\Psi|^{4}}{\rho_{0}}+\beta(\gamma+1)\left(\frac{N|\Psi|^{2}}{\rho_{0}}\right)^{\gamma}\Big{]}\Big{)}\Psi(\vec{x}).$

For hetero-nuclear quantum droplets, a generatization of Eq. (52) has been reported in Ref. [19] for Rb-K mixtures and particular values of the $\mathrm{a}_{\alpha\beta}$ coupling parameters.