Many-body and temperature effects in two-dimensional quantum droplets in Bose–Bose mixtures

We study the equilibrium properties of self-bound droplets in two-dimensional Bose mixtures employing the time-dependent Hartree–Fock–Bogoliubov theory. This theory allows one to understand both the many-body and temperature effects beyond the Lee–Huang–Yang description. We calculate higher-order corrections to the excitations, the sound velocity, and the energy of the droplet. Our results for the ground-state energy are compared with the diffusion Monte Carlo data and good agreement is found. The behavior of the depletion and anomalous density of the droplet is also discussed. At finite temperature, we show that the droplet emerges at temperatures well below the Berezinskii–Kosterlitz–Thouless transition temperature. The critical temperature strongly depends on the interspecies interactions. Our study is extended to the finite size droplet by numerically solving the generalized finite-temperature Gross-Pitaevskii equation which is obtained self-consistently from our formalism in the framework of the local density approximation.

www.nature.com/scientificreports/ DMC data and the previous theoretical results for the energy and the equilibrium density 31 . Regarding self-bound droplets of single dipolar BECs, the TDHFB provides also satisfactory explanations to experimental results and gives best match with the latest QMC simulation 23 .
In this paper, we investigate many-body effects and impacts of higher-order quantum fluctuations on the ground-state properties of self-bound droplets of 2D symmetric Bose mixtures at both zero and finite temperatures using our HFB theory. At finite temperature this exotic states of matter remains largely unexplored most likely due to the self-evaporation i.e. non existence of collective excitations below the particle-emission threshold 1 . We calculate analytically the contribution to the sound velocity, the ground-state energy and the free energy from higher-order quantum and thermal fluctuations. At zero temperature, the energy has a minimum at a finite density corresponding to a self-bound liquid-like droplet state. The obtained ground-state energy shows an excellent concordance with the DMC results of Ref. 13 , indicating the relevance of our model. We analyze also the behavior of the depletion and the anomalous correlations of the droplet in terms of the equilibrium density. At finite temperature, we find that the self-bound droplet may occur only at a certain critical temperature well below the Berezinskii-Kosterlitz-Thouless (BKT) transition due to the crucial role played by thermal fluctuations effects. Such a critical temperature decreases as the strength of interspecies interactions grows. Furthermore, we show that our formalism provides an extended finite-temperature GPE in which higher-order logarithmic factors are added to the nonlinear term of the condensate. We use this model and discuss in particular the role of the quantum fluctuations play in the density profiles and the width of the droplet. To the best of our knowledge this is the first theoretical investigation of 2D self-bound Bose mixtures at finite temperature in the presence of higher-order corrections.

Fluctuations and thermodynamics of 2D Bose mixtures
We consider a weakly interacting 2D Bose mixture with equal masses. The dynamics of this system including the effect of quantum and thermal fluctuations is governed by the coupled TDHDB equations which can be written in compact form as 29-32 : where ρ j (r, t) is the single particle density matrix of a thermal component defined as jψ j + g j n 2 j /2 + g 12 drn 1 n 2 + E LHY , is the energy of the system with E LHY = 2 j=1 (g j /2) dr 2ñ j n j −ñ 2 j + |m j | 2 +m * j � 2 j +m j � * j 2 being the LHY correction to the energy.
In Eqs.(1) h sp j = −( 2 /2m j )� − µ j is the single particle Hamiltonian, µ j is the chemical potential of each component, δµ jLHY (r)� j (r) = g j ñ j (r)� j (r) +m j (r)� * j (r) is the relevant LHY term which is obtained self-consistently, ψ j (r) =ψ j (r) − � j (r) is the noncondensed part of the field operator with � j (r) = �ψ j (r)� , n cj = | j | 2 is the condensed density, ñ j = �ψ † jψ j � is the noncondensed density, m j = �ψ jψj � is the anomalous correlation, and n j = n cj +ñ j is the total density of each species. In 2D Bose gases, the intra-and interspecies coupling strengths are given by g j = 4π 2 / m ln 4e −2γ /a 2 j κ 2 , and g 12 = g 21 = 4π 2 / m ln 4 2 e −2γ /a 2 12 κ 2 , where a j and a 12 being the 2D scattering lengths among the particles (see, e.g., 13,33,34 ), γ = 0.5772 is Euler's constant. An adequate value of the cutoff κ can be obtained in the weakly interacting regime. In such a case, attraction (repulsion) can be reached when the scattering lengths are exponentially large (small) compared to the mean interparticle separation 13 .
The presence of the noncondensed and anomalous densities in Eq. (1) enables us to derive higher-order quantum corrections without any ad-hoc assumptions in contrast to the standard GPE. In our formalism ñ j and m j are related with each other via This equation which steems from the conservation of the Von Neumann entropy, represents the variance of the number of noncondensed particles 34,35 . Equation (2) clearly shows that the anomalous density is not negligible even at zero temperature ( I → 1 ), contrary to what has been argued in the literature. Hence, m is crucial for the stability of Bose gases. Its involvement in such systems leads to a double counting of the interaction effects 30 .
In order to calculate the elementary excitations and fluctuations of a homogeneous Bose mixture, we linearize Eqs. (1) using the generalized random-phase approximation (RPA): � j = √ n cj + δ� j , ñ j =ñ j + δñ j , and m j =m j + δm j , where δ� j (r, t) = u jk e ik·r−iε k t/ + v jk e ik·r+iε k t/ ≪ √ n cj , δñ j ≪ñ j , and δm j ≪m j 29, 30 . Since www.nature.com/scientificreports/ we restrict ourselves to second-order in the coupling constants, we keep only the terms which describe the coupling to the condensate and neglect all terms associated with fluctuations δñ j and δm j (see Methods). The obtained second-order coupled TDHFB-de Gennes equations which are similar to the Beliaev's equations 27,36,37 , provide correction terms to the Bogoliubov formula for the energy spectrum: 2 12 , and α =ḡ 2 n c2 /ḡ 1 n c1 . Here the density-dependent coupling constants ḡ j = g j (1 +m j /n cj ) have been introduced in order to reinstate the gaplessness of the spectrum 30 .
The knowledge of the noncondensed and anomalous densities allows one to predict higher-order corrections to the free energy. In the frame of our formalism, it can be written as: where is the ground-state energy, and E 0 = 1 2 2 j=1 g j (n 2 cj + 4n cjñj + 2ñ 2 j +m 2 j + 2n cjmj ) + g 12 n 1 n 2 . The second term in Eq. (8) accounts for the LHY quantum corrections. It can be computed using the above dimensional regularization where only the bound modes that have energy lower than the magnitude of the binding energy are included in the integral 19 . The subleading term in Eq. (7) which represents the thermal effects is finite. Gathering quantum and thermal fluctuations contributions to the free energy (7), we get here we employed the identity is the Riemann zeta function. Expression (9) extends naturally the results of Petrov and Astrakharchik 13 since it takes into account both many-body and temperature effects.
It is worth stressing that Eqs. (3)-(9) are self-consistent and must be solved iteratively.

Self-bound droplets
Now, we consider 2D symmetric Bose mixture with repulsive intraspecies interaction and attractive interspecies interaction where a −1 12 ≪ √ n ≪ a −1 . The atoms are chosen to have equal intra-component scattering lengths a 1 = a 2 = a and equal atom densities n 1 = n 2 = n , ñ 1 =ñ 2 =ñ , and m 1 =m 2 =m . For the sake of simplicity we put = m = 1.
Zero-temperature case. At zero temperature, the properties of self-bound ultradilute Bose mixtures can be analyzed by minimizing the ground-state energy with respect to the density or equivalently using the zeropressure condition P = µn − E/S = 0 , where S is the surface area 13 . According to the method outlined in Ref. 13 , we introduce a new set of coupling constants given as: g = 4π/ ln 4e −2γ /a 2 ǫ 0 and g 12 = 4π/ ln 4e −2γ /a 2 12 ǫ 0 , www.nature.com/scientificreports/ where ǫ 0 = 4e −2γ /a 12 a has been choosed in such a way that the condition g 2 = g 2 12 must be fulfilled. This implies that c s− = 0 which means that ñ − =m − = 0 . Then, the corrected sound velocity can be obtained via c 2 s = 2gn 1 −ñ/2n +m/2n (here we set c s+ = c s , ñ + =ñ , and m + =m for convenience). For the purpose of analytical tractability, we keep only lowest order in ñ and m . This gives: where n 0 = ǫ c / 2ge 8π is the equilibrium density which can be obtained by minimizing the ground-state energy (9) with respect to the density. The sound velocity at the equilibrium is defined as c 2 s0 = gn 0 . Clearly, Eq. (10) predicts an imaginary sound velocity which may lead to a complex energy functional. Similar behavior has been reported in 22,27 for 3D droplets. In the Petrov's work 1, 13 , such a dynamically unstable phonon mode has been completely ignored under the assumption that its contribution is negligibly small. To stabilize the sound velocity and obtain the associated ground-state energy, we should include higher-order fluctuations (see below).
The noncondensed and anomalous densities of the droplet corresponding to the sound velocity (10) read: and Figure 1 shows that m is larger than ñ regardless of the value of ln(a 12 /a) as in the case of self-bound droplets in 3D Bose mixtures 31 . Both densities are increasing with decreasing ln(a 12 /a).
Let us now calculate the ground-state energy for 2D symmetric Bose mixtures by seeking the effect of higherorder fluctuations where a numerical method is used to treat the involved integration. The results are depicted in Fig. 2.
We see from Fig. 2a that the variation of the energy-cutoff which depends on interspecies interactions may strongly change the position of the local minimum of the energy leading to affect the stability and the existence of the droplet. For instance, for ǫ 0 ≤ 0.1 , the local minimum disappears and the energy becomes positive indicating that the droplet may turn into a soliton-like many-body bound state in good agreement with the predictions of Refs. 15,19,41 .
In Fig. 2b we compare our results for the ground-state energy up to second order in m and ñ of the iteration method with the DMC data and the Bogoliubov theory 13 . We see that when ln(a 12 /a) gets larger, our results excellently agree with the DMC simulations and improve the standard Bogoliubov findings. This implies that for large ln(a 12 /a) , the HFB predictions become increasingly accurate due to the considerable role of higher-order terms arising from the normal and anomalous fluctuations. Our results diverge from the DMC simulations only for very small values of interspecies interaction ln(a 12 /a) < 5 and higher densities.
n eq n = ln(n/n 0 ) − 1 ln (a 12 /a) , m eq n =ñ n ln n/n 0 (ln(n/n 0 ) − 1) 8πe 8π ,   46,47 . In such a quasicondensate, the phase coherence governs only regime of a size smaller than the size of the condensate, marked by the coherence length l φ 48 . Therefore, below the BKT transition temperature one can use the HFB theory to describe the true BEC 34, 49, 50 even though it cannot predict the critical fluctuations near the BKT region. At finite temperature, the free energy becomes divergent since c s− = 0 results in an unstable droplet in contrast to the zero-temperature case. Hence, to properly study the finite-temperature behavior of the 2D self-bound droplet, the sound velocity must be finite: where δg ± = g ± g 12 .
Minimizing the resulting free energy, we could observe the equilibrium emergence of the droplet, as visible in Fig. 3. The temperature is normalized to T BKT which is defined for a symmetric mixture according to T BKT = π N ln(380/g)/S 50 . It has been demonstrated that the interspecies interaction plays a minor role near the (13) c 2 s± = δg ± n 1 + δg ± 4π ln nδg ± ǫ c e ,  www.nature.com/scientificreports/ BKT critical temperature 51,52 . We see that the free energy F diverges like n −1 as the density goes to zero due to the presence of thermal fluctuations effects. Well below the BKT transition i.e. 0 < T 0.3 T BKT , F develops a local maximum which corresponds to an unstable droplet, and a local minimum supporting a higher density stable self-bound solution. In such a regime, thermally excited atoms that occupy continuum modes are unbound and leave the droplet result in a process of self-cooling predicted earlier by Petrov 1 . The two solutions disappear at the critical temperature ( T = T c ≃ 0.3 T BKT ) revealing that the liquid-like droplet start to evaporate. Increasing further the temperature ( T > T c ), the free energy increases without any special structure and thus, the self-bound state loses its peculiar self-evaporation phenomenon and entirely destroys eventually. The same situation takes place for dipolar droplets in a single BEC [53][54][55] and in dual condensates 29,49 . Note that T c strongly relies on ǫ 0 and hence, on the interspecies interactions as we shall see below. The critical temperature above which the BEC-droplet phase transition occurs can be determined by minimizing the free energy. For δg − ≪ δg + , one has As shown in Fig. 4 for fixed density n/n 0 , the droplet critical temperature decreases with the interspecies interaction ln(a 12 /a) regardless the presence or not of the higher-order effects. For example for ln(a 12 /a) = 20 , the droplet reaches its thermal equilibrium at ultralow temperature ( T c ≃ 0.06T BKT ). We see also that higherorder corrections may reduce the critical temperature.

Generalized finite-temperature Gross-Pitaevskii equation
In this section, we consider the finite size effects on equilibrium properties of the self-bound droplet. The basic idea behind finite size contributions to the droplet's energy is that the quantities , ñ , and m must vary slowly at the scale of the extended healing length. As a consequence, we can include higher-order corrections locally as nonlinear terms in the TDHFB equations and treat them classically. For simplicity, we will ignore the dynamics of the thermal cloud and the anomalous correlations. Therefore, the TDHFB Eq. (1a) leads directly to the generalized finite-temperature GPE where α ≃ ln (en/n0) ln(n/en0)/8πe 8π , and αT ≃ ζ(3)T 3 ln(a12/a)/ 2 ln |�| 2 /e 4π+1 n0 |�| 4 [1/ ln |�| 2 /e 4π+1 n0 + 1] . Importantly, the generalized finite-temperature GPE (15) extends naturally the GPE of Ref. 13 since it takes into account higherorder quantum and thermal corrections.
The stationary solutions of Eq. (15) can be found via the transformation �(r, t) = �(r) exp(−iµt) . We solve the resulting static equation numerically using the split-step Fourier transform 4 . In Fig.5, we plot the density profiles as a function of the radial distance at both zero and finite temperatures. As can be seen in Fig.5a, the density n is flattened in accordance to the liquid character of the condensate. The obtained density is compared with the predictions of the GPE-LHY theory 13 . Our results show a slight deviation downwards for distance (14) T c T BKT ≃ ln[(n/2e 4π−2 n 0 ) − 1] 1/3 (πζ (3)) 1/3 ln(a 12 /a) ln[ln(a 12 /a)95/π] , www.nature.com/scientificreports/ r < 8 with respect to the findings of Ref. 13 owing to the higher-order quantum fluctuations. At temperatures T T c , the droplet exhibits a weak-temperature dependence (see Fig.5b). Whereas, at T ≥ T c , the droplet has a Gaussian-like shape pointing out that the system experiences droplet-BEC phase transition.
To evaluate the width of the self-bound droplet, we first use the following trial wavefunction : where R is the self-bound droplet width. Then, we minimize the resulting functional energy with respect to R. In the absence of the higher-order corrections, the width takes the form This analytical prediction is reported in Fig. 6 and compared with the numerical results of our generalized GPE (15). We see that the width of the droplet decreases exponentially versus the number of particles. The interspecies interactions a 12 /a lead also to reduce the width. The comparison between our predictions and those of Petrov 13 indicates that the higher-order quantum effects may shift the droplet width. At finite temperature one  www.nature.com/scientificreports/ can expect that the droplet size increases significantly only at temperatures T T c . Above such a temperature the self-bound droplet is in its thermal stabilization.

Discussion
We studied the equilibrium properties of symmetric self-bound droplets of 2D binary BEC beyond the standard LHY treatment, at both zero and finite temperatures. We computed higher-order corrections to the excitations spectrum, the sound velocity, the normal and anomalous correlations, and the free energy. These corrections improve the ground-state energy obtained from the Bogoliubov approach 13 predicting an energy in good agreement with recent DMC simulations owing to the non-negligeable role of higher order terms. At finite temperature, we revealed that the droplet occurs at temperature well below the BKT transition and destroys when the temperature becomes slightly larger than the ground-state energy of the droplet due to the thermal fluctuations effects. We found that the interspecies interaction tends to lower the critical temperature. We analyzed in addition the finite-size droplets in the framework of our generalized finite-temperature GPE. As outlined above, one can infer that in 2D mixtures, the droplet survives only in an ultradilute regime and at ultralow temperatures.
Our results could be extended in weakly interacting quasi-2D Bose mixtures as long as the following condition is fulfilled 0 < −a 3D 12 < a 3D ≪ l 0 13 , where a 3D and a 3D 12 are the 3D intra and interspecies scattering lengths, and l 0 is the oscillator length in the confinement direction. The creation of such 2D mixture droplets in the experiment, still remains a challenging question.

Methods
Derivation of the condensate fluctuations. As we concluded in the main text, for a thermal distribution at equilibrium and by working in the momentum space, one has 34 where ρ mn (k) is the Fourier transform of ρ mn (r − r ′ ) . After some algebra, expression (2) turns out to be given as: From Eq. (18) we can straightforwardly derive the expressions (3) and (4) describing the normal and anomalous correlations. For an ideal Bose gas where the anomalous density vanishes, I k = coth 2 (E k /2T) 30 .

TDHFB-de Gennes equations.
We use the generalized RPA which consists of imposing small fluctuations of the condensates, the noncondensates, and the anomalous components, respectively, as: � j = √ n cj + δ� j , ñ j =ñ j + δñ j , and m j =m j + δm j , where δ� j ≪ √ n cj , δñ j ≪ñ j , and δm j ≪m j 30 . We then obtain the TDHFB-RPA equations: and Remarkably, this set of equations contains a class of terms beyond second order. Note that we keep in Eqs. (19) and (20) only the terms which describe the coupling to the condensate and neglect all terms associated with δñ and δm owing to the fact that we restrict ourselves to second-order in the coupling constants.
Inserting the transformation δ� j (r, t) = u jk e ik·r−iε k t/ + v jk e ik·r+iε k t/ into Eq. (19), we find the second-order coupled TDHFB-de Gennes equations for the quasiparticle amplitudes u kj and v kj : where dr[u 2 j (r) − v 2 j (r)] = 1 , L j = E k + 2ḡ j n cj + 2g jñj + g 12 n 3−j − µ j , M j =ḡ j n cj , and A = g 12 √ n c1 n c2 . Equations (21) are appealing since they enable us to calculate in a simpler manner corrections to the excitations spectrum ε k± of homogeneous Bose mixtures (see the main text).