Dilute dipolar quantum droplets beyond the extended Gross-Pitaevskii equation

Dipolar quantum droplets are exotic quantum objects that are self-bound due to the subtle balance of attraction, repulsion, and quantum correlations. Here we present a systematic study of the critical atom number of these self-bound droplets, comparing the experimental results with extended mean-field Gross-Pitaevskii equation and quantum Monte Carlo simulations of the dilute system. The respective theoretical predictions differ, questioning the validity of the current theoretical state-of-the-art description of quantum droplets within the extended GrossPitaevskii equation framework and indicating that correlations in the system are significant. Furthermore, we show that our system can serve as a sensitive testing ground for many-body theories in the near future.


I. INTRODUCTION
For systems with competing interactions, quantum fluctuations can stabilize an otherwise collapsing system [1]. The balance of attraction and repulsion in these systems means that they share properties with liquids, despite being orders of magnitude more dilute. These quantum droplets were experimentally discovered by driving a dipolar Bose-Einstein condensate (BEC) into the strongly dipolar regime. However, instead of collapsing, the system formed stable droplets [2]. Since their discovery, many properties of dipolar quantum droplets have been observed and compared to the predictions of their current state-of-the-art theoretical description, the extended Gross-Pitaevskii equation (EGPE). These predictions include the stabilization of the droplets typically explained by the Lee-Huang-Yang (LHY) correction of the mean-field energy [3,4], their self-bound nature [5], collective modes [4,6], the emergence of striped states in confined geometries [7], and the existence of arrays of phase coherent droplets with transient supersolid properties [8][9][10]. In addition to dipolar systems, quantum droplets have also been observed in Bose-Bose mixtures [11][12][13][14].
Here we systematically study the critical atom number that is necessary to form the liquidlike droplet state, extending previous results with 164 Dy [5] to an order of magnitude higher atom number. Comparing the measured critical atom number to results obtained from numerical simulations of the EGPE seems to indicate a systematic shift to lower values. Experimental discrepancies compared to the EGPE predictions have also been observed in other related systems [10,11,15].
Motivated by this observation, we present a theoretical approach that goes beyond the framework of the EGPE. We solve the dilute many-body system using quantum Monte Carlo (QMC) calculations, in particular the path-integral ground-state (PIGS) method [16,17]. With this we can extract the critical atom number of self-bound droplets, in good agreement with the experimental measurements. As a meanfield theory, the EGPE framework is limited to the usage of the local density approximation, as well as to the perturbative regime at small gas parameters [18]. In contrast, our PIGS calculations intrinsically include effects due to the finite system size and of the finite interaction range, as well as particle correlations and quantum fluctuations. With our method we directly have access to the correlations in the system, which we use to extract the spatial pair correlation function, as well as the condensate depletion, which is increased compared to the prediction of Bogoliubov theory. These results suggest that in the density regime relevant for quantum droplets, the state-of-the-art EGPE framework is not able to reproduce all observable properties of these quantum droplets.

II. EXPERIMENT
For the experiments we use 162 Dy with a magnetic dipole moment μ = 9.93μ B , where μ B is the Bohr magneton. To quantify the relative strength of the dipole-dipole interaction compared to the contact interaction, we define the relative dipolar strength ε dd = a dd /a s in terms of the scattering length a s and the dipolar length a dd = μ 0 μ 2 m/12πh 2 . Hereh is the reduced Planck constant, μ 0 is the vacuum permeability, and m is the atomic mass. The dipolar length is different for the two studied isotopes and has a value of 129a 0 for 162 Dy and 131a 0 for 164 Dy, with the Bohr radius a 0 . Like all lanthanide atoms, 162 Dy exhibits a rich spectrum of Feshbach resonances that can be used to control the strength of the short-range contact interaction [19][20][21][22][23]. Here we use a specific double resonance at around 5.1 G [20] that, together with the high background scattering length a bg = 140(20)a 0 of 162 Dy [24][25][26], allows us to tune the scattering length a s from a contact-dominated sample (away from the resonances), to a dipolar-dominated sample closer to the zero crossing of the scattering length. Since the background scattering length of 162 Dy has so far not been determined to high precision, all calculated scattering lengths exhibit a systematic uncertainty of 15%.
While a BEC is a gaseous state, which means that it freely expands in the absence of an external trap, the droplet state is self-bound due to its intrinsic interactions [5,[27][28][29]. This self-bound character has been experimentally observed for 164 Dy [5], as well as for the related quantum droplets in Bose-Bose mixtures [11,12]. In order to probe this feature experimentally, we initially prepare a quasipure BEC with 4.5(3) × 10 4 162 Dy atoms and then closely follow the procedure presented in [5] to create a single self-bound droplet. After a variable evolution time we intentionally evaporate the droplet and subsequently image it after an expansion of either 8 or 30 ms, depending on the atom number. With this approach we observe atom number decay curves that settle at a constant atom number, the critical atom number of a self-bound droplet [5].
We want to compare our results to numerical EGPE and QMC simulations, as well as previous measurements with 164 Dy published in [5]. The two bosonic isotopes of dysprosium, 162 Dy and 164 Dy, only differ in their mass, which manifests in a shift of (m 164 − m 162 )/m 162 = 1.2% of the dipolar length, and thus on the scattering length axis. This effect is smaller than our uncertainty in this quantity, and we can thus safely neglect it. For the comparison with the 164 Dy measurements we use the Feshbach resonances discussed in [2,5] together with the latest measurement of the background scattering length in the droplet state a bg,164 = 69(4)a 0 [6]. With this we observe systematically lower critical atom numbers compared to the EGPE simulations, which could be accounted for by the uncertainty of the scattering length calibration. For our 162 Dy results, we first use the literature value of the background scattering length a bg,162 = 140(20)a 0 [24][25][26], together with our measured parameters of the Feshbach resonances, leading to a similar systematic shift to lower atom numbers as in the 164 Dy data. Again, this shift could be explained by the even larger uncertainty in the 162 Dy scattering length.
Next we use the sensitive scaling of the critical atom number with respect to the scattering length in order to extract the ratio of the two background scattering lengths a bg,162 /a bg,164 , free of the systematic uncertainties of their respective measurements. To put this scaling into context, we note that the current uncertainty of the background scattering length of 162 Dy is about 20a 0 . The critical atom number changes by an order of magnitude over a comparable range of the FIG. 1. Critical atom number of a self-bound dipolar quantum droplet for 162 Dy (blue circles) and 164 Dy [5] (black squares). We extract the critical atom number by analyzing the atom number decay curves. The theoretical boundary of the phases is obtained from numerical EGPE simulations. The red dashed and dash-dotted lines show the corresponding boundary as expected from an increased effective dipolar length due to a finite collision energy of 30-50 nK and 100 nK [30,31]. The red triangles show the results obtained by QMC simulations, with the error bars chosen to cover the uncertainties of both the statistical error and the nonuniversality. See the text for more information.
scattering length, showing that our measurement constitutes a very sensitive probe of the scattering length. Here we use the critical atom numbers derived from our EGPE simulations to extract the ratio a bg,162 /a bg,164 . Naturally, the same procedure can be done using any quantum many-body theory able to predict critical atom numbers for a self-bound droplet.
Starting from the latest measurement of the background scattering length of 164 Dy [6], we first shift the EGPE critical number curve 1 in order to minimize the difference from the experimental data for 164 Dy. From this point on we optimize the background scattering length a bg,162 to minimize the residual of our measurements with respect to this shifted theory curve. Taking into account the residual systematic uncertainty of the background scattering length of 164 Dy [6], we end up with a bg,162 = 140(7)a 0 , in agreement with the literature value [24][25][26], but with a significantly reduced uncertainty. More importantly, since it is independent of this residual systematic uncertainty, we can extract the ratio of the two background scattering lengths. Comparing the two isotopes, we find a ratio of a bg,162 /a bg,164 = 2.03(6) of their respective background scattering lengths.
This way we calibrate the scattering length of 162 Dy and show a summary of all measured critical atom numbers in Fig. 1. The experimental atom number uncertainties are chosen to cover the determination of N crit , with an additional 10% uncertainty of the imaging. The uncertainty of the scattering lengths is based on our experimental magnetic field stability of ∼2 mG, the knowledge of the used Feshbach resonances, and the uncertainty of the respective background scattering length.
Experimentally, we seem to find systematically lower critical atom numbers compared to EGPE predictions. The straightforward experimental source would be a shift of approximately 6a 0 in the scattering length axis. This is at the edge of the error bars for 164 Dy and within the error bars of the independent calibration for 162 Dy. Note that our calibration of the 162 Dy scattering length is consistent with our recent experimental observation of supersolid properties in an array of 162 Dy quantum droplets [8], which agrees with theoretical calculations in a narrow window of the scattering length. In order to make the presented measurements a sensitive benchmark for quantum many-body theories, an independent and precise measurement of the scattering length a s would be required for the two isotopes. Another possible experimental explanation would be a systematic uncertainty in the atom number determination, which we rule out by performing measurements with independent imaging techniques, all resulting in similar values within the quoted error bars.

III. THEORY
We now focus on possible theoretical explanations of systematically lower critical atom numbers of the self-bound state. Within the EGPE framework, it was shown in [30] that it is possible to account for finite-temperature effects in the twobody scattering problem by effectively enhancing the dipolar length a dd . In Fig. 1 we also show two theoretical curves for an enhanced dipolar length corresponding to a collision energy of 30-50 nK [31] (dashed red line), and 100 nK (dash-dotted red line). The resulting shift for the critical atom number is in good agreement with the experimental results, suggesting that such an effective enhancement of the dipolar length might play an important role in dipolar scattering at finite collision energies. In the future, one may resort to spectroscopic measurements on embedded impurities [32] in order to directly measure the temperature of these quantum droplets and therefore determine the actual strengths of this proposed correction.
Next we present an approach that goes beyond the current state-of-the-art EGPE framework. For this we solve the dilute many-body system at zero temperature using QMC simulations. For weakly interacting systems, the results obtained in QMC are in good agreement with mean-field predictions [33][34][35][36]. However, for strongly correlated systems, e.g., liquid helium [37,38], mean-field fails, while QMC still provides reliable predictions. To see whether correlations influence the properties of quantum droplets, we now focus on QMC simulations [34,39,40]. In particular, we use the PIGS method to determine the ground-state properties of ensembles of 162 Dy atoms at zero temperature. Although computationally extensive, this method intrinsically includes finite-range effects present in a more realistic description of the atom-atom interaction, the finite system size, and a correct description of correlations and quantum fluctuations. Compared to this, the EGPE approach relies on the local density approximation for the quantum fluctuations described by the LHY term [41,42]. Our approach is however limited to the usage of an effective Hamiltonian without bound states to describe the interaction between particles. This effective Hamiltonian includes the dipolar interaction and a central two-body interaction, with a repulsive core and a realistic C 6 coefficient [43]. In order to study whether the system is universal, we use different model potentials, with the respective parameters fixed to adjust the zero-energy s-wave scattering length to the range of the experimentally measured values. This is accomplished by solving the Lippmann-Schwinger equation associated with the T matrix of the full interaction.
First, we study the influence of correlations compared to the Bogoliubov prediction used to derive the LHY term [41,42] in the EGPE framework. Therefore, we simulate the equivalent homogeneous bulk system with a density of n = 5.88 × 10 21 m −3 for two different scattering lengths a s = 60a 0 and 90a 0 and extract the pair correlation function g(r). The density is chosen to correspond to the equilibrium density of a saturated droplet at a s = 60a 0 . Due to the anisotropy of the dipolar interaction, the correlation function g(r) shown in Fig. 2(a) depends on the direction with respect to the polarization axis. Perpendicular to the polarization axis, the pair correlation function is a monotonic function of the distance that resembles the one of a weakly interacting system. On the other hand, along the polarization direction it shows signatures of local ordering, as highlighted by a broad peak at short distances. In both directions, the length scale of the repulsion at short distances is caused by the repulsive core of the used two-body model potential.

033088-3
Another property that directly measures the strength of correlations is the quantum depletion 1 − n c /n, with the condensate fraction n c /n. The condensate fraction equals one for an ideal Bose gas at zero temperature and decreases when correlations are enhanced. In strongly interacting liquid helium the condensate fraction is typically below 10% [44]. For typical weakly interacting ultracold-atom experiments it is around 99%, while by largely increasing the scattering length, condensate fractions as low as about 90% [45] have been realized. In Fig. 2(b) we show the comparison between the PIGS prediction and the derived quantum depletion 1 − n c /n within the Bogoliubov theory for the weakly interacting Bose gas with no dipolar interaction, as well as with dipolar interaction [41,42], for the corresponding parameters of our quantum droplets. As it can be seen, the dipolar interaction leads to stronger correlations, as well as to a larger overall depletion of the condensate in the range of densities relevant for saturated quantum droplets. Compared to this, the PIGS results show that the effect of correlations is even stronger.
Now we turn to the study of the full finite system. We analyze realizations of the system with a different number of particles for a given s-wave scattering length and compute the ground-state energy in each case. Like for the EGPE simulations, we identify a self-bound droplet as a system with negative energy in the absence of an external trapping confinement. The self-bound droplets predicted by our PIGS calculations differ from those obtained in the EGPE approximation in the overall density profiles. This difference can be attributed to the presence of correlations, which can be quantified by looking at the two-body properties shown in Fig. 2. Distinguishing the two proposed theories experimentally may be possible in the future by looking at the density distributions of self-bound droplets. As in the experiments and the EGPE simulations, we find that there is a critical atom number below which the system ceases to be self-bound. Close to the critical atom number our PIGS calculations result in a lower peak density than predicted within the EGPE framework. Note that there is a non-negligible dependence of the critical atom number on the exact model of the two-body potential we employ, which indicates that the problem is nonuniversal in terms of the scattering length. The resulting critical atom numbers for several values of the scattering length are shown with red triangles in Fig. 1, with the error bars chosen to take into account the effect of the nonuniversality according to the analyzed model potentials, as well as the statistical errors of the simulations. The obtained critical atom numbers are always below the EGPE predictions and in good agreement with the experimental measurements. The improvement of the PIGS predictions with respect to the EGPE results points to the relevance of finite-range effects which enhance quantum correlations, similar to dilute Bose mixtures [46].

IV. CONCLUSION
We have systematically studied the critical atom number for a self-bound dipolar quantum droplet experimentally and have used these measurements to establish a comparison between the current state-of-the-art EGPE description and QMC simulations. Compared to EGPE results, we observe indications of a systematic shift of the experimentally measured critical atom numbers to lower values. Those values are nevertheless well reproduced by a zero-temperature QMC simulation based on the PIGS algorithm. In contrast to the EGPE, our PIGS calculations include finite-range effects in the interaction as well as finite-size effects, together with correlations and quantum fluctuations. This is used to extract the spatial pair correlation function and the condensate depletion, showing that in the relevant density regime of quantum droplets, correlations are enhanced and need to be included in a realistic description of the problem. In this way, our PIGS results indicate that correlations are beyond what a zero-temperature modified mean-field theory can capture. Alternatively, the inclusion of finite-temperature effects in the EGPE framework (through an effective renormalization of the dipolar interaction strength associated to finite collision energies in the two-body scattering problem) can also reproduce the critical number data, although the prediction is strongly dependent on the temperature used in the calculation. All in all, our results call into question the validity of the EGPE framework to fully describe the quantum droplets. After this we typically change the trap and/or the magnetic field for the actual experiments, as described in further detail for the performed experiments below.
Along theẑ axis we have a microscope objective allowing for in situ imaging with 1-μm resolution. We can use this microscope for far-detuned phase-contrast imaging as well as resonant absorption imaging. Both techniques, as well as an independent time-of-flight imaging along theŷ direction, result in similar atom numbers. The microscope can also be used to focus an additional optical dipole trap (λ = 532 nm) that has a calculated beam waist of ∼22 μm, to change the trap aspect ratio from the oblate CODT to a spherical or even prolate trap. 033088-4

Self-bound droplet measurements
To measure the critical atom number for self-bound droplets shown in Fig. 1 of the main text, as well as the measurement of the expansion velocity in Fig. 5, we apply a magnetic field gradient along theẑ direction after the preparation of the BEC. This applied gradient exactly compensates gravitational forces and leads to a shift of the magnetic field by −428 mG, which we compensate by ramping up the amplitude of the constant magnetic offset field at the same time. After this we reshape the trap within 20 ms into a spherical trap with a mean trap frequency ofω = 96(4) Hz. To do this we change the trap aspect ratio from λ = ω z /ω r = 96(2) Hz/25(2) Hz = 3.8 to λ = 98(2) Hz/94(4) Hz = 1.05, by applying an additional optical dipole trap along theẑ direction. Here ω z (ω r ) is the trapping frequency along (perpendicular to) the magnetic field direction. For such a spherical trap geometry the regular BEC at large scattering lengths and droplet phase at low scattering lengths are connected by a continuous crossover [3,28,47], as can also be seen in the phase diagram shown in Fig. 4. To reach the droplet state we ramp the magnetic field amplitude within 20 ms closer to the zero-crossing Feshbach resonance to reduce the scattering length. Subsequently, we let the atoms evolve for 10 ms to allow the droplets to form. Then we slowly switch off the trap by ramping down the intensities of the laser beams to ∼5% of their initial values, keeping a nearly constant trap aspect ratio. We then suddenly turn off the trap and levitate the cloud for various times before subsequently imaging the density distribution.
For the measurements of the expansion velocity shown in Fig. 5, we follow the expansion up to t TOF = 20 ms and subsequently image the atoms using far-detuned phase-contrast imaging.
For the measurements of the critical atom number, as shown in Fig. 1, we levitate the atomic cloud for a variable evolution time and then intentionally evaporate the droplets [5] by ramping up the magnetic field within 100 μs to B ≈ 6.0 G. At this field the ground state of the system is an expanding BEC and therefore the droplet evaporates back to the gaseous phase and expands freely. After an expansion of either 8 or 30 ms, depending on the atom number, we image the atomic cloud with resonant absorption imaging. With this we can observe atom number decay curves, which settle to a constant value, the critical atom number of the self-bound droplet. The shown scattering length range in Fig. 1 for 162 Dy corresponds to magnetic fields in the range from B = 5.293 G to B = 5.249 G.

APPENDIX B: FESHBACH RESONANCES AND THREE-BODY LOSS COEFFICIENT
The complicated spectrum of Feshbach resonances for lanthanide atoms allows one to control the short-range contact interaction by tuning the amplitude of the magnetic field close to one of the many resonances. In this work we use a specific double resonance (see Fig. 3) at a field of B 1 = 5.126(1) G and B 2 = 5.209(1) G with widths of B 1 = 35(1) mG and B 2 = 12(1) mG, respectively. In order to calibrate the magnetic field amplitude we use radio-frequency spectroscopy to get the absolute value of the magnetic field amplitude from the The two mentioned resonances, as well as a third resonance, can be seen in Fig. 3. The third resonance at 5.273 G has a width of only ∼1 mG. This resonance seems to vanish at even lower temperatures and therefore does not influence the atoms in a BEC. Because of this we do not include this narrow resonance in our consideration of the scattering length. Additionally, there is a broader resonance at B 3 = 21.95(5) G with a width of B 3 = 2.4(8) G [48], which still has a small effect on the scattering length in the magnetic field range considered in this work. Using the mentioned resonances, we can calculate the scattering length as a function of the magnetic field, with only the background scattering length as a free parameter [ Fig. 3(c)].
In addition to the atom-loss spectroscopy used to extract the position and the width of the Feshbach resonances, we also need to measure the three-body loss rate L 3 and check whether the observed shorter lifetimes of the 162 Dy droplets can be 033088-5 understood from the theory. To measure L 3 we prepare a thermal cloud at about 200 nK and then ramp up the magnetic gradient and again compensate the magnetic field shift by ramping up the offset field at the same time. After this we compress the atomic cloud by ramping up the powers in the CODT within 25 ms such that we reach a trap with trap frequencies of 83(4), 299(3), and 434 (2) Hz. Then we change the magnetic field to its variable final value within 3 ms and subsequently let the atoms evolve for up to 1 s. We then image the atoms after a 10-ms time of flight and fit the atom number N and the temperature T in order to extract the three-body loss rate to the differential equations, similarly to the methods described in [49], In these equations α is the two-body loss rate, which we measured in a dilute thermal cloud to be ∼20 s, and γ is connected to the three-body coefficient Hereω is the mean trap frequency, m is the mass of the atoms, and k B is the Boltzmann constant. The temperature of the sample increases due to the losses, because two of the colliding atoms can form a molecule and the third can gain the binding energy.
What we see in Fig. 3(d) is a three-body loss coefficient L 3 = 8 × 10 −41 m 6 /s away from the Feshbach resonances, which increases the closer we get to the resonance. Compared to the thermal cloud, the three-body loss coefficient in the BEC is decreased by a factor of 6 [50], leading to L 3 = 1.33 × 10 −41 m 6 /s away from the resonances. In the droplet region (indicated by the blue area in Fig. 3) we have an L 3 that is a factor of 4 (a s ∼ 105a 0 ) or even 15 (a s ∼ 80a 0 ) times larger than the value far away from the resonance. In this range we also see a small peak due to the very narrow Feshbach resonance located there, which however we do not observe anymore for lower temperatures. The observed large increase of L 3 explains the shorter lifetime that we observe for the self-bound droplets for lower scattering lengths.
For the comparison to the 164 Dy data from [5], we convert the given magnetic field amplitudes using the Feshbach resonances at B 164,1 = 7.117(3) G with a width of B 164,1 = 51(15) mG and a second resonance at B 164,2 = 5.1(1) G with a width of B 164,2 = 0.1(1) G.

APPENDIX C: EXTENDED GROSS-PITAEVSKII SIMULATIONS
We compare our experimental results to theory [27,28], by performing numerical simulations of the extended Gross-Pitaevskii equation using a simple interaction potential and taking the quantum fluctuations and three-body losses into account within a local density approximation. In this equation g = 4πh 2 a s /m is the contact interaction parameter and is the dipole-dipole interaction of polarized particles, with θ the angle between the polarization direction and the relative orientation of the dipoles, and μ = 9.93μ B the magnetic dipole moment of 162 Dy. We change the scattering length a s in the range of 60a 0 -115a 0 . Furthermore, we use the measured L 3 parameter for the respective scattering lengths. The EGPE (C1) uses two assumptions: the Born approximation for the interaction potential and the local density approximation. The second approximation is supported by QMC simulations [34] and a comparison of theory and experiment with erbium [4]. The first assumption was studied in [30] and needs to be adjusted at finite temperature. The departure from the Born approximation can be taken into account by an effective dipolar length a dd that is shifted by a few percent compared to the dipolar length within the Born approximation.
In order to get the theory curve in Fig. 1 we used two different methods leading to the same result. For both methods we choose V ext = 0 and start with an atom number N > N crit , initially prepared with an elongated Gaussian density distribution. Then we find the ground state by imaginary-time evolution of the EGPE. Next we can either simulate atom losses like in the experiment or repeat this process of finding the ground state with ever lower atom number to start with until we do not find a stable solution anymore. In the second method we get an uncertainty due to the step size that we choose for the atom number. For the first method we do realtime evolution of the EGPE in order to simulate the dynamics of three-body losses. Due to the losses, the density and the effective two-body attraction reduces with time, until we reach N = N crit , where the contributions of the effective two-body attraction and the quantum pressure cancel each other. This leads to the evaporation of the droplet into the gaseous phase, which leads to a steep decrease of the density by more than an order of magnitude. This suppresses further losses and the atom number stays almost constant as soon as the droplet has evaporated, validating the interpretation of the experimentally observed loss curves.

APPENDIX D: DROPLET PHASE DIAGRAM
We can calculate the phase diagram [28,47] of trapped dipolar atoms as a function of the trap aspect ratio λ and the scattering length a s using the EGPE and then either apply a Gaussian ansatz to solve it analytically or resort to full numerical simulations. Here we are only interested in a qualitative discussion and therefore restrict ourself to the Gaussian ansatz; the full simulations together with a measurement of the critical point λ c can be found in [51]. The calculated phase diagram for a cylindrically trapped BEC, with mean trap frequencyω = √ ω x ω y ω z = 30 Hz and containing 20 000 162 Dy atoms, is shown in Fig. 4 of the EGPE exists (white region in Fig. 4), which has a cloud aspect ratio close to that of the trap only weakly altered by magnetostriction [52]. At low scattering lengths, also only a single solution exists, with a cloud aspect ratio much less than 1 that is more or less independent of λ. This is the quantum droplet solution (blue region) that is only stable due to beyond mean-field corrections [3,4,28,47]. This solution has a peak density that is a factor of 10 higher than that of the BEC. For trap aspect ratios larger than λ c there exists a bistable region (gray region) where both solutions can be stable. Crossing the boundary leads to a modulational instability and therefore higher number of droplets than what is expected for the ground state of the system. In the case of λ < λ c the two phases are connected through a continuous crossover instead of the phase transition in the bistable region. In Fig. 4 we also show the way in which we prepare the single droplet ground state in the experiment (indicated by the red circles and the black arrows). We start with the atoms in our CODT with a trap aspect ratio λ = 3.8 and a scattering length a s ≈ 140a 0 , deep in the BEC regime. We then change the trap aspect ratio to λ ≈ 1 and then change the magnetic field amplitude to probe the crossover to the droplet phase. In the droplet regime we can then turn off the trap completely and observe the self-bound state in free space.

APPENDIX E: EXPANSION VELOCITY
While a BEC is a gaseous state, which means that it freely expands in the absence of an external trap, the droplet state is self-bound due to its intrinsic interactions [5,27,28] and therefore does not expand. To map out the range of self-bound droplets, we determine the expansion velocity v exp [4] by FIG. 5. Expansion velocity across the crossover from a BEC to a quantum droplet. We extract the expansion velocity v exp from the evolution of the widths of the atomic cloud in time of flight for up to 20 ms, averaged over five realizations. This procedure can be applied in both the x and y directions, for which we find comparable results. We thus plot here the average over both directions. As error bars we show the quadratic sum of the uncertainty of the determination of v exp along the two directions and for the scattering length the uncertainty due to the experimental field stability, the knowledge of the Feshbach resonances, and the background scattering length. The two insets show example single-shot images for a self-bound droplet (top) and an expanding BEC (bottom) after a 20-ms time of flight.
fitting the evolution of the observed widths of the atomic cloud σ TOF up to t TOF = 20 ms to σ TOF = √ σ 2 0 + v 2 exp t 2 TOF . In this σ 0 corresponds to the size at zero time of flight. The extracted expansion velocity across the complete crossover from stable BEC to single droplet state is shown in Fig. 5. Entering the regime with ε dd > 1, we observe a small decrease of the expansion velocity with higher relative dipolar strength, until at around a s ≈ 110a 0 , where we observe a sharp decrease. For a s 107a 0 we enter the self-bound regime, where we do not observe an increase in the size within our 1-μm imaging resolution. In this regime we also observe aberrations in the images (see the top inset in Fig. 5) due to the small radial size compared to our imaging resolution. After some time, depending on the magnetic field and on the initial atom number, we find that the cloud has expanded even in the droplet state. This is understood in terms of three-body decay until it reaches the critical number below which the droplet is not self-bound anymore. The short lifetime due to enhanced threebody losses closer to the resonances leads to the observed increase of the expansion velocity for lower scattering lengths.

APPENDIX F: EXPERIMENTAL DETERMINATION OF THE CRITICAL ATOM NUMBER
In order to measure the critical atom number of a selfbound droplet we look at atom number decay curves [5], as exemplarily shown in Fig. 6 for a scattering length of a s = 99(7)a 0 . We use a sequence of intentional evaporation and subsequent expansion that allows us to determine the atom number precisely without being limited by the finite resolution of our imaging system or by the high density of the droplets. This intentional evaporation is done by ramping up the magnetic field within 100 μs to B evap ≈ 6.0 G. We 033088-7 FIG. 6. Exemplary atom number decay curve of a self-bound droplet. The measured decay of the atom number is plotted as a function of the levitation time, averaged over ten realizations, and the error bars denote the respective standard deviations. After a fast decay for short times (gray circles), we observe a constant atom number (blue circles). To extract the critical atom number N crit (horizontal blue line) we analyze the atom number distribution (histogram on the right-hand side) and fit our convolution model to the blue data (black line).
know that the ground state of the system is a BEC for this magnetic field. Therefore, by quickly ramping up the field, we force the droplet back into a gaseous state, where it then expands freely. We let the atomic cloud expand for either 8 or 30 ms, depending on the atom number, and then we image the atomic cloud with resonant absorption imaging. We use this sequence for different magnetic fields in the range from B = 5.293 G to B = 5.249 G, corresponding to a s = 107a 0 and 78a 0 , respectively.
With this sequence we observe that the atom number decays very fast in the beginning, but then settles to a constant value. This behavior results from an initial fast three-body decay of the high-density droplet state, followed by a quick expansion of the droplet as it crosses the phase boundary to the gaseous state. The crossing of the phase boundary leads to a fast drop in density and thus suppresses further loss. To extract the critical number we employ a statistical evaluation procedure, because the critical atom number is reached at different times during these ten realizations due to the stochastic preparation of the initial conditions in our experiment. We show histograms of the atom number distribution on the right side of Fig. 6. The first histogram includes all data points (gray bars), while the second histogram only includes the points for which the atom number has settled to its final value (blue bars).
We use a phenomenological model (black line in the histogram of Fig. 6) in order to extract the critical atom number of a self-bound droplet [5]. This model consists of a convolution of a Gaussian and a Maxwell-Boltzmann distribution. We use a symmetric Gaussian distribution in order to represent broadening effects that result from statistical errors, e.g., noise in the imaging of the atomic cloud. To cover the possibility that residual collective excitations in the droplet lead to an early evaporation at atom numbers higher than the critical number, we use an asymmetric Maxwell-Boltzmann distribution. By fitting this phenomenological model we extract the critical atom number and two different widths, one for each distribution. As an error bar of the critical atom numbers in Fig. 1 we use the quadratic mean of these two widths, together with an overall 10% uncertainty due to the uncertainty in calibration of the imaging system.

APPENDIX G: EXPERIMENTAL DETERMINATION OF THE BACKGROUND SCATTERING LENGTH
Using the latest measurement of the background scattering length of 164 Dy [6] as a starting point, we use the sensitive scaling of the critical atom number to extract the background scattering length of 162 Dy. This procedure, which is described in detail below, assumes that the critical atom number does not depend on other parameters, e.g., the actual density distribution of the droplet.
Since we observe a systematic shift of our measured critical numbers compared to the results obtained from numerical simulations, we first shift the theoretical curve in order to minimize the difference (res 164 ) between the experimental data for 164 Dy [5] and the shifted theory. This shift of the simulated boundary between self-bound droplet and expanding BEC can be done along the atom number axis or along the scattering length axis, resulting in slightly different results. In Fig. 7(a) we show the obtained residual for a shift along the atom number axis. The curve presents a clear minimum for a shift of 30% of the corresponding simulated critical atom numbers to  lower values (marked by the vertical red line). As uncertainty we use the range where this residual from the shifted theory has doubled (shown as light red).
From this point on, we then optimize the background scattering length a bg,162 such that we get the smallest residual of the measured critical atom numbers for 162 Dy with respect to the shifted theory curve. This can be done since with our knowledge of the parameters of the used Feshbach resonances the background scattering length is the only free parameter not known precisely. In Fig. 7(b) we show the residual res 162 for the case of a shifted theory along the atom number axis. From this we get two uncertainties, one again extracted from a doubling of the corresponding residual (light red area) and the second from the propagation of the uncertainty from the shifted theory curve (darker red area). With this procedure with the shift along the atom number axis we get a value of a bg,162 = 139(135, 141)a 0 for the background scattering length of 162 Dy.
Similarly, this can be done for a shift along the scattering length axis resulting in a bg,162 = 141(136, 145)a 0 . Taking the average of these two procedures, we end up with a bg,162 = 140(4)a 0 , in good agreement with the quoted literature value [24][25][26]. Note, however, that the systematic uncertainties for the measurement of the background scattering length of 164 Dy [6] apply here as well, leading to a final value with the respective uncertainty of a bg,162 = 140(7)a 0 .
Free from the systematic shift of the measurement, we can extract the ratio of the respective background scattering lengths of the two isotopes 162 Dy and 164 Dy. With this we find a ratio of a bg,162 /a bg,164 = 2.03 (6).

APPENDIX H: EFFECTIVE RENORMALIZATION OF THE DIPOLAR LENGTH
One possible explanation of lower critical atom numbers of a self-bound droplet arises due to the complexity of the scattering problem for dipolar lanthanide atoms such as dysprosium. As it was pointed out in [30], taking into account the full scattering amplitude leads to alterations compared to the Born approximation. This can be accounted for within the EGPE framework by an effective shift in the dipolar length a dd that depends on the collision energy. This effective shift is on the order of 2% for a temperature of 10 nK and 10% for 100 nK. As our BEC has an initial temperature of about 15 nK and the preparation process likely leads to additional heating, it may seem reasonable to include this effect in our considerations. Since the droplets are self-bound, we cannot use standard time-of-flight expansion to measure the temperature of these states. Future methods may measure the temperature inside the droplets, e.g. with embedded impurities [32] similar to liquid helium droplets [53]. With this, we can also clarify the role played by thermal fluctuations at larger atom numbers [54].
In Fig. 1 we also show two theory curves for an enhanced dipolar length: a dd + 5% enhancement (dashed red line) corresponding to a temperature of 30-50 nK [31] and a dd + 10% enhancement (dash-dotted red line) corresponding to 100 nK. The observed theoretical shifts are in good agreement with the experimental shifts, suggesting that such an effective enhancement of the dipolar length might play an important role in dipolar scattering at finite temperature.

APPENDIX I: PATH-INTEGRAL GROUND STATE
Given a Hamiltonian, the PIGS method can be used to evaluate exactly many-body properties of a correlated Bose system, beyond the mean-field plus Lee-Huang-Yang approximation. Designed as a reduction of the Feynman path-integral formalism to zero temperature, particle coordinates at different (but close) imaginary times are sampled in chains, starting from a variational wave function T (R) that is located at the end points. Since propagation in imaginary time removes any component that is orthogonal to the true ground-state wave function φ 0 (R) in the asymptotic limit, samples of φ 0 (R) are realized at the center of the chains when the total propagation time is long enough In this expression, G(R i+1 , R i ) is the imaginary-time propagator between positions R i+1 and R i in a time step δτ , which is related to the actionŜ through the expressionĜ = e −Ŝ . In general, the actionŜ is not known, but since δτ is small, a low-order series expansion in powers of δτ can be successfully employed. In this work we use one of the fourth-order propagators of Ref. [55], which improves convergence when compared with other, simpler schemes based on a secondorder Trotter expansion. The main ingredient required to perform a PIGS simulation is the Hamiltonian. Here we use a model that includes both the dipolar interaction and an effective potential V HC with a repulsive core that prevents the system from collapsing. Assuming that all the dipoles are polarized along the Z axis, the Hamiltonian readŝ where r i, j and θ i, j are the relative polar coordinates between the atoms, m is the atomic mass, and C dd sets the strength of the dipolar interaction. In order to study whether there are universality properties in the system at the given conditions, we solve the Hamiltonian for three different V HC models The coefficient C 6 in these equations is known for dysprosium [43]. The other coefficients, σ 9 and σ 12 , are fixed such that the complete interaction (V HC plus dipolar) has the experimental s-wave scattering lengths. This is accomplished by solving the low-momentum limit of the scattering T matrix, as briefly described in Appendix J.
033088-9 One of the fundamental quantities that can be obtained from the PIGS simulations is the ground-state energy, which is negative for a self-bound droplet state. As in the experiments and for a given scattering length, we find that there is a critical number N c below which the system ceases to be self-bound. Figure 8 shows, for a s = 60a 0 , the ground-state energy obtained from the Hamiltonian in Eq. (I2) for the three V HC models, as a function of the total number of particles, together with a linear fit that determines the point N c where the energy is zero. As it can be seen, different models lead to slightly different predictions, and that defines, together with the statistical noise in the simulation, the error bar in the critical number N c . As the scattering length increases, higher values of both N c and its error are found, but they are all compatible with the experimental error bars. Unfortunately, the computational cost of the simulation grows very fast with the number of particles and we cannot reliably determine N c for scattering lengths larger than a s = 90a 0 , approximately.
The droplets predicted by PIGS differ from those obtained in the EGPE approximation not only in the critical number, but also in the overall density profiles. Figure 9 shows the integrated density profiles along the axial directions of the droplet, obtained from both methods, for 1000 [ Fig. 9(a)] and 2000 [ Fig. 9(b)] atoms and for a scattering length of 60a 0 . As it can be seen, for these (low) particle numbers, the profiles are different, the PIGS one being more spread but with a lower central density. Still, the difference reduces when the number of atoms is increased from 1000 to 2000. Increasing the atom number even more, we expect the differences in the density profile to decrease.

APPENDIX J: CALCULATION OF THE s-WAVE SCATTERING LENGTH FOR A TWO-BODY POTENTIAL
The s-wave scattering length of the combined two-body plus dipole-dipole interaction is obtained from the on-shell T matrix, in the limit of vanishing momentum transfer. The T matrix can be obtained solving the Lippmann-Schwinger with V l,m l ,m the matrix elements of the complete interaction and M the reduced mass of two atoms. Due to the anisotropy of the dipolar potential, the matrix elements of T , for different values of the quantum number l and l , are coupled. Moreover, the long-range character of the combined potential makes all partial waves contribute significantly, even at low scattering energies [56]. Due to the nature of the dipolar interaction, different scattering lengths corresponding to different (coupled) channels appear and read a l,m l ,m ≡ lim k→0 π T l,m l ,m (k, k) k , with l = |l ± 2|. Still, the dominant one is the s-wave scattering length, corresponding to l = l = m = 0. In practice, the low-momentum matrix elements T l,m l ,m (k, k) can be efficiently determined using the Johnson algorithm [57], which solves the Schrödinger equation and finds the logarithmic derivative of the wave function. Table I shows