Quantum Bose-Fermi droplets

We study stability of a zero temperature mixture of attractively interacting degenerate bosons and spin-polarized fermions in the absence of confinement. We demonstrate that higher order corrections to the standard mean field energy of the system can lead to a formation of Bose-Fermi liquid droplets -- self-bound incompressible systems in a three-dimensional space. The stability analysis of the homogeneous system is supported by numerical simulations of finite systems by explicit inclusion of surface effects. Our results indicate that Bose-Fermi droplets can be realized experimentally.


Introduction
Self-bound systems are quite common in nature. They appear at different scales. Atomic nuclei, Helium droplets, or astronomical objects like white dwarfs or neutron stars are some of the prominent examples. Their stabilization mechanism is due to a subtle balance of attractive forces and repulsive interactions.
Yet another, not known so far, self-bound systems -extremely dilute quantum liquid droplets of ultracold atoms in a mixture of two Bose-Einstein condensates of different species 2 Uniform mixture The mean-field energy E 0 of uniform system in a volume V , having N B = n B V and N F = n F V bosons and fermions, respectively, can be written in the form: where n B and n F are atomic densities, and the consecutive terms, energy densities, correspond to: (1) a kinetic energy of fermions with ε F = 2 k 2 F /2m F = 5 κ k n 2/3 F /3 being the Fermi energy, and k F = (6π 2 n F ) 1/3 the Fermi wave number, (2) a boson-boson interaction energy, and finally (3) a boson-fermion contact interaction energy. For a convenience we introduced the following notation: κ k = (3/10) (6π 2 ) 2/3 2 /m F , g B = 4π 2 a B /m B , and g BF = 2π 2 a BF /µ, where a B (a BF ) is the scattering length corresponding to the boson-boson (boson-fermion) interactions and m B , m F , and µ = m B m F /(m B + m F ) are the bosonic, fermionic, and reduced masses, respectively.
We assume weak interaction limit, i.e. the gas parameters are small: n 1/3 B a B ≪ 1 and n 1/3 F a BF ≪ 1. At low Fermi density, the kinetic energy of fermions, being proportional to n 5/3 F , is the largest contribution to the mean-field energy. This term tends to spread the particles all the way to infinity. Similar is the effect of repulsive boson-boson interactions, g B > 0. However, a sufficiently strong attraction, g BF < 0, can suppress this expansion, but the equilibrium reached is unstable. Higher order terms must come into a play to ensure stability. Perturbative approach suggests the Lee-Huang-Yang term (LHY), the zero-point energy of Bogoliubov vacuum of the Bose system: B . However, this is not enough. Contribution to the mutual boson-fermion interaction resulting from the higher order term in the Bose-Fermi coupling turns out to be the most important. These effects were considered on a theoretical ground [28][29][30]. Results for Bose-Fermi system across a broad Feshbach resonance as given in [30] are applicable only when fermionic density is much larger than the density of bosons. As it will be shown later, it is not favorable for droplets formation. Contribution to the Bose-Fermi interaction energy obtained in a frame of the second order perturbation theory in [29], is more general than results of [28] based on the renormalized T-matrix expansion. This energy has the form where w = m B /m F and α = 2w(g B n B /ε F ) are the dimensionless parameters, C BF = (6π 2 ) 2/3 2 a 2 BF /2m F , and the function A(w, α) is given in a form of integral [29] A(w, α) = 2(1 + w) 3w where Θ(x) is the step theta-function. The arguments of A(w, α) function are: w = m B /m F and α = 2w(g B n B /ε F ) = 16π n B a 3 B /(6π 2 n F a 3 B ) 2/3 . The above formula, Eq. (3), for α ≪ 1, i.e. in the limit when the Fermi energy is much larger than the chemical potential of bosons (assuming the mass ratio, w, is of the order of one) coincides with the results of [28].
Summarizing the above discussion, in the regime where both gas parameters are small, a B n 1/3 B ≪ 1 and a BF n 1/3 F ≪ 1, we approximate the energy of the dilute uniform system by the following expression: Because the energy is an extensive quantity, it reaches extremal values for infinite system. Therefore some physical constrains should be introduced. For a trapped system they are set by number of atoms in every component. Here, we are looking for such configurations minimizing the energy, Eq. (5), which are stable in the absence of any external confinement. This means that neither the volume V , nor the number of bosons N B , nor fermions N F , are controlled. The important role plays a pressure, p(n B , n F ) = −dE/dV instead. Outside of the droplet it is equal to zero. The same must be inside since we look for a stable droplet and ignore here a surface tension. On the other hand, the energy related to a surface is proportional to its area ∝ V 2/3 , thus it represents only a fraction ∝ V −1/3 of the internal energy of a droplet. This fraction vanishes in the limit of infinite system, reducing the importance of the surface effects. This will be demonstrated in the next section. With all these limitations, vanishing pressure is a necessary condition for the mechanical stability of the system: where we introduced the chemical potentials µ B(F ) = dE/dN B(F ) = ∂ε/∂n B(F ) of both species. Numerical solutions of Eq. (6) are shown graphically in Fig. 1 for several values of a BF /a B and in the case of 133 Cs -6 Li (w = 22.095) mixture (contours for 41 K -40 K (w = 1.025) system look similar). Pressure vanishes on the closed contours forming a kind of loops in n B − n F plane. If |a BF |/a B is less than some critical value, |a BF |/a B < η 0 , then Eq. (6) has no solutions. The contours marked by solid lines support negative energy solutions. They shrink with increasing |a BF |/a B . All points on a single contour define densities of mechanically stable droplets for a given values of interaction parameters. Volume of the droplet is not specified, it is a scale parameter, if fixed, allows to determine the number of particles. What is the most important, the energy of the system (for a fixed volume) varies along the contour, and reaches the minimal value, if: Eq. (7) originates in a necessary condition for the extremum of the energy density ε(n B , n F ) constrained to the zero-pressure line. The energy minima are marked in Fig. 1 by dots. Only these particular spots define systems which are stable with respect to evaporation process. If the initial ratio of number of particles of the Bose and Fermi components is different than corresponding to the above mentioned solution, some particles, mostly the excess ones, are simply evaporated from the droplet. With our choice of the zero energy, the stable self-bound droplets should be characterized by negative energy. For some range of the parameter, η 0 < |a BF |/a B < η c energies are positive. These are metastable states, marked by broken lines in Fig. 1. Only if |a BF |/a B exceeds some  critical value, |a BF |/a B > η c , the energy of droplets becomes negative. These are the stable droplets. Values of η 0 and η c can be found numerically. In general, except of a small region of the parameter |a BF |/a B , the larger the attraction the smaller the equilibrium densities (see inset of Fig. 1).
In order to reach desired ratio of |a BF |/a B two approaches are possible. One is to increase the |a BF |, the second is to tune the Bose-Bose scattering length, a B , to small values. Both scenarios assume utilizing of appropriate Feshbach resonances. Simultaneously one should avoid large values of densities, as it would lead to a relatively large atom number decay, due to the three-body recombination. This is why a life-time of a droplet is typically of the order of 10 ms [7,8]. We did estimate the life-time of the Bose-Fermi droplet based on the measured losses of Cs atoms due to all possible three-body recombinations (see Fig. 4b in [27]). Extracting the three-body loss rate, averaged over the whole system, at the critical ratio |a BF |/a B = 2.8 we estimate the life-time of a droplet predicted here, to be about 50 ms.
In the first column of Tab. 1 we list values of the critical ratio η c supporting stable liquid droplets, while in the second and third columns we list corresponding densities of Bose and Fermi species. We present the results for three different mixtures of different mass ratio, 41 K-40 K, 87 Rb-40 K, and 133 Cs- 6 Li. For all these mixtures bosons are in vast majority. Therefore, fermions can be treated as impurity immersed in bosonic cloud bringing analogy to polaron. Boson-fermion attraction mediates an effective attraction between fermionic atoms. It prevents expansion of fermions due to quantum pressure. Similar mechanism, effective attraction between distinct electrons mediated by interaction with phonons is responsible for a formation of Cooper pairs. The question of fermionic superfluidity of Bose-Fermi droplets seems to be legitimate. The same interaction can also induce an effective attraction between bosons. For a large enough number of bosons this might result in a collapse of bosonic component. Then, fermions start to play an important role. They are able to counteract, due to quantum pressure, the collapse of bosonic cloud -an analogy with white dwarf and neutron stars becomes immediate. Hence, studying of the 'atomic white dwarfs' in the laboratory seems to be possible with Bose-Fermi droplets.

Finite system analysis
Now we address the properties of finite Bose-Fermi droplets including the surface effects. Evidently, additional energy terms related to the density gradients has to be considered. Therefore, we apply the local density approximation and add the kinetic energy  [32,33]. We neglect the contribution due to higher order gradient terms. The total energy of a finite Bose-Fermi droplet is then given by E[n B (r), n F (r)] = dr(ε + ε B k + ε F k,W ). The time evolution of the Bose-Fermi system can be conveniently treated within the quantum hydrodynamics [34]. For that both bosonic and fermionic clouds are viewed as fluids characterized by the density and the velocity fields. Since we assume that bosons occupy a single quantum state, their evolution is governed just by the Schrödinger-like equation of motion which includes the mean-field, the LHY, and boson-fermion interaction terms. A special care regarding the fermions should be taken. It has been already discovered many years ago that oscillations of electrons in a many-electron atom can be described by the hydrodynamic equations [35]. We follow this proposal and write the hydrodynamic equations for fermions in the following way where n F (r, t) and v F (r, t)) denote the density and velocity fields of fermionic component, respectively. T is the intrinsic kinetic energy of a fermionic gas and is calculated including the lowest order gradient correction [31][32][33] δT with ξ = 1/9. The convenient way to further treat Eqs. (8) is to put them in a form of the Schrödingerlike equation by using the inverse Madelung transformation [36][37][38]. This is just a mathematical transformation which introduces the single complex function instead of density and velocity fields used in a hydrodynamic description. Both treatments are equivalent provided the velocity field is irrotational (vanishing vorticity). After the inverse Madelung transformation is applied, the equations of motion describing the Bose-Fermi mixture are turned into coupled Schrödinger-like equations for a condensed Bose field ψ B (n B = |ψ B | 2 ) and a pseudowavefunction for fermions ψ F = √ n F exp (iφ) (n F = |ψ F | 2 and v F = ( /m F )∇φ) Here, ξ ′ = 1 − ξ = 8/9. The bosonic wave function and the fermionic pseudo-wave function are normalized as N B,F = dr |ψ B,F | 2 .
To find the ground state of the Bose-Fermi droplet we solve Eqs. (10) using the imaginary time propagation technique [39]. The resulting ground state densities for 133 Cs -6 Li mixture for three different numbers of bosons and fermions are shown in Fig. 2. For the smallest droplets the peak densities for both bosons and fermions are slightly higher than predicted by the analysis based on a uniform mixture. For bigger droplets the peak density approaches the uniform mixture solution as expected, because the surface effects become less important. The stability of the droplets is verified by the real time propagation of Eqs. (10). The conclusion is that the droplet reaches the state of the minimal energy by evaporating mostly surplus atoms. . For the ratio |a BF |/a B equal to 2.0, which is below the critical value, both the bosonic and fermionic clouds spread out. In the case when |a BF |/a B = 3.6, i.e. larger than the critical one, the breathing droplet is formed. The width of fermionic part of a droplet behaves similarly and is not shown.
Finally, we check whether Bose-Fermi droplets can be achieved by cooling the bosonic and fermionic gases in a harmonic trap and then just by opening the trap. For that we choose the 133 Cs -6 Li mixture and set a B = 250 a 0 , where a 0 is the Bohr radius and a BF /a B = −2.0 and −3.6. Initial numbers of atoms are: N B = 100000, and N F = 10000. They are confined in the spherically symmetric harmonic traps with frequencies ω B /(2π) = 200 Hz for bosons and ω F /(2π) = 940 Hz for fermions. The trap parameters are chosen to match a radius of a droplet to be formed. We find the ground state of such a Bose-Fermi mixture by solving Eqs. (10) in imaginary time. Next, the confinement is removed in 1ms (2730 m B a 2 B / ). We monitor the evolution of the system by a period of about 8 ms. The width of bosonic component, dr r |ψ B | 2 /N B , is shown in Fig. 3. Fermionic width looks similar. Vertical line indicates the moment of time when the confinement is completely switched off. The densities preserve the spherical symmetry during the evolution. Clearly, the system stays bound and a droplet is formed when the ratio |a BF |/a B is above the critical value equal to about 2.8 (see dashed line in Fig. 3). Contrary, for low enough ratio |a BF |/a B the bosonic and fermionic clouds just spread out (solid line in Fig. 3).

Conclusion
Presented here analysis of stability of a mixture of ultracold Bose-Fermi atoms indicate that stable liquid self-bound droplets can be spontaneously formed when interspecies attraction is appropriately tuned. Droplets are stabilized by the higher order term in the Bose-Fermi coupling. We predict the values of interaction strengths as well as atomic densities corresponding to droplets of three different mixtures, suitable for experimental realization, 41 K-40 K, 87 Rb-40 K, and 133 Cs-6 Li. The LHY quantum correction for bosons, which is positive, makes bosonic and fermionic densities low enough to decrease the impact of three-body losses. We demonstrate by time dependent calculations that a Bose-Fermi droplet should be achievable by preparing the mixture of bosonic and fermionic atoms in a trap and then by slowly removing the confinement.
Quantum Bose-Fermi droplets bring into a play the higher order term in Bose-Fermi coupling. The role of this term was not studied extensively in experiments so far. Moreover the effect of this 'correction' seems to be somewhat elusive as reported in [40]. We think that this fact should not discourage future experiments towards investigation of this higher order effects in Bose-Fermi systems. Liquid droplets such as studied here, seem to be the best systems to this end.
The ultradilute self-bound Bose-Fermi droplets are novel, not known before form of the matter organization. Their composition, involving not only bosonic, but also fermionic component, might bring into play a rich variety of physical phenomena related to polaron physics, Cooper pairing mediated by bosons, as well as to fermionic superfluidity. Stabilization mechanism involving Fermi pressure brings some analogies to astronomical objects like white dwarfs or neutron stars. Dynamics of droplets' collisions, their merging and collisions can simulate some astronomical processes as well.