Self-bound ultra dilute Bose mixtures within Local Density Approximation

We have investigated self-bound binary ultra dilute bosonic mixtures at zero temperature within Density Functional Theory using a Local Density Approximation. We provide the explicit expression of the Lee-Huang-Yang correction in the general case of heteronuclear mixtures, and investigate the general thermodynamic conditions which lead to the formation of self-bound systems. We have determined the conditions for stability against the evaporation of one component, as well as the mechanical and diffusive spinodal lines. We have also calculated the surface tension of the self-bound state as a function of the inter-species interaction strength. We find that relatively modest changes of the latter result in order-of-magnitude changes in the calculated surface tension. We suggest experimental realizations which might display the metastability and phase separation of the mixture when entering regions of the phase diagram characterized by negative pressures. Finally, we show that these droplets may sustain stable vortex and vortex dimers.


I. INTRODUCTION
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 very recently [2,3]. The stability of these droplets can be explained as deriving from a subtle balance between the intra-species repulsive atom-atom interaction, and a tunable inter-species attractive interaction.
When the attraction between the two atomic species becomes larger than the single-species average repulsion, the mixture is expected to collapse according to meanfield (MF) theory. Possible stabilization mechanisms preventing the collapse of the mixture triggered by an attractive interaction have been proposed, as e.g. three-body correlations [4], spin-orbit coupling [5] and quantum fluctuations (Lee-Huang-Yang (LHY) mechanism [6]). In the latter case, the effective repulsion provided by the first beyond-mean-field (BMF) correction to the energy is enough to prevent the collapse and to stabilize the system, i.e. a density exists where these contributions balance each other and the droplets become self-bound, stable systems. The stabilization mechanism resulting from the inclusion of the LHY correction is also responsible for the existence of self-bound aggregates in one component dipolar systems [7][8][9][10] (where the anisotropic character of the dipole-dipole forces leads to the formation of filament-like self-bound droplets with highly anisotropic properties), in Rabi-coupled Bose-Bose mixtures [11] and in low-dimensional mixtures [12].
The formation of liquid drops in a Bose-Bose mixture has been recently addressed using the diffusion Monte Carlo (DMC) method [13], confirming the prediction for the stability of self-bound bosonic mixtures [1]. More recently, the properties of uniform Bose mixtures have been analyzed using the variational hypernetted-chain Euler Lagrange (HNC-EL) method [14], which includes pair correlations non-perturbatively and turns out to be computationally very fast as compared to the DMC method. In particular, the conditions for having a self-bound, stable mixture of 39 K atoms in two different internal states have been studied within the HNC-EL approach. Deviations from a universal dependence on the s-wave scattering lengths are found in spite of the low density of the systems [14].
These ultra dilute and weakly interacting quantum liquids -whose densities are orders of magnitude lower than that of the prototypical quantum fluid, namely liquid helium-can be ideal platforms to benchmark quantum many-body theories in actual experiments, and to study processes more difficult to address in the case of liquid helium. For this reason, determining the properties of the underlying uniform system is of special interest in itself and also represents a first step towards a better understanding of self-bound droplets.
In this work we study, within a density functional theory (DFT) approach in a Local Density Approximation (LDA), Bose-Bose mixtures at zero temperature (T ); in particular, we address the thermodynamic conditions which lead to the formation of self-bound states of the system. We study the general case of heteronuclear mixtures, which has not been considered before within the beyond-mean-field approach to Bose mixtures. As a case of study we discuss in detail a 23 Na-87 Rb mixture. In order to make a comparison with previous theoretical works, we also address a Bose-Bose mixture with equal masses, made of two hyperfine states of 39 K.
FIG. 1: Solid line is the P = 0 curve in the (n1, n2) plane for the 23 Na-87 Rb mixture with a12 = −80 a0. The big dot is the point of minimum energy per particle of the P = 0 curve. The dotted lines correspond to the configurations in the (n1, n2) plane with chemical potential of one species equal to zero, µ1 = 0 and µ2 = 0, respectively. In the closed narrow region delimited by these lines the system has both chemical potentials negative.
This work is organized as follows. In Sec. II we present the DFT approach to Bose-Bose mixtures. The method is applied in Sec. III, within LDA, to the description of the uniform system, allowing us to determine the main characteristics of the stability phase diagram, in particular the mechanical and diffusive spinodal lines obtained as outlined in the Appendix. The surface tension of the mixture is calculated in Sec. IV as a function of the interspecies attraction, and the structure of selected mixed droplets is presented in Sec. V. We show that these droplets may sustain stable vortex and vortex pairs in Sec. VI. Finally, a summary and outlook are given in Sec. VII.

II. DENSITY FUNCTIONAL APPROACH
Let us consider a uniform Bose-Bose mixture with two components (with masses m 1 and m 2 ) in a volume V , interacting with coupling constants g 11 = 4πa 11h 2 /m 1 , g 22 = 4πa 22h 2 /m 2 , and g 12 = 2πa 12h 2 /m r , where m r = m 1 m 2 /(m 1 + m 2 ) is the reduced mass. The intra-species s-wave scattering lengths a 11 and a 22 are both positive, while the inter-species a 12 is negative. The total number of bosons is N = N 1 + N 2 .
In the following n 1 , n 2 are the number densities, nor-malized such that V n 1 dr = N 1 and V n 2 dr = N 2 . Equivalently, we will also characterize the homogeneous mixture with the total number density, n = n 1 + n 2 , and concentration of the species 2, Y = n 2 /n. Although a precise knowledge of the finite-range details of the interatomic potential might be necessary for a more accurate description of such system [14], we consider the simpler description where the s-wave scattering lengths are assumed to be enough to fully characterize the inter-particle interactions, and leave to subsequent studies attempts to go beyond the contact interaction approximation. This strategy is similar to that successfully followed within the DFT approach to liquid helium and droplets [15][16][17].
Within the DFT framework, the total energy of the system is given by where n i = |Ψ i | 2 and i = 1, 2. E LHY is the BMF Lee-Huang-Yang term, which is necessary in order to yield self-bound configurations [1]. Functional minimization of the above functional leads to the Euler-Lagrange (EL) equations where µ i is the chemical potential of the i-species and Eq.(2) is the two-components version of the well-known Gross-Pitaevskii equation [18] with the addition of the BMF correction. The LHY correction to the mean-field theory of the mixture can be expressed as [1] At the mean-field level, the condensed Bose-Bose mixture collapses when the inter-species attraction becomes stronger than the geometrical average of the intra-species repulsions, |g 12 | > √ g 11 g 22 . Quantum fluctuations, embodied within the LHY energy term, stabilize the mixture. As shown in Ref. [1], the instability manifests itself in the fact that some of the energy contributions in E LHY acquire an imaginary component at small momenta. However, in the region mostly contributing to the LHY term, these modes are found to be insensitive to small variations of δg ≡ g 12 + √ g 11 g 22 , and also to its sign [1]. When evaluated at δg = 0, E LHY in Eq. (4) is well defined and free from imaginary contributions. We will use here the same approximation as in Ref. [1] to evaluate E LHY , namely we set g 12 + √ g 11 g 22 = 0. The explicit expression for f in Eq.(4) was given in Ref. [1] only for the particular case of equal masses. In the more general case m 1 = m 2 , which is addressed here -and considering δg = 0, as discussed above-one finds: where Here z ≡ m 2 /m 1 , x ≡ (g 22 n 2 )/(g 11 n 1 ), and k is dimensionless.
The integral (5) converges in spite of the presence of individually diverging terms in the integrand due to mutual cancellation of the singular terms. For the numerical evaluation of Eq. (4) we found convenient to calculate the improper integral (5) using the following transformation The right-hand side integral has been computed numerically using a Second Euler-McLaurin summation formula refined until some specified degree of accuracy is achieved [19]. In the particular case N 1 = N 2 ≡ N/2, m 1 = m 2 ≡ m, and g 11 = g 22 = g 12 ≡ g = 4πah 2 /m, Eq. (4) yields the well-known LHY correction for a system of N identical bosons in a volume V (i.e. with density n): where we have used Within the Local Density Approximation of DFT one can write where the energy density E LHY is evaluated at the local densities n 1 (r), n 2 (r): (9) The terms appearing in Eq. (3) are As a case of study we consider in the following a uniform mixture of 23 Na (a 11 = 54.5 a 0 ) and 87 Rb (a 22 = 100.4 a 0 ) atoms [20,21], where a 0 is the Bohr radius. We will take the inter-atomic scattering length a 12 as tunable at will. This mixture has been recently studied [22] and proven to be a good candidate to investigate interactiondriven effects in a superfluid Bose mixture with a largely tunable inter-species interactions (both repulsive and attractive). We are interested in the regime where selfbound states appear, i.e. when g 11 , g 22 > 0, g 12 < 0 and |g 12 | > √ g 11 g 22 . The energy per unit volume of the uniform system is and the total energy is E = EV . At T = 0, especially relevant stable states of the homogeneous mixture are those that correspond to zero pressure since isolate self-bound droplets must be at equilibrium with vacuum. Recalling that P (n 1 , n 2 ) = −E + µ 1 n 1 + µ 2 n 2 , and that one obtains for the pressure P (n 1 , n 2 ) = 1 2 where the BMF terms are evaluated according to Eqs. (10) and (11). Figure 1 shows the P = 0 curve in the (n 1 , n 2 ) plane computed with a 12 = −80 a 0 . The big dot represents the stable state configuration that, for this chosen a 12 value, has the minimum energy per atom E/N = E/n.
The densities associated to that minimum energy state are shown in Fig. 2 as a function of the interatomic scattering length a 12 . They have been computed as described above, i.e. selecting the minimum energy states among those satisfying the condition P = 0.
From the results of Fig. 2 it follows that self-bound states appear for a 12 < a c ∼ −62 a 0 . The critical value a c is consistent with that obtained from the condition g 12 = − √ g 11 g 22 (which is the instability condition at the mean-field level, i.e. with E LHY = 0), that is a 12 ∼ −60 a 0 .
In order for an atom of the i-th species to be bound in the mixture the chemical potential must be negative, µ i < 0. If it is not, the energy will be lowered by removing atoms from the system (evaporation). We have thus computed the limiting curves in the (n 1 , n 2 ) plane where the conditions µ 1 = 0, µ 2 = 0 are fulfilled, by using Eqs. (14) for the same value of a 12 considered above. The results are shown in Fig. 1 where the dotted lines represent the configurations with µ 1 = 0 and µ 2 = 0. Only the points within the closed narrow region shown in the figure are stable against evaporation, i.e. only inside this region the system verifies µ 1 < 0 and µ 2 < 0 simultaneously. Similarly, a very narrow stability region against evaporation has been found for the 39 K-39 K mixture within the HNC-EL approach [14].

B. Spinodal lines
Binary mixtures such as those described here are not thermodynamically stable at all densities n = n 1 + n 2 , temperatures and relative concentrations Y = n 2 /n. At T = 0, necessary and sufficient conditions for thermodynamic stability are expressed by the following inequalities [23]: where κ is the inverse compressibility -incompressibilityof the system. A positive incompressibility guarantees mechanical stability, and the condition on the chemical potential derivative guarantees diffusive stability. If one of these conditions is violated the mixture cannot exist as a single phase and must undergo phase separation.
The coexisting phases that appear may have different densities, different concentrations or both. The lines obtained setting to zero the above inequalities are called mechanical and diffusive spinodal lines, respectively. Systems such as nucleonic matter, 3 He-4 He liquid mixtures, and partially polarized liquid 3 He are examples of highly correlated systems for which the spinodal lines were determined in the past by solving similar equations [24][25][26]. The big dot represents the stable, minimum energy per particle mixture at P = 0. The dotted line is the P = 0 curve in the (n, Y ) plane. The two dash-dot curves correspond to µ1 = 0 and µ2 = 0 (also shown in Fig. 1 as a function of n1 and n2).
In our system, a straightforward calculation outlined in the Appendix yields: From this expression the mechanical spinodal line κ = 0 can be readily calculated. Similarly, the diffusive spinodal line is obtained by solving the equation: As outlined in the Appendix, this condition on the (∂µ 1 /∂Y ) P partial derivative can be cast in the following more convenient expression The mechanical and diffusive spinodal lines in the (n, Y ) plane are shown in Fig. 3 for the 23 Na -87 Rb mixture with a 12 = −80 a 0 . It can be seen from this figure that the stability region against evaporation -the narrow region where the chemical potentials are negative-is considerably reduced by thermodynamic stability conditions, and it is represented in the figure by the small, triangularshaped region delimited by the two dash-dot lines and the solid line. In particular, the point representing the stable minimum energy mixture is rather close to the diffusive spinodal line. This point is at P = 0; reducing n at constant Y amounts to decreasing P . Hence, from that point down to the diffusive spinodal line, the mixture is in a metastable state at negative pressure. Under these conditions, bubbles might appear in the mixture (cavitation phenomenon), eventually leading to a first order phase transition as thoroughly studied, both theoretically and experimentally, in liquid helium [27][28][29]. The stability of the 39 K-39 K mixture has been studied using the HNC-EL method [14], and it was found that the condition that limits the thermodynamic stability of such mixture arises from the mechanical and not from the diffusive spinodal, at variance with the Na-Rb mixture just described. Within the present DFT approach, we have also studied the 39 K-39 K mixture under similar conditions as those described in Ref. [14], where finite range interactions were used rather than contact interactions as done here. In particular, we have taken δa = −0.156 (in the same units used in Ref. [14]). The results for the equilibrium densities, chemical potentials and energy per particle are in agreement with their HNC-EL results. In fact we find an equilibrium density ratio n 2 /n 1 = 1.385 to be compared with the HNC-EL result, n 2 /n 1 = 1.380. We find that the total energy per atom is ǫ = −3.74 (in the energy units of Ref. [14]), whereas they find (using an effective range r ef f = 43.2) ǫ = −3.364. Figure 4 shows the phase diagram of the homogeneous 39 K-39 K mixture. As found above for the Na-Rb mixture, the stability region is considerably reduced by thermodynamic conditions, and it is represented in the figure by the triangular-shaped region delimited by the two dashdot lines and the dashed line. At variance with the Na-Rb case, the point representing the stable minimum energy per particle mixture is now close to the mechanical instead to the diffusive spinodal line. This is in agreement with the HNC-EL results [14]. We conclude that, not surprisingly, the stability phase diagram is very sensitive to the parameters defining the mixture. As in the Na-Rb case, bubbles are expected to appear in the 39 K-39 K mixture when the density (thus P ) is decreased from the stable minimum energy per particle point.
The existence of mechanical and diffusive spinodal lines in self-bound Bose-Bose mixtures might cause dynamic instabilities similar to those characterizing the expansion phase of a highly compressed nuclear spot created in the course of an energetic nucleus-nucleus collision, which triggered an enormous activity in the Nuclear Physics field in the 1980's, see e.g. Refs. [30][31][32][33] and references therein. In the BEC case, it is plausible that self-bound mixed droplets compressed by an external trap will expand upon release of the trap, bringing a large portion of the expanding droplet into the unstable region of the phase diagram (e.g. Figs. 3 and 4).
Related effects could be observed in experiments leading to cavitation, similarly to what found in liquid helium [27][28][29]. Cavitation bubbles could be created, e.g., by sweeping a large droplet with a laser beam: the pressure difference due to fore-to-aft asymmetry in the fluid structure around the laser spot could trigger the appearance of cavitation bubbles in the wake of the moving laser. A similar geometry was recently investigated [34] in numerical simulations of a moving thin wire in superfluid 4 He, where vortex dipoles shedding occurred, and where cavitation bubbles formed in the wake of the moving wire, which were found to be responsible for large part of the dissipation accompanying the wire motion.

IV. SURFACE TENSION OF SELF-BOUND BOSE-BOSE MIXTURES
The appearance of self-bound droplets implies the existence of a surface energy, and a surface tension associated to it. Cikojevič et al. [13] have fitted their DMC energies for self-bound droplets of 39 K atoms in two different internal states to a liquid droplet expression [16] in order to determine the surface tension of the mixture where E v , E s and E c are volume, surface and curvature energies. The surface tension of the fluid is estimated as σ = E s /(4πr 2 0 ), where the bulk radius r 0 is related to the equilibrium density of the liquid n 0 as 4πr 3 0 n 0 /3 = 1, implicitly assuming that the radii of the density profiles for both species are sensibly the same. Within DFT, one may address the surface tension of the mixture avoiding the fit procedure to a series of calculated droplets. The surface tension σ -actually the grand potential per unit surface [23]-of a fluid planar free surface is determined along the saturation line of the liquid-vapor (or liquid-liquid) two-phase equilibrium. In the present case, the line reduces to the P = 0 point, as the mixture is at T = 0. If the z-axis is taken perpendicular to the free surface, one has where A is the free surface area.
To avoid the complication of imposing different boundary conditions for the density profiles at the opposite ends of the simulation cell, it is more convenient to use a "slab" geometry characterized by a uniform density in the (x,y) plane and two "liquid"-vacuum planar interfaces perpendicular to the z-axis. Here "liquid" means a self-bound mixture of species 1 and 2, whose densities n 1 , n 2 in the bulk region of the slab are determined by the equilibrium conditions discussed before for the uniform system case.
The energy density of the inhomogeneous system with densities n 1 (r), n 2 (r) is, from Eq. (1): As a case of study, we consider again a mixture of 23 Na and 87 Rb atoms. We will look for self-bound states (i.e. a 12 < −62 a 0 , as determined from the uniform system calculation in the previous Sec.) of a number of atoms N 1 , N 2 contained in a cell of sides (L x , L y , 2L). The size 2L of the cell along z is chosen in such a way to guarantee a wide enough region outside the slab, where the densities n 1 , n 2 are essentially zero. We have obtained the equilibrium densities n 1 (r) = |Ψ 1 (r)| 2 and n 2 (r) = |Ψ 2 (r)| 2 from the solution of the coupled EL Eqs. (2) in the slab geometry for different values of the interaction strength a 12 . Several total density profiles, n(z) = n 1 (z) + n 2 (z), for the calculated equilibrium configurations are shown in Fig. 5. Only one-half of the simulation cell containing the slab is shown for clarity. Note the very different shapes of the interface separating a bulk region in the left part of the figure and the vacuum region to the right, and that the more negative is a 12 , the narrower the interface.
We have calculated the surface tension using these The results are shown in Fig. 6 on a logarithmic scale, to underline the huge variation of σ -which spans almost four orders of magnitude-as the inter-species interaction strength is varied.

V. DROPLETS.
We describe here numerical calculations of isolate, spherical self-bound droplets made of N 1 atoms of 23 Na and N 2 atoms of 87 Rb. To this end, we have solved the coupled EL Eqs. (2) to obtain the densities, n 1 (r) and n 2 (r), for different values of the inter-particle interaction strength a 12 .
In our calculations we arbitrarily fix the radius R of the droplet to be computed. The values of N 1 and N 2 thus depend upon the chosen value of a 12 , and must be such that in the central part of the droplet, where the density profiles are sensibly constant, the associated densities n 1 , n 2 are those of the lowest energy per particle state of the mixture in the uniform system (see Fig. 2). In practice, we started our calculation with a density profile which reproduces, at the center of the droplet, the bulk equilibrium values (n b 1 , n b 2 ) predicted for the uniform system. Then N i is fixed so that N i = 4 πR 3 n b i /3. If we choose instead values of the densities which are far from the equilibrium ones, during the minimization process the excess atoms move towards the outer droplet surface and form a background of excess species. The experimental counterpart of such behavior is the evaporation which accompanies any excess species in the forming droplet.
We show in Fig. 7 the density profiles for one such self-bound droplet corresponding to R = 10 5 a 0 and a 12 = −75 a 0 . The total number of atoms contained in the droplet is N ∼ 3.5 × 10 6 . Figure 8 shows the droplet energy per atom E/(N 1 + N 2 ) for different a 12 values. The crossing with the E = 0 line marks the critical value of a 12 for the formation of a self-bound droplet.
As discussed previously, self-bound mixtures, when subject to a tensile stress, might enter the metastable negative pressure region and eventually reach the mechanical or diffusive instability line. A way to achieve this experimentally would be to compress a fairly big selfbound droplet by applying an external harmonic trap and then let it expand upon releasing the trap, thus bringing a large portion of the bulk of the expanding droplet into the negative pressure region of the phase diagram (e.g. Figs. 3 and 4).
We have simulated this process numerically, using the time-dependent version of the EL equations governing the dynamics of the mixed droplets. We consider here a self-bound droplet of radius R = 10 5 a 0 , made of 2.7 × 10 6 Na atoms and 3.6 × 10 6 Rb atoms, interacting with a 12 = −80 a 0 and subject to the compression exerted by an external isotropic harmonic potential with frequency ω = 2π × 400 Hz. As a result we observe, after the trap release, a sudden expansion of the compressed droplet accompanied by the formation of well separated radial shells and the fragmentation of the latters into smaller radially distributed clusters, each characterized by the same relative composition of the original droplet. After such initial expansion, the collection of fragments contracts, eventually leading to a radial oscillation of the whole structure. Several snapshots taken during such evolution are shown in Fig. 9. Note that symmetry breaking of the densities of the emerging fragments occurs in spite of the spherical symmetry of the initial configuration and contact atom-atom interactions, which is a common manifestation of modulation instability against azimuthal perturbations.
The details of the fragmentation process depend on the amount of the initial compression. A gentle squeezing of the droplet leads instead to the excitation of an intrinsic mode of the droplet in the form of a breathing oscillation, whose frequency depends upon the incompressibility of the mixture, Eq. (17).

VI. VORTICES
Finally, we briefly address vortical states in mixed self-bound droplets. In particular, we study the stationary states where a singly quantized vortex and a doubly-quantized vortex are nucleated in the center of the droplet. Vortical states in self-bound droplets have been studied recently in dipolar Bose droplets [35], and found to be unstable as a consequence of the very anisotropic nature of such droplets. The spherical mixed Bose droplets studied in our work, however, might sustain stable vortices. Swirling self-bound droplets made of Bose mixtures have been studied recently in Ref. [36], where it was found that self-trapped vortex "tori" with double vorticity are stable topological defects when the droplet exceeds a certain critical size.
We consider first a 23 Na-87 Rb droplet of radius R = 2 × 10 5 a 0 for a 12 = −75 a 0 , and imprint on each component a vortex with quantization number m by multiplying the pure droplet wave function by a phase factor e imφ , φ being the azimuthal angle, and evolve this initial configuration in imaginary time until a stable structure is found. A singly (doubly) quantized vortex is shown in the top left(right) panel of Fig. 10.
We have studied the dynamical stability of these two configurations by evolving them in real time after a quadrupolar perturbation has been added to the droplet by multiplying its wave function by the phase e iǫxy (note that this adds kinetic energy but not angular momentum to the system). We have chosen the small constant ǫ such that the applied perturbation increases the kinetic energy of the system by a few percent.
While the singly-quantized vortex is robust against quadrupolar perturbation, we have found that the doubly-quantized vortex rapidly decays into a pair of singly-quantized vortices. Such vortices are shown in the bottom panels of Fig. 10. As a result of the angular momentum associated with the two vortices, the vortex dimer rotates as a rigid body around the center of mass of the droplet.
The velocity field associated with the added quadrupo- lar phase, together with the angular momentum stored in the doubly-quantized vortex result in surface capillary waves, which distort the droplet surface and are responsible for the apparent rotation of the droplet as a whole [37], as shown in the bottom panels of Fig. 10. In this particular case, the vortex dimer appears to rotate with a frequency ω = 3 × 10 15 a.u. It is worth mentioning that, at variance with vortices and vortex arrays in expanding unbound condensates, these vortical configurations are stable and similar to those recently found in rotating 4 He droplets [38,39]. A more systematic study of vortex arrays in self-bound droplets and the merging of droplets hosting vortices is currently underway [40].
A qualitative similar behavior is observed for the 39 K-39 K mixed droplet under similar conditions. Fig. 11 shows the fate of a doubly-quantized vortex in a K-K droplet subject to the same perturbation described previously. The vortex decays in a close pair of singlyquantized vortices, as shown in the lower panels of Fig.  11 resembling a partially fused vortex dimer. Yet, these results seem to indicate that, under suitable conditions, the m = 2 vortex may represent a robust, stable topological defect in mixed Bose droplets, as indeed found in Ref. [36].

VII. SUMMARY AND OUTLOOK
We have investigated the zero temperature phase diagram of self-bound ultra dilute bosonic mixtures made of two different species within the DFT-LDA approach, providing an explicit expression for the Lee-Huang-Yang correction in this general case. We determined the general thermodynamic conditions which permit the formation of self-bound systems. To this end, we have obtained simple expressions to calculate the mechanical and diffusive spinodal lines. We have shown that, depending on the mixture, the thermodynamic condition that limits its stability may be either of the spinodal lines, and found, in agreement with previous work on equal species mix-  (iii) and (iv) shows the distorted-core structure resulting from the dynamical evolution of (ii). (iii) and (iv) are taken at different times during the evolution initiated from (ii), and shows the apparent rotation (counter clockwise) of the droplet/core system. Color map represents the total density (in a −3 0 ). Lengths are in units of a0.
tures, that the region of stability in the (n 1 , n 2 ) plane of compositions is extremely narrow. The appearance of self-bound droplets implies the existence of a surface tension, at variance with most of cold gases, which are metastable, unbound systems. We have thus calculated the surface tension of the mixture freesurface and the density profile of some selected droplets. In particular, our results show that the surface tension changes by orders of magnitude when the inter-species interaction changes by only a factor of two. The realization of stable, self-bound ultra dilute mixtures opens the possibility of studying phenomena that are otherwise restricted to high-densities, strongly correlated superfluids like liquid helium and liquid helium mixtures, such as e.g. cavitation [27][28][29], free-droplet merging [40], 1D or 3D droplet collisions and merging [41,42] or rotating free-droplets [36,38,39]. In a similar context, the possibility of creating self-bound Bose-Fermi droplets [43] opens the possibility to extend to these ultra dilute systems the study of, e.g., cavitation [25] and swirling properties of the prototypical Bose-Fermi quantum mixture, namely the 3 He-4 He fluid mixture [16].
Using Eq. (13)   Noticing that the terms containing ∂f /∂x and ∂ 2 f /∂x 2 in the second derivatives of E LHY cancel out, i.e. can be easily computed by transforming the above equation using the method of the Jacobians [23]. At constant temperature, From the definition of the chemical potentials and pressure at zero T , it is easy to show that ∂(µ 1 , P ) ∂(Y, P ) = n 2 n 2 where E is the energy density for the uniform system. Thus the diffusive spinodal line can be obtained by solving the equation 11) or, equivalently, Eq. (19). Notice that, from the computing viewpoint, both spinodals involve the same ingredients.