Quantum Phases of Self-Bound Droplets of Bose-Bose Mixtures

Introduction.—Quantum droplets in both dipolar Kadau et al. (2016); Schmitt et al. (2016); Chomaz et al. (2016); Böttcher et al. (2019a) and binary condensates Cabrera et al. (2018); Cheiney et al. (2018); Semeghini et al. (2018); D’Errico et al. (2019); Burchianti et al. (2020) represent a new state of matter emerging in the mean-field unstable regime. They have attracted great interests in studying their fascinating properties, such as self-bound feature Schmitt et al. (2016); Chomaz et al. (2016); Cabrera et al. (2018); Semeghini et al. (2018), collective excitations Natale et al. (2019); Tanzi et al. (2019a), soliton to droplet transition Cheiney et al. (2018), droplet-droplet collision Ferioli et al. (2019), supersolid states Böttcher et al. (2019b); Tanzi et al. (2019b); Chomaz et al. (2019), Goldstone mode Guo et al. (2019), and so forth (see recent reviews Luo et al. (2021); Böettcher et al. (2021); Guo and Pfau (2021) and references therein). A widely adopted treatment for droplet states is the extend Gross-Pitaevskii equation (EGPE) which perturbatively incorporates quantum fluctuation into the Gross-Pitaevskii equation in terms of the Lee-Huang-Yang correction Lee et al. (1957); Schützhold et al. (2006); Lima and Pelster (2011, 2012); Petrov (2015). Although EGPE has provided satisfactory explanations to experimental observations, there are still discrepancies with experimental measurements Cabrera et al. (2018); Böttcher et al. (2019a); Böettcher et al. (2021).

Theoretical methods beyond EGPE have also been employed to study different aspects of the droplet states Saito (2016); Macia et al. (2016); Cinti et al. (2017); Bombin et al. (2017); Cikojević et al. (2019); Ota and Astrakharchik (2020); Boudjemâa and Guebli (2020); Hu et al. (2020); Hu and Liu (2020); Gu and Yin (2020). Among them, we studied the dipolar droplets using the Gaussian-state theory (GST) in which quantum fluctuation is included in a self-consistent manner Wang et al. (2020). Compared to other theoretical approaches beyond the EGPE theory, GST is computationally efficient and can be applied to realistic systems with large number of atoms. Interestingly, we found that, besides coherent state, dipolar droplets also contained two new macroscopic squeezed phases which were characterized by large second-order correlation and asymmetric atom-number distribution (AND), in striking contrast with those of the coherent-state-based droplets Wang et al. (2020); Shi et al. (2019). Whereas the asymmetric AND was experimentally observed Schmitt et al. (2016); Böttcher et al. (2019a), verifying the predication on the second-order correlation function requires further experiments beyond simple density measurements.

In this Letter, we study the quantum phases of self-bound droplets in quasi-two-dimensional (quasi-2D) binary condensates. Although this system has same quantum phases as those in dipolar droplets, here we discover three readily accessible signatures for the determination of the quantum states and the stablization mechanisms. Specifically, due to the multiple quantum phases containing in the droplet states, the radial size ($\sigma$) versus atom number ($N$) curve is of W shape, as opposed to the V-shaped curve expected for single quantum phase. Therefore, the observation of double dips on the $\sigma$-$N$ curve may rule out the pure coherent explanation of the droplet state. Moreover, the critical atom number (CAN), i.e., the minimal number of atoms required for forming a self-bound state, is determined by the quantum states of the condensate and its precise measurement of the CAN allows us to distinguish squeezed state from coherent one. Finally, the dip atom number (DAN), i.e., the atom number at the dip of the $\sigma$-$N$ curve, depends on both quantum phase and stabilization mechanism, which makes it a quantitative criterion for stability mechanism.

Formulation.—We consider a ultracold gas of $N$ ${}^{39}\text{K}$ atoms with $N_{\uparrow}$ atoms being in the hyperfine state $\left|\uparrow\right\rangle\equiv\left|F,m_{F}\right\rangle=\left|1,-1\right\rangle$ and $N_{\downarrow}$ in $\left|\downarrow\right\rangle\equiv\left|1,0\right\rangle$. The total Hamiltonian of the system,

consists of the single-, two-, and three-particle terms. In second-quantized form, the single-particle part reads

where $\hat{\psi}_{\alpha}\left(\mathbf{r}\right)$ ($\alpha=\uparrow,\downarrow$) are the field operators and $h_{\alpha}=-\hbar^{2}\nabla^{2}/(2M)-\mu_{\alpha}$ with $M$ being the mass of the atoms and $\mu_{\alpha}$ the chemical potential introduced to fixed the atom number in the $\alpha$th component. The two-body (2B) interaction Hamiltonian takes the form

where $g_{\alpha\beta}=4\pi\hbar^{2}a_{\alpha\beta}/M$ characterize the 2B interaction strengths with $a_{\alpha\beta}$ being scattering length between components $\alpha$ and $\beta$. The scattering lengths are tunable through Feshbach resonance and for scenario of interest to quantum droplets, the intra- and inter-species scattering lengths satisfy $a_{\uparrow\uparrow},a_{\downarrow\downarrow}>0$ and $a_{\uparrow\downarrow}<0$, respectively. Finally, the three-body (3B) interaction Hamiltonian can be expressed as

where $g_{3}$ is the 3B coupling constant which, for simplicity, is assumed to be independent of the spin components. Moreover, since we are only interested in the ground states of the system, we assume that $g_{3}$ is real and positive. The value of $g_{3}$ shall be determined by fitting the experimental data.

We study self-bound droplet states using the GST. Specifically, we assume that the many-body wave function of a binary gas adopts the variational ansatz Shi et al. (2018, 2019); Guaita et al. (2019)

where $\hat{\Psi}(\mathbf{r})\equiv\begin{pmatrix}\hat{\psi}(\mathbf{r})\\ \hat{\psi}^{\dagger}(\mathbf{r})\end{pmatrix}$ with $\hat{\psi}(\mathbf{r})=\begin{pmatrix}\hat{\psi}_{\uparrow}(\mathbf{r})\\ \hat{\psi}_{\downarrow}(\mathbf{r})\end{pmatrix}$ is the field operator expressed in the Nambu basis,

with $I_{2}$ being a $2\times 2$ identity matrix in the spin space, $\Phi(\mathbf{r})=\langle\hat{\Psi}(\mathbf{r})\rangle=\begin{pmatrix}\phi^{(c)}\\ \phi^{(c)*}\end{pmatrix}$ with $\phi^{(c)}=\begin{pmatrix}\phi^{(c)}_{\uparrow}\\ \phi^{(c)}_{\downarrow}\end{pmatrix}$ is the variational parameter describing the coherent part of the condensates, and $\xi(\mathbf{r},\mathbf{r}^{\prime})$ is the variational parameter that characterizes the squeezed part of the condensate. The ground state wave function can be found by minimizing energy $E_{0}=\left\langle\Psi_{\rm GS}\right|H\left|\Psi_{\rm GS}\right\rangle$ through imaginary-time evolution Shi et al. (2018, 2019); SM .

Physically, it is more convenient to factorize the Gaussian state wave function into a multimode squeezed coherent state Wang et al. (2020); SM , i.e.,

where $N^{(c)}=\int d\mathbf{r}\left[\phi^{(c)}(\mathbf{r})\right]^{\dagger}\phi^{(c)}(\mathbf{r})$ is the atom number in the coherent mode and $\hat{c}=\int d\mathbf{r}\left[\bar{\phi}^{(c)}(\mathbf{r})\right]^{\dagger}\hat{\psi}(\mathbf{r})$ with $\bar{\phi}^{(c)}(\mathbf{r})=\phi^{(c)}(\mathbf{r})/\sqrt{N^{(c)}}$ being the normalized mode function for the coherent component. In the squeezing operators, $\hat{s}_{i}=\int d\mathbf{r}\left[\bar{\phi}_{i}^{(s)}(\mathbf{r})\right]^{\dagger}\hat{\psi}(\mathbf{r})$ is the annihilation operator for the $i$th squeezed mode $\bar{\phi}^{(s)}_{i}\equiv\begin{pmatrix}{\phi}^{(s)}_{i,\uparrow}\\ {\phi}^{(s)}_{i,\downarrow}\end{pmatrix}$ that satisfy $\int d\mathbf{r}\left[\bar{\phi}_{i}^{(s)}({\mathbf{r}})\right]^{\dagger}\bar{\phi}^{(s)}_{j}({\mathbf{r}})=\delta_{ij}$. Furthermore, $\sinh\xi_{i}=\sqrt{N^{(s)}_{i}}$ with $N_{i}^{(s)}$ being the occupation number in the $i$th squeezed mode. The total number of atoms in the squeezed modes is $N^{(s)}=\sum_{i}N^{(s)}_{i}$ and the total squeezed density is $n^{(s)}({\mathbf{r}})=\sum_{i}N_{i}^{(s)}|\bar{\phi}_{i}^{(s)}({\mathbf{r}})|^{2}$. For convenience, it is always assumed that $N^{(s)}_{i}$ are sorted in descending order with respect to the index $i$. Mode $j$ is notably populated if $N^{(s)}_{j}/N^{(s)}\geq 0.1\%$ Wang et al. (2020).

In the main text, we focus on the self-bound droplets in quasi-2D geometry achieved by imposing a harmonic confinement, $V_{z}(z)=M\omega_{z}^{2}z^{2}/2$, along the $z$ axis Cabrera et al. (2018). When $\omega_{z}$ is sufficiently large, the motion of atoms along the $z$ direction is frozen to the ground state of $V_{z}(z)$. In additioin, we always fix the atom number ratio between two components at $N_{\uparrow}/N_{\downarrow}=\sqrt{a_{\downarrow\downarrow}/a_{\uparrow\uparrow}}$ by following the experiments Cabrera et al. (2018); Semeghini et al. (2018). Numerical results can be conveniently checked using the virial relation, $E_{k}+E_{2B}+2E_{3B}=0$, where $E_{k}$, $E_{2B}$, and $E_{3B}$ are kinetic, 2B, and 3B energies, respectively SM . We point out that results for 3D droplets are presented in the Supplemental Material SM .

Three-body coupling strength.—Let us first fix the value of $g_{3}$ for K atom. Previously, $g_{3}$ of Dy atom was determined by fitting the AND of the droplet states Wang et al. (2020). Here we instead determine $g_{3}$ by fitting the radial size $\sigma$ of the quasi-2D droplets. Following the experiment Cabrera et al. (2018), we extract $\sigma$ by fitting the total density with a two-dimensional Gaussian. In Fig. 1(a), we plot the radial size $\sigma$ as a function of the effective scattering length $\delta a\equiv a_{\uparrow\downarrow}+\sqrt{a_{\uparrow\uparrow}a_{\downarrow\downarrow}}$ (or, equivalently, the magnetic field $B$) computed with $g_{3}=6.65\times 10^{-39}\hbar{\rm m}^{6}/{\rm s}$. Good agreement between the numerical and the experimental results is achieved. We point out that varying the value of $g_{3}$ does not change the shape of the red line; it only shifts the red line vertically. Therefore, to achieve better agreement for the fitting, one may have to make $g_{3}$ spin dependent.

To further demonstrate the validity of the fitted $g_{3}$, we show, in Fig. 1(b) and (c), the $N$ dependence of $\sigma_{\rm GST}$ under different magnetic fields. As can be seen, although discrepancies exist at small $N$, good agreements are still achieved for large $N$. Particularly, the main feature of the experimental results is captured by the numerical simulations (see below for details). Hereafter, we shall always use the fitted $g_{3}$ for all numerical results presented below.

Phase diagram.—Given that $N_{\uparrow}/N_{\downarrow}=\sqrt{a_{\downarrow\downarrow}/a_{\uparrow\uparrow}}$, it is found numerically that the coherent and squeezed modes satisfy $\phi^{(c)}_{\uparrow}({\mathbf{r}})/\phi^{(c)}_{\downarrow}({\mathbf{r}})=\sqrt{N_{\uparrow}/N_{\downarrow}}$ and $\phi^{(s)}_{i,\uparrow}({\mathbf{r}})/\phi^{(s)}_{i,\downarrow}({\mathbf{r}})=\sqrt{N_{\uparrow}/N_{\downarrow}}$ for any $i$, respectively. As a result, $N^{(c)}_{\alpha}/N_{\alpha}$ is independent of the spin component $\alpha$, which allows us to define the coherent fraction as $f_{c}\equiv N^{(c)}/N$ to characterize the property of the condensate, instead of considering the coherent fraction of individual spin component.

Figure 2(a) shows the distribution of $f_{c}$ on the $(\delta a,N)$ parameter plane. Similar to dipolar droplets Wang et al. (2020), there are four distinct regions: expanding gas (EG), squeezed-vacuum state (SVS), squeezed-coherent state (SCS), and coherent state (CS). The EG region lies below the CAN (solid line or $\triangle$), in which the attractive two-body interaction is insufficient for the gas to form a self-bound state. Both SVS and SCS phases contain a single squeezed mode except SCS phase also consists of a coherent component. The SVS-to-SCS transition (dash-dotted line or $\diamondsuit$) breaks the Z₂ symmetry of the condensate and is a third-order phase transition. The CS phase is dominated by the coherent component with $f_{c}>0.75$ and the SCS-to-CS transition is of first order. Interestingly, as shown by the solid and the dash-dotted lines, the boundaries for the self-bound state and the SVS-to-SCS transition can be found nearly analytically with high accuracy SM . We may also determine the boundary of SCS-to-CS transition by analyzing the stability of Bogliubov excitation in the CS phase SM .

The phase diagram is most conveniently reproduced by plotting the number of the notably-populated squeezed modes $\mathcal{S}$. As shown in Fig. 2(b), both SVS and SCS phases are featured by single squeezed mode; while CS phase may contains a large number of squeezed modes even if the squeezed component becomes negligibly small. Physically, squeezing in SVS and SCS phases are macroscopic quantum state originating from two-body attraction Shi et al. (2019); Wang et al. (2020); while squeezing in CS phase represents quantum depletion induced by 3B repulsion Wang et al. (2020). Similar to dipolar droplets, the presence of squeezing also significantly modifies the AND of these quantum states SM . The properties of these quantum phases can be further explored by examining the peak density $n_{p}$ and the radial size $\sigma_{\rm GST}$ as plotted in Fig. 2(c) and (d), respectively. Among them, $\sigma_{\rm GST}$ exhibits visible difference as compared to the distributions of other quantities.

Radial size.—To explore the properties of $\sigma_{\rm GST}$, we plot the $N$ dependence of $\sigma_{\rm GST}$ for $\delta a=-3.062a_{0}$ in Fig. 3. As a comparison, the radial size numerically computed via EGPE, $\sigma_{\rm EGPE}$, is also plotted. Surprisingly, $\sigma_{\rm GST}$ is W-shaped function of $N$, in striking contrast to the V-shaped $\sigma_{\rm EGPE}$. To understand its origin, we note that the total energy per atom can be expressed as SM

where $E_{0}=E_{k}+E_{2B}+E_{3B}$ are the total energy. On the right hand side, the first term is the kinetic energy, the second term is contributed by the 2B interaction, and the last term is the energy associated with the stabilization force, i.e., $\nu=3$ for 3B repulsion and $5/2$ for quantum fluctuation. Furthermore, $\tilde{g}_{2}\,(>0)$ and $\tilde{g}_{\nu}\,(>0)$ are reduced interaction parameters (RIPs). From Eq. (7), the equilibrium size can be immediately found to be

which is a V-shaped function of $N$. More specifically, $\sigma$ diverges as $N$ approaches the CAN

grows asymptotically as $\sigma\approx\sqrt{N}\left[(\nu-1)\tilde{g}_{\nu}/\tilde{g}_{2}\right]^{1/[2(\nu-2)]}$ for large $N$, and reaches its minimal value at the DAN

In contrast, we note that the peak density, $n_{p}\propto N/\sigma^{2}=\left\{(\tilde{g}_{2}N-1)/[(\nu-1)\tilde{g}_{\nu}N]\right\}^{1/(\nu-2)}$, is always a monotonically increasing function of $N$ and converges to a constant proportional to $\left\{\tilde{g}_{2}/[(\nu-1)\tilde{g}_{\nu}]\right\}^{1/(\nu-2)}$ for large $N$. It is worthwhile to point out that $\sigma(N)$ depends not only on the stabilization mechanism through $\nu$ but also on the many-body wave function via the RIPs $(\tilde{g}_{2},\tilde{g}_{\nu})$ SM .

To carry out quantitative comparisons, we assume that all density profiles are Gaussian functions. Then for pure squeezed and coherent states with 3B interaction, we obtain the RIPs $(\tilde{g}_{2}^{(s)},\tilde{g}_{3}^{(s)})$ and $(\tilde{g}_{2}^{(c)},\tilde{g}_{3}^{(c)})$, respectively, satisfying $\tilde{g}_{2}^{(s)}=3\tilde{g}_{2}^{(c)}$ and $\tilde{g}_{3}^{(s)}=15\tilde{g}_{3}^{(c)}$ SM ; Wang et al. (2020). In addition, for the coherent-state-based EGPE theory, the RIPs are $(\tilde{g}_{2}^{(c)},\tilde{g}_{5/2}^{(c)})$ SM . Correspondingly, we plot $\sigma_{\nu}^{(c,s)}(N)$ using Eq. (8) in Fig. 3, where the subscript denotes the stabilization mechanism and the superscript denotes the quantum states. As can be seen, $\sigma_{3}^{(s)}$ agrees very well with $\sigma_{\rm GST}$ around the left dip where the condensate is a SVS. This agreement also indicates that the density profile of a SVS is well approximated by a Gaussian function. Around the right dip, the quantum state of the condensate is more complicated than a pure coherent state. Therefore, $\sigma_{3}^{(c)}$ only exhibits rough agreement with $\sigma_{\rm GST}$ for large $N$ where the quantum state is dominated by coherent component. Still, the remaining discrepancy can be explained as the density profile of the coherent state is significantly deviated from a Gaussian function SM . As to $\sigma_{5/2}^{({c})}$, it agree well with $\sigma_{\rm EGPE}$, which confirms the variational calculation. One should note that, although $\sigma_{5/2}^{(c)}$ and $\sigma_{3}^{(c)}$ have the same CANs, they define different DANs.

It is now clear that the W-shaped $\sigma_{\rm GST}$ stems from the multiple quantum phases in the self-bound droplets. Remarkably, this feature can be vaguely seen in the experimental data Cabrera et al. (2018). In fact, as shown in Fig. 1(b) and (c), after the completion of the left dip, the experimental data of the radial size starts to decrease again as one further increases $N$. Since radial size will eventually grow as $\sqrt{N}$ in the liquid phase, the drop after the first dip then signals the existence of the second dip and, consequently, the W shape. As the coherent-state-based EGPE theory only predicts a V-shaped radial size curve, a clear observation of the double dips may qualitatively exclude the EGPE explanation of the droplet states.

More interestingly, variational analyses can even lead to quantitative criteria for quantum phases and stabilization mechanisms. To see this, let us examine $N_{\rm cri}$, the CAN extracted analytically from Eq. (8). At first sight, it is surprised that $N_{\rm cri}$ is independent of $\tilde{g}_{\nu}$. This can be understood as follows. Because $\varepsilon_{k}$ and $\varepsilon_{2B}$ have the same power dependence on $\sigma$, the 2B attraction can be canceled out exactly by the kinetic energy when $N=N_{\rm cri}$. Then, in the absence of stabilization force, a self-bound state can be of arbitrary radial size. Once the stabilization force is turned on, in order to have a self-bound state of the same atom number, the radial size of the droplet has to be infinite to make the repulsive interaction energy vanish. This universality of $N_{\rm cri}$ is highly nontrivial as it allows our theory to be examined by experimental data regardless of the accuracy of the fitted $g_{3}$. We point out that the above analysis is not applicable in 3D geometry where a self-bound droplet always has a finite size SM .

Still, the CAN depends on the quantum state of the condensate through the RIP $\tilde{g}_{2}$. Then corresponding to the pure squeezed and coherent states, we have distinct CANs $N_{\rm cri}^{(s)}$ and $N_{\rm cri}^{(c)}$ which satisfy $N_{\rm cri}^{(c)}=3N_{\rm cri}^{(s)}$. The precise measurements of CAN can therefore be used to distinguish the quantum state of the condensates. In Fig. 2(a), we show the experimentally measured CAN versus the reduced scattering length. From the rightmost three data points, it seems that the quantum state at the self-bound boundary is coherent state; while the quantum state for the other three data point is unclear. However, after carefully examining these data, we find that high precision measurements are still needed to accurately determine the quantum phase SM . In fact, to precisely determine CAN, the experimental data should be fitted using proper curves corresponding to appropriate stabilization mechanisms. In particular, the experimental data should be of V shape close to the self-bound boundary.

Finally, we note that even the stabilization mechanism can be identified using DAN $N_{\rm dip}$, a quantity depending on $\nu$. From Eq. (10), it is clear that, corresponding to three $\sigma_{\nu}^{(c,s)}(N)$, there are three distinct DANs that satisfy $N_{{\rm dip},3}^{(c)}=3N_{{\rm dip},3}^{(s)}$ and $N_{{\rm dip},5/2}^{(c)}=(9/2)N_{{\rm dip},3}^{(s)}$, where $N_{{\rm dip},3}^{(s)}=2/\tilde{g}_{2}^{(s)}$. Therefore, precise measurement of DAN allows us to distinguish not only the quantum phase but also the stabilization mechanism. Moreover, $N_{\rm dip}$ is also independent of $\tilde{g}_{\nu}$, which makes the determination of the stabilization mechanism independent of the fitting of $g_{3}$.

Conclusion.—We have provided a detailed phase diagram for self-bound droplets of quasi-2D binary condensates beyond the coherent-state based mean field theory. Among all quantum phases, SVS and SCS are two macroscopic quantum squeezed states that demand experimental verification. More importantly, other than measuring AND or second-order correlation function, we have found three easily accessible experimental signatures: i) the observation of the double dips on the $\sigma(N)$ curve rule out the explanations of the liquid droplets based on single quantum state; ii) the precise measurements of $N_{\rm cri}$ or $N_{\rm dip}$ further allow us to distinguish the squeezed and the coherent states; iii) we may even determine the stabilization mechanism from $N_{\rm dip}$. More remarkably, both $N_{\rm cri}$ and $N_{\rm dip}$ are independent of the strength of the stabilization force, which makes the comparison between the theoretical and the experimental results immune to any unknown parameters.

Acknowledgement.—We thank enlightening discussions with H. Hu, X.-J. Liu, C. R. Cabrera, C. Navarrete-Benlloch, T. Pfau, and H. P. Buechler. This work was supported by National Key Research and Development Program of China (Grant No. 2017YFA0304501), by the NSFC (Grants No. 11974363, and No. 12047503), by the Strategic Priority Research Program of CAS (Grant No. XDB28000000), and by the Open Project of Shenzhen Institute of Quantum Science and Engineering (Grant No. SIQSE202006).

Supplemental Material

In the Supplemental Material, we present the derivation for the imaginary-time evolution in Gaussian-state theory (Sec. SM-I), the atom number distribution of a Gaussian state (Sec. SM-II), the virial relation (Sec. SM-III), the density profile (Sec. SM-IV), the boundary of the SVS-to-SCS transition (Sec. SM-V), the boundary of the SCS-to-CS transition (Sec. SM-VI), the variational analysis for radial size (Sec. SM-VII), and the main properties of the 3D droplets (Sec. SM-VIII).