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.
We thank both referees for their comments on the manuscript. Both referees agree that the topic we consider is timely and interesting. Unfortunately they express worries that our approach is not valid because of applying a mean field-approach and pseudo-wave-function formalism which, as they suspect, cannot account for a correct treatment of fermionic systems in a droplet.
The Referee 1 points two major defects of the manuscript:
1-The methods used are not valid to describe fermion dynamics.
2-The methods used are not valid to describe fermions in self-bound droplets.
Referee’s 2 criticism is based upon the same arguments. The referee says:
1. The use of a single complex function to represent all of the fermions in this system is not physically reasonable. There are significant issues both for the ground state and for dynamics.
The arguments on the necessity of accounting for a discrete spectrum of bound fermions are given in the two reports. In addition the referees point out to problems of the unphysical phase of the pseudo-wavefunction for fermions.
The Referee 2 requests the following changes:
1. Use a methodology that allows for the small number of discrete modes available to the fermions.
2. With that methodology show the bounds on the number of bosons and fermions for a droplet to exist and how that varies with interaction strength.
3. For reasonable interaction strengths, then give corresponding peak densities in SI units to allow comment on three-body loss.
In the response below we show that the above mentioned worries are not justified and we substantiate our results by additional arguments and results based on numerical calculations. Because both reports address similar issues we will write a single answer where we carefully address all points raised by the two referees.
Resubmitted version of our manuscript is substantially modified. In particular:
1. We added a new section no. 4, entitled: Finite system analysis – atomic-orbital approach, where we show results based on the Hartree-Fock formalism for fermions which, contrary to the mean field description, accounts for a discreet character of the fermionic spectrum
a. In the first part we present dynamical calculations for the ground state with a discreet single-particle basis for fermionic atoms (the size of the basis equals to the number of fermions). The base functions are modified adiabatically, starting from non-interacting trap system evolving towards interacting one without any trapping potential eventually. We consider both bound and unbound states. New figure 5 is added.
b. In the second part we present results of extensive numerical diagonalization of the system Hamiltonian in basis of 2600 oscillatory functions. We investigate two cases: unbound and self-bound systems. In both cases we show single-particle spectrum of fermionic subsystem. Figures 6 and 7 are added.
2. We added discussion of losses and new Fig.4 illustrating how number of bosons in droplets decreases due to the three-body collisions. We give physical values of parameters for the Cesium-Lithium mixture as used in the Chin Cheng’s group experiment .
We did not included in the resubmitted manuscript the request #2 of the second Referee: `With that methodology show the bounds on the number of bosons and fermions for a droplet to exist and how that varies with interaction strength.’
We have to admit that this is indeed a very interesting problem. We have shown that droplets exist even for as little as 35 fermions and 350 bosons. However systematic studies using the Hartree-Fock method require extensive, time consuming numerical work. In our opinion comprehensive answer to the referee’s question deserves a separate publication.
Below we present more elaborate discussion in reply to the referees’ concerns.
1. Comment on the referees’ criticism of using the pseudo-wavefunction formalism and on the mean-field description of fermionic systems.
We want to stress that the paper can be divided into two-parts: a semi-analytic (to some extend) approach where we specify conditions for equilibrium of a self-bound system with a free surface. This method, valid in the limit of infinite system, is general, and contrary to the approximate approach of D. Petrov does not base upon diagonalization of a quadratic form. It looks as if the referees are ignoring this part. In our opinion this part is very important. Results of this part were supported by numerical calculations based on the hydrodynamic approach.
In a stationary situation, both infinite system and finite system case, our approach is simply the standard Thomas-Fermi method with the Weizsacker correction included.
Evidently, the Thomas-Fermi model, introduced to describe electronic cloud and binding energies of multielectron atoms is far from the accuracy of sophisticated quantum-chemist approaches to multielectron systems. We do not claim that we are such accurate. But the TF approach in not totally wrong. It gives quite reasonable estimation (with 10% accuracy) of the binding energies of atoms, quite small systems though. Our predictions prove that the quantum fluctuations are able to stabilize the Bose-Fermi mixtures and can lead to formation of liquid droplets. We are convinced that in a case of large systems, the approach used by us is quite accurate, and we are not aware of any better than mean-field approach for systems having about 1000, 10000, or 100000 fermions or more, as we show in Fig. 2. No discreet treatment is possible for the system this big!
The static approach is generalized then to a dynamical situation by introducing a velocity field. The corresponding hydrodynamic equations are brought to the form of the Schoedinger-like equation for the pseudo-wavefunction which results from a kind of “complexification” of the density and the velocity “potential”. The transformation assumes that the velocity field is irrotational, or more precisely that it is defined on a simply connected support, and a velocity potential – a phase, can be used instead.
We explain the issue of complex pseudo-wavefunction for fermions in details in the comments to the reports and we cannot add anything substantial to this discussion without repeating the same arguments. So we believe that our approach is correct as long as Thomas-Fermi model is justified.
2. Elaborated discussion of the arguments against our approach bringing the controversy on continuous versus discreet approach for many fermion system.
The referees say, that self-bound Bose-Bose droplets have only few, or none, bound excited states. It’s true, all bosons can occupy the same state so one bound singe-particle state can support Bose-Bose-droplets, this is not enough in the case of fermionic component because of the Pauli principle. For fermions number of bound singe particle states gives an estimation of number of fermions which can be “trapped” in the droplet.
The referees, based on observation of  claim that no more than 10 one-particle states can be trapped in bosonic cloud, thus discreetness of the fermionic spectrum is crucial for bound systems. The referee 2 gives also example of Bose-Bose droplets, which have only few excited states.
The referee’s observation is correct. Number of bound states in effective potential formed by bosonic atoms must not be smaller than the number of fermions. We want to stress that we deal here with 3D situation, contrary to the effective 1D case of . Note that in 3D the energy states are highly degenerate, every angular momentum state L, is (2L+1)-fold degenerate. Our approach, leading to the mean-field energy of fermionic component is based on estimation of the number of bound states in a uniform potential. And this estimation is correct up to the leading order in the Fermi energy. For a spherically symmetric harmonic oscillator, the number of bound states of energies not larger them $m \hbar \omega$ grows as $m^3$. The TF model recovers the same scaling.
Because this issue seems to be controversial, what is expressed also in the comment of Mathew Davis, in the corrected version of the manuscript we included the subsection entitled “Finite system analysis – atomic orbital approach” devoted exclusively to justification of the mean field method. We support the hydrodynamic results by results obtained in the Hartree-Fock method accounting for discreetness of fermionic orbitals.
In the first part of this section we defined the Hartree-Fock formalism equivalent to the energy functional used by us. Then we find densities of droplet by adiabatic following of the ground state of 35 fermions and 350 bosons starting from decoupled system and gradually increasing the mutual coupling. We used the basis of 35 fermionic states, which were dynamically modified according to the Hartree-Fock equations coupled to the extended Gross-Pitaevskii equation. Resulting densities of droplets agree very well with those obtained by the pseudo-wavefunction formalism. In addition we showed that for too weak Bose-Fermi attraction the system is not bound and its radius spreads in time after releasing from the trap.
In the second part we used a huge basis of oscillatory wave-functions to find fermionic single-particle states in the effective bosonic potential for a small system of 35 fermions and 350 bosons. We want to stress that these are several-month-lasting calculations. Within this approach we show that number of states bound by bosonic cloud is exactly equal to the number of fermions.
We believe that this extensive and time-consuming calculations are convincing for the two referees, moreover they justify usage of the mean field approach as well as the pseudo-wave-function formalism.
Finally, as the referee 2 requested, we assumed realistic values of densities of Cesium-Lithium systems , for which we calculated the peak atomic densities, estimated the loss rate and showed the results of numerical simulations of droplet dynamics with losses included. We showed that the lifetime of droplet is sufficiently long there. To illustrate this analysis we included new figure.
We think that we gave answered to all the concerns of the both referees. We modified the manuscript to account for all their criticisms and we hope that our manuscript, in the present form, will be accepted for publication.
 B.J. DeSalvo et al. Phys. Rev. Lett. 119, 233401 (2017)
 M. Wentzel et al. Physica Scripta 93, 104004 (2018)