Microscopic pairing theory of a binary Bose mixture with interspecies attractions: bosonic BEC-BCS crossover and ultradilute low-dimensional quantum droplets

In the weakly interacting regime, quantum phase of ultracold atomic Bose gases is typically determined by their mean-field interactions (Dalfovo1999, ). Attractive mean-field interactions can induce mechanical instability towards collapse (Donley2001, ). This common viewpoint, however, is radically changed due to the seminal work by Petrov (Petrov2015, ), who proposed that the mean-field collapse could be prevented by the repulsive force provided by quantum fluctuations, i.e., the celebrated Lee-Huang-Yang (LHY) correction to the energy functional (LeeHuangYang1957, ). Although the beyond-mean-field LHY correction is usually small, it can be made comparable to the mean-field energy by experimentally tuning the interatomic interactions with the Feshbach resonance technique. As a result, self-bound liquid-like quantum droplets may form, even in free space without container (Petrov2018, ; FerrierBarbut2019, ; Kartashov2019, ). Petrov’s ground-breaking proposal has now been surprisingly confirmed in single-component Bose gases with anisotropic dipolar forces (FerrierBarbut2016, ; Schmitt2016, ; Chomaz2016, ; Bottcher2019PRA, ; Bottcher2019PRX, ) and in two-component Bose-Bose mixtures with contact inter-particle interactions (Cabrera2018, ; Cheiney2018, ; Semeghini2018, ; Ferioli2019, ; DErrico2019, ). It opens a new rapidly developing research field (Bottcher2020, ), where the beyond-mean-field many-body effect could be systematically explored, both experimentally (FerrierBarbut2016, ; Schmitt2016, ; Chomaz2016, ; Bottcher2019PRA, ; Bottcher2019PRX, ; Cabrera2018, ; Cheiney2018, ; Semeghini2018, ; Ferioli2019, ; DErrico2019, ) and theoretically (Petrov2016, ; Baillie2016, ; Wachtler2016, ; Baillie2017, ; Li2017, ; Cappellaro2017, ; Jorgensen2018, ; Cui2018, ; Staudinger2018, ; Parisi2019, ; Cikojevic2019, ; Tylutki2020, ; Shi2019, ; Cikojevic2020, ; Wang2020, ; Hu2020a, ; Ota2020, ).

Despite the great success of Petrov’s proposal, strictly speaking, it is not a consistent microscopic theory. This is particularly clear for three-dimensional Bose-Bose mixtures, with which the Petrov prototype theory of quantum droplets was constructed, within the Bogoliubov approximation (Petrov2015, ; Larsen1963, ). As the mean-field solution is not stable towards collapse, one of the two gapless Bogoliubov modes becomes softened and acquires a small imaginary component. This results in a complex LHY energy functional (Petrov2015, ; Cikojevic2019, ; Hu2020a, ; Ota2020, ), which is not physical. To circumvent this technical issue, Petrov assumed a weak dependence of the LHY energy functional on the interspecies interaction strength and fixed the LHY energy to the value on the verge of the collapse where the mechanical instability first sets in (Petrov2015, ). Hereafter, we will refer to such an approximation as Petrov’s prescription. Due to the intrinsic inconsistency in the Petrov prototype theory, the resulting energy of quantum droplets shows an appreciable deviation from the numerically accurate diffusion Monte Carlo (DMC) predictions (Cikojevic2019, ). The predicted critical number for the droplet formation also seems to be larger than the one measured experimentally, both in dipolar Bose gases (Bottcher2019PRA, ) and in Bose-Bose mixtures (Cabrera2018, ).

Recently, we developed a consistent microscopic theory to remove the annoying loophole in the Petrov theory of quantum droplets (Hu2020a, ). The crucial ingredient of the theory is the inclusion of a weak bosonic pairing between different components due to the attractive interspecies interactions. The pairing explicitly removes the unstable softened Bogoliubov excitation and turns it into a stable gapped mode, as in the conventional Bardeen-Cooper-Schrieffer (BCS) theory for interacting fermions (BCS1957, ). An apparent advantage of the pairing theory is that, it is variational and therefore predicts an upper bound for the ground-state energy. Remarkably, in three dimensions our pairing theory leads to an improved agreement with the DMC simulation, for the energy and the equilibrium density of quantum droplets (Hu2020a, ).

In this work, we would like to provide more details of the microscopic pairing theory in three dimensions (Hu2020a, ), on the equation derivation and the numerical calculation. In particular, we show that the weakly interacting pairing in our previous study can be naturally generalized to the strongly interacting regime, in which the interspecies scattering length diverges across a Feshbach resonance. Therefore, we predict the existence of a strongly interacting Bose droplet, where a significant fraction of atoms turns into bound pairs. This is precisely analogous to the crossover physics from a Bose-Einstein condensate (BEC) to a BCS superfluid (NSR1985, ; Hu2006, ; Diener2008, ), recently observed with fermionic ⁴⁰K and ⁶Li atoms (Regal2004, ; Zwierlein2004, ). Across the bosonic BEC-BCS crossover, the interspecies scattering length becomes positive and small, and the strongly interacting Bose droplet eventually disappears and turns into a gas of tightly-bound molecules or dimers via a first-order phase transition.

We then focus on the new cases of low-dimensional quantum droplets and examine systematically their bulk properties. Our pairing theory turns out to work extremely well in one dimension. For nearly all the interaction strengths at which quantum droplets exist, we find an excellent agreement between the theory and the DMC simulation for the ground-state energy (Parisi2019, ). In two dimensions, quantum droplets emerges for an arbitrarily weak interspecies interaction strength. Due to the weakness of the interspecies attraction, our pairing theory does not differ too much with the Petrov theory (Petrov2016, ). Both theories fail to have a good agreement with the DMC simulation, presumably due to the beyond-LHY-correction that becomes increasingly important in two dimensions (Mora2009, ; Astrakharchik2009, ). Nevertheless, it is remarkable that with increasing interspecies attraction our pairing theory predict a critical attraction, above which quantum droplets cease to exist. This critical value is consistent with the threshold for zero-crossing in dimer-dimer scatterings from a four-body problem in two dimensions (Guijarro2020, ). Above the threshold, the effective interaction between dimers (i.e., tightly bound bosonic pairs) changes from weakly attractive to weakly repulsive, indicating the instability of quantum droplets in the few-body limit.

The rest of the paper is organized as follows. In the next section (Sec. II), we introduce the model Hamiltonian for two-component Bose-Bose mixtures and present the pairing theory, i.e, the Bogoliubov theory with bosonic pairing. We mention briefly the concept of a bosonic BEC-BCS crossover. In Sec. III, the connection of our pairing theory to the conventional Bogoliubov theory without pairing is discussed. In Sec. IV, Sec. V, and Sec. VI, we consider the three-, one-, and two-dimensional cases, respectively. In three dimensions, we show the existence of a strongly interacting Bose droplet and map out the phase diagram of the bosonic BEC-BCS crossover. In the low-dimensional cases, we discuss in detail the comparison of the pairing theory to the benchmark DMC simulations and the physics near the threshold, at which the Bose droplet starts to disappear. Finally, Sec. VII is devoted to conclusions and outlooks.

We consider a two-component Bose-Bose mixture as in the seminal work by Petrov (Petrov2015, ). To be specific, let us focus on homonuclear mixtures such as the ³⁹K-³⁹K mixture, with which the masses of the two components are the same, i.e., $m_{1}=m_{2}=m$. In the presence of the intraspecies interactions $g_{11}$ and $g_{22}$, and the interspecies interactions $g_{12}=g_{21}$, the system in free space can be described the following Hamiltonian density,

	$\displaystyle\mathscr{H}\left(\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle\mathscr{H}_{0}\left(\mathbf{x}\right)+\mathscr{H}_{\textrm{int}}\left(\mathbf{x}\right),$		(1)
	$\displaystyle\mathscr{H}_{0}\left(\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle\sum_{i=1,2}\phi_{i}^{\dagger}\left(\mathbf{x}\right)\left[-\frac{\hbar^{2}\nabla^{2}}{2m}-\mu_{i}\right]\phi_{i}\left(\mathbf{x}\right),$		(2)
	$\displaystyle\mathscr{H}_{\textrm{int}}\left(\mathbf{x}\right)$	$\displaystyle=$	$\displaystyle\sum_{i,j=1,2}\frac{g_{ij}}{2}\phi_{i}^{\dagger}\left(\mathbf{x}\right)\phi_{j}^{\dagger}\left(\mathbf{x}\right)\phi_{j}\left(\mathbf{x}\right)\phi_{i}\left(\mathbf{x}\right),$		(3)

where $\phi_{i}(\mathbf{x})$ ($i=1,2$) is the annihilation field operator of the $i$-species bosons and $\mu_{i}$ is the chemical potential. In three or two dimensions, the use of contact inter-particle interactions leads to the well-known ultraviolet divergence, so the bare interaction strengths $g_{ij}$ need regularization and are to be replaced by the $s$-wave scattering lengths $a_{ij}$ or the binding energies $E_{B}^{ij}$. For example, in three dimensions we may write,

$$\frac{1}{g_{ij}}=\frac{m}{4\pi\hbar^{2}a_{ij}}-\frac{1}{\mathcal{V}}\sum_{k}\frac{m}{\hbar^{2}\mathbf{k}^{2}},$$

(4)

where the volume $\mathcal{V}$ (or area in two dimensions and length in one dimension) will be set to unity hereafter.

We use the conventional path-integral formalism to describe our bosonic pairing theory, following its fermionic counterpart (He2015, ; Hu2020b, ). We are interested in calculating the thermodynamic potential $\varOmega$ from the partition function,

$$\mathcal{Z}=\int\mathcal{D}[\phi_{1},\phi_{2}]e^{-\mathcal{S}},$$

(5)

where the action is given by,

$$\mathcal{S}=\int dx\left[\sum_{i=1,2}\bar{\phi}_{i}\left(x\right)\partial_{\tau}\phi_{i}\left(x\right)+\mathscr{H}\left(x\right)\right].$$

(6)

Here, we have used the standard notations $x\equiv(\mathbf{x},\tau)$ and $\int dx\equiv\int d\mathbf{x}\int_{0}^{\beta}d\tau$, and $\beta\equiv 1/(k_{B}T)$.

Due to the attractive interspecies interaction (i.e, $g_{12}<0$), we anticipate the pairing between different species. To make it evident, we explicitly introduce an auxiliary pairing field $\Delta(x)$ and take the Hubbard–Stratonovich transformation to decouple the Hamiltonian density for interspecies interactions, i.e.,

$$\exp\left[-g_{12}\int dx\bar{\phi}_{1}\bar{\phi}_{2}\phi_{2}\phi_{1}\right]=\int\mathcal{D}\left[\Delta\left(x\right)\right]\exp\left\{\int dx\left[\frac{\left|\Delta\left(x\right)\right|^{2}}{g_{12}}+\left(\bar{\Delta}\phi_{2}\phi_{1}+\bar{\phi}_{1}\bar{\phi}_{2}\Delta\right)\right]\right\}.$$

(7)

The action then becomes,

$$\mathcal{S}=\int dx\left[-\frac{\left|\Delta\left(x\right)\right|^{2}}{g_{12}}-\left(\bar{\Delta}\phi_{2}\phi_{1}+\bar{\phi}_{1}\bar{\phi}_{2}\Delta\right)\right]+\sum_{i=1,2}\int dx\left[\bar{\phi}_{i}\left(\partial_{\tau}-\frac{\hbar^{2}\nabla^{2}}{2m}-\mu_{i}\right)\phi_{i}+\frac{g_{ii}}{2}\bar{\phi}_{i}^{2}\phi_{i}^{2}\right].$$

(8)

For the pairing field $\Delta(x)$, it suffices to take a uniform saddle-point solution $\Delta(x)=\Delta>0$. At the same level of the Bogoliubov approximation, at zero temperature we assume the bosons condensate into the zero-momentum states, i.e.,

$\displaystyle\phi_{i}\left(x\right)$

$\displaystyle=$

$\displaystyle\phi_{ic}+\delta\phi_{i}\left(x\right),$

(9)

with a positive wave-function $\phi_{ic}>0$. Following the Bogoliubov approximation, the intraspecies interaction terms may be approximated as,

$$\frac{g_{ii}}{2}\bar{\phi}_{i}^{2}\phi_{i}^{2}\simeq\frac{C_{i}^{2}}{2g_{ii}}+\frac{C_{i}}{2}\left(4\delta\bar{\phi}_{i}\delta\phi_{i}+\delta\bar{\phi}_{i}\delta\bar{\phi}_{i}+\delta\phi_{i}\delta\phi_{i}\right),$$

(10)

where $C_{i}=g_{ii}\phi_{ic}^{2}$. As a consequence, we find the effective action $\mathcal{S}\simeq\beta\mathcal{V}\varOmega_{0}+\mathcal{S}_{B}$, where the condensate thermodynamic potential $\varOmega_{0}$ is given by,

$$\varOmega_{0}=\sum_{i=1,2}\left(-\mu_{i}\phi_{ic}^{2}+\frac{C_{i}^{2}}{2g_{ii}}\right)-\frac{\Delta^{2}}{g_{12}}-2\Delta\phi_{1c}\phi_{2c},$$

(11)

and the quantum fluctuations around the condensates have the contribution,

	$\displaystyle\mathcal{S}_{B}$	$\displaystyle=-\int dx\Delta\left(\delta\bar{\phi}_{1}\delta\bar{\phi}_{2}+\delta\phi_{2}\delta\phi_{1}\right)+\int dx\sum_{i=1,2}$
		$\displaystyle\left[\delta\bar{\phi}_{i}\left(\partial_{\tau}+\hat{B}_{i}\right)\delta\phi_{i}+\frac{C_{i}}{2}\left(\delta\bar{\phi}_{i}\delta\bar{\phi}_{i}+\delta\phi_{i}\delta\phi_{i}\right)\right],$		(12)

with $\hat{B}_{i}(x)\equiv-\hbar^{2}\nabla^{2}/(2m)-\mu_{i}+2C_{i}$. By introducing a Nambu spinor $\Phi(x)=[\delta\phi_{1}(x),\delta\bar{\phi}_{1}(x),\delta\phi_{2}(x),\delta\bar{\phi}_{2}(x)]^{T}$, we may recast $\mathcal{S}_{B}$ into a compact form,

$$\mathcal{S}_{B}=\int dxdx^{\prime}\bar{\Phi}\left(x\right)\left[-\mathscr{D}^{-1}\left(x,x^{\prime}\right)\right]\Phi\left(x^{\prime}\right),$$

(13)

where the inverse Green function of bosons is given by,

$$\mathscr{D}^{-1}=\left[\begin{array}[]{cccc}-\partial_{\tau}-\hat{B}_{1}&-C_{1}&0&\Delta\\ -C_{1}&\partial_{\tau}-\hat{B}_{1}&\Delta&0\\ 0&\Delta&-\partial_{\tau}-\hat{B}_{2}&-C_{2}\\ \Delta&0&-C_{2}&\partial_{\tau}-\hat{B}_{2}\end{array}\right].$$

(14)

Due to the delta function $\delta\left(x-x^{\prime}\right)$ in $\mathscr{D}^{-1}(x,x^{\prime})$, which we do not explicitly show in the above equation, it is convenient to work in momentum space by performing a Fourier transform. After replacing $-\partial_{\tau}$ with the bosonic Matasubara frequencies $i\omega_{m}$ (i.e., $\omega_{m}=2\pi mk_{B}T$ with $m\subseteq\mathbb{Z}$) and performing the analytic continuation $i\omega_{m}\rightarrow\omega+i0^{+}$, i.e.,

$$-\partial_{\tau}\rightarrow\omega+i0^{+},$$

(15)

and taking the replacement

$$\hat{B}_{i}\rightarrow B_{i\mathbf{k}}=\varepsilon_{\mathbf{k}}-\mu_{i}+2C_{i}$$

(16)

with $\varepsilon_{\mathbf{k}}=\hbar^{2}\mathbf{k}^{2}/(2m)$, it is straightforward to explicitly write down the expression of $\mathscr{D}^{-1}(\mathbf{k},\omega)$. By solving the poles of the bosonic Green function, i.e., $\det[\mathscr{D}^{-1}(\mathbf{k},\omega\rightarrow E(\mathbf{k}))]=0$, or more explicitly,

$$\omega^{4}-\omega^{2}\left[\left(B_{1\mathbf{k}}^{2}-C_{1}^{2}\right)+\left(B_{2\mathbf{k}}^{2}-C_{2}^{2}\right)-2\Delta^{2}\right]+\left[\left(B_{1\mathbf{k}}^{2}-C_{1}^{2}\right)\left(B_{2\mathbf{k}}^{2}-C_{2}^{2}\right)-2\left(B_{1\mathbf{k}}B_{2\mathbf{k}}+C_{1}C_{2}\right)\Delta^{2}+\Delta^{4}\right]=0,$$

(17)

we obtain the two Bogoliubov spectra,

$$E_{\pm}^{2}\left(\mathbf{k}\right)=\left[\mathcal{A}_{+}\left(\mathbf{k}\right)-\Delta^{2}\right]\pm\sqrt{\mathcal{A}_{-}^{2}\left(\mathbf{k}\right)+\Delta^{2}\left[\left(C_{1}+C_{2}\right)^{2}-\left(B_{1\mathbf{k}}-B_{2\mathbf{k}}\right)^{2}\right]},$$

(18)

where we have defined,

$$\mathcal{A}_{\pm}\left(\mathbf{k}\right)=\frac{\left(B_{1\mathbf{k}}^{2}-C_{1}^{2}\right)\pm\left(B_{2\mathbf{k}}^{2}-C_{2}^{2}\right)}{2}.$$

(19)

It is useful to explicitly compare the structure of our pairing theory with that of the widely-used Petrov theory. For this purpose, here we briefly review the Petrov’s prototype theory of quantum droplets. Within the Bogoliubov approximation (Larsen1963, ), we decouple the interspecies interaction Hamiltonian density ($\phi_{1c}=\phi_{2c}=\phi_{c}$),

	$\displaystyle g_{12}\bar{\phi}_{1}\bar{\phi}_{2}\phi_{2}\phi_{1}$	$\displaystyle\simeq$	$\displaystyle\frac{D^{2}}{g_{12}}-D\left(\delta\bar{\phi}_{1}\delta\phi_{1}+\delta\bar{\phi}_{2}\delta\phi_{2}\right)$		(33)
			$\displaystyle-D\left(\delta\bar{\phi}_{1}\delta\bar{\phi}_{2}+\delta\bar{\phi}_{1}\delta\phi_{2}+\textrm{H.c.}\right),$

where H.c. stands for taking the Hermitian conjugate and $D\equiv-g_{12}\phi_{c}^{2}>0$ seems to play the role of the pairing field $\Delta$ in our pairing theory. However, there is a slight difference. The Bogoliubov decoupling shown in the above generates two additional terms in the quantum fluctuation action $\mathcal{S}_{B}$, i.e., $-\int dxD(\delta\bar{\phi}_{1}\delta\phi_{1}+\delta\bar{\phi}_{2}\delta\phi_{2})$ and $-\int dxD(\delta\bar{\phi}_{1}\delta\phi_{2}+\delta\bar{\phi}_{2}\delta\phi_{1})$. In momentum space, the inverse Green function of bosons $\mathscr{D}^{-1}(\mathbf{k},\omega)$ then becomes,

$$\mathscr{D}^{-1}=\left[\begin{array}[]{cccc}\omega-B_{\mathbf{k}}&-C&D&D\\ -C&-\omega-B_{\mathbf{k}}&D&D\\ D&D&\omega-B_{\mathbf{k}}&-C\\ D&D&-C&-\omega-B_{\mathbf{k}}\end{array}\right],$$

(34)

where $B_{\mathbf{k}}\equiv\varepsilon_{k}-\mu+2C-D$. The existence of the two additional terms leads to the two Bogoliubov spectra,

$$\tilde{E}_{\pm}\left(\mathbf{k}\right)=\sqrt{\left(B_{\mathbf{k}}-C\right)\left(B_{\mathbf{k}}+C\mp 2D\right)}.$$

(35)

Here, we have used the tilde to distinguish the dispersion relations from those of the pairing theory. In this case, it is easy to see that the condensate thermodynamic potential takes the form,

$$\varOmega_{0}=-2\mu\phi_{c}^{2}+g\phi_{c}^{4}+g_{12}\phi_{c}^{4}.$$

(36)

By minimizing $\varOmega_{0}$ with respect to $\phi_{c}^{2}$, we obtain the restriction,

$$\mu=g\phi_{c}^{2}+g_{12}\phi_{c}^{2}=C-D,$$

(37)

or equivalently $B_{\mathbf{k}=0}=C$, which ensures the gapless Bogoliubov spectra. Therefore, we may rewrite down the condensate thermodynamic potential,

$$\varOmega_{0}=-\frac{C^{2}}{g}-\frac{D^{2}}{g_{12}},$$

(38)

the two dispersion relations,

	$\displaystyle\tilde{E}_{-}(\mathbf{k})$	$\displaystyle=$	$\displaystyle\sqrt{\varepsilon_{\mathbf{k}}\left(\varepsilon_{\mathbf{k}}+2C+2D\right)},$		(39)
	$\displaystyle\tilde{E}_{+}(\mathbf{k})$	$\displaystyle=$	$\displaystyle\sqrt{\varepsilon_{\mathbf{k}}\left(\varepsilon_{\mathbf{k}}+2C-2D\right)},$		(40)

and also the LHY thermodynamic potential,

$$\varOmega_{\textrm{LHY}}=\frac{1}{2}\sum_{\mathbf{k}}\left[\tilde{E}_{+}\left(\mathbf{k}\right)+\tilde{E}_{-}\left(\mathbf{k}\right)-2\left(\varepsilon_{\mathbf{k}}+C\right)\right].$$

(41)

In three dimensions, using Eq. (4) we replace the bare interaction strengths $g$ and $g_{12}$ with the $s$-wave scattering lengths $a$ and $a_{12}$, respectively. Therefore, we obtain $\varOmega_{\textrm{3D}}=\varOmega_{0}+\varOmega_{\textrm{LHY}}$,

	$\displaystyle\varOmega_{\textrm{3D}}$	$\displaystyle=$	$\displaystyle-\frac{m}{4\pi\hbar^{2}}\left[\frac{C^{2}}{a}+\frac{D^{2}}{a_{12}}\right]+\frac{1}{2}\sum_{\mathbf{k}}\left[\tilde{E}_{+}\left(\mathbf{k}\right)+\tilde{E}_{-}\left(\mathbf{k}\right)\right.$		(42)
			$\displaystyle\left.-2\left(\varepsilon_{\mathbf{k}}+C\right)+\frac{2\left(C^{2}+D^{2}\right)}{\hbar^{2}\mathbf{k}^{2}/m}\right].$

The integration over the momentum can be easily calculated, by using the identity

$$\sum_{\mathbf{k}}\left[\sqrt{\varepsilon_{\mathbf{k}}\left(\varepsilon_{\mathbf{k}}+\alpha\right)}-\varepsilon_{\mathbf{k}}-\frac{\alpha}{2}+\frac{\alpha^{2}}{8\varepsilon_{\mathbf{k}}}\right]=\frac{\left(2m\right)^{3/2}\alpha^{5/2}}{15\pi^{2}\hbar^{3}}$$

(43)

in three dimensions. We arrive at,

$$\varOmega_{\textrm{3D}}=-\frac{m}{4\pi\hbar^{2}}\left[\frac{C^{2}}{a}+\frac{D^{2}}{a_{12}}\right]+\frac{8m^{3/2}C^{5/2}}{15\pi^{2}\hbar^{3}}\mathcal{F}_{3}\left(\frac{D}{C}\right),$$

(44)

where $\mathcal{F}_{3}(\alpha)\equiv(1+\alpha)^{5/2}+(1-\alpha)^{5/2}$. To calculate the total energy, we note that unlike the pairing gap $\Delta$ in our pairing theory, the variable $D$ is not a variational parameter. Therefore, there is an ambiguity to determine the variables $D$ and then $C=\mu+D$ for a given chemical potential $\mu$. Nevertheless, we may assume that

$$\frac{D}{C}=-\frac{g_{12}}{g}\simeq-\frac{a_{12}}{a},$$

(45)

so that $C=\mu a/(a+a_{12})$ and $D=-\mu a_{12}/(a+a_{12})$. As quantum droplets emerge when $D/C=-a_{12}/a>1$ in three dimensions (Petrov2015, ), we immediately find that the Bogoliubov spectrum $\tilde{E}_{+}(\mathbf{k})=\sqrt{\varepsilon_{\mathbf{k}}(\varepsilon_{\mathbf{k}}+2C-2D})$ becomes complex and consequently the function $\mathcal{F}_{3}(\alpha)$ in Eq. (44) is ill-defined (Petrov2015, ; Cikojevic2019, ; Hu2020a, ). To cure this problem, we may simply set $\alpha=D/C=1$ in the function $\mathcal{F}_{3}(\alpha)$ (Petrov2015, ) and hence the LHY term become independent on the interspecies interaction strength. This Petrov’s prescription is now widely taken in the theoretical studies of quantum droplets. We note also that, to calculate the total energy, we may further approximate $C=(2\pi\hbar^{2}a/m)n$ and $\mu=[2\pi\hbar^{2}(a+a_{12})/m]n$, which leads to the total energy,

$$\frac{E_{\textrm{3D}}}{N}=\frac{\pi\hbar^{2}\left(a+a_{12}\right)}{m}n+\frac{256\sqrt{\pi}}{15}\frac{\hbar^{2}a^{5/2}}{m}n^{3/2}.$$

(46)

In other words, to determine the density $n$, we take the derivative of the first term in Eq. (44) only with respect to the chemical potential $\mu$. This approximation is well justified for a conventional weakly-interacting Bose gas in three dimensions (Salasnich2016, ). However, it may not be convincing for quantum droplets, where the second term in Eq. (44) may become comparable to the first term.

Let us now consider the pairing theory in three dimensions, showing some technical details behind our previous work (Hu2020a, ) and describing the strongly interacting Bose droplet at the BEC-BCS crossover. By regularizing the bare interaction strengths $g$ and $g_{12}$ in terms of the $s$-wave scattering lengths $a$ and $a_{12}$, we rewrite the thermodynamic potential Eq. (30) into the form,

	$\displaystyle\varOmega_{\textrm{3D}}$	$\displaystyle=$	$\displaystyle-\frac{m}{4\pi\hbar^{2}}\left[\frac{C^{2}}{a}+\frac{\Delta^{2}}{a_{12}}\right]+\frac{1}{2}\left(\mathcal{I}_{+}+\mathcal{I}_{-}\right),$		(47)
	$\displaystyle\mathcal{I}_{\pm}$	$\displaystyle=$	$\displaystyle\sum_{\mathbf{k}}\left[E_{\pm}\left(\mathbf{k}\right)-\left(\varepsilon_{\mathbf{k}}+C+\Delta\right)+\frac{\left(C\pm D\right)^{2}}{2\varepsilon_{\mathbf{k}}}\right].$		(48)

$\mathcal{I}_{-}$ can be directly calculated, with the help of the identity Eq. (43),

$$\mathcal{I}_{-}=\frac{16m^{3/2}}{15\pi^{2}\hbar^{3}}C^{5/2}\left(1+\frac{\Delta}{C}\right)^{5/2}.$$

(49)

To calculate $\mathcal{I}_{+}$, we introduce a new variable $t\equiv[\hbar^{2}k^{2}/(2m)]/(2C)$ and $\alpha\equiv\Delta/C$ and rewrite

$$\mathcal{I}_{+}=\frac{16m^{3/2}}{15\pi^{2}\hbar^{3}}C^{5/2}h_{3}\left(\alpha\right),$$

(50)

where the function

	$\displaystyle h_{3}\left(\alpha\right)$	$\displaystyle\equiv$	$\displaystyle\frac{15}{4}\int_{0}^{\infty}dt\sqrt{t}\left[\sqrt{\left(t+1\right)\left(t+\alpha\right)}-t\right.$		(51)
			$\displaystyle\left.-\frac{1+\alpha}{2}+\frac{\left(1-\alpha\right)^{2}}{8t}\right].$

By combining $\mathcal{I}_{+}$ and $\mathcal{I}_{-}$, we obtain the (regularized) LHY thermodynamic potential ($C=\mu+\Delta$),

$$\varOmega_{\textrm{LHY}}=\frac{8m^{3/2}}{15\pi^{2}\hbar^{3}}\left(\mu+\Delta\right)^{5/2}\mathcal{G}_{3}\left(\frac{\Delta}{\mu+\Delta}\right),$$

(52)

where $\mathcal{G}_{3}(\alpha)\equiv(1+\alpha)^{5/2}+h_{3}(\alpha)$. Compared with the function $\mathcal{F}_{3}(\alpha)\equiv(1+\alpha)^{5/2}+(1-\alpha)^{5/2}$ in the last section, we find interestingly that the role of $(1-\alpha)^{5/2}$, which is not well-defined for $\alpha>1$, is now taken by the new function $h_{3}(\alpha)$. Therefore, we obtain the total thermodynamic potential,

	$\displaystyle\varOmega_{\textrm{3D}}$	$\displaystyle=$	$\displaystyle-\frac{m}{4\pi\hbar^{2}}\left[\frac{\left(\mu+\Delta\right)^{2}}{a}+\frac{\Delta^{2}}{a_{12}}\right]$		(53)
			$\displaystyle+\frac{8m^{3/2}}{15\pi^{2}\hbar^{3}}\left(\mu+\Delta\right)^{5/2}\mathcal{G}_{3}\left(\frac{\Delta}{\mu+\Delta}\right).$

It takes essentially the same form as the thermodynamic potential Eq. (44) in the standard Bogoliubov theory, except that the ill-defined function $\mathcal{F}_{3}(\alpha)$ is now replaced by $\mathcal{G}_{3}(\alpha)$, and the pairing gap $\Delta$ is variational and should be determined by minimizing $\varOmega_{\textrm{3D}}(\Delta)$. For a given chemical potential $\mu$, we therefore minimize $\varOmega_{\textrm{3D}}$ to find the saddle-point value of the pairing order parameter $\Delta_{0}$. For nonzero $\Delta_{0}\neq 0$, we obtain $\varOmega_{\textrm{3D}}(\mu,\Delta_{0})$ and calculate $n=-\partial\varOmega_{\textrm{3D}}(\mu,\Delta_{0})/\partial\mu$.

In the weakly interacting regime, where $a+a_{12}\simeq 0$, we typically find that the chemical potential is much smaller than either the parameter $C$ or the pairing gap $\Delta_{0}$ (Hu2020a, ). This could be easily understood from the $\Delta$-dependence of $\varOmega_{0}$ and $\varOmega_{\textrm{LHY}}$, as shown in Eq. (53). We note that two terms in $\varOmega_{0}$ are large and have opposite sign. Each of them (i.e., absolute value) is much larger than $\varOmega_{\textrm{LHY}}$. Therefore, when we minimize $\varOmega$ with respect to $\Delta$, we only need to minimize $\Omega_{0}$. This leads to the condition,

$$\frac{\mu+\Delta_{0}}{a}+\frac{\Delta_{0}}{a_{12}}\simeq 0.$$

(54)

Hence, as a result of $a_{12}\sim-a$, we obtain

$$\mu\simeq-\left(1+\frac{a}{a_{12}}\right)\Delta_{0}\ll\Delta_{0},C.$$

(55)

Due to the smallness of $\left|\mu\right|$, it is reasonable to neglect the $\mu$-dependence in $\varOmega_{\textrm{LHY}}$ and also the term $\mu^{2}$ in $\varOmega_{0}$. Therefore, we obtain,

$$\varOmega_{\textrm{3D}}\simeq-\frac{m}{4\pi\hbar^{2}}\left[\frac{2\mu\Delta+\Delta^{2}}{a}+\frac{\Delta^{2}}{a_{12}}\right]+\frac{32\sqrt{2}m^{3/2}}{15\pi^{2}\hbar^{3}}\Delta^{5/2}.$$

(56)

By taking the derivative with respect to $\mu$ and taking the saddle-point value $\Delta=\Delta_{0}$, we find

$$n=-\frac{\partial\varOmega_{\textrm{3D}}}{\partial\mu}\simeq\frac{m}{2\pi\hbar^{2}a}\Delta_{0}.$$

(57)

Replacing $\Delta_{0}$ by the density $n$ everywhere in $\varOmega_{\textrm{3D}}$ and calculate $E_{\textrm{3D}}=\varOmega_{\textrm{3D}}+\mu n$, we finally arrive at an approximate energy for small densities,

$$\frac{E_{\textrm{3D}}}{N}=-\frac{\pi\hbar^{2}}{m}\left(a+\frac{a^{2}}{a_{12}}\right)n+\frac{256\sqrt{\pi}}{15}\frac{\hbar^{2}a^{5/2}}{m}n^{3/2}.$$

(58)

It is useful to compare this analytic expression with the energy functional obtained by Petrov using his prescription (Petrov2015, ), i.e., Eq. (46). It is interesting to see that these two energy functionals have the exactly same LHY term. The reproduce of the Petrov’s approximate LHY term in our pairing theory suggests that Petrov’s prescription is actually very reasonable, at least for the case of equal intraspecies interactions considered here (Hu2020a, ). However, it is worth emphasizing that, quite unexpectedly, the mean-field energy term in our pairing theory (i.e., the first term in Eq. (58)) changes a lot. It is weakened by a factor of $-a/a_{12}<1$, compared with the conventional mean-field expression used by Petrov (Petrov2015, ), i.e.,

$$\left(g+g_{12}\right)\frac{n}{4}\rightarrow\frac{\pi\hbar^{2}(a+a_{12})}{m}n.$$

(59)

This difference partly comes from our regularization of the bare interaction strengths, which is rigorously treated in the pairing theory. As the beyond-mean-field LHY effect becomes dominant in quantum droplets, a consistent treatment of the potential regularization is necessary. Therefore, it is not a surprise why our regularized mean-field energy becomes different from the widely-accepted conventional expression.

Let us now consider the strongly interacting regime, where interspecies scattering length $\left|a_{12}\right|$ is significantly larger than the intraspecies scattering length $a$. In Fig. 1, we report the density dependence of the energy per particle $E/N$, when the interspecies attraction increases from $a_{12}=-1.5a$ to $a_{12}=-\infty$, crosses the Feshbach resonance (i.e., $a_{12}$ jumps from $-\infty$ to $+\infty$), and further increases towards the formation of tightly bound molecules ($a_{12}\rightarrow 0^{+}$). Remarkably, with increasing interspecies attraction, the energy per particle decreases steadily from the weakly interacting limit (i.e., $a_{12}<-1.5a$, see the top left of the figure, which has been studied in Ref. (Hu2020a, )) to the strongly interacting regime ($a_{12}\rightarrow\pm\infty$). This is very similar to what happens at the fermionic BEC-BCS crossover (NSR1985, ; Hu2006, ; Diener2008, ), which is realized by tuning the attractive interaction between two unlike fermions across a Feshbach resonance (Regal2004, ; Zwierlein2004, ). Here, we always find a minimum in the energy per particle (see the dashed line), until the interspecies scattering length becomes positive and smaller than a threshold $a_{12,\textrm{crit}}\simeq 3.6a$. The existence of a minimum in the energy per particle, i.e., $\partial(E/N)/\partial n=(\mu-E/N)/n=0$, precisely implies the appearance of a self-bound Bose droplet in a vacuum, which has a zero pressure $P=n\mu-E/\mathcal{V}=0$. Therefore, we observe that a Bose droplet may exist in the strongly interacting regime, where the interspecies scattering length diverges. Importantly, the equilibrium density of the droplet, at which the minimum energy per particle is reached, does not significantly increase as the interspecies scattering length diverges. We find that the equilibrium density attains its largest value near the Feshbach resonance and is about a few times larger than the equilibrium density in the weakly coupling regime (i.e., at $a_{12}=-1.5a$).

The experimental confirmation of a strongly interacting Bose droplet would be a highly non-trivial task. In a single-component Bose gas, it is well-known that a strongly interacting Bose gas is difficult to reach due to severe three-body loss. Indeed, in the recent experiments the lifetime of a ³⁹K-³⁹K Bose droplet is limited to about ten milliseconds, due to the intraspecies three-body loss in a particular component (i.e., the $\left|F=1,m_{F}=0\right\rangle$ hyperfine state) (Cabrera2018, ; Semeghini2018, ). This is not a serious problem, since we can choose a binary Bose mixture with small intraspecies three-body loss rate, for example, the heteronuclear ⁴¹K-⁸⁷Rb mixture, where the lifetime of the droplet was found to be much longer (DErrico2019, ). As the equilibrium density of a strongly interacting Bose droplet is at the same order as a weakly interacting droplet in magnitude, a large intraspecies three-body loss could be avoided. The fundamental difficulty then may come from the interspecies three-body loss, which remains unknown at the moment. Theoretically, the solution of the three-body problem of two like bosons and one unlike boson and the calculation of the related interspecies three-body loss rate would be an interesting research topic to consider in future works.

On the BEC side of the Feshbach resonance, the Bose droplet ceases to exist below the threshold $a_{12,\textrm{crit}}\simeq 3.6a$, as highlighted in Fig. 2. This is closely related to the formation of a gas of dimers towards the BEC limit, which also happens at the fermionic BEC-BCS crossover (NSR1985, ; Hu2006, ; Diener2008, ). To show this, in Fig. 3(a) we present the minimum energy per particle $E_{\min}/N$ and the half of the bound-state energy $\epsilon_{B}/2$ of a dimer, as a function of the ratio $a/a_{12}$. It is readily seen that a gas of non-interacting dimers becomes energetically favorable over a Bose droplet once we pass the crossing point of the two energy curves at $a/a_{12}\simeq 0.23$ or $a_{12}\simeq 4.35a$. Interestingly, the pairing order parameter $\Delta_{0}$ appears to have a maximum around the crossing point, as shown in Fig. 3(b).

The above observation leads us to sketch a general phase diagram for a binary Bose mixture across the Feshbach resonance for the attractive interspecies interactions, as illustrated in Fig. 4. When we increase the interspecies attraction from $a_{12}=0^{-}$ (i.e., the far-left part of the figure), the mixture is initially in the miscible phase. Upon reaching the mean-field collapsing point at $a_{12}=-a$, it turns into a weakly interacting Bose droplet as experimentally observed (Cabrera2018, ; Semeghini2018, ). Across the Feshbach resonance ($a/a_{12}=0$), the Bose droplet becomes strongly interacting with a significant portion of Cooper pairs, whose size is comparable to the mean-distance between bosons. By further increasing the interspecies attraction to the threshold $a_{12,\textrm{crit}}\simeq 3.6a$, the Bose droplet disappears and changes into a gas of dimers via a first-order transition.

It is worth noting that, since the intraspecies scattering length $a$ is small, the interspecies scattering length $a_{12}$ at the threshold would also be small and positive. This means that the gas of dimers would be extremely difficult to reach adiabatically in the time scale of actual experiments, due to its deep energy level. Instead, we would observe a metastable state of the mixture, which is an immiscible phase (i.e., a phase-separation phase where the two components do not overlap in real space) at the threshold $a_{12,\textrm{crit}}\simeq 3.6a$. By further increasing the interspecies attraction towards a vanishing scattering length ($a_{12}\rightarrow 0^{+}$), the mixture will again enter a miscible phase at $a_{12}=a$ and connect smoothly to the miscible phase at $a_{12}=0^{-}$.

We now turn to discuss low-dimensional quantum droplets, starting from the one-dimensional case. In one dimension, the contact interaction is well defined and does not need regularization. The interaction strength can be characterized by using the dimensionless interaction parameter, such as $\gamma=mg/(n\hbar^{2})=-2/(na)$ and $\eta=mg_{12}/(n\hbar^{2})=-2/(na_{12})$, where

	$\displaystyle a$	$\displaystyle=$	$\displaystyle-\frac{2\hbar^{2}}{mg}<0$		(60)
	$\displaystyle a_{12}$	$\displaystyle=$	$\displaystyle-\frac{2\hbar^{2}}{mg_{12}}>0$		(61)

are the one-dimensional $s$-wave scattering lengths. Following Ref. (Parisi2019, ), we choose the binding energy of a dimer of two bosons in different species, i.e.,

$$\varepsilon_{B}=\frac{\hbar^{2}}{ma_{12}^{2}},$$

(62)

as the units of energy, and the inverse scattering length $\left|a\right|^{-1}$ as the units of density.

In the Bogoliubov theory, the LHY thermodynamic potential Eq. (41) in one dimension is easy to calculate. By performing the one-dimensional integration,

$$\sum_{\mathbf{k}}\left[\sqrt{\varepsilon_{\mathbf{k}}\left(\varepsilon_{\mathbf{k}}+\alpha\right)}-\varepsilon_{\mathbf{k}}-\frac{\alpha}{2}\right]=-\frac{\left(2m\right)^{1/2}\alpha^{3/2}}{3\pi\hbar},$$

(63)

we obtain

$$\varOmega_{\textrm{LHY}}=-\frac{2m^{1/2}}{3\pi\hbar}\left[\left(C+D\right)^{3/2}+\left(C-D\right)^{3/2}\right].$$

(64)

By taking $C=gn/2$ and $D=g_{12}n/2$ as before in $\varOmega_{\textrm{LHY}}$ and adding the mean-field energy $(g+g_{12})n^{2}/4$, we obtain the energy per particle predicted by the Bogoliubov theory,

$$\frac{E_{\textrm{1D}}}{N}=\left(g+g_{12}\right)\frac{n}{4}-\frac{\left(2m\right)^{1/2}}{6\pi\hbar}g^{3/2}\mathcal{F}_{1}\left(\frac{g_{12}}{g}\right)n^{1/2},$$

(65)

where $\mathcal{F}_{1}(\alpha)\equiv(1+\alpha)^{3/2}+(1-\alpha)^{3/2}$. It is interesting to note that, in one dimension the LHY energy functional is negative so the force provided by quantum fluctuations is attractive. It is to be balanced by the repulsive mean-field force at $g>-g_{12}>0$. Therefore, somehow counterintuitively, the formation of quantum droplets is driven by quantum fluctuations (Petrov2016, ). As the mean-field solution is stable, the energy in Eq. (65) does not suffer from the issue of complex number as we encounter earlier in three dimensions. Nevertheless, following Ref. (Petrov2016, ) we may still use the Petrov’s prescription and take $g_{12}=-g$ in Eq. (65) to define an energy per particle,

$$\frac{E_{\textrm{1D}}}{N}=\left(g+g_{12}\right)\frac{n}{4}-\frac{2m^{1/2}}{3\pi\hbar}g^{3/2}n^{1/2}.$$

(66)

In the pairing theory, we calculate the thermodynamic potential Eq. (30) in one dimension. As in the three-dimensional case, we similarly separate the LHY thermodynamic potential into two parts, $\mathcal{I}_{-}$ and $\mathcal{I}_{+}$. It is straightforward to obtain $\mathcal{I}_{-}$ with the help of the identity Eq. (63) and we obtain,

$$\mathcal{I}_{-}=-\frac{4m^{1/2}}{3\pi\hbar}C^{3/2}\left(1+\frac{\Delta}{C}\right)^{3/2}.$$

(67)

For $\mathcal{I}_{+}$, we instead find

$$\mathcal{I}_{+}=-\frac{4m^{1/2}}{3\pi\hbar}C^{3/2}h_{1}\left(\alpha\right),$$

(68)

where the function $h_{1}(\alpha)$ is given by,

$$h_{1}\equiv 3\int_{0}^{\infty}dt\left[t^{2}+\frac{1+\alpha}{2}-\sqrt{\left(t^{2}+1\right)\left(t^{2}+\alpha\right)}\right],$$

(69)

and is plotted in Fig. 5. Therefore, we obtain the total thermodynamic potential ($C=\mu+\Delta$),

$$\varOmega_{\textrm{1D}}=-\frac{C^{2}}{g}-\frac{\Delta^{2}}{g_{12}}-\frac{2m^{1/2}}{3\pi\hbar}C^{3/2}\mathcal{G}_{1}\left(\frac{\Delta}{C}\right),$$

(70)

where $\mathcal{G}_{1}(\alpha)\equiv(1+\alpha)^{3/2}+h_{1}(\alpha)$.