Quantum-fluctuation-driven crossover from a dilute Bose-Einstein condensate to a macro-droplet in a dipolar quantum fluid

In a joint experimental and theoretical effort, we report on the formation of a macro-droplet state in an ultracold bosonic gas of erbium atoms with strong dipolar interactions. By precise tuning of the s-wave scattering length below the so-called dipolar length, we observe a smooth crossover of the ground state from a dilute Bose-Einstein condensate (BEC) to a dense macro-droplet state of more than $10^4$ atoms. Based on the study of collective excitations and loss features, we quantitative prove that quantum fluctuations stabilize the ultracold gas far beyond the instability threshold imposed by mean-field interactions. Finally, we perform expansion measurements, showing the evolution of the normal BEC towards a three-dimensional self-bound state and show that the interplay between quantum stabilization and three-body losses gives rise to a minimal expansion velocity at a finite scattering length.


I. INTRODUCTION
The extraordinary success of ultracold quantum gases largely relies on their simple description. Notably the actual inter-particle interactions, which might be complex or even unknown at short range, are very well captured in terms of simple mean-field (MF) potentials, proportional to the local particle density, n, and accounting for the average mutual effect of all neighboring particles [1]. In particular, the MF contact interaction is solely parametrized by the s-wave scattering length, a s , which in turn can be widely tuned by means of Feshbach resonances (FRs) [2]. The MF treatment of a Bose gas leads to the famous Gross-Pitaevskii equation (GPE), which is very powerful in describing the underlying ground-state physics of a Bose-Einstein condensate (BEC), whereas the Bogoliubov-de Gennes (BdG) spectrum of collective modes thoughtfully accounts for excitations in the BEC [1].
Beyond the great achievements of dilute gases as testbed for MF theories, the quest for beyond-MF effects has triggered great interest in the ultracold community. The general question of how the many-body ground-state of bosons is modified by quantum fluctuations (QFs) of elementary excitations was first addressed by Lee, Huang, and Yang (LHY) in the 1950's [3]. The so-called LHY term, which accounts for the first-order correction to the condensate energy, scales for a contact-interacting gas as a s n na 3 s . While in the weakly-interacting regime the effect of QFs is negligible and difficult to isolate from MF contributions, it can be sufficiently amplified by increasing a s via a FR. Based on this concept, recent beautiful experiments with alkali have observed clear shifts of * Francesca.Ferlaino@uibk.ac.at the BdG spectrum and equation of state caused by the LHY term in strongly-interacting Fermi [4][5][6] and Bose gases [7,8].
While in these measurements the role of QFs remains ancillary, it has been recently pointed out [9] that in systems with various tunable interactions of different origin, the overall MF energy can be made small and the LHY term becomes dominant. In this extreme regime, a novel phase of matter is expected to appear, namely a liquid-like droplet state. For purely contact-interacting gases, this situation is hardly realizable since it would require, for instance, Bose-Bose mixtures with coincidental overlapping FRs [9]. In contrast, dipole-dipole interaction (DDI) offers genuinely this possibility in a singlecomponent atomic gas by competing with the isotropic MF contact interaction [10,11]. In the MF picture, a paradigm of the competition between DDI and contact interaction is embodied by the ability of quenching a dipolar BEC to collapse by varying ε dd = a dd /a s , where a dd = µ 0 µ 2 m/12π 2 is a characteristic length set by the DDI, with m the mass and µ the magnetic moment of the atoms [12,13]. Here stands for the reduced Planck constant and µ 0 for the vacuum permeability. In general, because of the special geometrical tunability of DDI with the external trapping potential and dipole orientation, the stability and phase diagram remarkably depend on λ = ν / / /ν ⊥ , where ν / / (ν ⊥ ) is the trapping frequency along (perpendicular to) the dipole orientation [11,12,14].
In parallel, recent breakthrough experiments with an oblate dysprosium (Dy) dipolar BEC (λ > 1) have shown that when quenching up ε dd , the system, instead of collapsing, forms a metastable state of several small droplets [15,16]. This observation has triggered an intense debate on the nature of such a state and its underlying stabilization mechanism [17][18][19][20][21][22][23]. These works have eventually converged in confirming a scenario based on arXiv:1607.06613v1 [cond-mat.quant-gas] 22 Jul 2016 QFs [16,20,21] and meanwhile raise fundamental questions on the parameter space of the effect, the singleto-multi droplet phase diagram, and the lifetime of the state. Aside from the Dy specificity, will other dipolar systems, lying in different parameter ranges -e.g. different dipolar strength (a dd ), geometry (λ) or atom number N -form a droplet state? When does QF bear stabilization and when does this picture break down? The fundamental questions of the actual ground-state and of the resulting phase diagram as a function of λ, ε dd and N , remain widely unanswered beyond the quench dynamics and metastablity.
The present work addresses these questions using a BEC of erbium (Er) atoms [24], confined in a cigarshaped harmonic trap with large elongation along the dipole orientation (λ 1). In this condition, the DDI is strongly attractive, imposing a stringent condition to the system stability: already at ε dd ≈ 1, a dipolar BEC of about 10 5 atoms should readily collapse within the MF description [10,12]. We explore the behavior of the dipolar gas by studying both the ground-state evolution and the dynamics when increasing ε dd . Our work, based on complementary observations on density distributions, elementary excitations of the BdG spectrum, expansion dynamics, and lifetime of the dipolar quantum gas, show the existence of a crossover from a dilute BEC to a dense single droplet solution and quantitatively identify the driven role of QFs in the system dynamics.

II. EXPERIMENTAL PROCEDURES
The atomic properties of Er offer a privileged platform to explore a variety of interaction scenarios. Beside its strongly magnetic character and its many FRs [26], Er has several stable isotopes. This feature adds an important flexibility in term of the choice of the background a s [27]. In our early work on Er BECs, we employed the 168 Er isotope, which has a background a s of about twice as large as the dipolar length, a dd = 65 a 0 [28,29].
In the work reported here, we produce and use a BEC of 166 Er in the lowest internal state. This isotope provides us with two major advantages. First, its background a s is comparable to its dipolar length, a dd = 65.5 a 0 , realizing ε dd = a dd /a s ≈ 1 without the need of Feshbach tuning. Second, 166 Er features a very convenient FR at ultralow magnetic-field values, B. To precisely map a s as a function of B, we use a spectroscopic technique based on the measurement of the energy gap of the Mott insulator state in a deep threedimensional optical lattice [29,30]. A detailed description is given in the Supplementary Material [25]. Between 0 and 3 G we observe a smooth variation of a s , which results from two low lying FRs whose centers are fitted to 0.05(5) G and 3.0(1) G respectively; see Fig. 1. This feature gives an easy access into the ε dd > 1 regime, allowing variation of ε dd from 0.70(2) to 1.58 (18)  The data points (circles) are extracted from spectroscopic measurements in a lattice-confined gas and the solid line is a fit to the data with its statistical uncertainty (grey shaded region [25]). (Upper inset) Zoom in of ε dd as a function B. The grey dashed line marks ε dd = 1; see also the other figures. The lower inset illustrates the geometry of our experimental setup, the relevant axes (x,y,z), the crossedoptical-dipole-trap beams (shaded regions), the magnetic field orientation (green arrow), and the//-and ⊥-imaging view axes (blue arrows).
data [2], we extract a s (B) valid for B in the [0.15, 2.5]−G range, which we will use along this manuscript [25]. We achieve Bose-Einstein condensation of 166 Er using an all-optical scheme very similar to Ref. [28] with cooling parameters optimized for 166 Er [25]. In brief, we drive forced evaporative cooling at a magnetic field B = 1.9 G, corresponding to a s = 81(2)a 0 (ε dd = 0.81(2)). In this phase B is oriented along the vertical z axis. At the end of the evaporation, we obtain a BEC of N = 1.2 × 10 5 atoms with a condensed fraction above 80%.
To reach the λ 1 regime, we slowly modify, in the last step of the evaporation, the confining potential to the final cigar shape, with typical frequencies (ν x , ν y , ν z ) = (156(1), 17.2(4), 198(2)) Hz. Simultaneously, we decrease B to 0.8 G (a s = 67(2) a 0 ) and then change the magneticfield orientation to the weak trapping axis (y) while keeping its amplitude constant [25]. Finally, we ramp B to the desired target value (and equivalently a s ) in t r [25], hold for a time t h and perform absorption imaging of our gas after a time-of-flight (ToF) of t ToF . Two imaging setups are used in order to measure either the density distribution integrated along the dipoles (//-imaging) or perpendicular to them (⊥-imaging) [25]. Figure 1 (lower inset) illustrates the final geometry of our system with ν / / = ν y , ν ⊥ = (ν 2 x + ν 2 z )/2 giving λ = 0.097 (3), and defines the relevant axes.
Here, we explore the properties of the system when the (color online) Density profiles in the BECto-droplet crossover : (a-c) 2D column density distributions probed with //-imaging and (d) corresponding central cuts along the x = 0 line (dots) after a ramp of tr = 10 ms and holding of t h = 6 ms to different as: 93 a0 ((a) and circles in (d)), 57 a0 ((b) and diamonds in (d)), 50 a0 ((c) and triangles in (d)). Each distribution is obtained by averaging 4 absorption images taken after tToF = 27 ms. In (d), the lines show the central cuts of the 2D bimodal fit results. The solid (resp. dashed-dotted) lines show the fit to a two-Gaussian (resp. MF TF plus Gaussian) distribution and the dashed (resp. dotted) lines show the broad thermal Gaussian of the corresponding fits. For visibility, not all MF TF plus Gaussian bimodal fits are shown.
DDI is made attractive enough to overcome the repulsive MF contact interaction (λ 1, ε dd > 1), after adiabatically changing (t r ≥ 45 ms) or quenching (t r = 10 ms) a s to its target value [25]. For t r ≥ 45 ms, the system evolves following its ground state and gives access to the slow dynamics, whereas for the t r = 10 ms case we can probe the fast dynamics and study the relaxation towards an equilibrium. The key question is whether QFs protect the system from collapsing. Indeed, in this regime, the MF treatment would imply that the attractive BEC becomes unstable, leading to a two-fold dramatic consequence [1]. First, some modes of the BdG spectrum acquire complex frequencies. Second, in a trap, the density distribution of the cloud undergoes a marked change on short time scales (≤ 1/ν ⊥ ), described as a "collapse", which can develop into a rapid loss of coherence [12,31], and pattern formations, such as anisotropic atom bursts ("Bosenova") and special d-wave-type structures, as observed in rubidium [32] and dipolar gases of chromium [12,13], respectively. This fast dynamics has been proved to be well encompassed by GPE simulation [13,14,33].

III. DENSITY DISTRIBUTION
In a first set of experiments, we study the evolution of the ToF density distribution of our dipolar Er BEC with a s . Figure 2 (a-c) shows the 2D column density profiles acquired with //-imaging for t r = 10 ms and t h = 6 ms ( 1/ν ⊥ ), and for different a s . To analyze our data, we fit the 2D density profiles to bimodal distributions. We illustrate their results on the central cuts (x = 0) in Fig. 2(d). For a s > a dd (a), the density distribution follows the one of a dilute BEC with the expected MF Thomas-Fermi (TF) profile on top of a broad Gaussian distribution, accounting for the thermal atoms. When lowering a s below a dd , we observe a drastic change of the density profile marked by the shaping of a narrower and denser central core (b-c), whose profile starts to deviate from the usual MF TF one. Indeed, the //-images are better described by the sum of two Gaussian functions already for a s < 70 a 0 as the fit to this distribution gives a smaller residue than the MF TF plus Gaussian bimodal fit; see Fig. 2 (d). In contrast with the evolution of the central peak, the distribution of the thermal atoms, encompassed by the broad Gaussian function, remains mainly unaffected by the change of a s . This remarkable behavior evidences an absence of a sizable heating and an apparent decoupling of the evolution of the coherent and thermal parts. The central peak remains visible down to a s ∼ 50 a 0 (ε dd ∼ 1.3) (c), and exhibits a long-lived character from several tens to hundreds of ms. As discussed in details later and suggested in Ref. [23], three-body (3B) collisions regulate the lifetime of the high-density component, as this loss mechanism exhibits a high-power dependence on the density. We note that we observe a similar qualitative behavior of the density distribution when using an adiabatic ramp of a s . However, in this case 3B losses already sets in during the ramp reducing the importance of the central peak.
Being a smoking-gun for long-range phase coherence, the survival of a bimodal profile in the ToF distribution suggests a persistent coherent behavior far beyond the MF instability threshold. The absence of a collapse advocates the outbreak of an additional stabilization mechanism.

IV. COLLECTIVE OSCILLATION
In a second set of experiments, we unveil the origin of stabilization mechanism by studying the elementary excitations of the coherent cloud. This is a very powerful probe of the fundamental properties in quantum degenerate gases [1,34]. In particular, collapse is intimately related to the softening of some collective modes at the MF-instability threshold. We focus here on the axial mode, which is the lowest-lying excitation in the system. It corresponds to a collective oscillation of the condensate length along y (R / / ) with frequency ν axial . The axial oscillation comes along with a smaller-amplitude oscilla- In the latter case, quenches to small as excite higher-energy modes and the system is not anymore in the linear response regime (open squares). The inset in (b) exemplifies a measurement of R / / (triangles) and its fit to a damped-sine (solid line) for as = 80 a0. We typically fit 4 to 5 oscillations for all our as.
tion of the radial sizes in phase-opposition; see Fig. 3(a). As a result, this mode has a mixed character between a compression and a surface mode [1]. The compression character is particularly relevant since it involves a change in the density and it is therefore sensitive to the LHY corrections [35].
We excite the axial mode either by ramping B during the final preparation stage or by transiently increasing the power of the vertical optical dipole trap (ODT) beam, after ramping B to B f . Here ν / / is abruptly changed from 17 Hz to typically 21 Hz, kept at this higher value for 8 ms and finally set back to 17 Hz. Following the excitation, we let the cloud evolve for a variable t h and image its ToF density distribution with ⊥-imaging. To extract ν axial , we probe the axial width R / / of the central coherent component of the gas [25] with t h and fit it to a damped-sine; see inset in Fig. 3(b). Figure 3 shows the observed ν axial normalized to the trapping frequency ν / / [36] as a function of a s for adiabatic (b) and non adiabatic (c) ramps. Both cases exhibit a similar qualitative behavior. For a s > a dd , the oscillations shows a smooth dependence on ε dd with ν axial increasing by about 5% around 1.70 ν / / . When lowering a s , the oscillation of the coherent part remains visible well below the ε dd = 1 threshold and ν axial exhibits a marked increase. ν axial /ν / / grows up to 2.6(1) at a s = 54 a 0 for t r = 100 ms (b). For t r = 10 ms (c), ν axial /ν / / first increases similarly to the adiabatic case (b), reaches a maximum of ∼ 2.13(7) at 57 a 0 (ε dd = 1.15), and finally decreases for even smaller a s (open squares). The latter behavior can be explained by the fact that the larger quenches in the interaction excites additional high-energy modes while it drives the system away from the linear response regime [37].

V. THEORY
To account for our observation and discern between the MF instability picture and QF mechanisms, we develop a beyond-MF treatment of our system at T = 0. The coherent gas is described here by means of the generalized non-local nonlinear-Schrödinger equation (gNLNLSE), which includes the first-order correction from QF, i. e. the LHY term, and 3B loss processes. The gNLNLSE reads as [20,23] i ∂ψ ∂t is the non-interacting Hamiltonian and V (r) = 2π 2 m η=x,y,z ν 2 η η 2 the harmonic confinement with η = x, y, z. The MF chemical potential, µ MF (n(r), dd ) = gn(r) + d 3 r V dd (r − r )n(r ), results from the competition between short-range interactions, controlled by the coupling constant g = 4π 2 a s /m, and the DDI term with V dd (r) = µ0µ 2 4πr 3 (1 − 3 cos 2 θ) and θ the angle sustained by r and the dipole moment µ. Here n(r) = |ψ(r)| 2 . The beyond-MF physics is encoded in the LHY term, leading to an additional repulsive term in the chemical potential , ∆µ(n, dd ) = 32 [38,39] using local-density approximation [40]. The last non-Hermitian term in Eq. (1) accounts for 3B loss processes [41]. In our calculations, we use the experimentally determined values of the 3B recombination rate of the condensate, L 3 (a s ) [25].
As discussed in Refs. [22,23], due to the repulsive LHY term, Eq. (1) sustains stable ground-state solutions for any a s and λ. For pancake traps (λ > 1), the solution of Eq. ((1)) is not unique. The phase diagram reveals three types of solutions: the one of a dilute BEC, a single droplet solution, and a third one, which separates the previous two phases, that corresponds to a metastable region of multi-droplet states. The latter has been observed in Dy experiments [15]. However, the single droplet solution appears difficult to access because of the overhead multi-droplet state and the stringent 3B loss mechanisms. Remarkably, in cigar-shaped traps (λ < 1) Eq. (1) has only one possible solution. In the ε dd parameter space, the corresponding wave function exhibits a smooth crossover from a dilute BEC to a single, high-density, macro droplet solution for increasing ε dd . It is worth to notice that the crossover physics, e. g. the formation and lifetime of the droplet state, is expected to crucially depend on the 3B collisional processes. In the following, we will concentrate on the λ < 1 case, which corresponds to our experimental setting.
The continuous and smooth change of the static properties of the system with increasing ε dd is consistent with our observations on the evolution from a dilute into an high-density state; see Fig. 2.
Based on Eq. (1), we theoretically investigate the dynamics of the coherent gas. In order to compare as close as possible to our experimental results, we precisely account for the experimental sequence by performing real time evolution (RTE) starting from the ground state of Eq. (1) at a s = 67 a 0 with N = 1.2 × 10 5 atoms. We simulate a linear ramp in a s from 67 a 0 to a variable final value of a s in t r , followed by a compression of the axial trap from ν / / = 17.3 Hz to 21 Hz during 8 ms. We record then the axial width from the standard deviation of n(r), σ y = y 2 , as a function of the subsequent holding time. The evolution of σ y is well fit by a sinusoidal function, whose frequency constitutes our theoretical prediction of ν axial .
In Fig. 3, we present our calculations with and without the LHY term. From a direct comparison with the experiments, one observes an excellent agreement of the theory including QFs, thus ruling out the MF scenario. In quench experiments (c), the effect of QFs is particularly evident since 3B losses do not have enough time to fully develop at the end of t r , leaving the system in the highdensity regime. In this condition, QFs play the crucial role of stabilizing the system, and drive the formation of a special coherent state, namely a single macro-droplet [20][21][22][23]. This provides a qualitatively different phase diagram than the one from the MF treatment, which predicts a collapse at a s ≈ 64 a 0 . For t r = 100 ms, a more stringent interplay between QFs and 3B losses arises, both mechanisms being able to drive the system out of the instability region. Although the system is expected to be stable even within the MF description down to a s ≈ 57a 0 , the effect of QFs is visible in a sizable shift of ν axial , which better matches the experimental data; see Fig. 3

VI. LIFETIME OF THE HIGH-DENSITY STATE
To further investigate the respective roles of 3B losses and QFs, we study the time evolution of the atom number of both the central core (N core ) and thermal (N th ) components along the BEC-to-droplet crossover. Since in the droplet regime the core density n core (r) dramatically increases, 3B losses are expected to play an important role even for moderate and low values of L 3 [23]. Notwithstanding, 3B losses and QFs exhibit different power dependencies on n(r) (see Eq. (1)) and thus the atom-loss dynamics should disclose the competition between diluteness and a droplet character. We fit a double exponential function to the data (solid lines) [25]. (d) From the fit, we deduct the mean in-situ density of the core,nc (see text), for t h = 4 ms (triangles) and 16 ms (squares) and as a function of as. The error bars include the statistical errors on the fit and on L3. The solid lines show results of the RTE including the LHY correction for the t h = 0 ms (red), t h = 25 ms (blue). Figure 4 (a-b) shows N core and N th , extracted from the double Gaussian fit as a function of a s after a nonadiabatic (a) and adiabatic (b) change of a s . Both cases show a similar evolution. When lowering a s , N core is first constant for a s > a dd , then shows a sharp drop starting around a s ∼ a dd and finally saturates for lower a s . We note that in the non-adiabatic case, N core decreases slower as compared to the adiabatic case because of the shorter timing involved. Remarkably, N th remains mainly unaltered over the whole range of a s and the whole system does not show any appreciable heating. This suggests that the condensed atoms, which are ejected from in the core, leave the trap instead of being transfer to the thermal component, confirming a picture in which the thermal and the condensed component have uncoupled dynamics.
We now compare the experiment with the theory, which, as previously, precisely accounts for the experimental sequence and its timing by performing RTE along Eq. (1). We compute here the final atom number N = n(r)d 3 r as a function of a s with and without the LHY term. Remarkably, the observed evolution of N core is very well reproduced by our beyond-MF calculations (solid lines), whereas in absence of the LHY stabilization, the calculations predict losses in the condensed core to occur at values of a s too large compared to the measured ones; see Fig. 4(a,b,d).
Finally, we investigate the in-trap time evolution of N core after quenching a s in the droplet regime; see Fig. 4 (c).
Our measurements reveal three different timescales for the losses. At very short t h (≈ 0 − 8 ms), N core is constant and the droplet state preserves its high density. It follows a fast decay dominated by 3B (≈ 8 − 25 ms), having the effect of ejecting atoms from the core until N core decreases slightly below 10 4 atoms. At this point, the loss dynamics substantially slows down (≈ 25 − 1000 ms).
Using the relation 1 Ncore dNcore dt h = −L 3n 2 , we extract the mean in-situ densityn = n(r) 2 of the high-density component in the BEC-to-droplet crossover. Figure 4 (d) shows our results for t h = 4, 16 ms. We observe a prominent increase ofn across the threshold a s ∼ a dd , and a surviving high density state deep into the MF instability regime. The formation of the droplet state is particularly visible for the t h = 4-ms case. Here,n grows from 6.2(9) × 10 20 m −3 at a s = 67 a 0 to a maximum of 35(7)×10 20 m −3 at a s = 57 a 0 , while it is slightly reduced to ∼ 24 × 10 20 m −3 at a s ∼ 46 a 0 . These observations are qualitatively reproduced by our simulations including the LHY correction and 3B losses.
Our results together with the good agreement between theory and experiments provide an alternative confirmation of the central role of beyond-mean-field physics. The lifetime of the high-density core reveals, on the one hand, the activation of the LHY term and the crossover toward a dense droplet state, and on the other hand the counteracting role of 3B losses in regulating the maximum density in the droplet regime.

VII. EXPANSION DYNAMICS
Besides their unlike stability diagram, collective excitations, and density distribution, a dilute BEC in the MF regime and a quantum droplet are also expected to exhibit a markedly different expansion dynamics. While the first is confined by an external trapping potential and thus freely expands in its absence, a droplet state is selfbound by its underlying interaction in analogy with the He-droplet case [20][21][22][23]. As in our previous discussions, the evolution from a trap-bound to a self-bound solution is expected to be regulated by the interplay between QFs and 3B loss processes.
We investigate the expansion dynamics of our system for various a s . To preserve the high density of the coherent component, our measurements focus on short timescales with t r = 10 ms and t h = 5 ms. After preparing the system at the desired a s , we abruptly switch off the ODT, let the gas expand for a variable t ToF , and probe the cloud width using the //-imaging. We fit the observed density distribution to a double-Gaussian function, as previously described. To extract the width σ η of the high-density core (n core ), we compute the second moments σ 2 η = η 2 n hG (r)dr along η = x, z. Figure 5 (a) exemplifies the ToF evolution of σ η=x at a s = 93a 0 , 64a 0 , and 55.5a 0 . When entering the ε dd > 1 regime, atoms in the high-density core exhibits a marked slowing down of the expansion dynamics, which can not be explained within the MF approach.
To systematically study this effect, we repeat the above measurements for different values of a s (i. e. ε dd ). From σ x (t ToF ), we extract the value of the expansion velocity, v x , by fitting the data to σ x (t ToF ) = σ 2 x,0 + v 2 x t 2 ToF . Figure 5 (b) shows v x in an ε dd range from 0.7 to 1.5. When the system approach the droplet regime with decreasing scattering length (a s < a dd ), v η undergoes a strong reduction and drops to a minimum equal to v x = 0.40 (2) µm/ms at about 56 a 0 (ε dd ∼ 1.17). For fur-ther lowering of a s , v x starts to increase again. A similar behavior is observed for v y . We note that only the highdensity component reveals this intriguing dependency on a s , whereas the thermal part shows a mainly constant expansion velocity.
Our observations agree well with the theory including the LHY term; see solid line in Fig. 5(b). The ToF evolution is calculated using a multi-grid numerical scheme [25]. We record the evolution of σ η with t ToF and extract the corresponding expansion velocities from the asymptotic behavior of dσ η /dt ToF . Our simulations show a slowing down with a minimum of v x = 0.32 µm/ms at a s ∼ 56 a 0 (ε dd ∼ 1.17), followed by an increase at lower a s [42]. In contrast, calculations in absence of beyond-MF corrections fail to reproduce the experimental data. Here, the velocity is first more drastically reduced above the MF instability threshold ε dd ∼ 1 than it is expected with LHY corrections, it then already increases at this threshold. The first point relies on the trivial slowing down of a BEC whose mean repulsion energy is weakened (by reducing a s or decreasing its population N core ). The second point reveals a collapsing behavior that gives rise to an explosive evolution of the density profile. The minimal velocity is here found to be v x = 0.56 µm/ms at a s = 68 a 0 , which is a much higher value than both our experimental results and our theory predictions including the LHY correction.
The expansion behavior can be qualitatively well understood considering the so-called released, or internal energy, E R . This is the energy of the system when substracting the energy related with the confinement [1]. In the MF scenario, E R > 0 as long as the ground state is stable. The BEC expands ballistically and v η decreases for decreasing a s and N . In the unstable regime, the expansion velocity depends crucially on the value of t h at which the trap is switched off due to the occurrence of an in-situ collapse dynamics. Contrary, in the presence of QF, a stable ground state always exists. The sharp variability in t h is expected to be suppressed. Assuming a fixed N core (i.e. no 3B losses), one can show that E R decreases with decreasing a s and eventually reaches E R < 0 for a s < a SB , marking the onset of the self-bound (SB) solution (e.g. a SB = 56a 0 for N = 1.2 × 10 5 ) [25]. However, in stark contrast to the MF case, E R increases with decreasing N core in the droplet regime. We note that a SB is then shifted to lower values when N core gets reduced by 3B losses, thus affecting the lifetime of the self-bound solution.
The existence of a minimal expansion velocity is thus a direct consequence of the competition between the decrease of E R for decreasing a s at a fixed N core , and the increase of E R for decreasing N core in the droplet regime. In the crossover regime, the system smoothly evolves towards a fully self-bound state (v η = 0) until 3B losses, occurring in trap or in the initial phase of the expansion, set in to unbind the system and to reduce the lifetime of the droplet state.

VIII. CONCLUSION
In summary, we have demonstrated the existence of the crossover from a dilute BEC to a quantum droplet state driven by QFs. Our experiments not only demonstrate that LHY stabilization is a general feature of strongly dipolar gases, but show as well an excellent quantitative agreement with our parameter-free theory, which is based on a generalized GP equation with LHY correction. The agreement concerns the lowest-lying excitation mode, the evolution of the atom losses, and the ToF expansion. Concerning the latter, we observe a prominent reduction of the expansion velocity to a minimum value, which provides a fingerprint of the competition between 3B losses and LHY stabilization. ACKNOWLEDGEMENT We thank M. Baranov We prepare an ultracold gas of the 166 Er isotope following a similar trapping and cooling scheme as the one employed for 168 Er [28,43]. We load a crossed-ODT from a narrow-line MOT of 3×10 7 166 Er atoms at about 10µK. At the end of the MOT sequence, the atoms are automatically spin-polarized in their lowest Zeeman sub-level [43]. The dipole orientation follows the one of the external applied magnetic field, B. In our experiment, the latter is controlled by independent tuning of the components B x , B y , B z along the experimental coordinate system (x, y, z), as defined in Fig. 1 (lower inset).
The ODT results from the crossing at their foci of two red-detuned laser beams at a wavelenght of 1064 nm. One beam propagates horizontally along the y-axis, and the other propagates vertically and nearly collinear to the zaxis. The z-beam has a maximum power of 10 W and an elliptical profile defined by its waists of (w  z by time averaging the frequency of the first-order deflection of an Acousto-Optic Modulator (AOM). This scanning scheme enable both an efficient loading of the MOT into the ODT (> 30% of the atoms are loaded) and an adiabatic and controlled tuning of the trap aspect ratio λ over a broad range. We achieve Bose-Einstein condensation of 166 Er atoms by means of evaporative cooling in the crossed ODT at |B| = B z = 1.9 G (a s = 80(2) a 0 ). Typically, we first rapidly (in 600 ms) reduce the power and aspect ratio w (y) x /w (y) z of the ybeam from 24 W to 4 W and 10 to 1.6, respectively. We further decrease the power of the y-beam from 4 W to 0.3 W in 3 s in an exponential manner and then exponentially increase the aspect ratio w (y) x /w (y) z from 1.6 to 8 in 2.5 s. The final trap frequencies are typically of (ν x ; ν y ; ν z ) ∼ (40; 40; 180) Hz. We finally obtain BECs of typically N = 1.2 × 10 5 atoms with more than 80% condensed fraction and a temperature T ∼ 70 nK. We typically measure N and the condensed fraction from a bimodal fit of the 2D column density distribution measured along //-imaging with t ToF = 27 ms. T is extracted from the evolution of the thermal size of the bimodal fit with t ToF varying from 14 to 28 ms.

Experimental setup and axes
In our setup we define the orthonormal (x, y, z)coordinate system in the following way: the vertical axis z is aligned with gravity and the y axis with the horizontal ODT beam; see Fig. 1. The polarizing magnetic field is created by three orthogonal pairs of coils. These pairs of coils define an orthonormal (X, Y, Z)-coordinate system with Z = z and (X, Y ) rotated by a small angle θ as compared to (x, y). The magnetic field components B X , B Y , B Z , each created by each pair of coils, can be controlled independently. We estimate θ to be about 15 o by sensing directly the atomic cloud, as its dipolar character makes it sensitive both to the trap geometry and to the magnetic field direction.
The small tilt θ between the dipoles and y causes a small reduction of the mean DDI energy and corresponding small shifts of the expected characteristics compared to the one predicted for θ = 0: the MF collapse threshold should appear at a lower a s and, for a given a s , ν axial /ν / / and v η are shifted respectively down and up. We have experimentally evaluated the shift of ν axial deep in the stable BEC regime (|B| = B 2 X + B 2 Y = 2 G) to be of the order of 2% and in the droplet regime to be about 10 to 15%, as confirmed also by our theoretical predictions.
Finally, we also note a tilt between the//-imaging beam and our reference frame, corresponding to an angle of θ im / /,0 ∼ 28 o compared to y and θ im / / = 13 o compared to Y in the xy-plane. The ⊥-imaging axis is tilted by ∼ 15 o , mainly in the xz-plane. Such tilts shift the observed size compared to the ideal case of imaging along and perpendicular to the dipoles. Such an effect is not expected to impact the measurement of the collective frequencies, whereas it might bring a systematic shift of v x because of a mixed projection of v X and v Y , the two first velocity components in the (X, Y, Z)-coordinate system, which are respectively perpendicular and along the dipoles.
In the theoretical calculations presented in the main text (Eq. (1), RTE and Gaussian ansatz), for simplicity, we do not account for these angles, whose effects are estimated to be smaller than our systematics.
Precise determination of the a s -to-B conversion is a delicate issue, especially in the case of complex species, such as Er, for which comprehensive multi-channel calculations are still out of reach and the knowledge of a s should thus rely on experimental investigations. We perform lattice modulation spectroscopy in a threedimensional optical lattice. From the measurements of the energy gap in the Mott insulator state we extract a s (B). Our lattice experiments focuses in the region of low B-field [0, 2.5 G] and are based on a lattice setup and procedure similar to the one described in Ref. [29]. In brief, after producing the BEC, we load the atoms in a three-dimensional optical lattice by exponentially increasing the lattice-beam power in 150 ms. The typical final depths are (s x , s y , s z ) = (20, 20, 100), given in unit of the respective recoil energies h × (4.2, 4.2, 1.05) kHz. At these lattice dephts, the gas is in the Mott insulator state. We then vary B = (0, 0, B z ) to the desired value by rapidly changing B z in 2 ms, either just before or just after loading the lattice. We use the latter option for the smallest B values at which L 3 is enhanced because of its proximity to the near-zero-field resonance. In this case, we further hold 20 ms to make sure the magnetic field is fully established before performing the modulation.
To perform spectroscopy measurements, we sinusoidally modulate s y at a variable frequency ν m for 90 ms with a peak-to-peak amplitude of about 30%. Finally, we ramp down the lattice depths to zero in 150 ms, and measure the recovered condensed fraction as a function of ν m from //-imaging ToF picture. For the smallest B values considered, we also ramp B back to 2 G in 2 ms before switching off the lattice-beams in order to again minimize 3B loss effects.
When varying ν m , we observe a resonant depletion of the condensate due to particle-hole excitations. The resonance condition in the Mott-insulator regime is given by Here U s , U dd and V dd are respectively the on-site contact interaction, the on-site dipolar interaction and the nearest-neighbor dipolar interaction along y in the corresponding extended Bose-Hubbard model. U dd and V dd can be accurately predicted from the knowledge of the lattice depths and dipole orientation and in our typical experimental condition, they are equal to h × −344.8 Hz and h × 31.5 Hz respectively. By subtracting the theoretical dipolar contributions to the measured frequency, we extract U s , which is directly proportional to the scattering length a s . A precise mapping of a s in the ultralow B-field region is then obtained by repeating the above measurements at various B values; see Fig. 1.
In the region [0, 2.5] G, our lattice spectroscopy reveals the presence of two FRs, one located at about zero B field and the other one at about 3 G. The existence and position of these two FRs agree with our Feshbach spectroscopy measurements performed in an harmonically trapped thermal cloud, where the maxima in 3B losses approximately pinpoint the resonance positions. In this measurement, further FRs are identified at 4.1 G and 5 G.
The scattering length of 166 Er can be parametrized by the following simple expression [2] a s (B) = a bg (B) 1 + in which the specific positions (B i ) and widths (∆B i ) of the two first resonances as well as the background scattering length are obtained from a fit to our lattice spectroscopy measurement. From the fit, we obtain B 1 = 48(45) mG, ∆B 1 = 39(20) mG, B 2 = 3.0(1) G, ∆B 2 = 110(35) mG. The B-dependent background scattering length a bg (B) accounts for the overlapping resonances and reads as a bg (B) = 62(4) + kB with k = 5.8(1.2) a 0 /G. We also account for the small effect of the two next resonances, whose positions B 3 , B 4 and widths ∆B 3 , ∆B 4 are fixed to their estimates from the loss-spectroscopy measurements to 10 mG. We check that the precise values of this parameters has little effect on our empirical description along Eq. A.2.

Ramps in scattering length
Our measurements rely on controlled variations of the scattering length a s (B). In our experiments, we either adiabatically change a s using t r = 45 ms or we quench it using t r = 10 ms. The adiabatic condition for a s reads as As shown in Fig. 6 For our theoretical description, we use a linear change of a s (t), similar to case (ii).

Determination of the three-body recombination rate coefficient
Since three-body inelastic losses play a crucial role in the many-body dynamics and lifetime of the droplet state, we run a dedicated set of measurements to determine L 3 . We first prepare a non-degenerate thermal sample of Er atoms at T = 490(10) nK in an harmonic trap. We then record the decay of the atom number, N th , as a function of t h in a range from 0 to 1 s. N th is obtained from a Gaussian fit to the measured density distribution.
To fit the time evolution of N th , we use the integrated 3B rate equation, which reads as 1 Here, n 2 is the mean square density on the cloud. To describe the scaling of n 2 with N th , we use its prediction for an ideal gas at thermal equilibrium at T whose state occupancies follow Boltzmann law and take into account the anti-evaporation effect [44]. Then L th 3 is extracted from a fit of N th (t h ) along: where N 0 is the atom number at t h = 0 ms and γ 0 is with k B the Boltzmann constant andν = (ν x ν y ν z ) 1/3 . We account for the a s -dependence of L th 3 near a FR by repeating the measurement at different B. We check that the measured L th 3 does not depend on the B orientation and measure its B-dependency using |B| = B y and for B varying in 0.1 to 1.9 G. Note that, due to the bosonization effect, the L 3 in a quantum degenerate bosonic gas is equal to L th 3 /3. Figure 7 shows L 3 in a quantum degenerate bosonic gas of 166 Er as a function of a s using the measured a s -to-B conversion. Despite the existence of many coupled molecular potential in Er, we measure a comparable low L 3 (a s ) as compared to the typical values reported in alkali atoms. L 3 (a s ) varies between 1.7(3) × 10 −40 m 6 /s at a s = 18(17) and 3.2(3) × 10 −42 m 6 /s at a s = 80(2) a 0 .

Time-of-Flight measurements
For our ToF measurements, we abruptly extinguish the ODT in about 100 µs. In order to accurately image our gas while minimally influencing its dynamics during the expansion, we split the ToF in two parts. During a first part, lasting t ToF − t B , B remains unchanged and the dynamics occur at the original a s and dipole orientation. At time t B before the image is taken, B is modified both in amplitude and in direction, in order to correctly set the quantization axis for the imaging light to be σ − polarized. The amplitude |B| of the imaging field is chosen constant and equal to 0.31 G for all a s considered. t B is set to be 12 ms for //-imaging and 15 ms for ⊥-imaging where the change in B is more drastic for the typical dipole orientation (Y ) used in this experiment. We note that our resolution limit for both //-and ⊥imaging is estimated to be 3 µm. Moreover the effective pixel sizes are set to 8.4 µm and 3.1 µm in our setup. These limit the size of the structure we are able to observe as well as the minimal t ToF we can use, which is typically t ToF ≥ 16 ms.

Averaging experimental density distribution
In order to obtain a better image quality and resolution, we typically average four experimental absorption images taken in the same condition and with the same experimental series (i.e. within less than a few hours interval). In order to average the images we first define a region of interest (ROI) of typically 300 µm × 300 µm containing the cloud shadow and translate the ROI to superimpose the cloud centers. To estimate the translation amplitudes for each individual image, we use the center from a simple Gaussian fit to the 2D density distribution ROIs. In this averaging process, we use a sub-grid resolution of 1/10 of a pixel to more accurately superimpose the centers. We fit the averaged density distribution after binning back to the original resolution. We note that fits on the sub-pixel resolved images give similar results.

Extracting the frequency of the collective modes
As stated in the main text, we focus on the axial mode, which reveals itself by a collective oscillation of the axial size R / / of the BEC, along with smaller amplitude oscillations of the radial size in phase opposition. We extract its frequency ν axial by studying the larger amplitude oscillation of R / / . For this, we probe the ToF density distribution of the gas with ⊥-imaging after a ToF of 30 ms. We focus on the sizes of the central, high-density component of the cloud, which we study as a function of t h for different a s . We note that the precise shape of the column density profile is expected to change as a function of a s and this in a different way for the two axis x and y under observation in ⊥-imaging. This complicates the analysis, in particular compared to the //-imaging where both axis are nearly equivalent. Here, we extract the sizes of the central component, using various methods and select the most reliable method according to a s . Typically, we use a bimodal MF TF plus Gauss fit for a s ≥ 57 a 0 . For a s ≤ 57 a 0 we select a central region in the cloud and perform a simple Gaussian fit on it. Such a determination should give access to the variations, if not to its absolute value, of R / / with t h at fixed a s and thus enable to determine ν axial .
To fit ν axial , we use a damped-sine function of t h . Typically we fit the evolution of R / / for t h varying from 0 to few hundreds of ms, depending on the damping rate observed. The upper value of t h used is never less than 150 ms such that at least 4 to 5 oscillations are observed and fitted. We also note that for our shortest t r = 10 ms we typically discard the first 4 ms of the evolution in order to ensure that the magnetic field is safely stabilized at its target value. 9. Bimodal fits of the density distribution in //-imaging.
To quantitatively analyse the experimental column density distribution imaged along //-imaging, we perform bimodal fits on the 2D averaged profiles n / / (x, z). The bimodal fits are made of the sum of two peak distributions, describing respectively a high-density, coherent part and a thermal incoherent background. To account for the change of the profile of the density distribution across the BEC-to-droplet crossover, we use two types of fitting functions f fit (x, z): (i) A sum of a MF TF and a Gaussian distribution, which account respectively for the coherent and the thermal part. For the integrated column density, the MF TF distribution The sum of two Gaussian distributions. Typically the central Gaussian can be anisotropic with any orientation axis.
As expected, the thermal background is broader than the central coherent component and for the t ToF considered its width is typically at least 1.6 times larger. It is first fitted by taking out the central part of the density distribution at radius typically smaller than 1.3 times its width.

Describing the atom number decay
In Figure 4 (c) and in the main text we have briefly described the evolution of the atom number in the coherent part N core with t h . There our aim was to extract the mean density and we did not expand on the mere description of the N core (t h ).
We note that N core (t h ) is well accounted for by a double exponential decay evolution of respective amplitude N 0 (1 − p) and N 0 p. N 0 is the initial atom number. Each decay corresponds to a different time constant, respectively t fast and t slow , and starts after a different delay time, respectively t d and t D . We fix t D = t d + 2t slow . For t r = 10 ms t d is approximately constant and equal to 3.65 ms. The absence of evolution for t h < t d indicates that the magnetic field may not have reached the target value in the first ms. For t r = 45 ms, t d can be set to 0. t fast , t slow and p (and N 0 ) depend on a s , typically de-creasing with it. The evolution is illustrated in Fig. 8 for t r = 10 ms.
Due to its simplicity, Eqs. (A.11) and (A.12) permit a much more flexible simulation of the exact experimental conditions and sequences compared to the obviously more exact but numerically much more cumbersome simulation of the gNLNLSE. We have checked that the results of the Gaussian ansatz approach are in excellent agreement both qualitative and to a large extent also quantitative to full simulations of the gNLNLSE, in what concerns lowest-lying excitations, atom losses, and expansion velocities.

Self-bound droplets
The Gaussian ansatz approach allows for an intuitive understanding of the degree of self-bound (SB) character of the system. As mentioned in the main text, a SB solution is characterized by a negative released energy, E R < 0. In absence of losses, we may evaluate E R by means of the Gaussian Ansatz for the ground-state of a trapped BEC with scattering length a s and fixed N . Figure 9 shows the results for E R for different N values. Whereas E R increases with growing N for large a s , for small a s in the droplet regime E R increases with decreasing N . For each N there is a finite scattering length a SB such that if a s ≤ a SB the droplet will be fully self-bound (v η = 0). Given its N -dependence, a SB shifts to lower values with decreasing N . 3B losses add a time dependence on N and thus on a SB that governs the lifetime of the droplet state. We note that for small (t r , t h ) E R may change its sign during the ToF due to atom losses, i.e. a SB solution may unbind during the ToF. In our experiments, the interplay of losses and LHY stabilization leads to a minimal released energy, that translates into a minimal expansion velocity as shown in the main text.

Simulation of the ToF expansion using the gNLNLSE
ToF expansion is simulated using the gNLNLSE by means of a multi-grid method, i.e. dynamically enlarging the numerical grid following the expansion. This is necessary due to the obvious difference in length scales at the beginning and at the end of the ToF. We note that the precise description of the ToF dynamics is very relevant, since in contrast to standard cases, nonlinear dynamics and losses here may play an important role during the