Self-bound quantum droplets in atomic mixtures

Self-bound quantum droplets are a newly discovered phase in the context of ultracold atoms. In this work we report their experimental realization following the original proposal by Petrov [Phys. Rev. Lett. 115, 155302 (2015)], using an attractive bosonic mixture. In this system spherical droplets form due to the balance of competing attractive and repulsive forces, provided by the mean-field energy close to the collapse threshold and the first-order correction due to quantum fluctuations. Thanks to an optical levitating potential with negligible residual confinement we observe self-bound droplets in free space and we characterize the conditions for their formation as well as their equilibrium properties. This work sets the stage for future studies on quantum droplets, from the measurement of their peculiar excitation spectrum, to the exploration of their superfluid nature.

Very recently the stabilization mechanism generated by the LHY correction has been recognized responsible for the formation of a different class of self-bound quantum systems, i.e. dipolar droplets [4][5][6][7][8][9][10].While attractive mixtures create spherical droplets, in dipolar gases droplets are elongated along the dipole direction and thus strongly anisotropic.The different geometry, together with the different kind of interactions governing the stabilization, lead to important differences in the properties of the two objects and enriches the range of phenomena that can be explored.A primary example of this concerns the excitation spectrum, where self evaporation is expected to occur only in the mixture case, being strongly related to the droplet isotropic shape [11].
In this work we study experimentally the formation of spherical quantum droplets in a homonuclear bosonic mixture.Exploiting magnetic Feshbach resonances we tune the interatomic scattering lengths to reach the interaction regime where the mixture is predicted to be self-bound.We implement an optical levitating potential with negligible residual confinement along all directions, which allows to have long interrogation times and access the droplet properties in free space.We probe the mixture phase-diagram, proving the existence of a self-bound phase and identifying the critical conditions for its formation.We characterize the equilibrium properties of the droplet and we find a good agreement with the theoretical predictions in [1].Finally, we discuss the dynamics observed in the droplet formation and compare it to the result of a numerical simulation.
We create self-bound droplets using two hyperfine states of 39 K, namely |F = 1, m F = 0 (state 1) and |F = 1, m F = −1 (state 2).Feshbach resonances allow to tune the mutual contact interactions as represented in Fig. 1a as a function of the magnetic field B [12].The intra-species scattering lengths a 11 and a 22 are both positive, while the inter-species a 12 is negative.We define an effective scattering length for the mixture δa = a 12 + √ a 11 a 22 , which becomes negative for B < B c , setting the threshold for collapse in the usual MF picture [2].The stabilization effect of the LHY correction arXiv:1710.10890v1[cond-mat.quant-gas]30 Oct 2017 predicted in [1] appears exactly here.Contrary to the case of a single species [13], in a mixture of BECs the MF and LHY terms have a different dependence on the interparticle scattering lengths.While the MF energy, E M F , is proportional to |δa| and thus vanishes close to B c , the LHY correction, E LHY , scales with the average of a 11 and a 22 , thus becoming comparable to the MF term in this regime.Moreover, the two terms have a different dependence on the density n, since E M F ∝ n 2 while E LHY ∝ n 5/2 .This means that when the MF contribution becomes negative leading to an uncontrolled increase of density and eventually to collapse, the positive LHY term, having a steeper dependence on n, arrests the collapse and stabilizes the system.In this regime the mixture can be found in two different phases depending on the total atom number N .When N is larger than a critical number N c [1], the mixture forms a self-bound liquid-like droplet.Below that threshold, the kinetic energy overcomes the MF attraction and the system goes back into an expanding gas phase, labeled as LHY gas in the phase diagram of Fig. 1b.
We prepare a Bose-Einstein condensate (BEC) of 39 K atoms in state 2, in a crossed dipole trap, created by three red-detuned laser beams, with trapping frequencies ω x = 2π × 250(10) Hz, ω y = 2π × 240 (10) Hz and ω z = 2π × 280(10) Hz, along the axes sketched in Fig. 1c.A homogeneous magnetic field is used to tune the scattering lengths as in Fig. 1a.Starting with a BEC with up to 4 × 10 5 atoms, we ramp linearly the magnetic field in 20 ms to a desired target value and then we apply a radio-frequency pulse of 10 µs to transfer ∼ 50% of the atoms in state 1.In order to observe the subsequent evolution for sufficiently long times, remaining within the field of view of our imaging system, gravity compensation is required.The vertical position of a red-detuned elliptical laser beam is modulated in time with an acoustooptical modulator at a frequency of 4 kHz, such that the averaged potential experienced by the atoms provides a gradient opposite to gravity (red beam in the inset of Fig. 1c).A large waist on the horizontal direction (y) and a suitable time-modulation along the vertical direction (z) guarantee negligible residual curvatures on all directions [14].At the end of the RF pulse, we switch off the dipole traps and switch on the levitating potential to observe the expansion of the mixture.After a variable waiting time we record the density profile of the cloud via absorption imaging along the y direction [14].We fit it with a two-dimensional gaussian and measure the size along x and z as the half width at 1/e 2 .In order to probe the different phases expected for the mixture, we study the expansion of the cloud in three different points of the phase diagram of Fig. 1b.In Fig. 1c we report the average size σ = √ σ x σ z as a function of the expansion time t.The mixtures prepared with δa > 0 (orange diamonds) or δa < 0 and N < N c (light blue triangles) show the typical gas behavior: when released from the dipole trap they expand at a finite velocity.For δa < 0 and N > N c (blue circles), instead, the size of the cloud remains constant for 10 ms, proving the formation of a self-bound droplet.In order to characterize the droplet phase, we also perform measurements of the total atom number and of the relative population in state 1 and 2. After a variable expansion time, we perform a Stern-Gerlach separation of the two components, by applying a magnetic field gradient along the vertical direction z, so that we can count separately N 1 and N 2 .In Fig. 2 we report the evolution of the size σ, measured as in Fig. 1c, of the total atom number N = N 1 + N 2 and of the ratio N 1 /N 2 , for B = 56.45(1)G.The data in Fig. 2a show that the mixture is in the self-bound regime only up to a critical time t c , while afterwards it expands as a gas.We estimate the value of t c by fitting our data with a heuristic piecewise function composed by a flat line at short times, plus a linear growth after t c (dashed line in Fig. 2a).The behavior of σ(t) is explained by the evolution of N (t) in Fig. 2b.The atom number drops quite rapidly in the first 10 ms, so that the system follows the phase-diagram trajectory sketched in the inset of Fig. 2b: at a given expansion time the atom number goes below N c = N (t c ) and the system undergoes a droplet-to-gas transition.In order to un-derstand the loss dynamics observed in Fig. 2b, we have to consider two different effects.The first one is threebody recombination, which causes strong losses, mainly in state 1 [14].The second one is related to a stabilization mechanism of the mixture: the droplet forms with a specific population imbalance, N 1 /N 2 = a 22 /a 11 [1].It can only bear a deviation from that value which corresponds to an increase of the atom number in one of the two components δN i /N i ∼ |δa|/a ii , i = 1, 2. Any excess of atoms beyond this threshold is not bound to the droplet and expands as a gas.Combining these two mechanisms, we can interpret the behavior observed in Fig. 2b-c.In the first 4 ms the ratio N 1 /N 2 decreases, meaning that the mixture is mainly losing atoms in state 1.This is due both to the stronger three-body losses in that component and to the initial ratio N 1 /N 2 ∼ 1 being significantly larger than its equilibrium value (green area in Fig. 2c).The droplet thus expels the exceeding atoms in state 1, which expand away from it and get effectively lost during the Stern-Gerlach sequence.When the mixture reaches the equilibrium population imbalance, the ratio N 1 /N 2 stabilizes, meaning that three-body losses in state 1 are accompanied by a release of atoms in state 2 into the unbound fraction.This is also compatible with the observed bimodal density profiles in the measurements of Fig. 2a, with a dense central part constant in size, corresponding to the droplet, surrounded by low density tails corresponding to the unbound atoms.We fit the profiles with the sum of two gaussians and associate the size of the droplet to the width of the central one (inset of Fig. 2a).
We repeat the measurements of Fig. 2 for different values of the magnetic field B. In Fig. 3 we report the results for the critical number N c and for the size σ and ratio N 1 /N 2 also measured at t c .The size reported in Fig. 3b corresponds to the geometrical average of the measured σ z and σ x , plotted together with the aspect ratio σ x /σ z .We compare these experimental results with the theoretical predictions obtained from [1] and we find a good agreement within our error bars over the whole magnetic field range we explored.
It is important to stress that the procedure we use to form the droplet starting from a single component BEC is a non-adiabatic process.Moreover, due to significant three-body losses in the system, the droplet has a limited lifetime.For this reason we decided to start from a BEC in an almost spherical confinement, with in-trap size close to the droplet one for B < 56.5 G, i.e. ∼ 2 µm.Minimizing the initial dynamics, we could observe the equilibrium configuration of the droplet, as visible from the measurements in Fig. 2, where the size of the cloud is essentially constant up to t c .However, for the two largest values of the magnetic field we considered, the evolution of σ(t) is qualitatively different (Fig. 4).In these cases the equilibrium size of the droplet, significantly larger than the BEC size, is reached after an ini-0.7 0.5 FIG.3: Critical atom number Nc (a) and equilibrium properties σ (b) and N1/N2 (c) of the droplet as a function of the magnetic field B. In b) we also report the aspect ratio σx/σz (diamonds).The curves in a) and b) correspond to the theoretical predictions for the metastable (dashed) and stable (solid) self-bound solutions [1].In c) the theoretical curves are obtained as in Fig. 2c, considering the equilibrium value of N1/N2 (solid line) and the allowed deviations (green area).The vertical error bars include the uncertainty coming from the determination of tc and the statistical uncertainty.The horizontal ones are due to the uncertainty in the magnetic All error bars correspond one standard deviation.
tial expansion.After transient, the size is constant for a given time and then it increases again when the atom number drops below N c .In this case we perform a different heuristic fit to deduce the critical time t c , that includes an additional linear growth of σ at early times (blue dashed line in Fig. 4).In order to verify that the observation timescale is sufficient to reach equilibrium, we perform a simulation of the droplet dynamics in our experimental conditions [14].We integrate numerically a system of two generalized Gross-Pitaevskii equations, which include first order quantum corrections in the local chemical potential [9] via the two-species LHY term discussed in [1].In Fig. 4 we compare the experimental data for σ(t) at B = 56.45(1)G and B = 56.64(1)G with the result of the simulations.In order to capture the stabilization of the size at long times, three-body losses are not included in the numerical simulations.For the smallest value of B, the size slightly oscillates around its equilibrium value and finally stabilizes there, consistently with the plateau observed in the experiment.At large B, the equilibration time of the simulation is compatible with the initial transient of the experimental data, justifying the assumption that the central plateau in the evolution of σ corresponds to the equilibrium size of the droplet.The study of this dynamics allows to add a further consideration.For B = 56.64(1)G, even if the sample is prepared out of equilibrium, no oscillation of the size is visible, i.e. there is no signature of excitations.In the simulation this is associated with an expulsion of atoms from the droplet to the expanding low density tails.This is compatible with the predicted self evaporation, since discrete excitation modes are allowed only for N > 10 6 at this magnetic field [1].In our experiment it was difficult to distinguish a self-evaporation effect in the losses dynamics discussed above, but the observed stabilization of the size is a first hint of the occurrence of this phenomenon.
FIG. 4: Evolution of the size σ for B = 56.45(1)G (circles) and B = 56.64(1)G (diamonds).We compare the experimental data with the result of the numerical simulations described in the text (solid lines).The dashed lines correspond to the heuristic fits performed to estimate the critical time tc in the two regimes.

In conclusion, our work represents the first experimental observation of self-bound droplets in an atomic
Bose-Bose mixture in free space.The measurement of the size and of the critical atom number as a function of the MF attraction provides a confirmation of the theoretical model proposed by Petrov [1].The dynamics observed in the droplet formation indicates the main perspective of this work, which is the study of self evaporation, providing a quantitative evidence for it and probing how it disappears by increasing the atom number, or by adding an external confinement to induce an anisotropic geometry.Further studies may include the effect of a coherent coupling between the two atomic species in the droplet, exploring how its equilibrium properties are modified [15].It will also be interesting to investigate differences and analogies with dipolar droplets and different self-bound quantum fluids like helium clusters [16,17], probing the superfluid behavior [18,19] and comparing the excitation spectra [20].Finally our work triggers the possibility to investigate the formation of self-bound droplets in different atomic mixtures, possibly with reduced three body loss rates and longer lifetimes [21][22][23][24].
While completing this work, we became aware of a related work, reporting the observation of droplets in atomic mixtures in a one-dimensional lattice: C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell, arXiv:1708.07806(2017).
We acknowledge insightful discussions with D. Petrov, A. Simoni, S. Stringari, L. Tarruell, L. Tanzi and our colleagues of the Quantum Degenerate Gases group at LENS.We thank G. Spagnolli for assistance in the early stages of the experiment.This work was supported by the ERC Starting Grant AISENS No. 258325 and by EC-H2020 Grant QUIC No. 641122.

Optical levitating potential
The optical potential used to levitate the atoms against gravity is created by a single elliptical beam, whose vertical position is modulated in time.The laser wavelength is λ = 1064 nm and the waists of the beam along z and y (axes defined as in Fig. 1c in the paper) are w z = 23(3)µm and w y = 985(8)µm.While the intensity of the beam is fixed, its vertical position z 0 changes periodically in time.The modulation function is defined on a single period T as where A = 50 µm is the modulation amplitude and C = 7 µm is an offset chosen to place the zero-curvature point at z = 0, that corresponds to the position of the atoms.The beam is moved along the vertical direction by means of an acousto-optic modulator whose radiofrequency is modulated in time, thus changing the diffraction angle.The time-averaged potential along z is then where V 0 is the potential amplitude.The profile of V (z) is reported in Fig. 5a).By adjusting the power in the laser beam to 2.4 W, the gradient produced in the central part of V (z) is opposite to gravity (Fig. 5b) and thus allows to levitate the atoms and access long observation times.We use a modulation frequency ν mod = 1/T = 4 kHz, which is fast enough that the atoms are only affected by the averaged potential.The residual curvatures created by the levitating potential are weak along all directions.Along z, the calculated curvature is zero at z = 0, is positive above and negative below, reaching |ω z | ∼ 2π × 20 Hz at a distance of ∼ 15µm from the center (Fig. 5c).Along the other two directions, the curvature is positive but very weak: ω x = 2π × 2.2 Hz and ω y = 2π × 7 Hz.While the confinement along these two axes is well controlled by measuring the power and waists of the laser beam, we would like to perform an in-situ calibration of the curvature along z, which is more sensitive to misalignments.To perform a precise experimental measurement of ω z is actually quite difficult, since the shape of the potential is far from being harmonic and the estimated curvatures are very weak.We perform two different measurements: in the first one we provide an upper bound to the local curvature in z = 0, while in the second one we estimate the global effect produced by the levitating potential on a larger scale.The measurement of the local curvature is performed as follows: we measure the trapping frequency of a dipole trap aligned on the z = 0 position, first in presence of a magnetic levitating potential (ω magn ) and then with the optical levitation (ω opt ).The curvature generated by the optical levitation is then estimated as opt − ω 2 magn .We find ω magn = 2π × 129(1) Hz and ω opt = 2π × 128(1) Hz, which are compatible within the error bars.From this we can only deduce an upper bound to ω z coming from the error bars of the measurement, i.e. |ω z | < 2π × 16 Hz.In the second test, we compare the evolution of a Bose-Einstein condensate in |1, −1 at 7.5 a 0 in the levitating potential and in free space.In Fig. 6 we report the measured Thomas-Fermi radius of the cloud along z, as a function of the expansion time.The two datasets correspond to the expansion in free-fall (purple) and in the levitating potential (orange).We compare them to the theoretical curves corresponding to the expansion of the BEC in free space and in a harmonic confining potential along z.When the vertical size R y increases above ∼ 30µm, we see a deviation from the free-space expansion, which is compatible with the evolution of the BEC with ω z = 2π × 12 Hz.This measurement is mostly sensitive to the curvature at large distances, so that it also provides only an upper bound to the effective confinement experienced by the droplets around z = 0. Finally, note that the sizes of the droplets we observe are always isotropic along the two measured directions, being σ x and σ z compatible within their error bars.Along x the estimated confinement is very weak, ω x = 2π × 2.2 Hz.The fact that we do not observe any deformation of the droplet from its predicted spherical geometry let us conclude that the effect of the levitating potential is negligible also on the vertical direction.μ FIG.6: Expansion of a BEC in the levitating potential vs free space.We measure the Thomas-Fermi radius along the vertical direction Rz after releasing a BEC from its dipole traps.The purple points correspond to the free-space evolution, while the orange ones to the expansion in the levitating potential.The theoretical curves are the result of the evolution calculated from the Gross-Pitaevskii equation in free-space (solid purple) and with a vertical confinement ωz = 2π × 12 Hz (dashed orange).

Absorption imaging of dense clouds
The droplets we observe have typical densities n ∼ 10 15 cm −3 .Such high density clouds are difficult to probe with standard absorption imaging, due to saturation effects.In order to image the atoms in the F=1 state of 39 K, we need to use light at two different wavelengths: the first one, called "repumper" light, excites the atoms in the F = 2 state; the second one makes the atoms cycle on the F = 2 → F = 3 transition, which provides the needed absorption.In order to avoid saturation effects and properly record the density profile of the droplet we use the following technique.We detune the repumper light from resonance by ∼ 100 MHz.The optical density of the droplet for such detuned laser is smaller than one.In this way we transfer in F = 2 a cloud with a smaller atom number, but with the same spatial density profile of the droplet.We can then properly measure the size of the cloud, avoiding saturation effects.
The resolution of our imaging system is δ res = 1.5(1) µm, measured on a non-interacting BEC loaded in a dipole trap with a trapping frequency of 300 Hz.Droplets are typically few µm large, so we need to take into account the finite imaging resolution to properly estimate their size.The data for σ reported in the paper are all deduced by deconvolving the measured size with the imaging resolution as σ = σ 2 meas − δ 2 res .During the imaging sequence all magnetic fields need to be off to avoid that the atomic transitions are shifted by the Zeeman effect.The minimum time to switch off the magnetic field B is ∼ 500µs.During this time, while the field goes to zero, the interatomic scattering lengths change.This could affect the measurement of the size of the cloud, especially when densities are very high, i.e. in the droplet or in the gas phase at short expansion times.We want to estimate this effect, by calibrating it on a mixture whose size can be predicted by a wellestablished theory, i.e. in the stable mean-field regime.Considering the data reported in Fig. 1c of the paper, we can compare the size measured at t = 0 for the mixture at δa > 0 (orange diamonds) with the theoretical prediction obtained with the coupled Gross-Pitaevskii equations for the double BEC at the same value of B. We find that the two results are compatible within the experimental error bars: σ meas = 2.4(1)µm and σ theor = 2.3µm.We conclude that the switching off of the field has a negligible effect.

Three-body losses
We perform measurements of three-body loss rates on single-component BECs in state 1 and 2 on the magnetic field range probed in the experiment.We analyze the data according to the model in [25] and we find K 111 /3! = 7 × 10 −28 cm 6 /s and K 222 /3! = 1 × 10 −29 cm 6 /s, with an uncertainty on the absolute values of a factor 2. By performing tests on thermal mixtures, we find that three-body recombination in state 1 is the dominant source of losses.Quantitative estimations for all the different loss rates in the mixture are quite demanding and go beyond the purposes of this work.Detailed measurements would be anyway interesting for future studies, in order to quantitatively characterize the loss dynamics observed in Fig. 2b-c of the paper.

Numerical simulations
The Bose-Bose mixture is described by two coupled Gross-Pitaevskii (GP) equations which include the twospecies LHY term discussed in [1].This is achieved by means of a local density approximation, as recently considered for the description of the formation of macroscopic quantum droplets in dipolar condensates [7,9,26].Namely, the evolution of the system is governed by two equations of the form (i = 1, 2) where µ i = δE/δn i is the local equilibrium chemical potential for the species i (see e.g.[27]).We consider the case of equal masses (m 1 = m 2 = m) and define g ij = 4πa ij /m the coupling constants and n i the density of species i.In the present case, the total energy E consists of two contributions: the usual mean-field GP term .
When the two species are prepared with the same density distribution, modulo a scale factor α (corresponding to different atom numbers, namely n 2 = αn 1 ≡ αn) and the system is close to the onset of the mean-field collapse regime, g 12 + √ g 11 g 22 < 0 (with g ii > 0), one can further simplify the above expression by setting g 2 12 = g 11 g 22 and neglecting small finite-δg corrections [1] E LHY = 8m Then, Eq. (3) (with ( 7), ( 8)) is solved by initially loading the condensate in state 2 in the ground state of the external trapping V ext , corresponding to the crossed dipole trap of the experiment (by means of a standard imaginary time evolution [28]), and then transferring half of the population to state 1, at t = 0.The subsequent evolution is obtained by solving Eq. ( 3) by means of a split-step method which uses fast Fourier transforms (see e.g.[29]).

FIG. 1
FIG. 1: a) Intra and inter-species scattering lengths between the hyperfine states |1, 0 (state 1) and |1, −1 (state 2) of 39 K, tuned by an external magnetic field B via Feshbach resonances.The resulting MF energy of the mixture is proportional to the effective scattering length δa, which becomes negative at Bc = 56.85G. b) Phase diagram for the mixture as a function of the atom number N and of the magnetic field B. c) Evolution of the cloud in free space for three different points of the phase diagram in b).The upper rows show the difference between the evolution of the density profiles in the gas and droplet phases.Inset: Schematic representation of the geometry of the experiment.

μFIG. 2 :
FIG. 2: Time evolution of σ (a), N (b) and N1/N2 (c) in the droplet phase at B = 56.45(1)G.The inset in a) reports the density profile of the droplet close to tc, together with the fitted bimodal function.In the inset in b) we draw a sketch of the trajectory followed by the system in the mixture phasediagram during the time evolution, due to losses.The dashed line in a) represents the heuristic fit described in the text and used to identify the critical time tc.In c) the solid line represents the theoretical equilibrium value N1/N2 = a22/a11 and the green area between the dashed lines includes the allowed deviations δNi/Ni ∼ |δa|/aii.The error bars represent the statistical uncertainty and correspond to one standard deviation.

FIG. 5 :
FIG.5: a) Profile of the levitating potential V (z).b) Profile of V (z) around z = 0 compared with gravity.c) Curvature produced by the levitating potential along the vertical direction.