Ultradilute Quantum Droplets in the Presence of Higher-Order Quantum Fluctuations

Abstract

We investigate the effects of higher-order quantum fluctuations on the bulk properties of self-bound droplets in three-, two- and one-dimensional binary Bose mixtures using the Hartree–Fock–Bogoliubov theory. We calculate higher-order corrections to the equation of state of the droplet at both zero and finite temperatures. We show that our results for the ground-state energy are in a good agreement with recent quantum Monte Carlo simulations in any dimension. Our study extends to the finite temperature case where it is found that thermal fluctuations may destabilize the droplet state and eventually destroy it. In two dimensions, we reveal that the droplet occurs at temperatures well below the Berezinskii–Kosterlitz–Thouless transition temperature.

1. Introduction

Quantum fluctuations which are crucially tied to Heisenberg’s uncertainty principle are ubiquitous in nature. This fascinating phenomenon is often associated with a change in the energy of the system. Quantum fluctuations can alter virtually all aspects of matter leading to the emergence of novel phase transitions such as quantum liquid droplets, supersolids and superglasses.

Ultracold quantum gases offer new vistas for the quantum simulation of interacting many-body systems. It has long been common that the presence of quantum fluctuations in a Bose gas can modify its equation of state, which at first order is given by the seminal Lee–Huang–Yang (LHY) correction [1]. On the theoretical side, the LHY correction to the ground-state energy has been calculated for the Bose–Einstein condensate (BEC) with a pure contact interaction [1,2,3,4,5,6] and for dipolar BECs [7,8,9,10] using different approaches. Experimental evidence of the LHY for a single BEC has been reported in [11,12,13,14] for both bosonic and fermionic systems.

On the other hand, the LHY corrections to the ground-state energy for homogeneous Bose–Bose mixtures were addressed in refs. [15,16] using the Bogoliubov theory. At finite temperatures, such LHY fluctuations have been computed using different models (see e.g., [17,18,19]). In d-dimensional binary BECs, quantum stabilization stems from the interplay of the attractive mean–field interaction, proportional to the density of the system n, and the repulsive LHY correction proportional to $n^{d/2}$ [20,21]. Quite recently, quantum corrections of self-bound Bose mixtures beyond LHY have been discussed in several works [22,23,24,25].

A similar stabilization scenario has been observed in single and binary dipolar Bose gases, where the competition between the mean-field energy associated with short and long-range interactions and the LHY corrections arrests the dipolar instability at high density, forming ultradilute quantum droplets [26,27,28,29,30]. Furthermore, quantum fluctuations play an intriguing role in low-dimensional dipolar BECs. They change their nature from repulsive to attractive due to the peculiar momentum dependence of the dipole–dipole interactions [8,31,32]. The dipolar instability is halted by such unconventional beyond mean-field corrections, leading to the formation of a droplet phase. The LHY quantum corrections lead also to shift the equation of state of disordered BECs [33,34].

Recent numerical results based on Quantum and Diffusion Monte Carlo simulations (QMC and DMC) [21,35,36,37,38,39,40] have verified the essential role played by the quantum fluctuations for the description of the self-bound droplets of both dipolar and nondipolar BECs realized in recent experiments. One year after the theoretical prediction, the experiment performed by Pfau et al. [26,27,28] confirmed the existence of quantum droplets in dipolar Bose gases. This exotic state of matter has also been successfully realized in experiments on the homonuclear mixture [41,42,43] and the heteronuclear mixture [44].

The aim of this paper is to investigate the influence of higher-order quantum fluctuations upon quantum droplets of symmetric Bose–Bose mixtures in three (3D), two (2D) and one (1D) dimensions using the self-consistent Hartree–Fock–Bogoliubov (HFB) theory. This model which, we developed recently to remove the handicaps of the standard Bogoliubov prescription [20,21], can self-consistently evaluate the quantum and thermal fluctuations exhibited by weakly interacting binary BECs [18,24,25]. It has been shown that the HFB theory is able to produce excellent predictions for the condensate, depletion, pair correlation function, critical number of particles and ground-state energy that have been measured experimentally and with QMC [18,24,25,30].

In 3D geometry, we show that at zero temperature, our theory captures genuine higher-order quantum effects predicted from recent DMC simulations for the ground-state energy [37]. At finite temperature, we calculate the free energy and find that the droplet destabilizes when the temperature becomes slightly larger than the ground-state energy of the droplet. Our study is extended to the regime of relatively strong interactions. We reveal that such higher-order quantum corrections become more pronounced for large values of the interspecies interaction. The latter may also considerably shift the critical temperature.

In the 2D case, the obtained ground-state energy shows a good agreement with the DMC results of ref. [21], indicating the relevance of our model. At finite temperature, we find that the self-bound droplet can survive only at a certain critical temperature well below the Berezinskii–Kosterlitz–Thouless (BKT) transition due to the crucial role played by thermal fluctuation effects. Such a critical temperature decreases as the strength of interspecies interactions grows.

We then confine ourselves to analyze the bulk properties of 1D quantum droplets. Beyond the LHY, quantum corrections to the ground-state energy are accurately worked out and an excellent agreement is found with the QMC of ref. [38]. Similarly to the 3D and 2D cases, we point out that at finite temperature, the thermal fluctuations are able to destroy the droplet.

2. HFB Theory

We consider weakly interacting two-component BECs with equal masses, $m_1=m_2=m$, in d-dimensions. The dynamics of this system, including the effect of quantum and thermal fluctuations, is modeled by the coupled HFB equations, which can be written as [18,24,25,29,30]:

$$i\hslash \frac{d{\mathsf{\Phi }}_j}{dt}=\left(h_j^{sp}+g_j{n}_j+g_{12}{n}_{3-j}+\delta {\mu }_{j\mathrm{LHY}}\right){\mathsf{\Phi }}_j,$$

(1)

where $h_j^{sp}=-({\hslash }^2/2m)\mathsf{\Delta }+U\left(\mathbf{r}\right)-{\mu }_j$ is the single particle Hamiltonian, $U\left(\mathbf{r}\right)$ is an external potential, ${\mu }_j$ is the chemical potential of each component and $g_j$ and $g_{12}=g_{21}$ are the intraspecies and the interspecies coupling constants, respectively. In order to ensure the gaplessness of the spectrum, we renormalize the coupling constants as: ${\overline{g}}_j=g_j(1+{\tilde{m}}_j/n_{cj})$ [18]. Such a renormalization scheme was used for a single-component Bose gas a long time ago (see for reviews [45,46]) and provides excellent results for the density profiles and the collective modes, and the agreement with experiments is better than other approaches. In Equation (1), ${\widehat{\overline{\psi }}}_j\left(\mathbf{r}\right)={\widehat{\psi }}_j\left(\mathbf{r}\right)-{\mathsf{\Phi }}_j\left(\mathbf{r}\right)$ is the noncondensed part of the field operator with ${\mathsf{\Phi }}_j\left(\mathbf{r}\right)=\langle {\widehat{\psi }}_j\left(\mathbf{r}\right)\rangle$, $n_{cj}={\left|{\mathsf{\Phi }}_j\right|}^2$ is the condensed density, ${\tilde{n}}_j=\langle {\widehat{\overline{\psi }}}_j^{†}{\widehat{\overline{\psi }}}_j\rangle$ is the noncondensed density, ${\tilde{m}}_j=\langle {\widehat{\overline{\psi }}}_j{\widehat{\overline{\psi }}}_j\rangle$ is the anomalous correlation and $n_j=n_{cj}+{\tilde{n}}_j$ is the total density of each species. The presence of these quantities enables us to derive the LHY term without any ad hoc assumptions, in contrast to the standard generalized GPE. The relevant LHY term is given by

$$\delta {\mu }_{j\mathrm{LHY}}\left(\mathbf{r}\right){\mathsf{\Phi }}_j\left(\mathbf{r}\right)=g_j\left[{\tilde{n}}_j\left(\mathbf{r}\right){\mathsf{\Phi }}_j\left(\mathbf{r}\right)+{\tilde{m}}_j\left(\mathbf{r}\right){\mathsf{\Phi }}_j^{*}\left(\mathbf{r}\right)\right],$$

(2)

which can be calculated self-consistently.

Let us now assume a uniform binary Bose mixture. The Bogoliubov excitations energy can be obtained by linearizing Equation (1), employing the generalized random-phase approximation [18]: ${\mathsf{\Phi }}_j=\sqrt{n_{cj}}+\delta {\mathsf{\Phi }}_j$, ${\tilde{n}}_j={\tilde{n}}_j+\delta {\tilde{n}}_j$, and ${\tilde{m}}_j={\tilde{m}}_j+\delta {\tilde{m}}_j$, where $\delta {\mathsf{\Phi }}_j(\mathbf{r},t)=u_{jk}{e}^{i\mathbf{k}·\mathbf{r}-i{\epsilon }_{k}t/\hslash }+v_{jk}{e}^{i\mathbf{k}·\mathbf{r}+i{\epsilon }_{k}t/\hslash }\ll \sqrt{n_{cj}}$, $\delta {\tilde{n}}_j\ll {\tilde{n}}_j$, and $\delta {\tilde{m}}_j\ll {\tilde{m}}_j$ [8,9]. Keeping only terms up to second order in the coupling constants, we obtain for symmetric mixtures with $n_1=n_2=n/2$, ${\tilde{n}}_1={\tilde{n}}_2=\tilde{n}/2$, ${\tilde{m}}_1={\tilde{m}}_2=\tilde{m}/2$ and $g_1=g_2=g$ the following formula of the Bogoliubov energy spectrum:

$${\epsilon }_{k\pm }=\sqrt{E_k^2+2E_k{\mu }_{\pm }},$$

(3)

where $E_k={\hslash }^2{k}^2/2m$, and ${\mu }_{\pm }=n_c\overline{g}(1\pm g_{12}/\overline{g})$. For small momenta $k\to 0$, the Bogoliubov dispersion relation is phononlike ${\epsilon }_{k\pm }=\hslash c_{s\pm }k$, where $c_{s+}=c_s\sqrt{1+g_{12}/\overline{g}}$, and $c_{s-}=c_s\sqrt{1-g_{12}/\overline{g}}$ are the sound velocities in the density and spin channels [19], respectively, with $c_s=\sqrt{\overline{g}{n}_c/m}$ being the corrected sound velocity of a single BEC [6]. In the opposite limit $k\to \infty$, the excitations spectrum (3) reduces to the free particle law.

The condensed depletion and the anomalous density are defined in d-dimensions as [6]:

$${\tilde{n}}_{\pm }=\frac{1}{2}\int \frac{d^{d}k}{{\left(2\pi \right)}^d}\left[\frac{E_k+{\mu }_{\pm }}{{\epsilon }_{k\pm }}\sqrt{I_{k\pm }}-1\right],$$

(4)

and

$${\tilde{m}}_{\pm }=-\frac{1}{2}\int \frac{d^{d}k}{{\left(2\pi \right)}^d}\frac{{\mu }_{\pm }}{{\epsilon }_{k\pm }}\sqrt{I_{k\pm }},$$

(5)

where $I_{k\pm }={\cot \mathrm{h}}^2\left({\epsilon }_{k\pm }/2T\right)$ [18]. Equations (4) and (5) are self-consistent since they rely on the density-dependent coupling $\overline{g}$ and must be solved iteratively.

In the spirit of our HFB theory, the energy of the system including the LHY corrections is defined as:

$$\begin{array}{c}\hfill E=E_0+\frac{1}{2}\sum _{\pm }\int \frac{d^{d}k}{{\left(2\pi \right)}^d}\left({\epsilon }_{k\pm }-2E_k-{\mu }_{\pm }\right),\end{array}$$

(6)

where $E_0=g(n_c^2+4n_c\tilde{n}+2{\tilde{n}}^2+{\tilde{m}}^2+2n_c\tilde{m})+g_{12}{n}^2$ is the ground-state energy. The second term in Equation (6) accounts for the LHY quantum corrections. It can be computed using the dimensional regularization. The subleading term in Equation (7), which represents the thermal effects, is finite.

It is worth stressing that integrals associated with anomalous density (5) and the ground-state energy (6) are ultraviolet divergent in both 3D and 2D cases and necessitate to be regularized [6,18]. We use the dimensional regularization that is asymptotically accurate for weak interactions [24,25].

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:

$$F=E+T\int \frac{d^{d}k}{{\left(2\pi \right)}^d}ln\left(\frac{2}{\sqrt{I_{k\pm }}+1}\right).$$

(7)

The minimization of the free energy with respect to the density permits us to check the existence of the self-bound droplet at finite temperatures, as we see below.

3. Three-Dimensional Droplets

The regime of interest corresponds to repulsive intraspecies interaction $g_j=4\pi {\hslash }^2{a}_j/m>0$, and attractive interspecies interaction $g_{12}=g_{21}=4\pi {\hslash }^2{a}_{12}/m<0$ with $a_j$ and $a_{12}$ being the intraspecies and the interspecies scattering lengths, respectively.

After regularization, the ground-state energy (6) takes the following dimensionless form [24]:

$$\begin{array}{cc}\hfill \frac{E_{3D}}{N{E}_{03}}=& \pi {\left(\frac{\delta a}{a}\right)}_{+}\left(n{a}^3\right)+\frac{32\sqrt{2\pi }}{15}{\left(n{a}^3\right)}^{3/2}\hfill \\ \hfill & \times \sum _{\pm }{\left(\frac{\delta a}{a}\right)}_{\pm }^{5/2}{\left(1+\frac{{\tilde{m}}_{\pm }-{\tilde{n}}_{\pm }}{n}\right)}^{5/2},\hfill \end{array}$$

(8)

where $E_{03}={\hslash }^2/m{a}^2$, and ${(\delta a/a)}_{\pm }=1\pm (a_{12}/a)$. Here, the noncondensed and the anomalous fractions can be obtained from Equations (4) and (5), respectively, as:

$$\begin{array}{cc}\hfill \frac{{\tilde{n}}_{\pm }}{n}=& \frac{5}{96\sqrt{2}}{\left(\frac{n}{n_{\mathrm{eq}}}\right)}^{1/2}{\left(\frac{\delta a}{a}\right)}_{+}\hfill \\ \hfill & \times \sum _{\pm }{\left(\frac{\delta a}{a}\right)}_{\pm }^{3/2}{\left(1+\frac{{\tilde{m}}_{\pm }-{\tilde{n}}_{\pm }}{n}\right)}^{3/2},\hfill \end{array}$$

(9)

and

$$\begin{array}{cc}\hfill \frac{{\tilde{m}}_{\pm }}{n}=& \frac{5}{32\sqrt{2}}{\left(\frac{n}{n_{\mathrm{eq}}}\right)}^{1/2}{\left(\frac{\delta a}{a}\right)}_{+}\hfill \\ \hfill & \times \sum _{\pm }{\left(\frac{\delta a}{a}\right)}_{\pm }^{3/2}{\left(1+\frac{{\tilde{m}}_{\pm }-{\tilde{n}}_{\pm }}{n}\right)}^{3/2}.\hfill \end{array}$$

(10)

Obviously, for $\tilde{m}=\tilde{n}=0$, Equation (8) reduces to Petrov’s equation [20]. The energy (8) can be rewitten in terms of the equilibrium density $n_{\mathrm{eq}}{a}^3=\left(25\pi /16384\right){\left(\delta a_{+}/a\right)}^2$ predicted by Petrov [20] as:

$$\begin{array}{cc}\hfill \frac{E_{3D}}{E_{03}}& =-3\left(\frac{n}{n_{\mathrm{eq}}}\right)+\frac{1}{2\sqrt{2}}{\left(\frac{n}{n_{\mathrm{eq}}}\right)}^{3/2}\hfill \\ \hfill & \times \sum _{\pm }{\left(\frac{\delta a}{a}\right)}_{\pm }^{5/2}{\left(1+\frac{{\tilde{m}}_{\pm }-{\tilde{n}}_{\pm }}{n}\right)}^{5/2}.\hfill \end{array}$$

(11)

It is useful to compare our analytical expression which we solve iteratively with the energy obtained from the theory of [22] based on the inclusion of higher-order terms in the sound velocities, calculated perturbatively, and with the QMC simulations [37].

Figure 1 shows that our findings are in excellent agreement with the QMC simulations [37], practically in the whole range of density, and improve upon the existing theoretical results. This reveals the importance of higher-order fluctuations in the stabilization and in the robustness of the droplet. We see also as $|a_{12}/a|$ become larger, both the HFB theory and Petrov’s approach [20] diverge from the QMC data [37] in the high-density limit $n/n_{\mathrm{eq}}\gtrsim 1$. Notice that for $a_{12}/a<-1.1$, the effects of fluctuations higher than first order are not important. However, such quantum fluctuations become more prominent for relatively large interspecies interactions, as shown in Figure 2a. This is indeed understandable since both the normal and anomalous fluctuations originate themselves from the interaction.

To gain further insights into the stabilization mechanism, we analyze the density dependence of the energy per particle in the strongly interacting regime, i.e., $|a_{12}|\gg a$. Figure 2b depicts that as interspecies attraction becomes larger the energy decreases (depth of the local minimum increases), indicating that a stable droplet still survives even for a strongly interacting regime. One can expect that this behavior persists even when $a_{12}$ crosses the Feshbach resonance [23]. It is noteworthy that strong interactions and higher-order fluctuations do not considerably affect the equilibrium density of the self-bound droplet. The latter can be obtained by minimizing the energy with respect to the density.

At finite temperature, the 3D free energy can be evaluated from Equation (7). Gathering quantum and thermal fluctuation contributions, we obtain [24]:

$$\begin{array}{cc}\hfill \frac{F_{3D}}{N{E}_{03}}=& \frac{E_{3D}}{N{E}_{03}}-\frac{\sqrt{2\pi }}{180}{\left(n{a}^3\right)}^{-5/2}{\left(\frac{T}{E_{03}}\right)}^4\hfill \\ \hfill & \times \sum _{\pm }{\left(\frac{\delta a}{a}\right)}_{\pm }^{-3/2}{\left(1+\frac{{\tilde{m}}_{\pm }-{\tilde{n}}_{\pm }}{n}\right)}^{-3/2},\hfill \end{array}$$

(12)

where $E_{3D}$ is given in Equation (8). The lowest thermal approximation ($\tilde{m}=\tilde{n}=0$) implies that the free energy (12) reduces to its standard form (i.e., without higher-order corrections). The self-consistent solutions of Equation (12) up to second-order for several values of temperature and interspecies interactions are shown in Figure 3. We see that at temperatures higher than the critical temperature ($T>T_c\approx 5\times {10}^{-4}{E}_0$) at which the thermal fluctuations compensate the repulsive quantum fluctuations, i.e., disapearence of the local minimum, the self-bound droplet is completely evaporated due to the strong thermal effects.

Importantly, the interspecies interactions may significantly increase the critical temperature. For instance, our HFB predicts a critical temeprature $T_c\approx 1.3\times {10}^{-4}{E}_0$ for $a_{12}=-1.1a$ [24], which is around four times lower than $T_c$ predicted above for $a_{12}=-3a$.

4. Two-Dimensional Droplets

Let us now consider homogeneous 2D Bose–Bose mixtures. It has been shown that below the BKT transition temperature [47,48], one can use the HFB theory to describe the true BEC [6,25,49] even though it cannot predict the critical fluctuations near the BKT region.

In 2D Bose gases, the intra- and interspecies coupling strengths are defined as: $g_j=4\pi {\hslash }^2/\left[mln\left(4e^{-2\gamma }/a_j^2{\kappa }^2\right)\right]$ and $g_{12}=g_{21}=4\pi {\hslash }^2/\left[mln\left(4^2{e}^{-2\gamma }/a_{12}^2{\kappa }^2\right)\right]$, where $a_j$ and $a_{12}$ are the 2D scattering lengths among the particles (see, e.g., [2,21,25]), and $\gamma =0.5772$ is Euler’s constant. An adequate value of the cutoff $\kappa$ 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 with the mean interparticle separation [21].

Using again the dimensional regularization to evaluate the ground-state energy (6) [25], one then analytically continues the result to finite coupling including a low-energy cutoff ${ϵ}_c={\hslash }^2{\kappa }^2/m\gg {\mu }_{\pm }$ [25]. We thus obtain

$$\begin{array}{c}\hfill E_{2D}=E_0+\frac{m}{8\pi {\hslash }^2}\sum _{\pm }{\mu }_{\pm }^{2}ln\left(\frac{\sqrt{e}{\mu }_{\pm }}{{ϵ}_c}\right).\end{array}$$

(13)

Following Petrov’s method [21], we introduce a new set of coupling constants given as: $g=4\pi /ln\left(4e^{-2\gamma }/a^2{ϵ}_0\right)$ and $g_{12}=4\pi /ln\left(4e^{-2\gamma }/a_{12}^2{ϵ}_0\right)$, where ${ϵ}_0=4e^{-2\gamma }/a_{12}a$ was chosen in such a way that the condition $g^2=g_{12}^2$ must be fulfilled. In Figure 4, we compare the resulting ground-state energy up to third order in $\tilde{m}$ and $\tilde{n}$ of the iteration method with the DMC data [21] and the Bogoliubov theory [21]. We observe that our results excellently agree with the DMC simulations and improve both the Bogoliubov theory and our recent second-order HFB calculation [25]. This indicates that the HFB predictions become increasingly accurate due to the crucial role of higher-order terms.

Collecting both quantum and thermal fluctuations contributions to the free energy (7), we obtain [25]

$$\begin{array}{c}\hfill F_{2D}=E_{2D}-\sum _{\pm }\frac{m\zeta \left(3\right)}{{\hslash }^2{\mu }_{\pm }}{T}^3,\end{array}$$

(14)

where $\zeta \left(3\right)$ is the Riemann zeta function.

At finite temperature, the minimization of the free energy (7) reveals that a higher-density stable self-bound solution supported by the local minimum can be formed only at temperatures well below the BKT transition, as seen in Figure 5a. At higher temperatures, the self-bound state loses its peculiar self-evaporation phenomenon and eventually entirely destroys owing to the thermal fluctuations effects. Remarkably, the variation of the energy-cutoff, ${ϵ}_0$, which depends on interspecies interactions, may strongly change the position of the local minimum of the free energy, giving rise to shift the critical temperature. For instance, for ${ϵ}_0=2.8$, the droplet emerges at very low temperatures $0<T\lesssim 0.18T_{\mathrm{BKT}}$, while for ${ϵ}_0=5.8$, the self-bound droplet occurs at $0<T\lesssim 0.35T_{\mathrm{BKT}}$ (see Figure 5b).

5. One-Dimensional Droplets

It is well known that in 1D geometry, the fluctuations are strong enough to destroy the Bose condensate. However, in the weak-coupling regime in which the ratio between the interaction energy and the kinetic energy is much smaller than unity, the mean-field theory can be safely used at both zero and finite temperatures.

Let us now consider a 1D symmetric homogeneous Bose mixture in a box of size L. It is worth stressing that the 1D droplets nucleate in the dominantly repulsive regime $\delta g=g+g_{12}>0$, where $g_j=-2{\hslash }^2/\left(m{a}_j\right)$, is opposite to the 3D case.

The ground-state energy can be computed straightforwardly via Equation (6), yielding:

$$\frac{E_{1D}}{L}=\frac{n^2}{4}\delta g-\frac{2\sqrt{m}}{3\pi \hslash }\sum _{\pm }{\mu }_{\pm }^{3/2}.$$

(15)

For $\tilde{n}=\tilde{m}=0$, the energy (15) reduces to that obtained in ref. [21] using the Bogoliubov prescription.

In Figure 6a, we show the ground-state energy predicted by the HFB theory up to first-order in $\tilde{n}$ and $\tilde{m}$ for $g_{12}/g=-0.7$ (i.e., $\delta g/g=0.3$) and compare it with the Bogoliubov prescription [21] and QMC simulation [38]. We find an excellent agreement between our HFB theory and the QMC data [38]. However, the discrepancy between the HFB results and the Bogoliubov theory of ref. [21] can be attributed to the quantum fluctuation corrections. In the regime of strong interactions, one can expect that the Bogoliubov theory significantly deviates from the QMC data, since it is valid only in the dilute regime.

At finite temperature, the free energy reads:

$$\frac{F_{1D}}{L}=\frac{E_{1D}}{L}-\frac{\pi \sqrt{m}{T}^2}{12\hslash }({\mu }_{+}^{-1/2}+{\mu }_{-}^{-1/2}),$$

(16)

where $E_{1D}$ is given in Equation (15). In the absence of higher-order quantum and thermal fluctuations, we recover the results of [50]. At sufficiently low temperature, the free energy has a local maximum which corresponds to an unstable droplet and a local minimum supporting a higher-density stable self-bound solution. Similarly to the 2D and 3D cases, by augmenting further the temperature, the liquid is predicted to destroy due to the thermal fluctuations.

6. Conclusions

We systematically studied the effects of higher-order quantum and thermal fluctuations on the bulk properties of utradilute quantum droplets in a Bose–Bose mixture utilizing our HFB formalism. These higher-order quantum corrections turn to improve the ground-state energy, bringing our model in excellent agreement with recent QMC and DMC findings in any dimension. Relatively strongly interacting 3D quantum droplets were analyzed. We also revealed that the self-evaporation phenomenon is affected by the higher-order correlations and interspecies interactions. The results of the present work are not only significant to understand the fundamental properties of the ultradilute quantum liquids but also offer the possibility to experimentally explore them in low dimensions.