Self-consistent Description of Bose-Bose Droplets: Modified Gapless Hartree-Fock-Bogoliubov Method

We define a formalism of a self-consistent description of the ground state of a weakly interacting Bose system, accounting for higher order terms in expansion of energy in the diluteness parameter. The approach is designed to be applied to a Bose-Bose mixture in a regime of weak collapse where quantum fluctuations lead to stabilization of the system and formation of quantum liquid droplets. The approach is based on the Generalized Gross -- Pitaevskii equation accounting for quantum depletion and anomalous density terms. The equation is self-consistently coupled to modified Bogoliubov equations. The modification we introduce resolves the longstanding issue of missing phonon-branch excitations when higher order terms are included. Our method ensures a gapless phononic low-energy excitation spectrum, crucial to correctly account for quantum fluctuations. We pay particular attention to the case of droplets harmonically confined in some directions. The method allows to determine the Lee-Huang-Yang-type contribution to the chemical potential of inhomogeneous droplets when the local density approximation fails.


I. INTRODUCTION
Standard theoretical description of atomic ultracold weakly interacting Bose gas is based on the Gross-Pitaevskii (GP) equation [1,2] for a Bose-Einstein condensate wavefunction, and on Bogoliubov equations [3] giving low energy excitation spectrum. This simple approach accounts for interactions in the lowest order perturbation in a small parameter, √ na 3 ≪ 1, where n is the density of atoms and a is the s-wave scattering length characterizing interaction potential. The formalism predicts a gapless phononic excitation spectrum in agreement with the Hugenholtz -Pines theorem [4]. This is a very important feature of any Bose condensate, heralding its superfluid nature according to the Landau criterion [5].
This simple description occurs to be insufficient in some situations. Recent discoveries of quantum droplets and dipolar supersolids formed by ultracold, dilute Bosonic gases [6][7][8][9][10] proved that stability of these exotic states of matter critically depends on higher order terms in expansion of energy in the diluteness parameter. Stabilizing role of these terms was pointed out by D. Petrov [11] predicting possibility of formation of quantum droplets. These theoretical and experimental discoveries triggered revival of theoretical interest in studies of ground state energy and excitation of weakly interacting Bose systems.
These studies have a long, lasting for many decades, history [12][13][14][15][16][17][18][19][20][21][22]. In [23] the two-leading terms in the expansion of diluteness parameters of the energy of hardsphere Bosons were found. This important result showed that energy of a dilute system depends only on the swave scattering length a, not on the details of interaction potentials. Soon later, S. T. Beliaev [24] formulated approach accounting for higher order terms in perturbation in √ na 3 . This achievement paved a way towards exploration of such systems like liquid 4 He where interactions are not weak at all. It is commonly expected that effects of higher order terms can be significant only in strongly interacting systems.
In the pioneering paper of D. Petrov [11], another possibility is noticed. Namely, if one arranges such conditions that dominant mean field interaction energy vanishes, higher order terms dominate. This might happen in a Bose-Bose mixture with repulsive intra-species interaction and sufficiently strong inter-species attraction equalizing the repulsion. Mean-field description predicts a collapse of a wave-function but the first order correction to the interaction energy, the repulsive Lee-Huang-Yang (LHY) energy term [23], might stabilize the system at finite density.
Self-bound quantum droplets are formed if interactions are in a certain range [11]. The system remains dilute and weakly interacting. Its density is of the order of 10 15 cm −3 . Moreover, droplets, if sufficiently large, have a constant bulk density independent on number of particles. This property is commonly attributed to a liquid state, although droplet's density is by orders of magnitude smaller than other liquids.
Calculations of LHY-like contributions to the energy, for systems in various configurations became at the center of interests of many research groups. The LHY energy is found for homogeneous Bose-Bose mixtures in 3D [11], as well as 1D or 2D geometry [7] or at dimensional crossover [34,35]. Similar studies give LHY energy of dipolar gases in 3D [36][37][38] as well as in lower dimensions, at dimensional crossover [35,[39][40][41], or in a lattice potential [42].
LHY energy, once obtained, becomes a crucial ingredient of the extended GP equation as suggested in the pioneering paper of D. Petrov [11], where theoretical description of droplets is formulated in two steps. In the first step, the LHY energy is found for a stable uniform system. In the second step, this energy is added to the energy functional and extended GP equation is obtained. The equation is applied then in a region of interactions where the mean-field description predicts instability, manifested by imaginary LHY energy and sound velocity of the soft-mode.
Evidently the approach has some drawbacks -the LHY energy is imaginary and, in addition, corresponds to a homogeneous system. Omitting a small imaginary contribution to the LHY energy might be excused close to the droplet formation point. The LHY energy term is constructed using two kinds of excitations -so called hard and soft modes. The imaginary part of the LHY energy comes from the soft mode contribution. It is very reasonable to assume that the soft-mode excitation energy of a nonuniform stable droplet, if determined selfconsistently, is real and very small as compared to the hard mode energy. The efforts to cure imaginary excitation energies focused on accounting for another, in addition to LHY term, contributions to the energy functional. In [43,44] a pairing energy was included while in [45], next to the LHY energy term was added. These terms could be essential at stronger interactions, but at the nearest proximity to the collapse, the LHY term alone is sufficient.
Using the LHY energy of an uniform system to describe a droplet has limited applicability. This approach is valid only if the local density approximation is justified, i.e. if density does not change significantly on a distance of a healing length. Droplets strongly squeezed in one or two directions do not meet this criterion. Their density profile changes abruptly in the confined direction. In such a case one should use a self-consistent approach accounting for quantum fluctuations of a nonuniform system. In particular, a stable droplet solution of generalized GP equation should be used as a source of fluctuation field in the Bogoliubov-de Gennes equations. Simultaneously the same fluctuations should enter GGP equations to give a stable droplet profile.
Here we develop the self-consistent approach capable to handle the issue -the method that permits to describe quantum droplets not using local density approximation. The only paper to date, considering non-homogeneous system beyond the local density approximation, assumes a priori the given density profile and cannot be generalized to account for a self-consistent density, [35]. Still there does not exist a method that enables to calculate the LHY energy using a density profile modified by that LHY energy. The need of such treatment arises in the systems with sufficiently strong external trapping. As a bonus, the approach should give real sound velocities.
The method defined in this paper relies on Generalized Gross-Pitaevskii equation (GGP) describing quantumdroplet wave-function coupled to Bogoliubov equations. The GGP equation accounts for quantum depletion and anomalous density which are to be obtained from solutions of Bogoliubov equations self-consistently. Unfortunately, accounting for quantum fluctuations both in GGP as well as in Bogoliubov equations meets a serious problem -an unphysical gap in energy spectrum appears, so phonon-excitation branch is missing [13,19,46]. Phonon excitations are crucial for low energy properties of the system.
The issue, in a case of single component Bose gas, is solved in [47] where Gapless Hartree-Fock-Bogoliubov approach based on the GGP equation is suggested. The approach is design to describe frequencies and damping of excitations at temperatures close to critical one.
Modified Gapless Hartree-Fock-Bogoliubov (MGHFB) method introduced here, tackles the gap problem in a different way than in [47], because for realistic parameters, the atoms in droplet phase are weakly interacting and a role of high energy contribution to quantum depletion and anomalous density is negligible. Situation is distinct from that at larger temperatures considered by S. A. Morgan [47]. We show that phononic spectrum can be recovered provided that Bogoliubov equations are modified by introducing a relatively small (controlled by the diluteness parameter) shift, δµ/µ ∼ √ na 3 , of the value of chemical potential. Instead of a value given by mean-field wavefunction solution of the GGP equation, the eigenenergy of the zero-mode, the one which restores broken U(1) gauge symmetry, is used, [48]. This modification makes a colossal change -it cures pathological behavior of low-energy modes.
Our self consistent approach accounts for terms responsible for formation of quantum droplets -quantum depletion and renormalized anomalous density. The first one describes density of noncondensed particles in ground state of a system, while the second measures correlations of pairs. Anomalous density is diverging with increasing cutoff in momentum space. We introduce a physical procedure of regularization which cancels the divergent part. The regularization procedure as well as its numerical implementation is yet another important result of our paper.
This paper describes the first part of our work devoted to self-consistent description of quantum droplets. Here we formulate MGHFB method -first, for a homogeneous case, next we focus on inhomogeneous system, and finally we discuss a two component inhomogeneous Bose-Bose mixture close to the transition to a droplet state.
Results of the present paper are applied in [49] where we use the MGHFB method to find contribution to the chemical potential originating from quantum fluctuations for a Bose-Bose droplet squeezed in one spatial dimension by a harmonic potential and not confined in two remaining dimensions. This geometrical settings allows to investigate quantum fluctuations in the entire range of geometric settings: from 3D to 2D arrangements. Local density approximation cannot be used for a tight direction in such a case.
This paper is organized as follows. In Sec. II we investigate a single component Bose gas to introduce main lines of our approach, first for a homogeneous system and next we generalize the method to introduce a beyond local density approach. We derive GGP equation introducing quantum depletion and renormalized anomalous density. We formulate MGHFB approach and discuss low energy excitations. Here we also discuss a numerical procedure to calculate renormalized anomalous-density. In Sec. III we use the method derived for a single component Bose gas to specify MGHFB method for the case of a Bose-Bose mixture. In Appendix A we prove that low energy excitations are phonon-like. In Appendix B, we show details of semiclassical calculations.

II. SINGLE COMPONENT BOSE GAS
To systematically introduce accounting for quantum fluctuations method of a self-consistent description of a Bose system, we investigate a single component case first.
Calculations are less involved than in a two-component mixture case. In a derivation we use effective two particle interaction potential U (r). The potential serves as a tool to formulate the method. It is chosen to give a value of the s-wave scattering length of the true interparticle potential.
For the reasons that shall become clear below, we choose the potential U (r) to be a positive bell-shapelike, extending over a distance of the order of σ, which is much larger than the s-wave scattering length, a, i.e. σ ≫ a. The system is characterized also by two other length scales. The first one is the so called healing length given by ξ = √ mang where n is the density of the gas, and g = 4π 2 a/m a . The second one, d, describes a characteristic length-scale of density variations. The density is not uniform if the gas is confined by some external potential. Below we assume that σ is the smallest of the two length scales, i.e. σ ≪ ξ, d.
The Hamiltonian of the system: generates dynamics of the bosonic fieldψ(r, t): where H 0 (r) = − 2 2ma △ + V (r). In the following we consider the ground state of the system and focus on a dilute gas limit i.e. na 3 ≪ 1. It is known that under such conditions most of particles populate a single mode of the system -a Bose-Einstein condensate is present. We use a standard procedure and divide the field operator into a mean field, ψ(r)e −iµt = ψ (r, t) , corresponding to a condensate, and fluctuations,δ(r, t)e −iµt , describing out-ofcondensate component:ψ(r, t) = ψ(r) +δ(r, t) e −iµt . In what follows we assume ψ to be a real function. Note that ∂ t ψ(r) = 0. Inserting the above into Eq. (2) we have: which after taking the ground-state mean value gives: In the above we introduced anomalous density: and quantum depletion: In writing Eq. (4) we notice that δ = 0, and we neglect the cubic term, ∼ δ †δδ . Moreover, to obtain a formalism depending only on the s-wave scattering length, we introduced some further approximations: whereũ(k) are Fourier components of the interaction potential,ũ(k) = dr e −ikr U (r), in particular,ũ(0) = dr U (r). All these approximations are justified since it is assumed that ψ 2 and quantum depletion δn(r) vary on the length scales, ∼ d, much larger than σ. Similarly, off-diagonal quantum depletion δn(r, r ′ ) varies on distances of the order of a healing length, ξ ≫ σ. Only anomalous density m(r ′ , r) depends on high momenta excitation, thus changes on a smaller scale. We will discuss this issue in the following part.
In the ground state, these quantities do not depend on time. Subtracting Eq. (4) from Eq. (3) and preserving only terms linear inδ, we get: Fluctuationsδ(r, t) can be expanded into eigenmodes u ν (r) and v ν (r): whereα ν ,α † ν are annihilation and creation operators of bosonic quasiparticles, bosonic commutation relation imply that: Finally Bogoliubov equations are obtained: Note, that terms m(r ′ , r) and δn(r ′ , r) which appear in GGP equation, Eq. (4), are not present in Bogoliubov equations, Eq. (15). This inconsistency may lead to a gap in the excitation spectrum or imaginary values of quasiparticle energies even in the ground state of the system. Both these phenomena are nonphysical artifacts of approximations made. We cure this problem by modification of Bogoliubov equations Eq. (15). The chemical potential µ (the one entering the GGP equation) is replaced by another (but close) value µ 0 : The analysis of the above approximation is discussed in details in Appendix A. This replacement Eq.(16) is a crucial point of our approach. It is done ad hoc to compensate for the inconsistent approximations. However, this small, as it will be shown later, modification of µ in Bogoliubov equations has important consequences. In case of systems which are not bound (at least in one direction) it ensures a correct, gapless, phonon-like excitation spectrum (see Appendix A for proof of that property). We postulate that the consistent formulation of MGHFB method should be based on the following, modified Bogoliubov equations (15): where µ 0 is the lowest energy solution of the following eigenproblem: where u 0 (r) and v 0 (r) = −u 0 (r) are zero energy eigenvectors recovering broken U (1) gauge symmetry, [48]. For completeness of the defined here approach, we remind that ψ(r) is a solution of the GGP equation, Eq. (4): The equations formulated above involve explicitly a particular form of the potential U (r), although their solutions depend on the low energy component of U (r) only. Anyway, the approach is impractical. Below we reformulate the method starting from a homogeneous case, to pinpoint some possible simplifications and further approximations which will grant solutions in terms of physical quantity, namely the s-wave scattering length. These reformulation allows to generalize the formalism to tackle inhomogeneous system.

A. Homogeneous system
In a homogeneous case we can expand quasiparticles' eigenfunctions, u ν (r) and v ν (r), into Fourier series: so solutions of Bogoliubov equations have the form: In the above n 0 = ψ 2 (r) is the condensate density, are energies of free particles, and coefficients A k and B k are defined as: The Fourier transform of the interaction potentialũ(k) = dr e −ikr U (r) is assumed to be non-negative. Inserting the above into definitions of m, and δn, and noticing that averaging should be performed with quasiparticle vacuum, we arrive at the expression giving δn(r ′ , r): as well as the pair correlation function m(r ′ , r): Trying to connect the above observables with the scattering length we use Born expansion of the T-matrix which takes the form: where g = 4π 2 ma a is the T-matrix and a is the s-wave scattering length.
Above, we assumed that width of potential U (r) is equal to σ. This implies that width of Fourier transform u(k) is ∼ 1/σ. This fact enables us to estimate the integral appearing above as: As a result from Eqs. (27) and (28) we obtain: According to our assumption, the range of the potential is much larger than the scattering length a ≪ σ, therefore the above can be estimated as: The above estimation gives The same arguments, based on the shape of interaction potentials, give thatũ(k) ≃ũ(0) for k ≪ 1/σ. This implies that the term E k +n0ũ(k) is proportional to ∼ 1/k 4 for k ≫ 1/ξ. As ξ ≫ σ thus the integral in Eq. (24) converges fast enough 'not to feel' the shape of the potentialũ(k) but only the valueũ(0) ≃ g. As a result δn(r, r ′ ) depends only on g and we have In case of dilute system, as we deal with, we have na 3 ≪ 1 thus δn ≪ n 0 . Additionally one may find that δn(r, r ′ ) changes on a length scale equal to ξ which again depends only on a (and not on the particular choice of u(k)). Eventually δn(r, r ′ ) depends only on the s-wave scattering length.
Similarly, one can find that n0ũ(k) 2ε k , present in Eq. (26), is proportional toũ(k)/k 2 for k ≫ 1/ξ. As a result the anomalous average, m(r, r), strongly depends on the shape of the potential U . However, the true quantity of interest is not the anomalous density but the chemical potential. From Eq. (4) it follows that in case of homogeneous system it takes the form: Substituting m(r ′ , r) from Eq. (26) into Eq. (32) we obtain: In the derivation above we approximateũ(0) ≃ g in the expressionũ(0)δn ≃ gδn. This approximation cannot be used however, when the term n 0ũ (0) is considered, because it involves large quantity, n 0 ≫ δn. In attempt to expressũ(0)n 0 by physical parameters, we wave to include higher order contribution to the scattering matrix, Eq. (27). As a result Eq. (33) takes the form: The expression entering integral above is very similar to the term discussed in case of quantum depletion. The integrated function of momenta (1/E k −1/ε k ) converges on the scale of 1/ξ and therefore the integral in Eq.(34) gives a finite value depending only on g. This well-behaved integral is often referred to as the renormalized anomalous density m R (up to the multiplicative constant g): The chemical potential, accounting for contribution originating in quantum fluctuations is therefore: Inserting into the above δn, m R from Eqs. (31) and (35) and using n = n 0 + δn we finally arrive at: where we additionally approximated n 0 a 3 ≃ na 3 . As we see the chemical potential depends only on the s-wave scattering length and not on any other detail of the U (r) potential. Eq. (37) gives chemical potential of the homogeneous system including contributions from quantum fluctuations ∆µ = gn 32 It follows from the above discussion that the anomalous density can be estimated to be of the order |m| ≃ n 0 a σ . As a ≪ σ we find that |m| ≪ n 0 . Moreover, δn ≪ n 0 . It means that in our case the condensate density ψ 2 (r) is much larger than m(r, r ′ ) and δn(r, r ′ ).

B. Inhomogeneous system
In case of nonuniform system we are not able to give analytic expressions linking directly the renormalized anomalous density and quantum depletion to the s-wave scattering length of the interparticle potential U (r). Our goal instead, is to formulate the approach in a way which allows to find the quantities in question numerically. In particular we will give a prescription of calculating the renormalized anomalous density avoiding all unphysical singularities. In what follow we assume that crucial from the physical point of view properties of δn and m remain valid also in the inhomogeneous case. This means that inhomogeneity is not too strong.
The interpartical potential U (r) enters explicitly the GGP equation, Eq. (4), through the term C : We proceed similarly as in the homogeneous case. We substituteũ(0) by an appropriate order of the T-matrix expansion, Eq. (27). We use the fist term of the series u(0) = g in terms involving small parameter δn, while expansion up to the second order is adapted in the dominant term proportional to a condensate density: In the above both quantities ψ 2 (r) 2E k and dr ′ U (r ′ − r)m(r ′ , r) depend on σ (i.e. on high energy modes). In full analogy with homogeneous case we define a renormalized anomalous density: H 0 + g ψ 2 (r) + 2δn(r) + m R (r) ψ(r) = µψ(r), (41) which should be supplemented by a normalization condition: where N denotes number of atoms. In the following part of this section we study GGP equation, Eq. (41), in depth. Our goal is to formulate an approach to effectively determine δn and m R . As discussed above, in homogeneous system the modes contributing to these quantities depend on low energy physics only, i.e.ũ(0) ≃ g. In Appendix B we show that the same takes place in inhomogeneous case. Here we shortly describe the main lines of the proof.
We divide the space of Bogoliubov modes into two parts: low and high energy one. We choose E c to be the energy dividing both sectors. This energy is chosen to be low enough for the low-energy modes to depend effectively only onũ(0) ≃ g. In this energy sector the Bogoliubov equations (17) can be approximated to read: Similarly, the zero-mode equation, Eq. (18), now is: These equations are to be solved numerically. On the other hand, in the high energy sector we solve analytically the Bogoliubov equations (17) using semiclassical approximation. Accordingly, both quantum depletion and anomalous density are divided into low and high energy components, i.e. n(r) = n L (r) + n H (r) and m(r, r ′ ) = m L (r, r ′ ) + m H (r, r ′ ), where low energy quantum depletion is: and analogously, anomalous density: The low energy contributions n L (r) and m L (r, r ′ ) can be calculated numerically while the high energy components can be obtained using semicalssical approximation. High energy part of anomalous density, m R H (r), requires renormalization, silmilarly as in a homogeneous case. According to Eq. (B24) : where k c (r) is the momentum dividing low and high energy sectors and is given by the equation ε(k c (r), r) = E c , and We define A(k, r) = 2 k 2 2ma + V (r) − µ 0 + 2gψ 2 (r) and B(k, r) = gψ 2 (r). One can clearly see that momenta k giving dominant contribution to the integral in Eq. (48), are of the order of k ∼ m a gψ 2 (r)/ 2 = 1/ξ(r). They are much smaller than 1/σ and the entire contribution depends only on low momenta components of the interaction potential,ũ(0). This observation has been already used to derive above formula where we explicitly wrote g ≃ũ(0) instead ofũ(k) (see Eqs. (B23) and (B24)). As a result, the renormalized anomalous density m R = m L + m H reads: Similarly the high-momenta contribution to the quantum depletion can be brought to the form (see Appendix B for details): The above discussion not only gives a prescription how to obtain the regularized anomalous density, m R , but it is also a direct proof that a value of m R (r) does not depend on a particular shape of the interaction potential and only on the s-wave scattering length, i.e. onũ(0) ≃ g.
In a view of this fact we might think about alternative approach to calculate the regularized anomalous density m R , directly from solutions of Bogoliubov equations (44) depending solely on g. In this approach we don't use the decomposition into low and high energy contributions, nor use the high energy semiclassical formulae given by Eqs. (50) and (51). In Appendix B we show that such a direct approach gives a closed expression (see Eq. (B27)) where m(R, ∆r) = m(R + ∆r 2 , R − ∆r 2 ), R = r+r ′ 2 and ∆r = r − r ′ , and lim ∆r→0 R = r. The above equation is equivalent to which is identical to the formula resulting from approximating the interaction potential by the regularized Fermi-Huang zero-range pseudopotential, U (r) = δ(r) ∂ ∂r r, [23,50]. In this approach we simply have In the above we described two possible schemes of implementing the method. The first one assumes calculation of m R from Eq. (52) what in practice is restricted to situations where solutions of Bogoliubov equations (44) are known analytically. In such case an analytical expression of m R , might be also accessible (see [34] as an example). On the other hand, if only numerical solutions of Bogoliubov equations (44) are in reach, semiclassical calculations become a solution of the problem and are the only possibility in practice. In case of numerical approach only low energy modes are available because every numerical approach uses a finite lattice spacing which introduces a high momenta cut-off.
Some comments are in order now. The first question is if both m R (r, r) and δn(r, r) are finite. It follows from Eqs. (50) and (51) that high energy components of m R and δn, are finite indeed. But it is not necessarily true for low energy components, m L and δn L . Obviously they are finite if we deal with discrete energy levels because of finite number of states at low energy sector. But if a system is not confined spatially in only one direction, contributions of low energy modes scale like 1/k and infrared problem appears. In such a case both δn and m are infinite and one needs to introduce Bogoliubov method using density-phase representation. This, however, is not a subject of the present paper. If a spatial confinement is missing in more than one direction the same kind of 1/k scaling does not lead to the infrared catastrophe. Now, we comment on obtaining the renormalized anomalous density using semiclassical method. The error of this method is rooted in replacing a sum over discreet Bogoliubov energies by an integral over wavevectors ν → dk. At the lower limit of integration, which is a sphere of radius k c , there is an uncertainty in a 'smooth' connection of low (discreet) and high (continuous) momenta sectors. In fact, the cutoff can k c be chosen anywhere between the two neighboring discreet energies E c , and E c + ∆E c , where 2 2ma k 2 c = E c . Evidently it leads to some uncertainty, ∆k c , in determination of the cutoff momentum: All this comes about to an error ε in the anomalous density to be of the order: In obtaining Eq. (56) we assumed, based on Eq. (48), that high-energy contribution to the anomalous density has the form, m H ∝ kc dk 1 E(k) . If ∆k c decreases with increasing the cutoff momentum k c , the error goes eventually to zero. However the problem appears in the case when ∆k c stays constant. Below we give examples of the both cases.
The first one it is a system confined in the z-direction by a harmonic potential of a frequency ω z , while it is not confined in the remaining directions. Such system is discussed in [49]. At high energies the discreet energy levels are equally spaced ∆E c = ω z = const. and from Eq. (55) we obtain ∆k c ∝ 1/k c , i.e. the error goes to zero.
In the second case we consider a system confined by a box-like potential of the z-edge size equal to L, with periodic boundary conditions in this direction (considered in [34]). The system is not confined in the two other directions. In such a case ∆E c ∝ k c and Eq. (55) gives a constant cutoff-independent error. In this case it turns out that to get a correct result, the value of k c has to be chosen exactly at the middle between the two energy levels i.e. k c = 2π L n c + 1 2 . We use the above method in the case of harmonic trapping and therefore we do bot need to worry about the above discussed error.

III. BOSE-BOSE MIXTURE
Having analyzed a single component Bose gas we now move to a case of Bose-Bose mixtures where quantum droplet state might exist.
The system is described by the Hamiltonian being a sum of a single particle Hamiltonian, H 0 , and two-body interaction energy. The first term contains kinetic energy as well as energy related to external trapping potential of the two species. The two-body interaction is assumed to be of the form: where U i (r) and U 12 (r) are interaction potentials. We assume that the two components, denoted as '1' and '2', have equal masses (to simplify the calculation). The standard mean field approach is based on the two coupled stationary Gross-Pitaevskii equations. This approach predicts a transition from a homogeneous solution to a state which is localized and eventually collapses (tends to infinite density). The transition occurs when interparticle interactions are properly tuned. In a simplest case of uniform mixture of species with equal masses the instability occurs when √ g 11 g 22 + g 12 ≤ 0. Just before the transition point, on its stable side, Bogoliubov equations [11,51] support two kinds of excitation. They are known as soft and hard modes. A velocity of soft modes goes to zero when approaching the transition point and becomes imaginary just passed the transition. A homogeneous solution becomes unstable in this region. On the other hand the sound velocity of hard modes is always real and much larger than the sound velocity of the soft modes. Thus, we consider in the following only hard-modes contribution to the LHY energy. This energy does not change significantly while crossing the transition point. It is enough though, to calculate the LHY energy directly at the transition and apply this expression also to the systems at a proximity to this point.
We need to add that in the approach presented below, one could obtain Bogoliubov equations of the soft modes as well. And then, the energies resulting from this equations would be real. Still this would lead to more complicated calculations, which almost does not change the properties of the system. That is why, in what follows, soft modes are not taken into account.
To further simplify calculations we assume that a twoparticle interaction potential is the same for both atomic species, U 11 = U 22 = U . Moreover we assume that both gases have equal number of atoms N and they are both placed in the same external potential V (r). The standard description of excitation of such system is given by the Bogoliubov method where response of the system to small perturbations is analyzed [52]: where ψ 1,2 = ψ 1,2 are condensate wave-functions, ψ 1 = ψ 2 . Due to the assumed symmetry of exchanging of the species 1 ↔ 2, it is convenient to introducê Inserting Eq. (58) into the above we find that: Introduction of fieldsψ ± allows to decoupleψ − modes fromψ + . The only non-vanishing mean field is the mean field of the soft mode ψ ≡ ψ + .
We have argued that contribution of fluctuations of soft modesδ + to the LHY energy can be neglected. The hard mode fluctuationsδ ≡δ − give the only important contribution to the LHY term. Justification of the introduced above notions of the 'soft' and 'hard' modes will become clear when we introduce GGP equation and Bogoliubov equations. Summarizing this discussion, at the proximity to the transition it is sufficient to consider only two fields: the mean field of the soft mode ψ(r, t) and fluctuations of the hard modeδ(r, t).
Bogoliubov equations giving quantum fluctuations,δ, can be obtained by linearization of the Heisenberg equation for this operator δ(r, t) = e −iµt ν u ν (r)e −iεν tα ν + v * ν (r)e iεν tα † ν : We assume U 12 = −U , which is justified close to the transition point and substitutedψ(r, t) = e −iµtψ (r) and ψ(r, t) = e −iµt ψ(r). We choose ψ(r) to be a real function. The above equations indicate that excitation energies ε ν are proportional to energy of repulsive interactions, dr ′ ψ(r)U (r − r ′ )ψ(r ′ ) ∝ gψ 2 (r). It does not vanish at the critical point. This justifies the name 'hard modes' coined for excitation triggered by perturbationδ. Now we turn our attention to the mean field ψ. Averaging the Heisenberg equation over the vacuum of hard modes and retaining quadratic terms inδ only, in particular neglecting all terms involvingδ + , we get: where δn(r ′ , r) = δ † (r ′ , t)δ(r, t) , m(r ′ , r) = δ (r ′ , t)δ(r, t) , U s = (U + U 12 )/2 and U d = (U − U 12 )/2. One can check that δn and m are indeed real and time independent functions.
The problem is analogical to a single component gas, Eq. (38) -we have to show that Eq. (64) does not depend on details of interaction potentials, but only on scattering lengths related to them. As previously we argue that both the density n(r) = ψ 2 (r) as well as δn(r, r ′ ) change on the length scale given by a healing length ξ being much larger than a range of the interaction potentials σ assumed to be the smallest length scale, ξ ≫ σ. There-fore the sum of all terms, C, in Eq. (64) involving integrals and interaction potentials U s and U d can be simplified: C ≃ũ s (0)ψ 2 (r)ψ(r) + gδn(r)ψ(r) In the above we used δn(r) = δn(r, r) and approximated dr ′ U (r−r ′ )m(r ′ , r)ψ(r ′ ) ≃ ψ(r) dr ′ U (r−r ′ )m(r ′ , r). This term depends on the range of the potential σ, thus on high energy modes.
Again, we notice that the dominant term in Eq. (65) is proportional to the density of atoms, ψ(r) 2 ≫ δn(r), thus while substitutingũ s (0) = (ũ(0)+ũ 12 (0))/2 by scattering length g and g 12 , we must use the second order Born expansion for both U and U 12 potentials, Eq. (27). These second order terms, depending on high energy modes, together with the anomalous density, give regularized anomalous density: depending only on low momenta k ≃ 0 part of the interaction potentials i.e. on their scattering length only. At the proximity to the critical point where U 12 ≃ −U the above expression simplifies This way we recover the problem of a single component inhomogeneous Bose gas and entire discussion of the previous section does apply.
Our self-consistent description of Bose-Bose mixture at the proximity to the transition to the droplet state is based on the GGP equation: and Bogoliubov equations Eq. (63) where we set U (r) = gδ(r). In Eq. (68) we introduced g 12 = −g + δg.
As before, we deal here with problem of chemical potential µ appearing in Bogoliubov equations. To have a consistent gapless approach we shall replace the chemical potential µ → µ 0 in these equations. The chemical potential µ 0 has to be found from Bogoliubov equation determining the zero-mode wavefunction, u 0 (r) = −v 0 (r), when excitation is set to zero, [48], ε 0 = 0: The replacement µ → µ 0 ensures that excitation spectrum is gapless and amplitudes of Bogoliubov modes have correct limit at low energies. However, a time dependence of fluctuations is the same as it was assumed, i.e. δ(t) = e −iµt ν u ν (r)e −iεν t + v ν (r)e iεν t . As a result of the above discussion, Bogoliubov equations of the hard mode take the form: Solutions of Bogoliubov equations allow to find quantum depletion δn(r) = ν =0 |v ν (r)| 2 and renormalized anomalous density, m R (r) = ∂ ∂|∆r| (|∆r|m (R, ∆r)) ∆r→0 , where m(r, r ′ ) = ν =0 u ν (r)v * ν (r ′ ) and R = r+r ′ 2 and ∆r = r−r ′ . To formally finish formulation of the method we supply it by a normalization condition: Here we need to add that the above method used in the case of the systems studied in [34] (gas placed between a two infinite plates -uniform system and periodic boundary conditions are used) gives exactly the same results as obtained in [34].
We now consider a system where numerical solution of the Bogoliubov equations (70-70) is a necessity. Then, as in the single component case, using semiclassical method, we obtain where k c (r) is given by equation ε ν (k c (r), r) = E c and A(k, r) = 2 k 2 2ma + V (r) − µ 0 + gψ 2 (r), B(r) = gψ 2 (r). Using the same method we find Note, that the problem has to be solved self-consistently because δn and m R depend on ψ which in turn is a solution of GGP equation, Eq. (68), which involves δn and m R as essential ingredients.

IV. CONCLUSIONS AND FUTURE OUTLOOK
In this paper we formulate the method allowing for a self-consistent treatment of coupled Generalized Gross-Pitaevskii and Bogoliubov de-Gennes equations to obtain ground state and low-energy excitations of a Bose-Bose mixture accounting for higher order terms in expansion of system energy in diluteness parameter. We argue, that while terms originating in quantum fluctuations are important at the level of Generalized Gross-Pitaevskii equations, they should be omitted in Bogoliubov-de Gennes equations. Consistently, a heuristic modification of the value of chemical potential by an amount being of the same order as other neglected terms, gives correct spectrum of excitation energies and corresponding eigenvectors of Bogoliubov-de Gennes equations.
In more details, our approach is based on three pillars which define closed, self-consistent system of equations. The pillars are: i) Generalized Gross-Pitaevskii equation allowing to find droplet's wavefunction ψ. This equation accounts for quantum fluctuations, i.e. quantum depletion and renormalized anomalous density given by solutions of Bogoliubov-de Gennes equation. ii) Modified Bogoliubov-de Gennes equations allowing to find quantum depletion and anomalous density. The Bogoliubovde Gennes equations are coupled to the Generalized Gross-Pitaevskii equation via the mean field, ψ. Modification of chemical potential entering Bogoliubov-de Gennes equations ensures gapless phononic excitation spectrum. iii) Regularization of anomalous density allowing to remove a cut-off-dependent high energy contribution.
The MGHFB method defined in this paper is applied by us in [49] to find quantum contribution to the LHY chemical potential of Bose-Bose droplets confined in one direction by a harmonic potential in the whole range of geometries from 2D to 3D case, including the crossover region.
Our approach should be also very usefull to describe small droplets squeezed by external potentials in some directions. In such case accounting for a back-action of a modified by the LHY term droplet density on the LHY energy in a self-consistent way is necessary.
Eq.(A9) and Eq.(A10) allow to split the 'quadratic in energy' Hermitian operator on the left hand side of Eq.(A8) into 'large', G(r), and 'small' δG(r) components: δG(r) = E k E k + 2(H (r) +ũ(0)ψ(r) 2 ) , (A12) and write Eq.(A8) as follows: Note that the lowest energy Bogoliubov mode is equal to g − 0,0 (r) = u 0 (r), thus G(r)g − 0,0 (r) = H(r)u 0 (r) = 0. Modes of the lowest excitation energy can be approximated by plane waves of momentum k ⊥ on top of g − 0,k (r) ≈ u 0 (r) profile, i.e. f − 0,0 = e ik ⊥ r ⊥ u 0 (r). The first order of perturbation in E k gives the low energy excitation spectrum: The second order perturbation gives corrections proportional to δG 2 which are of the order of E 2 k . This completes our proof that MGHFB method leads to phononlike gapless low energy excitation spectrum.
Finally we want to show that substitution of chemical potential µ by µ 0 is crucial for the above result. If the value µ which results form the GGPE is used in Bogoliubov equations, the zero mode eigenvector u 0 (r) and corresponding eigenenergy ε 0 are to be found from the equation: (H 0 (r) − µ 0 + (µ 0 − µ) +ũ(0)ψ 2 (r))u 0 (r) = ε 0 u 0 (r).
(A15) Comparing the above equation with Eq.(A7) one finds the excitation energy of the zero mode, ε 0 = µ 0 −µ ≡ δµ. It follows from Eq.(A14), that in in such a situation the excitation spectrum misses a phonon branch and has an energy gap: ε 2 k = δµ(δµ+2 u 0 |ũ(0)ψ 2 |u 0 )+2E k (δµ+ u 0 |ũ(0)ψ 2 |u 0 ). (A16) The relative difference between the two chemical potential can be estimated as µ − µ 0 ≪ gn, i.e. the modification is of the order quantities gδn or gm which are not accounted for in Bogoliubov equations. Therefore, replacement µ → µ 0 is of the same order of accuracy as other approximations made. Note however, that even if δµ is small, the energy gap, Eq.(A16), is much larger. It is of the order of √ gnδµ, which is not a small parameter.
where ̺ ν is the constant phase-space density which is found from normalization condition 1 = dr |C ν (r)| 2 = drdp (2π ) 3 δ (ε ν − ε ν (p/ , r)) ̺ ν . (B13) In the above u ν,en and v ν,en are still not fully defined as we have only modulus of these functions. A direction of k(r) is also not defined. For our purposes one can assume u ν,en (r), v ν,en (r) ≥ 0 and k(r) = k(r)e k where e k is the unit vector independent of r which we choose to have uniform distribution on the unit sphere.
We find: In the above the integrand is important only for such k, for which we may approximateũ(k) ≃ũ(0) ≃ g. This substitution gives: A(k, r) = E k + V (r) − µ 0 + 2gψ 2 (r), B(k, r) = gψ 2 (r) ε ν (k, r) = A 2 (k, r) − B 2 (k, r) As written in the main body of the paper we want to calculate m R (r) directly from solutions of Bogoliubov equations (44) where we use g instead of dr ′ U (r − r ′ ).
In such case we may also use the semiclassical method to solve these equations. Repeating the steps as above we obtain By noticing that 1 2εv = 1 2 k 2 /ma + (− 1 2 k 2 /ma + 1 2εv ) and having in mind that we are interested in the limit |∆r| → 0, we rewrite the above as .