Self-bound Bose-Fermi liquids in lower dimensions

We study weakly interacting mixtures of ultracold atoms composed of bosonic and fermionic species in 2D and 1D. When interactions between particles are appropriately tuned, self-bound quantum liquids can be formed. We show that while formation of these droplets in 2D is due to the higher order correction terms contributing to the total energy and originating in quantum fluctuations, in 1D geometry the quantum fluctuations have a negligible role on formation of the self-bound systems. The leading mean-field interactions are then sufficient for droplet formation in 1D. We analyse stability conditions for 2D and 1D systems and predict values of equilibrium densities of droplets.


Introduction
Although the physical space is three-dimensional, lower dimensional spaces are quite often considered on a theoretical ground. This is because they are simpler for analysis, nonetheless able to capture the basic physics of the three-dimensional space. On the other hand the Hall effect, fractional statistics, Berezinskii-Kosterlitz-Thoules transition, Tonks-Girardeau gas are only some examples of physical phenomena attributed exclusively to 2D or 1D spaces.
It is possible to impose constraints on real 3D systems such that from the kinematic point of view they behave as low dimensional. In a case of ultracold quantum gases such tight confinement in one or two spacial directions can be obtained by a proper shaping of external fields forming a trap, making them highly anisotropic. If the excitation energy in the tightly confined direction(s) is the largest energy scale of the problem, the lower dimensional physics comes into a play.
Recently quantum liquids attracted a lot of attention of both experimentalists [1][2][3][4][5][6][7] and theorists [8][9][10][11][12]. Quantum self bound droplets can be formed in dipolar condensates or twocomponent Bose-Bose mixtures [8]. Recently conditions for formation of droplets in a Bose-Fermi mixture were found [13]. The life-time of droplets is limited by three-body collisions and at present experiments it is of the order of tens milliseconds. In lower dimensions it is expected that this life-time can be extended because of reduced phase-space available to colliding atoms. Low dimensional systems are therefore of the key interest [14][15][16][17]. In this work we study a possibility of formation of the self-bound low dimensional liquid droplets in a two component mixture of ultracold bosonic and fermionic atoms.

Uniform 2D mixture
The mean-field energy density of an uniform Bose-Fermi mixture is given by where n B and n F are bosonic and fermionic densities, respectively, and β = 2π 2 /m F . The first term in Eq. (1) corresponds to kinetic energy of fermions. The next two terms account for energy of intra-species (boson-boson) and inter-species (boson-fermion) interactions. The coupling constant g σσ ′ is related to the 2D scattering lengths a σσ ′ by g σσ ′ = 2π 2 µ σσ ′ 1 ln(ǫ/(a 2 σσ ′ κ 2 )) , where κ is a cut-off momentum, ǫ = 4 exp(−2γ), γ is Euler's constant, and µ σσ ′ = m σ m σ ′ /(m σ + m ′ σ ) is the reduced mass [14,18]. Here σ ∈ {B, F }. 2D scattering has few key differences in comparison to the 3D case -the 2D scattering length, a 2d , is always positive valued; scattering amplitude is energy dependent for all possible values of a BB and a BF . The dependence of energy is weak, so it is always possible to use a fixed value of coupling strength provided that high energy contributions are eliminated by a proper choice of a value of a cut-off momentum κ. The weakly-interacting regime requires a σσ ′ ≪ r 0 (weak repulsion) or a σσ ′ ≫ r 0 (weak attraction). Here r 0 is inter-particle separation. We show that it is possible to find such energy range of colliding particles, and thus to choose a suitable value for the cut-off momentum, for which the intra-species interaction is weakly repulsive and interspecies interaction is weakly attractive.
In the case of a two-dimensional Bose gas we deal with a quasicondensate. In such a case the lowest order correction to the mean field energy is given by [19] where ε LHY is the Lee-Huang-Yang (LHY) [20] energy density, Ω is a volume entity, ω k = (g BB n B 2 k 2 /m B + ( 2 k 2 /2m B ) 2 ) 1/2 . The contribution to the bosonic repulsive energy, E LHY = Ωε LHY , due to the zero-point energy can be obtained via standard Bogoliubov theory.
The above expression, representing interacting Bose gas, was obtained by treating ultraviolet divergences using the discreet lattice model as a regularization scheme. Therefore integration in Eq. (2) over momentum is carried up to the cut-off momentum κ. The 2D integration gives the following contribution: The higher-order correction in the Bose-Fermi coupling in a three-dimensional mixture were estimated via several theoretical approaches [21][22][23], the results from which agree within appropriate validity regimes. In order to evaluate the correction to the Bose-Fermi interaction due to coupling between density fluctuations in the two species, we begin with the relevant term given in Ref. [22]: Similarly as for Bose system we use here a discreet lattice model as the regularization scheme, and thus integration over momenta are carried up to the cut-off momentum determined by the lattice. Therefore the term used in [22] that counters ultraviolet divergence has been dropped. Here The Fermi momentum is given by k F = √ 4πn F . The coefficients u k and v k follow from the relations: where Herek c = κ/k F and α = 2w(g BB n B /ǫ F ), where ǫ F is the Fermi energy, and w = m B /m F . Let ε B = g BB n 2 B /2 + ε LHY and ε BF = g BF n B n F + δε BF . The total energy density of the system: formally depends on a cut-off momentum κ via coupling constants and to specify the system the value of κ = κ c has to be set. If chosen correctly the total energy density of the system is almost independent of the choice. Guided by analogous analysis of D. Petrov and G. Astrakharchik [14] who considered Bose-Bose droplet in two dimensions, we observe that the energy density reaches a maximal value when varying the cut-off at fixed values of bosonic and fermionic densities, ∂ε 2D /∂κ 2 | κ=κc = 0. κ c is our choice for the cut-off momentum. Evidently, it depends on the densities of the species. Two comments are rather obvious: i) close to the maximum there is no dependence on κ, ii) for correctly chosen cut-off, the beyond-mean field terms should be small corrections to the leading mean-field energy.
To better justify this choice for the Bose-Fermi system we consider a limit of infinitely weak attraction between the two components of the mixture, i.e. n B a 2 BF → ∞. In this limit bosonic and fermionic atoms become independent. Energy of bosonic subsystem in the limit of vanishing density, n B → 0, should approach the energy of a two-dimensional single component Bose gas, given by the famous Schick formula [24]. In Appendix A we show that this is indeed the case. This is a very important suggestion that our choice of the cut-off shall give correct value of ground state energy of Bose-Fermi system.
The cut-off depends on densities, but on the other hand, the densities of species forming a droplet are fixed by interactions, thus depend (weakly) on the cut-off. The problem of choosing the cut-off must be solved in a self-consistent way. Our strategy is the following. If the cut-off is set then the interactions' couplings are fixed. Stability conditions of the selfbound system give equilibrium densities of a droplet. Having those, the cut-off corresponding the maximum of energy, for fixed atomic densities, can be determined. The self-consistency can be easily reached. The procedure converges very fast and gives both, the optimal cut-off and the densities of Bose and Fermi components. We focus here at the weakly interacting regime assuming Bose-Fermi attraction, i.e. a 2 BF n B,F ≫ 1, and Bose-Bose repulsion, i.e. a 2 BB n B,F ≪ 1. The initial guess for the cutoff momentum is chosen from the bare minimum stability condition of the mixture, similarly as in [14,25]. Determinant of the second derivatives of mean-field energies must vanish at the stability edge: This leads us to the following equation: which can be solved for κ c , and the abovementioned iterative procedure can be initiated. The initial value of κ c happens to be a quite good guess.
The necessary condition for the appearance of a liquid droplet is vanishing pressure i.e., where the chemical potentials µ σ = ∂ε 2D /∂n σ . Due to the quadratic form of the dominant mean-field energy of the Bose-Fermi mixture in 2D, at the edge of the mean-field instability, i.e., for a small and negative δg = g BF + √ βg BB , we can use arguments of [8] and show that the equilibrium densities are close to the line n eq B /n eq F ≈ β/g BB , which turns out to be In the following we consider a weakly interacting mixture of 133 Cs and 6 Li, corresponding to the mass ratio w = 22.09, with varied a BF /a BB . Fig. 1(a) demonstrates a graphical representation of the numerical solution from Eq. (9) for a particular case of 133 Cs-6 Li with a BF /a BB = 10 4 . The zero-pressure line forms a closed contour on n B − n F plane. The droplet forms at equilibrium densities, for which the energy density of the system constrained by the zero-pressure requirement (Eq. (9)) attains the minimal value. This implies In Fig. 1(a), the energy minimum is marked by a dot. The variation of the energy density, ε 2D , along the zero-pressure line is shown in the inset of Fig. 1(a). The particular density ratio n B /n F , for which the system's total energy density is minimal, corresponds to the dot at the zero-pressure contour. In Fig. 1(b), we examine the flatness of various components of the total energy density functional around the chosen cut-off momentum, a B κ c = 0.034 obtained via the selfconsistent approach. One can see that close to κ c , once the higher-order corrections are accounted for along with the mean-field contributions, the energy density functionals show only a slight variance with the cut-off momentum. This suggests existence of a plateau around κ = κ c , where droplet equilibrium densities are almost independent on the choice of a particular cut-off.
In Fig. 2, we show the equilibrium density ratios of liquid droplets in 133 Cs-6 Li mixture with a BF /a BB . The solid line shows the approximate results governed by Eq. (10), and the circles represent the full solutions obtained via direct inspections of energy minima on the zero-pressure contour as explained in Fig. 1. The agreement between the complete and approximate solutions remains very good. In the inset of Fig. 2, we show the bosonic equilibrium density as a function of a BF /a BB . The fermionic part behaves similarly.

Finite 2D mixture
Now we investigate the density profile of the Bose-Fermi droplets with finite number of particles. We follow the similar treatment adopted in context of the 3D Bose-Fermi droplets [13]. In order to incorporate the surface effects, it's required to consider additional density gradient terms in system energy. This is done within local density approximation via inclusion of the kinetic energy component and the Weizsäcker correction [26] to the kinetic energy, for the fermionic component. For a 3D Fermi gas ξ = 1/9 [27,28], for a 2D case (Eq. (7)), however, parameter ξ starts to depend on the number of fermions [29]. The dependence is rather weak and for the number of 6 Li atoms considered by us here it is fairly to put ξ ≈ 0.04. The total energy of a finite Bose-Fermi droplet is then given by E[n B (r), n F (r)] = dr(ε 2D + ε B k + ε F k,W ). Both of the system's constituents, bosonic and fermionic clouds, can be regarded as fluids that can be treated within standard quantum hydrodynamics [30,31] by introducing the density and the velocity fields. Hydrodynamic equations can be reworked into a form of coupled set of Schrödinger-like equations via inverse Madelung transformation [32][33][34]. In order to find the density profiles of the bosonic and fermionic species, the coupled set of Schrödinger equation of motion are solved by employing imaginary time propagation technique [35]. See Appendix B and Appendix C for details. Figure 3 shows the ground state densities for 133 Cs -6 Li mixture for three different numbers of bosons and fermions. With increasing number of particles, as the droplets grow in size, the surface effects diminish as expected, and consequently, the peak densities approach the ones predicted by the analysis based on a uniform mixture in thermodynamic limit. The stability of the droplets is further verified by the real time propagation of the coupled set of Schrödinger equation of motion.

1D mixture
We now look at the possibilities for droplet formation in the case of attractive inter-and repulsive intraspecies interactions of a 1D Bose-Fermi mixture. The mean-field energy  density of 1D uniform mixture is given by The 1D coupling constant can be re-expressed in terms of the scattering length g 1d σσ ′ . The weakly interacting regime requires |a 1d σσ ′ |n σ ≫ 1. The 1D LHY correction for the Bose-Fermi mixture can be obtained from Eq. (2). Unlike the 2D case we no longer require to introduce a cut-off. The integration over the entire momentum range can be performed analytically. We find ε 1d LHY = −2g BB n B (m B g BB n B ) 1/2 /(3π ). The correction to the Bose-Fermi interaction in 1D, as it can be followed from Eq. (4), turns out to be δε 1d Finally, combining all the contributions, the total energy of the 1D system The investigation for the possibility of a stable droplet formation in 133 Cs-6 Li follows the same route discussed already in context of the 2D systems. The red circles in Fig. 4 show the equilibrium densities of 133 Cs as function of a 1d BF /|a 1d BB |. The equilibrium densities of corresponding 6 Li are shown by blue circles in the inset of Fig. 4. One can see that the smaller the value a 1d BF /|a 1d BB |, the greater the quantity |a 1d σσ ′ |n σ , which must be much larger than unity in order to ensure the applicability of our theory.
It can be numerically checked that I 1d is always positive, which makes the Bose-Fermi correction term attractive in nature, similar to the LHY term. Interestingly, we find that an effectively attractive mixture satisfying all condition of droplet formation exist within the mean-field description itself. In contrast to the 2D and 3D Bose-Fermi droplets and also to the Bose-Bose droplets, the tiny correction terms, although present, do not play any crucial role in formation of a stable droplet in 1D.
In Fig. 5 we show the bosonic and fermionic densities of a finite one-dimensional 'meanfield' droplet for a 1d BF /|a 1d BB | = 0.032 (g 1d BF /g 1d BB = −360) and different numbers of atoms. In both cases N B /N F is of the order of one hundred. Since the number of fermions considered is small, we use here an approach which treats fermionic atoms individually. We use the Hartree-Fock formalism. Assigning a single-particle orbital to each fermionic atom, φ F j (r), and assuming all the bosons to occupy the same state φ B (r), we solve the following set of time-dependent Hartree-Fock equations where j = 1, ..., N F and the total fermionic and bosonic densities are n F = NF j |φ F j | 2 and n B = N B |φ B | 2 , as well as its time-independent version (see Ref. [13] for details). In this way we find self-consistently the single-fermion orbitals, i.e. one-particle eigenstates of the effective potential created by both bosonic and fermionic components and hence densities. Clearly, as Fig. 5 suggests, going to larger samples the solution changes its character from the soliton-like (black solid and dashed lines, left corner) reported by us already some time ago (Ref. [36]) to the one with fully developed flat part typical for droplets (red solid and dashed lines). Noticeably, appearance of a train of Bose-Fermi solitons for appropriately tuned Bose-Fermi interaction, as claimed in Ref. [36], has been confirmed in recent experiment [37]. In 1D the fermions can be considered as impurities -much smaller in fraction in comparison to 2D and 3D [13] cases, yet essential for the droplet formation.

Discussion
In conclusion, in this work we show that higher-order quantum corrections may lead an ultracold weakly interacting gas of Bose-Fermi mixture towards a self-bound liquid state of matter. We find that in contrast to the 3D case [13], where droplets emerge above certain critical value of a BF /a BB , 2D ultradilute liquids are formed for arbitrary intra-species and inter-species interactions. Interestingly in 1D geometry, no quantum fluctuations are needed to form self-bound droplets.
Bose-Fermi systems are exotic due to the non-trivial interplay between dimension dependent scaling of kinetic energy of fermions, mean-field interactions and their higherorder contribution. The role of quantum pressure of fermions decreases while going from 3D to 1D geometry. Unlike in 3D systems, where fermionic kinetic energy scales like n 5/3 F , in 2D space the kinetic energy scales as n 2 F , and in 1D the scaling is n 3 F . This is why self-bound droplets in 2D can be formed for almost vanishing mean-field energy what is not the case for the 3D system. While recent experiments with Bose-Bose droplets have vindicated the role of beyond-mean-field LHY correction, the quantum Bose-Fermi droplets promises to be an ideal platform for probing another higher order quantum correction in Bose-Fermi interaction originating from the density fluctuations of bosonic and fermionic species.
The quadratic scaling of kinetic energy on fermionic density in 2D makes these systems closely analogous to Bose-Bose systems and hence the most promising candidates for investigating the liquid phases in highly anisotropic pancake shaped geometry with tight harmonic trap in the transverse direction. In future, it will be interesting to understand the nature of higher-order quantum correction in Bose-Fermi interaction, and the role they play in liquid formation, throughout the entire dimensional crossovers [16]. For a two-dimensional Bose gas alone, in order to reproduce low-energy scattering properties of the exact potential, one has [14,19] 1 where κ is a cut-off momentum. The formula is equivalent to the one used in the main text. Energy of a two-dimensional Bose gas is then given by After integrating over momenta, the energy density becomes and the first (mean-field) and the second (LHY correction) terms can be already identified in Eqs. (1) and (3), respectively, in the main text. By introducing g BB from (A.1) into Eq.
where x = n B a 2 BB and y = (κ a BB /2e −γ ) 2 are dimensionless parameters. In the limit of small values of parameter x, the energy density is given by Schick formula [24] where a new, density dependent, cut-off z = y/x ∝ κ 2 /n B has been introduced. Obviously, the second term in (A.6) vanishes in the limit of small x, whereas the first one approaches one. The low density limit for the Bose gas is then recovered. The energy becomes flat in a wider range of the cut-off parameter z when x gets smaller. However, in this paper we consider a two-dimensional Bose-Fermi mixture and description of boson-fermion interaction in two dimensions also requires introducing an appropriate cut-off momentum. We equalize both cut-offs and propose an iterative procedure to find its value. We start with, in principle, arbitrary chosen initial value of the cut-off but the value obtained by the minimum stability condition, Eq. (8), could be a good guess. Then the bosonic and fermionic densities are obtained based on Eqs. (9) and (11). Next, we utilize the condition ∂ε 2D /∂κ 2 = 0 to get a new value of κ. This is because the total energy of the system should only weakly depend on the cut-off. Therefore, the value of the cut-off should be close to the position of the wide maximum of the total energy, ε 2D (κ) for fixed densities. Then the next iterations are performed until the bosonic and fermionic densities, as well as position of maximum do not change. The value of the cut-off found in this way is κ c = 0.0340 and the densities are: n B = 3.97 × 10 −5 a −2 BB and n F = 4.10 × 10 −6 a −2 BB , see Fig. 1. We expect that even after adding fermions to bosons, when mutual attraction becomes infinitely small, n B a 2 BF → ∞, the energy of bosonic component will still satisfy the limit of low densities, the Schick formula. It is indeed the case, as shown in Fig. A1. The bosonic energy becomes more and more flat when the bosonic density gets smaller and the value of cut-off approaches the position corresponding to the maximum of ε B /ε Schick , which in turn approaches the value of one as expected (red dot in Fig. A1). This way we show that the cut-off we choose allows to recover the standard expression of energy of weakly interacting Bose gas in 2D.

Appendix B. Hydrodynamic equations
Hydrodynamic equations for a gas of neutral fermionic or bosonic atoms can be derived based on quantum kinetic equations for reduced density matrices [38][39][40]. From the whole hierarchy of equations, the one for the one-particle density matrix is of most practical interest where V ( r 1 − r 2 ) is the two-particle interaction term and V ext ( r, t) is the external potential. Eq. (B.1) involves the two-particle density matrix ρ 2 . In the limit r 1 → r 2 (= r), Eq. (B.1) results in the continuity equation ∂n( r, t) ∂t where the density and velocity fields are defined as follows: n( r, t) = lim r1→ r2 and χ( r 1 , r 2 , t) is the phase of the one-particle density matrix Eq. (B.1) can be rewritten by introducing the center-of-mass, r = ( r 1 + r 2 )/2, and the relative position, s = r 1 − r 2 , coordinates. Then, by taking the derivative of Eq. (B.1) with respect to the coordinate s, the hydrodynamic Euler-type equation of motion is obtained in the limit s → 0 where the kinetic-energy stress tensor, T, is given by ∂ 2 σ( r, s, t) ∂s k ∂s l (B.5) and the force, F ( r, t), resulting from interactions between atoms is The kinetic-energy stress tensor for fermions can be easily calculated within the Thomas-Fermi approximation. The one-particle Wigner function within the Thomas-Fermi approximation in a three-dimensional space is given by where η() is the unit step function. The one-particle density matrix is calculated according to ρ 1 ( r, s) = d 3 p (2π ) 3 w( r, p)e i p s/ , which implies that where p F ( r) = (6π 2 n F ( r)) 1/3 is the local Fermi momentum. Then, from Eq. (B.5) one gets the kinetic-energy stress tensor T kl = [ 2 /(30π 2 m F )](6π 2 n F ) 5/3 δ kl . Assuming the case of non-interacting fermions (spin-polarized atoms at low temperature), the Eqs. (B.2) and (B.4) become a closed set of hydrodynamic equations There is, of course, a space for improvements here. Commonly, gradient corrections are locally added to the Thomas-Fermi kinetic energy density (diagonal part of the kineticenergy stress tensor), the simplest one being the Weizsäcker correction in the form of T W = ξ( 2 /8m F )( ∇n F ) 2 /n F , with ξ = 1/9 [26][27][28]. Then the last term in the second equation of Eqs. (B.6) is modified by adding the derivative δT W /δn F .
We further assume that the fermionic flow is irrotational, i.e. the vorticity vanishes, ∇ × v F = 0. By using one of the vector identities, The set of equations (B.6) was first used by Ball et al. [31] to study the oscillations of electrons in a many-electron atom induced by ultraviolet and soft x-ray photons. For a gas of neutral bosonic atoms, in the mean-field approximation when the system's many-body wave function is ϕ( r 1 ) ϕ( r 2 ) · ... · ϕ( r N ), the one-particle density matrix is given by The kinetic-energy stress tensor possesses now off-diagonal elements and leads to the socalled quantum pressure term in the equation of motion. Assuming only contact interactions between bosonic atoms (with a being the scattering length), the Eqs. (B.2) and (B.4) become the set of the following equations (B.8) Eqs. (B.8) are, in fact, the hydrodynamic representation of the Gross-Pitaevskii equation. They can be improved with respect to the boson-boson interactions by considering the beyond mean-field correction -the Lee-Huang-Yang correction. The mixture of bosonic and fermionic atoms is then described by equations (B.7) and (B.8) modified by a term representing the interactions between bosons and fermions. This term, of the form of δE BF /δn F for fermions and δE BF /δn B for bosons, where E BF is the total Bose-Fermi interaction energy, acts as an additional potential in the Euler-like equations of motion. What is crucial for the existence of Bose-Fermi droplets, E BF energy includes the beyond mean-field correction. Although the above derivation of hydrodynamics of the Bose-Fermi mixture was performed explicitly in a three-dimensional space, it can be repeated in two-and one-dimensional cases as well. The changes which should be done are related to the local fermionic kinetic energy (for instance, in a two-dimensional space T kl = ( 2 π/m F ) n 2 F δ kl and T W = ξ( 2 /8m F )( ∇n F ) 2 /n F , with ξ = 0.04) and beyond mean-field corrections for Bose-Bose and Bose-Fermi interaction terms.

Appendix C. Inverse Madelung transformation
As shown in the previous Section, the Bose-Fermi mixture can be treated by using equations like (B.7) and (B.8). The very convenient way to further tackle Eqs. (B.7) and (B.8) is to put them in a form of the Schrödinger-like equations by using the inverse Madelung transformation [30]. 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). This assumption is obviously fullfilled for bosonic component since we consider that bosons populate a single quantum state. It is also true for fermions in a stationary case (zero velocity field) and in the case of dynamics studied by us while the trap confining mixture is adiabatically removed (in this case we a posteriori check that vorticity is zero).
After the inverse Madelung transformation is applied: ∂ψ F ( r, t) = n F ( r, t)e −iχ( r,t) (C.1) where v( r, t) = ∇χ( r, t), Eqs. (B.7) and (B.8) are turned into where now the main objects to be solved are the bosonic wave function, ψ B , and the fermionic pseudo-function, ψ F . The coefficients are: ξ ′ = 1 − ξ = 0.96, A = . The bosonic wave function and the fermionic pseudo-wave function are normalized as N B,F = dr |ψ B,F | 2 . The integral I c and all the parameters ω, α, and β and those appearing in coefficients A, B, and C are defined in the main text.
Note that Eqs. (C.2) describe two-dimensional Bose-Fermi mixture with beyond meanfield Lee-Huang-Yang correction for boson-boson interactions and analogous correction for the interaction between bosons and fermions.