Observation of Quantum Droplets in a Heteronuclear Bosonic Mixture

We report on the formation of heteronuclear quantum droplets in an attractive bosonic mixture of 41K and 87Rb. We observe long-lived self-bound states, both in free space and in an optical waveguide. In the latter case, the dynamics under the effect of a species-dependent force confirms their bound nature. By tuning the interactions from the weakly to the strongly attractive regime, we study the transition from expanding to localized states, in both geometries. We compare the experimental results with beyond mean-field theory and we find a good agreement in the full range of explored interactions. Our findings open up the production of long-lived droplets with important implications for further research.


I. INTRODUCTION
Interactions are ubiquitous in nature and understanding their effects is a primary challenge in physics, especially in many-body quantum systems. Generally, meanfield (MF) theories well reproduce a plethora of phenomena related to interparticle interactions. However, this approach fails when quantum fluctuations are nonnegligible. The first-order correction to the MF energy, the so called Lee-Huang-Yang (LHY) term, was first calculated in 1957 [1] to describe bosons with hard-core repulsion. In ultracold gases the LHY term is normally negligible and only recently it has been experimentally investigated [2][3][4][5][6]. Notwithstanding, there are situations where the LHY and MF contributions can be of the same order. For instance, in attractive bosonic mixtures the repulsive LHY term may stabilize the system against collapse, leading to the formation of self-bound droplets. Although very dilute, these states have a liquid-like behaviour, characterized by a core with uniform density in the large atom number limit [7].
While originally predicted for Bose-Bose mixtures, quantum droplets have been first realized in dipolar gases, by exploiting the competition between contact repulsion and anisotropic dipole-dipole attraction [8][9][10][11][12][13]. Recently, tuning contact interactions through a Feshbach resonance, quantum droplets have been also observed in a spin mixture of 39 K, both in the presence of an external potential [14,15] and in free space [16]. By studying collisions between droplets it has been possible to observe the crossover between compressible and incompressible regimes [17]. Free space 39 K droplets have a short lifetime (only few ms), limited by three body losses [16,17], motivating the quest for stable long-lived droplets in different * derrico@lens.unifi.it † burchianti@lens.unifi.it atomic mixtures. Longer lifetimes could indeed allow the investigation of many peculiar features of these states, like the characterization of the incompressible regime and the observation of self-evaporation [7].
In this work, we study the formation of quantum droplets in a heteronuclear Bose-Bose mixture of 41 K and 87 Rb where the two species experience different trapping potentials. We exploit Feshbach resonances for tuning the interspecies interaction. In particular, we explore the transition from the weakly to the strongly attractive regime, beyond the MF threshold for collapse, where self-bound quantum droplets are predicted to exist. We observe that in the absence of any external trapping potential, the strongly attractive mixture remains localized on a timescale of several tens of milliseconds, consistently with the formation of long-lived self-bound droplets. Furthermore, we study the size evolution and the centerof-mass dynamics of the mixture in a horizontal optical waveguide, in the presence of a species-dependent magnetic force. Also in this geometry, where both the LHY correction and the radial confinement provide a mechanism to prevent the collapse [18], we observe the formation of a localized bound state, whose center-of-mass follows a trajectory intermediate between those of two weakly attractive 41 K and 87 Rb clouds. By studying the dynamics, we extract the number ratio of 41 K and 87 Rb atoms forming the bound state, and we find it consistent with the value predicted by the theory [7]. Our experimental findings are well reproduced by numerical simulations performed by using two coupled generalized Gross-Pitaevskii (GP) equations [19,20], including the LHY correction for heteronuclear mixtures [7,21].

II. EXPERIMENT
We produce a binary condensate of 41 K and 87 Rb both prepared in the |F = 2, m F = 2 state in a crossed opti-arXiv:1908.00761v2 [cond-mat.quant-gas] 6 Nov 2019 cal dipole trap [22]. To tune the interspecies interaction, two microwave pulses transfer both species in the |F = 1, m F = 1 state, and a vertical homogeneous magnetic field is increased to 71 G, lying between two Feshbach resonances [23]. Due to their different masses, 41 K and 87 Rb atoms have a differential vertical gravitational sag and their equilibrium positions are separated by ∼ 15 µm in the optical trap. We compensate for this offset by compressing the dipole trap and by adding a vertical magnetic field gradient b z = −16 G/cm [24,25], exploiting the different magnetic moments of the two species (µ K = −0.83 µ B and µ Rb = −0.52 µ B , at 70 G), as shown in Appendix A. At this stage, we normally have a binary condensate with N K 1 − 2 × 10 4 and N Rb 4 − 6 × 10 4 atoms, confined in an approximately harmonic potential with different frequencies for the two atomic species: (ν x , ν y , ν z ) (130, 90, 180) Hz for 41 K and (100, 70, 130) Hz for 87 Rb. The homonuclear and heteronuclear scattering lengths are respectively a K = 65a 0 [26], a Rb = 100.4a 0 [27], and a KRb 0a 0 [28,29]. To enter the attractive regime, we then change the magnetic field to the desired value B F in 30 ms, thus tuning a KRb , as detailed in Appendix A. The mixture is completely characterized in terms of the intraspecies g K = 4π 2 a K /m K , g Rb = 4π 2 a Rb /m Rb and interspecies g KRb = 2π 2 a KRb /m r coupling constants, with m K and m Rb the atomic masses and m r the reduced mass. The onset of the MF collapse regime corresponds to δg = g KRb + √ g K g Rb = 0 (a KRb = −75.4a 0 ) [30]. We explore the crossover from the weak (a KRb < 0, δg > 0) to the strong (δg < 0) attractive regime. For the measurements in free space, the trap beams and the magnetic field gradient b z are simultaneously switched off and the atoms are imaged in time-of-flight (TOF). For the expansion in the horizontal waveguide (along thê y direction) one of the two beams of the dipole trap (labeled as crossed beam in Fig. 1a) is slowly turned off by linearly reducing its intensity in 10 ms, and the atoms are imaged in-situ after a given evolution time. In the waveguide the atoms experience a radial trapping potential with average frequencies ν r 160 Hz (120 Hz) and an axial anti-trapping potential, produced by the magnetic field gradient, with ν y −1.8 Hz (−0.6 Hz) for 41 K ( 87 Rb).

A. Dynamics in free space
We first consider the dynamics in free space. In Fig. 1b we show typical absorption images of the two condensates in the xz plane and TOF = 24.5 ms for δg > 0 (top) and δg < 0 (bottom). In the former case both 41 K and 87 Rb samples freely expand, while in the latter case we observe that a fraction of both species bounds and forms a small and dense component. The different shape of 41 K and 87 Rb clouds for δg < 0 is a consequence of our imaging procedure. To detect both atomic species, in fact, we first take an absorption image of 41 K sample at the magnetic field B F , then we switch off B F and after further 2.5 ms of TOF we take the image of 87 Rb atoms. The actual size of the bound state can be measured only in the first image, after which the state is dissociated and 87 Rb cloud expands. In Fig. 1c we plot the average width (rms) σ = √ σ x σ z of 41 K sample as a function of TOF for both interaction regimes. Whereas for δg > 0 (open circles) σ increases in time, for δg < 0 (closed circles) it does not exceed our imaging resolution (5 µm) up to 28 ms of TOF, indicating the formation of a droplet state. In our mixture the droplet density is expected to be one order of magnitude smaller than in 39 K at the same δg, leading to a substantial reduction of the three-body losses. This is consistent with a two orders of magnitude longer lifetime. Such difference is mainly due to the larger scattering lengths of 41 K and 87 Rb, since the density scales as n ∼ (δg/g) 2 a −3 for approximately equal scattering lengths, as discussed in Appendix B. In the experiment, for TOF > 28 ms, we observe that the atomic clouds start to expand. We attribute this to the variation of a KRb , since for larger TOF the atoms fall off the region where the Feshbach field is spatially homogeneous. Nevertheless, this time is longer than the maximum observation time of 39 K droplets in free space [16]. In order to increase the observation time, one could exploit a magnetic gradient b z to levitate the droplet, with an average magnetic moment per particle: where N d K and N d Rb are the number of 41 K and 87 Rb atoms forming the bound state. However, we observed that, due to the species-dependent magnetic force, the droplet breaks up for b z 5 G/cm. Thus, we apply a lower magnetic field gradient b z = 1.1 G/cm, that, even if too small to compensate for gravity, reduces the variation of the Feshbach field during TOF [31]. In addition, it introduces a species-dependent acceleration sufficient to spatially discriminate the bound and the unbound atoms. In the images at δg < 0 ( Fig. 1b bottom), we note indeed the presence of a halo surrounding the localized clouds of both 41 K and 87 Rb. Due to the finite temperature, a fraction of atoms remains unbound. Moreover, for an efficient sympathetic cooling of the mixture, the atom number of 87 Rb is typically three or four times larger than that of 41 K, leading to an excess of unbound 87 Rb atoms. For 41 K the unbound halo is upshifted (|µ K | > |µ KRb |), whereas for 87 Rb is downshifted (|µ Rb | < |µ KRb |). By measuring the accelerations of the different components one could extract µ KRb and therefore the ratio of the atom numbers N d K /N d Rb in the droplet. However, in free space, the acceleration is dominated by gravity (for the droplet the ratio between magnetic and gravitational forces is only 0.035).

B. Waveguide dynamics
To minimize the effect of gravity, we hold the atoms in a horizontal waveguide and let them evolve under the dominant effect of a magnetic field gradient along theŷ direction. In Fig. 2 we show in-situ images of the two condensates as a function of the evolution time t in the waveguide. In this configuration we are able to detect the bound states for a longer time, since the center-of-mass displacement is less than 300 µm in 50 ms. For δg > 0 (top), both species expand as t increases, while for δg < 0 (bottom) we observe the formation of self-bound states. Also here, as well as in free space, a fraction of the atoms remains unbound and follows the same dynamics of the weakly-attractive gas. In Fig. 3 we plot the size σ y and the center-of-mass position y cm along the longitudinal direction. For δg < 0 we fit the atomic density distribution with two Gaussians to distinguish the bound and unbound components. For clarity, σ y and y cm are shown only for the bound components, since the unbound fraction behaviour is very similar to the one observed for δg > 0 (see Appendix C). The measured size for δg < 0 (closed symbols) is constant and is comparable with the imaging resolution in the waveguide (10 µm) [32], while for δg > 0 both species expand (open symbols). In the waveguide, the center-of-mass moves under the combined effect of a magnetic force and a residual component of the gravity, since the guide is tilted of θ 0.1 − 0.2 • with respect to the horizontal plane. Assuming a constant acceleration (see Appendix C), we extract the ratio of atom number of the two species in the bound state as where α K and α Rb are the accelerations of the unbound 41 K and 87 Rb clouds, and α d is the acceleration of the bound state. From α K and α Rb measured at δg > 0 , which is consistent with the predicted value of g Rb /g K = 0.85 [7]. We note that the center-of-mass motion allows a measurement of the atom number ratio which is not affected by the usual systematic uncertainty of N measurements.

C. Comparison with theory
Both in free space and in the waveguide, we study the full range of attractive interactions. In Fig. 4 we plot the width σ x of 41 K in free space after TOF = 26.5 ms (Fig. 4a) and σ y after t = 45 ms of expansion in the waveguide (Fig. 4b). In both cases we observe a transition from an expanding gas for a KRb −82a 0 to a localized self-bound state for strong attractive interactions (a KRb −82a 0 ), whose width is below the optical imaging resolution. In Fig. 4, we also show the MF threshold for collapse (dash-dot gray line) and the critical scattering length to enter the droplet regime with our experimental atom number (dashed red line) [7]. Both lines are calculated for a homogeneous system. In the presence of the radial confinement a stable bound state may exist also in-between the two vertical lines. This is due to the combined effect of the LHY term and the dispersion along the guide giving rise to the formation of "solitonic" solutions [18,33]. In order to compare the experimental results with theoretical predictions [7,21], we have performed numerical simulations by using two coupled generalized timedependent GP equations including the LHY term, taking into account the actual experimental preparation (see Appendix B for details). In Fig. 4 the colored areas correspond to the numerical predictions including the experimental uncertainty on the atom number. Considering the finite imaging resolution, we find that the theory well reproduces the observed behavior. From the simulations in free space (Fig. 4a) we find that, in the strongly attractive regime, almost spherical droplets form [7]. As expected, the critical scattering length to enter the droplet regime depends on the atom number. We find that this threshold actually occurs at larger values of |a KRb | with respect to the homogeneous case (see dashed red line in Fig. 4). This is due to the additional kinetic energy cost associated to the droplet surface, whose effects are more important for small droplets, as those formed in the experiment.
As for the waveguide dynamics, we observe expanding clouds in-between the two vertical lines instead of the predicted "solitonic" bound states. We attribute this to a twofold effect: (i) we adjust the initial size of the mixture to match the droplet size (smaller than the one of the solitonic ground state); (ii) the dynamics is triggered by a non adiabatic process. However, droplets formed for very large attractive scattering lengths (|a KRb | 95a 0 ), as already observed in previous experiments [14], eventually decay into "solitonic" solutions due to the decrease of the atom number induced by three body losses. In the waveguide, indeed, atomic losses do not affect the existence of self-bound states allowing their observation for a wider range of a KRb . In free space, instead, as confirmed by numerical simulations with three body losses (Appendix B), we are limited to |a KRb | 95a 0 . This corresponds to a broader range of δg with respect to the case of 39 K droplets, pushing the investigation far away from the MF collapse.

IV. CONCLUSIONS
In conclusion, we have observed quantum droplets in a mixture of two atomic species experiencing different potentials, providing a clear picture of their formation and evolution both in free space and under radial confinement. This proves that the presence of identical traps for the two components is not essential to the droplets production. We have found evidence that the lifetime of our droplets is substantially longer than those in 39 K, as expected. The droplet density obtained by our numerical simulations results in a lifetime up to 400 ms. This opens new perspectives for future studies of the droplets peculiar properties [7,21] and of their collective modes [33]. Indeed, the lower density of our sample allows for a larger ratio τ life /τ ∼ n −1 , with τ life and τ being respectively the lifetime and the characteristic time scale of the droplet dynamics (Appendix B). Other interesting directions could be the creation of droplets in reduced dimensionality [34][35][36][37][38], where quantum fluctuations become easily dominant, and the investigation of beyond LHY corrections [39] and other stabilization mechanisms [40,41]. Finally, the nucleation of vortex-antivortex ring pairs [42][43][44] or more exotic scenarios, like droplets clusters [45], could be studied with rotating droplets. In Fig. 5, we show the behaviour of the interspecies scattering lengths a KRb for 41 K and 87 Rb in the |F = 1, m F = 1 state as a function of the magnetic field B. Two Feshbach resonances around 40 G and 80 G can be seen. The values of a KRb predicted by the two collisional models in Ref. [28] and in Ref. [29] slightly differ by few a 0 , as shown by the thickness of the curve in Fig. 5 (b). All the values of a KRb given in this work are the average of the two models, with uncertainty equal to the half-deviation that dominates over the one due the magnetic field calibration.
In the experiment we produce the double condensate at a KRb 0 and then we tune the magnetic field in the range 65 G B 69 G to access the attractive interspecies interaction regime (−100 a 0 a KRb −50 a 0 ). In particular, we use two linear ramps: from 71 G to 68.5 G in 20 ms and then to the final value in other 10 ms. At these magnetic fields, the Zeeman effect is no longer linear (especially for 41 K) and the magnetic moments µ K and µ Rb of the hyperfine state of 87 Rb and 41 K in our mixture vary as shown in Fig. 6. In particular, in the range of B used in the experiment the two magnetic moments are quite different, enabling a different magnetic control of the two atomic species, while their relative variations ∆µ K /µ K and ∆µ Rb /µ Rb are below 4% and therefore neglected. The difference between the two magnetic moments is exploited both to superimpose the two condensates by means of a magnetic gradient b z , and to induce a different acceleration of the unbound and bound components. As mentioned in Section III C, we include in our analysis a loss term due to inelastic three-body collisions. In general, the three-body losses coefficient K 3 depends on the scattering length, scaling as a 4 in vicinity of a Feshbach resonance. In our case, we explored a very limited range of magnetic fields and scattering lengths a KRb , therefore we measured the atoms lifetime to obtain the rate of inelastic three-body collisions at a KRb = (−77.8 ± 1.6)a 0 . We measured the lifetime of the 41 K sample, since (i) the number of the minor-ity species is more sensitive to losses; (ii) the dominant three-body losses are due to K-Rb-Rb collisions [46], yielding a purely exponential decay for 41 K atoms. We measured a 1/e-lifetime τ K = 122(11) ms with an estimated 87 Rb peak density n 0,Rb = (6±2)·10 20 m −3 . From 1/τ K = (K 3 /2!)n 2 0,Rb ξ, with ξ = n 2 Rb n K dV /(n 2 0,Rb N K ), we obtain K 3 = (7 ± 4) · 10 −41 m 6 /s with an uncertainty due to that of n 0,Rb . This must be compared to published values for larger scattering lengths: K 3 10 −38 m 6 /s at a KRb −120a 0 [46] and K 3 3 · 10 −39 m 6 /s at a KRb −400a 0 [47].
Here f (z, u, x) > 0 is a dimensionless function, whose expression for z = 1 and u = 1 can be found in [21]. Following [7,21], we also consider this function at the mean-field collapse u = 1, f (z, 1, x). We note that the actual expression for f can be fitted very accurately with the same functional form of the homonuclear case (m 1 = m 2 ) [48] f (87/41, 1, where α and β are fitting parameters that in general depend on the value of the mass ratio m 1 /m 2 . For the present case we find α 1.554 and β 2.506, which are very close to the values of the approximated formula proposed in Ref. [48], which shows that only α is a function of the mass ratio, α = (m 2 /m 1 ) 3/5 = 1.571, whereas the value of β is universal, β = 5/2. Minimization of the action associated to Eq. (B1) leads to the following generalized GP equations (Euler-Lagrange equations) where and The above equations are solved by mapping the system (densities, wave functions, differential operators, etc.) on discrete equally spaced cartesian grids. The differential operators are represented by a 13-point discretization. The stationary GP equation is solved by imaginary time evolution. Dynamical equations have been solved by using Hamming's predictor-modifier-corrector method, initiated by a fourth-order Runge-Kutta-Gill algorithm [49]. The effect of three-body losses on the realtime dynamics of the mixture is simulated by adding to the energy functional in Eq. (B1) a dissipative term −(i/2) K 3 drn 1 (r, t)n 2 (r, t) 2 (with K 3 = 7 · 10 −41 m 6 /s, see Appendix A) which accounts for the dominant recombination channel, i.e. K-Rb-Rb.
We performed two different series of simulations of the experiment for the TOF in free space and for the expansion in the waveguide. In the first case, the ground state of the mixture is produced in two concentric harmonic potentials (for the two condensates) approximating the actual crossed dipole trap. Then the potential is switched off abruptly and an evolution in free space is followed taking into account three-body losses. For the simulation in the waveguide, we adopt a complete description of the potential experienced by the two species including the optical, magnetic and gravitational potentials. This allows us to reproduce the effect of a residual magnetic field and a differential gravitational sag of the two components. The results reported in Fig. 4 (b) are obtained neglecting three-body losses after checking that their contribution is negligible in the range of parameters investigated. All results of the simulations shown in Fig. 4 have been performed with N Rb /N K = g K /g Rb . Any unbalance from this ratio leads to an unbound expanding fraction which does not affect the bound fraction, but heavily affects computation time requiring a larger spatial grid. The simulations show the formation of self-bound states even when the trapping potentials of the two condensates are not identical and completely overlapped. A detailed investigation of this effect goes beyond the scope of the present work.
The equilibrium density of a droplet can be obtained by the condition P = − + i (∂ /∂n i ) n i = 0, where P is the pressure and is the energy density. For the heteronuclear case this condition gives: Here a 11 , a 22 and a 12 are the intra and interspecies scattering lengths and δa ≡ a 12 −a c 12 , with a c 12 the mean-field threshold for collapse. From Eq. (B8) we estimate that, in our mixture, the homogeneous equilibrium density n 0 1 varies from ∼ 2 · 10 14 atoms/cm 3 to ∼ 2 · 10 15 atoms/cm 3 for −6 a 0 < δa < −20 a 0 , corresponding to the range for which we observe a droplet state (Fig. 4 (a)). Furthermore, for 41 K and 87 Rb, the terms containing the masses are ∼ 1 and Eq. (B8) is not substantially different from the equation for the homonuclear case given in [7]: From Eq. (B9), it is evident that upon increasing all scattering lengths by a factor γ > 1, the densities decrease as γ −3 . To compare the equilibrium density of droplets formed either by 41 K-87 Rb mixture or with two spin states 39 K, we use Eq. (B8) and Eq. (B9), respectively. In the first case a 11 and a 22 are larger than those of 39 K (for which a 22 is tuned through a Feshbach resonance, while a 11 and a 12 are constant), resulting in a droplet density lower by approximately an order of magnitude at δa of a few −a 0 . As already mentioned in the main text, this also results in a longer lifetime (τ life ∼ n −2 ) and a larger ratio τ life /τ ∼ n −1 , with τ being the characteristic time scale of the droplet [7].
APPENDIX C: WAVEGUIDE DYNAMICS

Equations of motion
In this section we derive the equations of motion for both atomic species confined in the waveguide. We start by considering the case in which they are unbound. The anti-confining force experienced by the i-th atomic species along theŷ axis is mainly determined by the magnetic potential U y,i . This is given by the Feshbach field B F and the overlapping magnetic field gradient, which are produced by two different pairs of Helmholtz and anti-Helmholtz coils, along theẑ axis, respectively. Thus, we have where µ i is the magnetic moment and b y is theŷ component of the magnetic field gradient. For |b y y| B F , which is fulfilled in our experimental conditions, U y,i can be approximated as where we neglect the terms in (b y y) /B F higher than the second order. If we define ω 2 i = |µ i | /(m i B F )b 2 y , with m i the atomic mass, the equation of motion is reduced tö where y i is the center-of-mass position of the atomic cloud andÿ i its acceleration. The solution of equation Eq. (C3) is where the constant factors A 1 and A 2 are determined by the initial conditions. For short times (t min{ω −1 i }), Eq. (C4) can be approximated by its Taylor series expansion as which corresponds to a uniformly accelerated motion, with constant acceleration α i ≡ y i (0) ω 2 i . Eq. (C5) is a good approximation of the observed center-of-mass dynamics, and, indeed, it has been used for fitting the experimental data discussed in the main text. Here, y i (0) is the initial distance of the atomic cloud from the maximum of the magnetic potential andẏ i (0) is the initial velocity. The value of y i (0) is not zero since the position of the atoms in the dipole trap is displaced from the symmetry axis of the magnetic field coils. We found that settingẏ i (0) = 0 does not substantially affect the fitting results. Further we verified that, withẏ i = 0, for evolution times up to 60 ms the error in y i made using Eq. (C5) instead of Eq. (C4) is less than 1 %.
We now treat the case of the bound state assuming that it is formed by N d i atoms in each component (i = 1, 2), then its equation of motion is given bÿ where y d is the center-of-mass of the droplet, and F i is the total (constant) force acting on each component. This force can also be expressed in terms of the center-of-mass acceleration α i that the i−th component would have in the absence of the other species, namely F i = N d i m i α i . Then, Eq. (C6) can be rewritten as from which it is straightforward to get Eq. (2).

Experimental fit of the unbound state dynamics
In the droplet regime, we generally observe a portion of atoms left over the bound state. We ascribe this effect to the presence of a thermal fraction in both atomic species and an excess of 87 Rb atoms with respect to the droplet atom ratio. Thus, to distinguish the bound and unbound parts, we fit the atomic density distribution with a two component Gaussian function. In Fig. 4 we have shown only the results obtained for the bound portion. Here we report also the evolution of the center-of-mass position y cm of the unbound fraction for δg < 0 finding the same behaviour of the weakly attractive mixture (Fig. 7).
Following the same strategy described in the main text, from Eq. (2), where now α K and α Rb are the accelerations of the unbound 41 K and 87 Rb components at δg < 0, we obtain the atom number ratio of the two species in the bound state: N d K /N d Rb = (0.8 ± 0.7), again consistent with the expected value [7].