Lee-Huang-Yang effects in the ultracold mixture of $^{23}$Na and $^{87}$Rb with attractive interspecies interactions

The beyond-mean-field Lee-Huang-Yang (LHY) correction is ubiquitous in dilute ultracold quantum gases. However, its effects are often elusive due to the typically much larger influence of the mean-field energy. In this work, we study an ultracold mixture of $^{23}$Na and $^{87}$Rb with tunable attractive interspecies interactions. The LHY effects manifest in the formation of self-bound quantum liquid droplets and the expansion dynamics of the gas-phase sample. A liquid-to-gas phase diagram is obtained by measuring the critical atom numbers below which the self-bound behavior disappears. In stark contrast to trapped gas-phase condensates, the gas-phase mixture formed following the liquid-to-gas phase transition shows an anomalous expansion featuring a larger release energy for increasing mean-field attractions.

In this work, we study in two different ways the LHY effects in a double BEC of 23 Na and 87 Rb atoms with tunable attractive interspecies interactions. First, we observe the heteronuclear quantum liquid droplet in free space with more than 10 4 atoms when the interspecies Na-Rb scattering length is tuned into the MF collapse regime. Under optimized conditions, a low-number-density droplet with lifetime exceeding the observation time is observed. We also investigate the liquid-to-gas-phase transition and obtain the critical atom numbers at the phase boundary. Second, we measure the release energies of two types of gas-phase mixtures, the pure in-trap gas and the gas formed after a droplet crosses the liquid-to-gas transition, and observe their opposite dependence on the interaction strength. Coupled with calculations based on extended Gross-Pitaevskii equations (eGPEs), our results confirm the crucial contribution of E LHY and its effects in stabilizing the heteronuclear double BEC far into the MF collapse region.
Our experiment starts from an optically trapped double BEC of 23 Na and 87 Rb (denoted as species 1 and 2 hereafter, respectively) both prepared in their lowest-energy hyperfine Zeeman level |F = 1, m F = 1 [26]. To reveal the LHY effects, we focus on the near-zero MF region with δg = g 12 + √ g 11 g 22 ≈ 0. Here g i j = 2πh 2 a i j /M i j are the two-body interaction constants, with a i j the scattering lengths and M i j the reduced masses. To reach this regime, we use the Na-Rb Feshbach resonance at B 0 = 347.648 G to tune the Na-Rb scattering length. The scattering length is represented as a 12  lengths a 11 = 60.05a 0 [29,30] and a 22 = 100.13a 0 [30,31] remain almost constant.
For convenience, we define a dimensionless parameter δg ≡ δg/g to characterize the total E MF , with g = (g 11 + g 22 )/2. For the current system, the critical magnetic field for δg = 0 (and also δg = 0) is B c = 349.978 G, which corresponds to a critical interspecies scattering length a c 12 = −63.1a 0 . As shown in Fig. 1(a) (lower panel), for each value of δg < 0, the double BEC will undergo a transition from gas phase to droplet phase as the atom numbers are increased. With our atom numbers around 10 4 for both species, the theory predicts an observation window for the droplet phase for δg from −0.061 to −0.189 (349.910 G to 349.780 G). For δg from 0 to −0.061, droplet formation is still possible but atom numbers much larger than currently available are needed.
In the first experiment, we study the Na-Rb quantum droplet in free space by releasing the sample from the optical trap. We probe the system during the time of flight (TOF), as depicted in the inset of Fig. 1(a). To probe the atoms in situ at each TOF, we implement a two-species high-field detection protocol which consists of a partial optical pumping process followed by standard absorption imaging on the cycling transitions [32]. This partial imaging method is necessary as the droplets have very high optical depth and are difficult to detect directly [11,33]. Both the partial optical pumping and the absorption imaging are calibrated using standard methods [34,35]. The measured resolutions (1/ √ e half Gaussian widths) of our imaging system are 0.6 and 0.8 μm for 23 Na at 589 nm and 87 Rb at 780 nm, respectively, which are just enough to resolve the droplet. To maintain the imaging resolutions, we mount the objective and the probe light optics on high-precision vertical translational stages to keep the optical axis of the imaging system always aligned with the droplet.
While shutting down the optical trap provides a direct way of probing the droplet in free space, the nonadiabatic nature of the process can cause problems. Due to the confinement, the in-trap sample is smaller in size than the free-space droplet in its ground state, while the total energy of the system is much larger. In the worst scenario, the energy of the initial state on release from the trap is larger then the energy barrier of the expansion and a stable droplet can never be formed [30]. To mitigate this problem, the crossed optical dipole trap is configured in a nearly spherical shape with measured trap oscillation frequencies of 78 Hz (86 Hz) for 87 Rb ( 23 Na). This is in accordance with the understanding that a quantum droplet with isotropic short-range interactions should be spherical in free space. In addition, a carefully designed magnetic-field control sequence is used to improve the mode matching [30]. These steps allow us to observe droplet formation reliably for δg from −0.094 to −0.189. However, for δg from −0.061 to −0.094, due to the smaller binding energy, droplets can be observed only by a more sophisticated mode-matching method, aided by fast magnetic-field quenching at the instant of releasing the samples from the trap [30]. For this method to work, the starting magnetic field must be selected carefully so that the size of the in-trap sample matches that of the free-space droplet at the end of the magnetic-field quenching. Figure 1(b) shows a side-by-side comparison of the 23 Na and 87 Rb clouds in the gas phase and the droplet phase following the TOF. The signature self-bound behavior of the droplet can be observed clearly by comparing its dramatically different TOF expansion with that in the gas phase. For δg = −0.119 (349.849 G), where the droplet phase is expected, the sizes of both clouds stay nearly the same during the TOF. For δg = +0.344 (350.451 G), by contrast, the system remains in the gas phase, and the sizes increase steadily with TOF. To extract quantitative information, we fit the images with 2D Gaussian functions. In the droplet, the size should be the same for the two species; their densities and thus numbers should follow the ratio N 2 /N 1 = √ g 11 /g 22 = 1.51 [7]. When this ratio is not maintained, the droplet part shows up as a dense central peak, and the excess atoms of one species appear as a much larger-sized and expanding gaseous background surrounding the droplet [inset of Fig. 1(c)]. For this kind of "bimodal" distribution, the droplet parameters are extracted with a double 2D Gaussian fitting. Figure 1(c) shows the starkly different expansion behaviors in the droplet and the gas phases obtained from the images in Fig. 1 For the droplets formed by direct trap release from δg = −0.094 to −0.189, we observe two distinctively different dynamics in the evolution of the 23 Na and 87 Rb samples during the TOF. As an example, Figs. 2(a)-2(c) show the sizes, atom numbers, and number ratio of the droplet sample for δg = −0.116 (349.852 G). During the first 10 ms, the 23 Na and 87 Rb sizes increase very little, while atom losses are observed for both species. Bimodal distributions are observed from 0 to 5 ms as the initial number ratio is far from 1.51. After 5 ms expansion, the gaseous atoms surrounding the droplet are already too dilute to be detected.
After about 10 ms, the sample sizes start to increase while the atom numbers stay nearly constant. This is consistent with a phase transition from droplet to gas when the atom numbers are reduced to below the critical values [see Fig. 1(a)]. We fit the size evolution empirically with σ (t > t 0 ) = σ 2 0 + v 2 (t − t 0 ) 2 and σ (t < t 0 ) = σ 0 and obtain a lifetime in the droplet phase of about t 0 = 7 ms. Here v is the expansion velocity of the gas-phase sample. Similarly, as illustrated in Fig. 2(b), the critical atom numbers for the phase transition can be obtained by fitting the 23 Na and 87 Rb number evolution data with an exponential decay function with the critical number as the offset.
For several values of δg between −0.061 and −0.094, we have created longer-lived droplets by the magnetic-field quenching method. As shown in Figs. 2(d)-2(f), during the accessible TOF, the droplet size stays nearly the same although number losses are still observed. After 10 ms, the number ratio becomes close to 1.51 and the number loss slows down significantly. As the lifetime is longer than the usable TOF of 18 ms, we have not been able to measure the real lifetime and the critical numbers in this range of δg. More detailed investigations, e.g., by levitating the droplet, are warranted in the near future. Figure 3 shows the critical 23 Na and 87 Rb atom numbers for δg from −0.094 to −0.189. Although the range of δg is small, a fourfold change in the measured critical numbers is observed. The critical numbers calculated using coupled eGPEs with the LHY term included (solid curves) matches the measurements very well.
A major cause of the atom losses is three-body recombination as a result of the high number densities in the droplet [10][11][12]30]. For droplets with δg from −0.094 to −0.189, another possible loss mechanism is self-evaporation due to the imperfect matching between in-trap and free-space modes [36]. The short lifetimes in this region are the combined result of these loss mechanisms. In contrast, the droplet at δg = −0.089 [see Figs. 2(d)-2(f)] has relatively lower number densities and is formed with nearly optimal mode matching and is thus longer lived. Nevertheless, for both cases, the stabilized number ratios after 10 ms TOF are very close to 1.51 despite the different creation procedures and the large initial number mismatches [Figs. 2(c) and 2(f)].
We now turn to the gas-phase double BEC and release the sample directly from the trap (without the additional magnetic-field quenching) for δg from +0.344 to −0.094. To characterize the LHY effect, we measure the release energy The open squares represent E rel for gas-phase mixtures released from the trap directly, while the solid circles are E rel for the gas samples formed following the liquid-to-gas-phase transition. The red solid curve and the red dashed curve show E rel for these two cases, calculated with eGPEs. The magenta short dashed curve is E rel for the second case from variational theory. The blue solid curve shows the GPE calculation without the LHY term for comparison. E kin is the quantum pressure and E int = E MF + E LHY is the total interaction energy. When E MF is tuned to zero, E LHY becomes the only interaction term in E rel . For negative E MF , the fact that the sample expands rather than collapses is also a direct manifestation of the LHY effect.
We prepare the trapped double BEC and measure E rel from the TOF expansion velocities of the two species [30,37] for several different values of δg. The resulting E rel is shown as open squares in Fig. 4. In general, E rel depends on δg, N 1 , and N 2 and the density overlap between the two species. However, in this range of δg, the densities of the double BEC are still rather high even in the gas phase, and significant three-body losses are observed during the TOF. Thus, it is difficult to maintain constant N 1 and N 2 during the TOF for different δg, even though great efforts are made to keep the initial atom numbers in the trap fixed. Nevertheless, since E rel is a measure of the energy per atom, it is not affected severely by the loss.
We calculate E rel with the eGPEs for the corresponding atom numbers at the end of the TOF and find good agreement for δg 0, corresponding to E MF E LHY or E MF ≈ 0 (red solid curve in Fig. 4). We note that, without the LHY term, the calculation fails for δg < 0 due to collapse of the double BEC (blue solid curve in Fig. 4). At δg ≈ 0, where E MF almost cancels, the energy difference between the two calculations is about 1 nK. This is essentially the magnitude of E LHY , if we assume E kin does not change dramatically between the two cases.
The gas-phase mixture formed following the liquid-to-gasphase transition is also of great interest. As shown in Fig. 2(a), its expansion velocity can be measured and converted into another E rel even though there is no trap involved. Following the number losses, the total energy of the system changes from negative in the droplet phase to positive in the metastable and eventually gas phase. E rel here is thus just the sum of E kin , E MF , and E LHY at the phase-transition point [30]. The measured E rel is shown in Fig. 4 as solid circles for δg from −0.094 to −0.189. It has the opposite dependence on δg from that observed for the gas-phase mixture released from the trap. This is similar to previous observations on the expansion dynamics of the erbium dipolar system [17].
This behavior can be qualitatively understood from the competition between E LHY , E kin , and E MF for the droplet near the phase-transition point. As δg becomes more negative, the critical atom numbers at the phase-transition point become smaller, and the sample also shrinks. The positive E LHY and E kin terms increase faster than the negative E MF term decreases. The measured E rel thus increases for more negative δg. We note that, without E LHY , the magnitude of E MF is larger than E kin and the system collapses. This picture is confirmed semiquantitatively by the eGPE calculation, shown by the red dashed curve in Fig. 4. As a comparison, a variational calculation (magenta short dashed curve) using the experimentally measured atom numbers and the Gaussian ansatz is also shown; this actually agrees better with the experimental data. The deviation between the experiment and the eGPE calculation may be due to the inaccuracy in determining the transition point.
In conclusion, we have studied the heteronuclear BEC mixture of 87 Rb and 23 Na atoms with attractive intersperses interactions under conditions where the role of the LHY correction is amplified. The signatures of the LHY correction manifest clearly in both the self-bound behavior of the free-space quantum liquid droplet and the expansion of the gas-phase mixtures. Although the observation time is limited by the detection method, our result shows that heteronuclear droplets can be long lived. This will allow quantitative studies on the fundamental properties of the quantum droplet, such as the existence of a possible Bose pairing mechanism [39,40] with mass-imbalanced constituents. Another interesting topic is quantum droplets in lower dimensions and at the crossover between different dimensionalities [41][42][43]. The 23 Na and 87 Rb system here also holds great promise in creating large quantum liquid droplets with more than 10 6 atoms [44,45]. In this so far unreached regime, the droplet features a large core with constant density, and various peculiar collective excitation may be investigated [7].  lengths with the MOLSCAT package [46,47] using the molecular potentials in Ref. [29] for 23 Na and in Ref. [31] for 87 Rb. For 87 Rb, the scattering length a 22 near 350 G stays nearly the same as the well-calibrated low-field value of 100.13a 0 . On the other hand, for 23 Na, the scattering length a 11 varies smoothly as the magnetic field is changed from 0 to 350 G, as shown in Fig. 5. Near the experimentally used magnetic field for the current experiment, a 11 is 60.05a 0 . This is rather different from the low-field value of 54.5a 0 .
The variation of a 11 far away from the Feshbach resonance is due to the change of the singlet/triplet character of the 23 Na atomic states. This is the result of the competition between the hyperfine interaction and the Zeeman effect. Since the difference between the singlet and triplet 23 Na -23 Na scattering lengths is large [29], a 11 will be modified by the singlet/triplet character variation. For 87 Rb, since the singlet and triplet scattering lengths are very similar [31], a 22 is very insensitive to magnetic fields unless there is a Feshbach resonance.

Magnetic field gradient compensation
An important technical issue for observing the free-space droplet is the spatial varying magnetic field experienced by the droplet during the TOF. Due to the different magnetic dipole moments of 23 Na (μ b1 = 1.084 MHz/G) and 87 Rb (μ b2 = 0.837 MHz/G) near B 0 , a moderate gradient and curvature may generate a differential force large enough to tear the droplet apart. With the pair of Feshbach coils only, the measured magnetic field curvature near 347 G is 1.4 G/cm 2 and the averaged magnetic field gradient is about 0.6 G/cm. Experimentally, we observe the 23 Na and 87 Rb atomic clouds separating from each other following the TOF. To mitigate this problem, we placed another coil, with its current dynamically controlled following the TOF, on top of the Feshbach coil. This simple improvement reduces the magnetic field gradient to 0.11 G/cm and the magnetic field variation during the 18 ms TOF from ±10 mG to ±3 mG, which corresponds to a a 12 change of less than 0.26 a 0 . FIG. 6. Droplet formation by directly releasing the trapped sample. The blue dashed curves represent E tot = E kin + E pot + E int of the trapped double BEC, and the red solid curves are E tot in free space (i.e., E kin + E int ). (a) For very small | δg|, the free-space E tot (the red open circle) is higher than the barrier height on the large size side (the red solid dot). In this case, the sample will keep expanding and never forms the droplet. (b) For | δg| large enough, the sample does not have enough energy and will be confined near the energy minimum to form the droplet.
In the eGPEs, the magnetic field gradient is included as opposite direction linear potentials to the two species in the center-of-mass coordinate. According to the simulation, this small magnetic field gradient actually doubles the critical atom numbers.

In-trap and free-space mode matching
During the release of the double BEC sample from the optical trap, the trap energy suddenly is removed suddenly, but the size of the sample has no time to adapt to the new value. As shown in Fig. 6, in general the in-trap size is smaller than the size of the ground-state free-space droplet at the energy minimum. The sample size will thus start to expand after being released from the trap and its total energy (per particle) E tot = E kin + E int will also start to evolve following the red curves. For the cases with very small | δg| [ Fig. 6(a)], E tot is higher than the energy barrier at the right-hand side. In this case, the sample will just keep expanding after crossing the maximum point and never forms the droplet, while in Fig. 6(b), with | δg| large enough, the negative E tot of the released sample is not enough to overcome the energy barrier and the droplet can be formed, but its size will undergo a small-amplitude oscillation.
As described in the main text, to improve the mode matching, we used a nearly spherical shape trap for matching the sphere shape of the quantum droplet. In addition, as the in-trap size of the double BEC also depends strongly on δg, the magnetic field is controlled in a carefully designed sequence. After the double BEC is created, the magnetic field is first ramped up quickly across the Feshbach resonance to 358.4 G where a 12 has a small positive value of 46.1a 0 . After 350 ms for the system to stabilize, the magnetic field is ramped down to 350.451 G in 130 ms. This is just 0.424 G above B c and the corresponding a 12 is −39.5a 0 ( δg = +0.344). Here the double BEC is miscible and the attractive interaction increases the spatial overlap of the two condensates and also increases the in-trap densities substantially. Both effects significantly improve the mode matching and facilitate the free-space droplet formation. At this point, rapid atom losses are already observed; thus, the magnetic field holds here for only 10 ms. Afterward, the magnetic field is ramped to a target value within 5 ms before the droplet is detected in free space.
As depicted in Fig. 6(b), although the droplet can be observed after these improvements for | δg| large enough, the mode matching is not perfect. For δg from 0 to −0.094, with the current atom numbers, the droplet cannot be observed by simply releasing the double BEC from the trap. As shown in Fig. 7(a), even better mode matching can be achieved by finding a δg value where the in-trap size matches with the target droplet size and quenches δg at the same time when the trap is shut off. To achieve the fast magnetic-field quenching necessary for this method, we added a small magnetic coil that is capable of covering the magnetic field range in less than 10 μs. Figures 7(b) and 7(c) show two droplets formed with this method. For case (b) with δg = −0.089, the droplet lifetime is longer than the available observation time. For case (c) with δg = −0.061, the liquid-to-gas-phase transition is observed after about 7 ms. This short lifetime is probably limited by the atom numbers since at this δg, the critical atom numbers are much larger. Origin of E rel for the gas released from the trap. The blue dashed curve is the in-trap total energy E tot = E pot + E kin + E int , and the red solid curve is E tot without E pot when the sample is suddenly released from the trap. (b) E rel for the case following the liquid-to-gas-phase transition. The blue and orange dashed curves are E tot of the stable and meta-stable droplets, respectively. The red solid curve is E tot just on the phase-transition point, where the local minimum disappears. For both cases, following the size increase, the quantum pressure E kin and the interaction energy E int are converted to the kinetic energy of the atoms.

Three-body loss
As discussed in the main text, the major losses of atom number are from a three-body loss as a result of the high density of the droplet sample. We observe a faster loss of 87 Rb than that of 23 Na at almost all magnetic fields. We believe this is due to the larger three-body loss rate constant (7.2 × 10 −30 cm 6 /s) of the 87 Rb BEC than that of the 23 Na BEC (1.1 × 10 −30 cm 6 /s) [48,49].
While the three-body loss rate constants for condensed 87 Rb and 23 Na atoms are well known [48,49], those between 23 Na and 87 Rb with the weak interspecies interaction near B c are unknown. But from our previous few-body investigations [50], no heteronuclear Efimov resonances are expected near this region; thus, we can assume that the heteronuclear threebody loss should have a similar rate as single species 87 Rb or 23 Na samples.

Release energies
For the double BEC system studied here, the average release energy is defined as where v i is the expansion velocity. Figure 8 illustrates schematically the different origins of E rel for the in-trap gas and the gas formed from the droplet. For the first case [ Fig. 8(a)], E rel is dominated by E MF and decreases for smaller δg. For the gas from phase-transition case [ Fig. 8(b)], the total energy of the sample and thus E rel increases when δg gets more negative.
Experimentally, the RMS sizes are obtained from fitting to the absorption images. Different fitting functions are used for the two cases. For the sample directly released from the trap, the images are fitted with the parabola function in accordance with the Thomas-Fermi approximation for BECs with more than 10 4 atoms. For the gas mixtures formed following the liquid-to-gas transition, the Gaussian function, which is closer to the small atom number sample profiles, is used. 033247-6

Numerical simulation
For the condensate mixture of 23 Na and 87 Rb (denoted as species 1 and 2 hereafter, respectively) considered here, the coupled eGPEs including the LHY correction reads with m i the atomic mass, n i the number density, and μ i = g ii n i + g i j n j + δE LHY δn i , (i = j).
As introduced in Refs. [7,51], the LHY term for massimbalanced system can be described as with γ = m 2 /m 1 , x = g 2 12 /g 11 g 22 , and y = g 22 n 2 /g 11 n 1 . The dimensionless parameter f (γ , x, y) represents the renormalized integral with only numerical solution for the massimbalanced case.

Influence of the residue magnetic field gradient
As discussed in the Appendix, Sec. A 2, the magnetic field gradient compensation is not perfect in the current setup.
Although small, the residue gradient can in principle generate some influence to the droplet. The gradient ∂B/∂y is mainly along the vertical (y) direction. In the center-of-mass coordinate, the differential magnetic dipole moments lead to opposite direction linear potentials to the two species, where a c is the center-of-mass acceleration defined as Including V i in the eGPE, our simulation shows that the critical numbers needed for droplet formation can indeed be increased. The theoretical curve in Fig. 3 does not take this contribution into account. In future experiments, this should be investigated in more detail, for example, by studying the phase transition at various magnetic field gradients. For a better gradient compensation, additional coils will have to be installed.