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 superfluid12, and which was later extended to Fermi systems with spherical symmetry13 having in mind the liquid drop model of nuclei. This formalism is based on the random phase approximation (RPA) and the Thouless variational theorem20. RPA provides a strategy to decouple the collective motion of a many-body system from their individual motion21. In our case, RPA is equivalent to the linearization of the time-dependent equations that define the dynamics of the system in the vicinity of a variational minimum. A contemporary view of the theoretical basis that supports the Thouless theorem and its relevance for many-body theory is revisited in Ref.22. To study quantum droplets arising from vacuum fluctuations in a binary mixture of Bose gases, it is useful to implement Thouless variational principle to such a mixture.
Thouless variational theorem for a binary mixture of Bose gases
In an analysis of the collective modes of nuclei in RPA, Thouless20 showed that it is always possible to find a solution of the independent particle model that makes these collective modes as stable, i.e., they have real eigenfrequencies. This occurs because the solution that gives the lowest expectation value of the Hamiltonian also stabilizes the collective modes in RPA. The stability condition allows the construction of the variational principle and assures that the eigenvectors form a complete set in most cases; this assertion requires the proper definition of a scalar product. 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.
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} , \ldots ,{\vec{r}}_{{N^{{(a)}} }} ;{\vec{r}}_{1}^{\prime } , \ldots ,\vec{r}_{{N^{{(b)}} }}^{\prime } ) = \left[ {\prod\limits_{{i = 1}}^{{N^{{(a)}} }} {\chi ^{{(a)}} } (\vec{r}_{i} )} \right]\left[ {\prod\limits_{{i = 1}}^{{N^{{(b)}} }} {\chi ^{{(b)}} } (\vec{r}_{i}^{\prime } )} \right].{\text{ }} $$
(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
$$i\hbar \partial _t \psi ^{(\alpha )}(\vec{r},t)= {\hat{H}}_0^\alpha \psi ^{(\alpha )}(\vec{r},t) + \hat{ U}^\alpha \psi ^{(\alpha )}(\vec{r},t) $$
(3a)
$${\hat{H}}_0^\alpha = -\frac{\hbar ^2}{2m_\alpha }\nabla ^2 + V^\alpha _{ext}(\vec{r}),$$
(3b)
$$\hat{ U}^\alpha \psi ^{(\alpha )}(\vec{r},t)= \int dr^{\prime 3} \Big ({\mathcal {U}}_{\alpha \alpha }(\vec{r}, \vec{r}^\prime ;\rho ^{(\alpha )}(\vec{r}^\prime ,t)) + {\mathcal {V}}_{\alpha \beta }(\vec{r}, \vec{r}^\prime ;\rho ^{(\alpha )}(\vec{r}^\prime ,t),\rho ^{(\beta )}(\vec{r}^\prime ,t))\Big )\psi ^{(\alpha )}(\vec{r},t) $$
(3c)
In Eq. (3a), \(V^\alpha _{ext}(\vec{r})\) is the external potential and \(\rho ^{(\alpha )}=\psi ^{(\alpha )*}\psi ^{(\alpha )}\) the density of \(\alpha \)-atoms. The real functionals \({\mathcal {U}}_{\alpha \alpha }\) and \({\mathcal {V}}_{\alpha \beta }\) in the integral operator \(\hat{U}^\alpha \) may incorporate, in an effective way, interactions beyond the standard mean field approximation. The functional \({\mathcal {V}}_{\alpha \beta }\) represents the interaction between different species, i.e. \(\alpha \ne \beta \). In the case of EPGE with LHY, the density dependent interactions are superpositions of contact terms2,
$$\begin{aligned} {\mathcal {U}}_{\alpha \alpha }&= g_{\alpha \alpha }\rho ^{(\alpha )}(\vec{r}^\prime ) \delta ^3(\vec{r} – \vec{r}^\prime )\end{aligned}$$
(4a)
$$\begin{aligned} {\mathcal {V}}_{\alpha \beta }&= \left[ g_{\alpha \beta }\rho ^{(\beta )}(\vec{r}^\prime ) +\frac{4m_\alpha ^{3/5}g_{\alpha \alpha }}{3\pi ^2\hbar ^3}\left( m_\alpha ^{3/5}g_{\alpha \alpha }\rho ^{(\alpha )}(\vec{r}^\prime )+m_\beta ^{3/5}g_{\beta \beta }\rho ^{(\beta )}(\vec{r}^\prime )\right) ^{3/2}\right] \delta ^3(\vec{r} – \vec{r}^\prime ). \end{aligned}$$
(4b)
The stationary ground state solutions of Eq. (3) define the chemical potentials \(\mu _\alpha \),
$$\begin{aligned} \psi ^{(\alpha )}_0 (\vec{r},t)= e^{-i\mu _\alpha t/\hbar }\phi ^{(\alpha )}_0(\vec{r}),\quad \quad \alpha = a, b. \end{aligned}$$
(5)
Let us consider collective Bogolubov excitations with a harmonic time dependence,
$$\begin{aligned} \psi ^{(\alpha )}_{exc} (\vec{r},t)= e^{-i\mu _\alpha t/\hbar }\left( \phi ^{(\alpha )}_0(\vec{r}) + \sqrt{N^{(\alpha )}} \left( u_q^{(\alpha )}(\vec{r}) e^{-i\omega _q t} + v_q^{(\alpha )*}(\vec{r}) e^{i\omega _q t}\right) \right) . \end{aligned}$$
(6)
In this equation q is a label that determine the characteristics of the excitation derived, e. g., from its geometry. 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
$$\begin{aligned} \int d^3r \left[ u_q^{(\alpha )}(\vec{r})u_p^{(\alpha ) *}(\vec{r}) -v_q^{(\alpha )}(\vec{r})v_p^{(\alpha ) *}(\vec{r})\right] = \frac{\omega _q}{\vert \omega _q\vert } \delta _{qp}, \quad \omega _q\ne 0. \end{aligned}$$
(7)
Excitations with \(\omega _q = 0\) are assumed to be orthogonal to all other excitations with a different set of quantum numbers.
If \(\phi ^{(\alpha )}_0\) is a real function, the linear response approximation around the stationary function \(\psi ^{(\alpha )}_0 (\vec{r},t)\) is equivalent to the RPA approximation. Taking into account the real character of functionals \({\mathcal {U}}_{\alpha \alpha }\) and \({\mathcal {V}}_{\alpha \beta }\), and factorizing the terms that are either proportional to \(e^{ i\omega _q t}\) or \(e^{-i\omega _q t}\), Eq. (3) in the RPA approximation is equivalent to the equations
$$\begin{aligned}&\hbar \omega _q \zeta _q^{(\alpha )}= {\hat{\Delta }}_0^\alpha \eta _q^{(\alpha )} +2 \phi ^{(\alpha )}_0\int d^3r^\prime \left[ \left[ \frac{\delta {\mathcal {U}}_{\alpha \alpha }}{\delta \rho ^{(\alpha )}} + \frac{\delta {\mathcal {V}}_{\alpha \beta }}{\delta \rho ^{(\alpha )}}\right] \delta \rho ^{(\alpha )}_{q;eff}(\vec{r}^\prime ) \right] +2\phi ^{(\alpha )}_0\int d^3r^\prime \left[ \frac{\delta {\mathcal {V}}_{\alpha \beta }}{\delta \rho ^{(\beta )}}\delta \rho ^{(\beta )}_{q;eff}(\vec{r}^\prime ) \right] ,\end{aligned}$$
(8a)
$$\begin{aligned}&\hbar \omega _q \eta _q^{(\alpha )}={\hat{\Delta }}_0^\alpha \zeta _q^{(\alpha )} ,\end{aligned}$$
(8b)
$$\begin{aligned} \eta _q^{(\alpha )}&= u_q^{(\alpha )} + v_q^{(\alpha )},\,\, \zeta _q^{(\alpha )}=u_q^{(\alpha )}-v_q^{(\alpha )}\end{aligned}$$
(8c)
$$\begin{aligned}&{\hat{\Delta }}_0^\alpha = {\hat{H}}_0^\alpha + \int dr^{\prime 3}\big ({\mathcal {U}}_{\alpha \alpha }(\vec{r}, \vec{r}^\prime ;\rho ^{(\alpha )}_0(\vec{r}^\prime )) + {\mathcal {V}}_{\alpha \beta }(\vec{r}, \vec{r}^\prime ;\rho ^{(\alpha )}_0(\vec{r}^\prime ),\rho ^{(\beta )}_0(\vec{r}^\prime )\big ) – \mu _\alpha , \end{aligned}$$
(8d)
\(\delta \rho ^{(\alpha )}_{q;eff}= \phi ^{(\alpha )}_0\eta _q^{(\alpha )}\) could be directly interpreted as an the effective variation of the density for dynamical purposes since it multiplies the variational derivatives of \({\mathcal {U}}_{\alpha \alpha }\) and \({\mathcal {V}}_{\alpha \beta }\). In fact, using \({\hat{\Delta }}_0^\alpha \phi ^{(\alpha )}_0=0\), the continuity equation
$$\begin{aligned}&\frac{\partial }{\partial t}\left[ \delta \rho ^{(\alpha )}_{q;eff}e^{-i\omega _q t}\right] +\vec{\nabla }\cdot \vec{J}_{q;eff}^{(\alpha )} =0, \end{aligned}$$
(9a)
$$\begin{aligned}&\vec{J}_{q;eff}^{(\alpha )} =\frac{1}{2m_\alpha }\left[ \eta _q^{(\alpha )} \left( -i\hbar \vec{\nabla }\right) \phi ^{(\alpha )}_0 -\phi ^{(\alpha )}_0\left( -i\hbar \vec{\nabla }\right) \eta _q^{(\alpha )}\right] , \end{aligned}$$
(9b)
is derived; it implies a conservation relation for both the real and imaginary parts of \(\int d^3x \delta \rho ^{(\alpha )}_{q;eff}e^{-i\omega _q t}\).
Thouless variational theorem20 is now applied. In the system under consideration, the application of this theorem is based on the stability inherited to \(\psi ^{(\alpha )}_0 (\vec{r},t)\) from the stationary character of the ground state solutions of Eq. (3). The energy \(\hbar \omega _q\) is estimated after from Eq. (8) solving for \(\omega _q\), and averaging the estimates obtained from each atoms species, \(\alpha = a,b\). The variational condition is then written as
$$\begin{aligned}&\hbar \omega _{exact}\le \frac{1}{2} \frac{\sum _{\alpha ,\beta =a,b} \langle \eta _q^{(\beta )}\vert {\hat{\Xi }}^{\alpha \beta }\vert \eta _q^{(\alpha )}\rangle + \sum _{\alpha = a,b}\langle \zeta _q^{(\alpha )}\vert {\hat{\Delta }}_0^\alpha \vert \zeta _q^{(\alpha )}\rangle }{\sum _{\alpha =a,b} \langle \zeta _q^{(\alpha )}\vert \eta _q^{(\alpha )}\rangle } =:\hbar \omega _{\eta \zeta }\end{aligned}$$
(10a)
$$\begin{aligned} \langle \eta _q^{(\beta )}\vert {\hat{\Xi }}^{\alpha \beta }\vert\eta _q^{(\alpha )}\rangle &= \delta _{\alpha \beta } \langle \eta _q^{(\beta )}\vert {\hat{\Delta }}_0^\alpha \vert \eta _q^{(\alpha )}\rangle + 2 \int d^3r \eta _q^{(\beta )*}(\vec{r}) \phi _0^{(\alpha )}(\vec{r}) \int d^3r^\prime \left( \frac{\delta {\mathcal {U}}_{\alpha \alpha }}{\delta \rho ^{(\alpha )}} +\frac{\delta {\mathcal {V}}_{\alpha \beta }}{\delta \rho ^{(\alpha )}}\right) \eta _q^{(\alpha )}(\vec{r}^\prime ) \phi _0^{(\alpha )}(\vec{r}^{\prime }) \Big ] \nonumber \\&\quad +2 \int d^3r \eta _q^{(\beta )*}(\vec{r}) \phi _0^{(\alpha )}(\vec{r}) \int d^3r^\prime \left[ \frac{\delta {\mathcal {V}}_{\alpha \beta }}{\delta \rho ^{(\beta )}} \eta _q^{(\beta )}(\vec{r}^{\prime }) \phi _0^{(\beta )}(\vec{r}^{\prime })\right] . \end{aligned}$$
(10b)
The exact expressions for \(u_q\) and \(v_q\) would satisfy Eqs. (8) and they are unknown. In a variational calculation, an ansatz for the variational functions \(\eta _q^{(\alpha )}\) and \(\zeta _q^{(\alpha )}\), Eq. (8c), is done incorporating in as much as possible the expected physical characteristics of the particular system under consideration. Our ansatz for the excitation functions depends functionally on \(\phi _0\) and its derivatives and takes as variational parameters a relative weight between \(\eta _q^{(\alpha )}\) and \(\zeta _q^{(\alpha )}\) for each atomic species \(\alpha =a,b\),
$$\begin{aligned} \eta _q^{(\alpha )} = \eta _q^{(\alpha )}\left[ \phi _0^{(\alpha )},\partial \phi _0^{(\alpha )}\right] ,\quad \zeta _q^{(\alpha )} = \gamma _\alpha {\tilde{\zeta }}_q^{(\alpha )}\left[ \phi _0^{(\alpha )},\partial \phi _0^{(\alpha )}\right] . \end{aligned}$$
(11)
By requiring that \(\gamma _\alpha \) yields a minimum for \(\omega _{\eta \zeta }\) it results
$$ \gamma _a= C_a\sqrt{\frac{AB_b}{B_a\left( B_b C_a^2 + B_aC_b^2\right) }},\quad \gamma _b =\frac{B_aC_b}{B_bC_a}\gamma _a, $$
(12)
$$A= \sum _{\alpha ,\beta =a,b} \left\langle \eta _q^{(\alpha )}\vert {\hat{\Xi }}^{\alpha \beta }\vert \eta _q^{(\beta )}\right\rangle ,\quad B_\alpha = \left\langle {\tilde{\zeta }}^{(\alpha )}_q\vert {\hat{\Delta }}_0^\alpha \vert {\tilde{\zeta }}^{(\alpha )}_q\right\rangle ,\quad C_\alpha =\left\langle {\tilde{\zeta }}_q^{(\alpha )}\vert \eta _q^{(\alpha )}\right\rangle . $$
(13)
and the corresponding optimal variational excitation energy is then
$$\begin{aligned} \hbar ^2\omega _{opt}^2 = \frac{AB_aB_b}{B_aC_b^2 + B_bC_a^2}. \end{aligned}$$
(14)
A pair of variational ansatz for surface excitations
Since quantum droplets exhibit different regimes, e. g. compressible and incompressible, we decided to explore the relevance of two pairs of excitation functionals \(\eta ^{(\alpha )}_q \) and \({\tilde{\zeta }}^{(\alpha )}_q \), both of them with the simple structure,
$$\begin{aligned} \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 )}. \end{aligned}$$
(15)
The functions \(\vec{F}^\alpha (\vec{r})\) and \( \Lambda ^\alpha (\vec{r})\) incorporate geometric properties of the quantum fluid. The parameter \(\xi ^{(\alpha )}\) corresponds to a characteristic scale length for the \(\alpha \)-species. For these expressions for \(\eta ^{(\alpha )}_q\) and \(\zeta ^{(\alpha )}_q=\gamma _\alpha {\tilde{\zeta }}^{(\alpha )}_q\),
$$\begin{aligned} \delta \rho _{q;eff}^{(\alpha )} = \frac{\xi ^{(\alpha )}}{2} \vec{F}_q^\alpha \cdot \vec{\nabla }\left| \phi _0^{(\alpha )}\right| ^2, \quad \vec{J}_{q;eff}^{(\alpha )} =\left| \phi _0^{(\alpha )}\right| ^2\left[ \frac{1}{2} \gamma _\alpha \left( i\frac{\hbar }{m_\alpha }\vec{\nabla }\right) \Lambda ^\alpha _q\right] =:\vert \phi _0^{(\alpha )}\vert ^2 \vec{v}^{(\alpha )}. \end{aligned}$$
(16)
Thus, each component of \( \vec{F}_q^\alpha \) would determine the variations of the density \(\vert \phi _0^{(\alpha )}\vert ^2\) along different spatial directions, and the “velocity” field \(\vec{v}^{(\alpha )}(\vec{r})\) is determined by the gradient of \(\Lambda _q^\alpha \), i.e. \(\vec{v}^{(\alpha )}(\vec{r})\propto \vec{\nabla }\Lambda _q^\alpha \).
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). For a droplet \(\Lambda ^\alpha _q\) is taken as the solution of Laplace equation in spherical coordinates that is regular at the origin
$$\begin{aligned} \Lambda _{\ell m}^\alpha (\vec{r}) = {\mathcal {A}}_{\ell m} Y_{\ell m}(\theta ,\varphi ) \left( \frac{r}{\xi ^{(\alpha )}}\right) ^\ell , \end{aligned}$$
(17)
with \({\mathcal {A}}_{\ell m}\) a normalization constant to be evaluated according to Eq. (7). Even though in the present study we focus on quantum droplets, the bubble configuration could result useful for the description of the mixed-bubble regime in which bubbles of the mixed phase coexist with a pure phase of one of the components23. In such a case, one would take \(\Lambda _{\ell m}^\alpha (\vec{r}) ={\mathcal {A}}_{\ell m} Y_{\ell m}(\theta ,\varphi ) (r/\xi ^{(\alpha )})^{-(\ell +1)}\).
Concerning the function \(\vec{F}^\alpha _{\ell m}(\vec{r})\), two reasonable ansatz for spherical boundaries are identified,
Ansatz 1.
$$\begin{aligned} \vec{F}^\alpha _{\ell m}(\vec{r}) = {\mathcal {A}}^{(1)}_{\ell m} Y_{\ell m}(\theta ,\varphi ){\hat{r}} \Rightarrow \eta ^{(\alpha )}_{\ell m} ={\mathcal {A}}^{(1)}_{\ell m} Y_{\ell m}\xi ^{ (\alpha )}\partial _r \phi ^{ (\alpha )}_0\; \text {and}\;\zeta ^{(\alpha )}_{\ell m}=\gamma _\alpha {\mathcal {A}}^{(1)}_{\ell m} Y_{\ell m}(r/\xi ^{ (\alpha )})^\ell \phi ^{ (\alpha )}_0. \end{aligned}$$
(18)
It is equivalent to an initially infinitesimal movement of the droplet surface in the radial direction with an amplitude proportional to the angular solution of the wave equation, i.e., a spherical harmonic. This ansatz resembles that discussed by Berstch12 for a quantum fluid whose surface is defined by the \(z=0\) plane and for which \(\eta = \partial _z\phi _0 e^{ikx}\).
Ansatz 2.
$$\begin{aligned} \vec{F}^\alpha _{\ell m}(\vec{r})= \xi ^{(\alpha )}\vec{\nabla }\Lambda _{\ell m}^\alpha \Rightarrow \eta ^{(\alpha )}_{\ell m} ={\mathcal {A}}^{(2)}_{\ell m} Y_{\ell m}[(\xi ^{ (\alpha )}\partial _r) (r/\xi ^{ (\alpha )})^\ell ] [(\xi ^{ (\alpha )}\partial _r) \phi ^{ (\alpha )}_0]\; \text {and}\;\zeta ^{(\alpha )}_{\ell m}=\gamma _\alpha {\mathcal {A}}^{(2)}_{\ell m} Y_{\ell m}(r/\xi ^{ (\alpha )})^\ell \phi ^{ (\alpha )}_0. \end{aligned}$$
(19)
This ansatz is equivalent to \(\delta \rho _{q;eff}^{(\alpha )} \sim \vec{v}^{(\alpha )}\cdot \nabla \vert \phi _0^{(\alpha )}\vert ^2\sim d\vert \phi _0^{(\alpha )}\vert ^2/dt\), i. e., the surface of the droplet is displaced along the velocity field. Note that the second ansatz privileges an hydrodynamic interpretation, Eq. (19), while the first ansatz at \(t=0\) yields an infinitesimal deformation of the droplet perpendicular to its surface, Eq. (18).
A direct calculation shows that, in case the dynamics is governed by contact interactions as it occurs for the EGPE, the contribution of the interaction terms of these excitation functions to A and \(B_\alpha \), Eq. (13), is null. However, information about such interaction terms is still contained in the stationary solution \(\phi _0^{(\alpha )}\). For quantum droplets that can be generated without an external potential, the ground state order parameters depend only on the radial coordinate and are proportional to each other,
$$\begin{aligned} \phi ^{(a)}_0 = \sqrt{\frac{N^{(a)}}{N^{(b)}}}\phi ^{(b)}_0 \end{aligned}$$
(20)
the last equality being a consequence of Eq. (2). Under these conditions, the excitation energies become
Ansatz 1.
$$\begin{aligned} \hbar ^2 \omega ^2 \le – \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 } =\hbar ^2 \omega ^2_{A1}=\epsilon ^2_{A1} \end{aligned}$$
(21)
Ansatz 2.
$$\begin{aligned} \hbar ^2 \omega ^2 \le -\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} }=\hbar ^2 \omega ^2_{A2}=\epsilon ^2_{A2}, \end{aligned}$$
(22)
whenever an equal length scale
$$\begin{aligned} \xi ^{(a)} = \xi ^{(b)} =\xi , \end{aligned}$$
(23)
is chosen for both species. In these equations we notice the presence of the total mass of the droplet
$$\begin{aligned} M_T = N^{(a)}m_a + N^{(b)}m_b. \end{aligned}$$
(24)
Surface tension
In his seminal work24, 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
$$\begin{aligned} \omega _{exc}^2 = \ell (\ell -1)(\ell +2)\frac{\sigma }{m\rho _0R_0^3} \end{aligned}$$
(25)
and he 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
$$\begin{aligned} m\rho _0R_0^3 =\frac{3}{4\pi } M_T \end{aligned}$$
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 \ge 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:
-
1.
For ansatz 1, the \(\ell \)-dependence of the coefficient resembles that of Rayleigh and that of the liquid drop model of nuclei in the limit of \(N\rightarrow \infty \) nucleons25.
-
2.
For ansatz 2, the \(\ell \)-dependence of the coefficient is similar to that predicted for atomic nuclei by Bertsch13 and by Casas and Stringari14 using RPA and the density-density Green’s function formalism; in both cases, the structure of 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. Once adequate expressions for the order parameters \(\phi ^{(a)}_0\) are given, the approximate excitation energies are evaluated using both ansatz. The variational criteria establishes that the best approximation to the exact description of those excitations is that which provides the lowest excitation energies.
Following Lord Rayleigh ideas, the corresponding surface tension estimation would be given by the expressions
Ansatz 1.
$$\begin{aligned} \sigma ^{(A1)}_\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 } \end{aligned}$$
(26)
Ansatz 2.
$$\begin{aligned} \sigma ^{(A2)}_\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} } \end{aligned}$$
(27)
with
$$\begin{aligned} M = \frac{4\pi }{3}\frac{m_am_b}{N^{(a)}m_b + N^{(b)}m_a}. \end{aligned}$$
The density of a spherical liquid droplet conformed by N atoms is characterized by a radius \(R_0\), a skin width dR and a saturation density. The generic shape of such a droplet can be approximately described by a Boltzmann function,
$$\begin{aligned} \rho _B(R;N) = \frac{{\mathcal {B}}_1}{1 + \mathrm {exp}((R-R_0)/dR)}. \end{aligned}$$
(28)
For Bose–Bose mixtures of experimental interest and number of atoms ranging from \(10^4\)–\(10^7\), we have performed numerical calculations for the ground state derived from an EGPE for homo- and hetero-nuclear mixtures, and, for the homo-nuclear case, from the effective formalism that considers finite-range effects15,16. The resulting density profiles fits the Boltzmann density with a reliability above 0.998 in all considered cases. The resulting parameters are reported in the Supplementary Material. Using Eq. (28), the surface tension Eqs. (26–27) can be written in terms of the Fermi-Dirac integrals
$$\begin{aligned} F_s (z)= \frac{1}{\Gamma (s+1)}\int _0^\infty \frac{dx x^s}{1+e^{x-z}}. \end{aligned}$$
(29)
So that, for \(\ell > 1\)
$$ \sigma _\ell ^{(A1)}= \frac{\hbar ^2}{8M dR^4}\frac{(2\ell +1)!}{((\ell +2)!)^2}\left[ 1-\frac{1}{(1+ e^{R_0/dR})^2}\right] \frac{F_{2\ell }(R_0/dR)}{(F_{\ell +1}(R_0/dR))^2}, $$
(30)
$$\sigma _\ell ^{(A2)}= \frac{\hbar ^2}{16M dR^4}\frac{1}{\ell (4\ell ^2 -1)} \frac{F_{2\ell -3}(R_0/dR)+F_{2\ell -4}(R_0/dR)}{F_{2\ell }(R_0/dR)}, $$
(31)
expressions that correspond to ansatz 1 and 2 respectively.
Weber number for Boltzmann shaped droplets
Consider the collision of quantum droplets that initially were radially spaced at a distance \(2d_0\); they are assumed to be kicked off as described by a plane wave wavefunction factor,
$$\begin{aligned} \Psi (\vec{r}, t=0)&= \Psi _1(\vec{r}) + \Psi _2(\vec{r}) \nonumber \\&= \begin{pmatrix}\psi _{a_1}(\vec{r} + \vec{d}_0)\\ \psi _{b_1}(\vec{r}+ \vec{d}_0)\end{pmatrix} e^{i\vec{k}_0\cdot \vec{r}/2} +\begin{pmatrix}\psi _{a_2}(\vec{r} – \vec{d}_0)\\ \psi _{b_2}(\vec{r}- \vec{d}_0\end{pmatrix} e^{-i\vec{k}_0\cdot \vec{r}/2}. \end{aligned}$$
(32)
The mean value
$$\begin{aligned} {\mathcal {K}} = -\frac{\hbar ^2}{2}\sum _{i=1,2}\int d^3r \Psi _i^\dagger (\vec{r})\begin{pmatrix} \frac{\nabla ^2}{m_{a_i}} &{} 0\\ 0 &{} \frac{\nabla ^2}{m_{b_i}} \end{pmatrix} \Psi _i(\vec{r})\end{aligned}$$
(33)
is taken as a measure of the initial kinetic energy of the droplets. Assuming a non significant overlap between the droplets at \(t=0\), it can be decomposed as the sum of the translational kinetic energy of each droplet as a whole \({\mathcal {K}}_{trans}\), and an internal kinetic energy of the atoms within the droplets \({\mathcal {K}}_{int}\). Thus, \({\mathcal {K}}\approx {\mathcal {K}}_{trans} + {\mathcal {K}}_{int}\), with
$$\begin{aligned} {\mathcal {K}}_{trans}= \frac{\hbar ^2}{2}\left[ \frac{N^{(a_1)}k_0^2}{4m_{a_1}} + \frac{N^{(b_1)}k_0^2}{4m_{b_1}} + \frac{N^{(a_2)}k_0^2}{4m_{a_2}} + \frac{N^{(b_2)}k_0^2}{4m_{b_2}} \right] ,\quad {\mathcal {K}}_{int} = -\frac{\hbar ^2}{2}\sum _{i=1,2}\sum _{\alpha =a,b}\int d^3r\psi ^*_{\alpha _i}(\vec{r})\frac{\nabla ^2}{m_{\alpha _i}}\psi _{\alpha _i}(\vec{r}). \end{aligned}$$
(34)
The dimensionless quantity suitable for characterizing binary collisions of classical incompressible droplets is the Weber number19. It is a dimesionless quantity for each one of the droplets, defined as the ratio between the translational kinetic energy \({\mathcal {K}}_{trans}\) before the collision and the surface energy of excitation evaluated in terms of the surface tension of the droplet. This definition can be extended to quantum droplets as
$$\begin{aligned} \mathrm {We}_\ell =\frac{{\mathcal {K}}_{trans}}{R_0^2\sigma _\ell } \end{aligned}$$
(35)
where \(R_0\) is the initial radius of the droplet. For the collision of two identical quantum droplets like those described in this work, this expression becomes
Ansatz 1
$$\begin{aligned} \mathrm {We}_\ell ^{(A1)} = \frac{8\pi (dR k_0)^2}{3}\left( \frac{dR}{R_0}\right) ^2 \frac{((\ell +2)!)^2}{(2\ell +1)!}\left[ 1-\frac{1}{(1+ e^{R_0/dR})^2}\right] ^{-1} \left[ \frac{(F_{\ell +1}(R_0/dR))^2}{F_{2\ell }(R_0/dR)}\right] , \end{aligned}$$
(36)
Ansatz 2
$$\begin{aligned} \mathrm {We}_\ell ^{(A2)}=\frac{16\pi (dR k_0)^2}{3}\left( \frac{dR}{R_0}\right) ^2 (\ell (4\ell ^2 -1)) \frac{F_{2\ell }(R_0/dR)}{F_{2\ell -3}(R_0/dR)+F_{2\ell -4}(R_0/dR)}. \end{aligned}$$
(37)
Notice that the expressions of \(\mathrm {We}_\ell \) for the ground state of the quantum droplets given by Eqs. (36–37) make evident the relevance of the relation between the de Broglie wavelength and the skin width of the droplet, \(k_0 dR\), as a measure of the momentum that could be transfered to the droplet during the collision. The second relevant parameter is the ratio \(R_0/dR\) that determines the deformation of the droplet due to a multipole excitation.