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

## Abstract

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 -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 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.

###### pacs:

30.0, 67.85.-d, 03.75.Kk^{1}

## 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, , and accounting for the average mutual effect of all neighboring particles Pitaevskii and Stringari (2003). In particular, the MF contact interaction is solely parametrized by the -wave scattering length, , which in turn can be widely tuned by means of Feshbach resonances (FRs) Chin *et al.* (2010). 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 Pitaevskii and Stringari (2003).

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 Lee *et al.* (1957). The so-called LHY term, which accounts for the first-order correction to the condensate energy, scales for a contact-interacting gas as .
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 via a FR.
Based on this concept, recent beautiful experiments with alkali have observed clear shifts of the BdG spectrum and equation of state caused by the LHY term in strongly-interacting Fermi Altmeyer *et al.* (2007); Shin *et al.* (2008); Navon *et al.* (2010) and Bose gases Papp *et al.* (2008); Navon *et al.* (2011).

While in these measurements the role of QFs remains ancillary, it has been recently pointed out Petrov (2015) 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 Petrov (2015). In contrast, dipole-dipole interaction (DDI) offers genuinely this possibility in a single-component atomic gas by competing with the isotropic MF contact interaction Lahaye *et al.* (2009); Baranov *et al.* (2012).
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 , where is a characteristic length set by the DDI, with the mass and the magnetic moment of the atoms Koch *et al.* (2008); Lahaye *et al.* (2008). Here stands for the reduced Planck constant and 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 Baranov *et al.* (2012); Koch *et al.* (2008); Metz *et al.* (2009).

In parallel, recent breakthrough experiments with an oblate dysprosium (Dy) dipolar BEC () have shown that when quenching up , the system, instead of collapsing, forms a metastable state of several small droplets Kadau *et al.* (2016); Ferrier-Barbut *et al.* (2016).
This observation has triggered an intense debate on the nature of such a state and its underlying stabilization mechanism Xi and Saito (2016); Bisset and Blakie (2015); Blakie (2016); Wächtler and Santos (2016a); Bisset *et al.* (2016); Wächtler and Santos (2016b); Baillie *et al.* (2016). These works have eventually converged in confirming a scenario based on QFs Ferrier-Barbut *et al.* (2016); Wächtler and Santos (2016a); Bisset *et al.* (2016) and meanwhile raise fundamental questions on the parameter space of the effect, the single-to-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 (), geometry () or atom number - 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 , and , remain widely unanswered beyond the quench dynamics and metastablity.

The present work addresses these questions using a BEC of erbium (Er) atoms ^{2}*et al.* (2009); Koch *et al.* (2008).
We explore the behavior of the dipolar gas by studying both the ground-state evolution and the dynamics when increasing . 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 Frisch *et al.* (2014), Er has several stable isotopes. This feature adds an important flexibility in term of the choice of the background us: (). In our early work on Er BECs, we employed the isotope, which has a background of about twice as large as the dipolar length, Aikawa *et al.* (2012); Baier *et al.* (2016).

In the work reported here, we produce and use a BEC of in the lowest internal state. This isotope provides us with two major advantages. First, its background is comparable to its dipolar length, , realizing without the need of Feshbach tuning. Second, features a very convenient FR at ultralow magnetic-field values, .
To precisely map as a function of , we use a spectroscopic technique based on the measurement of the energy gap of the Mott insulator state in a deep three-dimensional optical lattice Mark *et al.* (2011); Baier *et al.* (2016). A detailed description is given in the Supplementary Material sup ().
Between and G we observe a smooth variation of , which results from two low lying FRs whose centers are fitted to G and G respectively; see Fig. 1. This feature gives an easy access into the regime, allowing variation of from to by changing from 2.5 G to 0.15 G; see Fig. 1 upper inset.
By fitting our data Chin *et al.* (2010), we extract valid for in the G range, which we will use along this manuscript sup ().

We achieve Bose-Einstein condensation of using an all-optical scheme very similar to Ref. Aikawa *et al.* (2012) with cooling parameters optimized for sup (). In brief, we drive forced evaporative cooling at a magnetic field , corresponding to (). In this phase is oriented along the vertical axis. At the end of the evaporation, we obtain a BEC of atoms with a condensed fraction above .

To reach the regime, we slowly modify, in the last step of the evaporation, the confining potential to the final cigar shape, with typical frequencies Hz. Simultaneously, we decrease to 0.8 G () and then change the magnetic-field orientation to the weak trapping axis () while keeping its amplitude constant sup (). Finally, we ramp to the desired target value (and equivalently ) in sup (), hold for a time and perform absorption imaging of our gas after a time-of-flight (ToF) of . Two imaging setups are used in order to measure either the density distribution integrated along the dipoles (-imaging) or perpendicular to them (-imaging) sup (). Figure 1 (lower inset) illustrates the final geometry of our system with , giving , and defines the relevant axes.

Here, we explore the properties of the system when the DDI is made attractive enough to overcome the repulsive MF contact interaction (, ), after adiabatically changing (ms) or quenching (ms) to its target value sup (). For ms, the system evolves following its ground state and gives access to the slow dynamics, whereas for the 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 Pitaevskii and Stringari (2003). 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 (), described as a “collapse”, which can develop into a rapid loss of coherence Roberts *et al.* (2001); Koch *et al.* (2008), and pattern formations, such as anisotropic atom bursts (”Bosenova”) and special -wave-type structures, as observed in rubidium Donley *et al.* (2001) and dipolar gases of chromium Koch *et al.* (2008); Lahaye *et al.* (2008), respectively. This fast dynamics has been proved to be well encompassed by GPE simulation Ueda and Saito (2003); Lahaye *et al.* (2008); Metz *et al.* (2009).

## Iii Density Distribution

In a first set of experiments, we study the evolution of the ToF density distribution of our dipolar Er BEC with . Figure 2 (a-c) shows the 2D column density profiles acquired with -imaging for ms and ms (), and for different . To analyze our data, we fit the 2D density profiles to bimodal distributions. We illustrate their results on the central cuts () in Fig. 2(d). For (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 below , 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 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 . 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 () (c), and exhibits a long-lived character from several tens to hundreds of ms.
As discussed in details later and suggested in Ref. Baillie *et al.* (2016), 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 . 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 Pitaevskii and Stringari (2003); Grimm (2007). 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 () with frequency . The axial oscillation comes along with a smaller-amplitude oscillation 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 Pitaevskii and Stringari (2003). The compression character is particularly relevant since it involves a change in the density and it is therefore sensitive to the LHY corrections Pitaevskii and Stringari (1998).

We excite the axial mode either by ramping during the final preparation stage or by transiently increasing the power of the vertical optical dipole trap (ODT) beam, after ramping to . 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 and image its ToF density distribution with -imaging. To extract , we probe the axial width of the central coherent component of the gas sup () with and fit it to a damped-sine; see inset in Fig. 3(b).

Figure 3 shows the observed normalized to the trapping frequency ^{3}^{4}

## 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 . 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 Wächtler and Santos (2016a); Baillie *et al.* (2016)

(1) |

where is the non-interacting Hamiltonian and the harmonic confinement with .
The MF chemical potential, , results from the competition between
short-range interactions, controlled by the coupling constant ,
and the DDI term with and
the angle sustained by and the dipole moment . Here . The beyond-MF physics is encoded in the LHY term, leading to an additional repulsive term in the chemical potential , . The function
is obtained from the LHY correction in homogeneous 3D dipolar BECs Lima and Pelster (2011, 2012) using local-density approximation ^{5}*et al.* (1998). In our calculations, we use the experimentally determined values of the 3B recombination rate of the condensate, sup ().

As discussed in Refs. Wächtler and Santos (2016b); Baillie *et al.* (2016), due to the repulsive LHY term, Eq. (1) sustains stable ground-state solutions for any and .
For pancake traps (), 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 Kadau *et al.* (2016). 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 () Eq. (1) has only one possible solution. In the parameter space, the corresponding wave function exhibits a smooth crossover from a dilute BEC to a single, high-density, macro droplet solution for increasing . 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 case, which corresponds to our experimental setting.

The continuous and smooth change of the static properties of the system with increasing 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 with atoms. We simulate a linear ramp in from to a variable final value of in , followed by a compression of the axial trap from Hz to Hz during ms. We record then the axial width from the standard deviation of , , as a function of the subsequent holding time. The evolution of is well fit by a sinusoidal function, whose frequency constitutes our theoretical prediction of .

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 , leaving the system in the high-density 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 Wächtler and Santos (2016a); Bisset *et al.* (2016); Wächtler and Santos (2016b); Baillie *et al.* (2016).
This provides a qualitatively different phase diagram than the one from the MF treatment, which predicts a collapse at .
For 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 , the effect of QFs is visible in a sizable shift of , which better matches the experimental data; see Fig. 3(b).

## 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 () and thermal () components along the BEC-to-droplet crossover.
Since in the droplet regime the core density dramatically increases, 3B losses are expected to play an important role even for moderate and low values of Baillie *et al.* (2016).
Notwithstanding, 3B losses and QFs exhibit different power dependencies on (see Eq. (1)) and thus the atom-loss dynamics should disclose the competition between diluteness and a droplet character.

Figure 4 (a-b) shows and , extracted from the double Gaussian fit as a function of after a non-adiabatic (a) and adiabatic (b) change of . Both cases show a similar evolution. When lowering , is first constant for , then shows a sharp drop starting around and finally saturates for lower . We note that in the non-adiabatic case, decreases slower as compared to the adiabatic case because of the shorter timing involved. Remarkably, remains mainly unaltered over the whole range of 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 as a function of with and without the LHY term. Remarkably, the observed evolution of 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 too large compared to the measured ones; see Fig. 4(a,b,d).

Finally, we investigate the in-trap time evolution of after quenching in the droplet regime; see Fig. 4 (c). Our measurements reveal three different timescales for the losses. At very short ( ms), is constant and the droplet state preserves its high density. It follows a fast decay dominated by 3B ( ms), having the effect of ejecting atoms from the core until decreases slightly below atoms. At this point, the loss dynamics substantially slows down ( ms).

Using the relation , we extract the mean in-situ density of the high-density component in the BEC-to-droplet crossover. Figure 4 (d) shows our results for ms. We observe a prominent increase of across the threshold , and a surviving high density state deep into the MF instability regime. The formation of the droplet state is particularly visible for the -ms case. Here, grows from at to a maximum of at , while it is slightly reduced to at . 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 self-bound by its underlying interaction in analogy with the He-droplet case Wächtler and Santos (2016a); Bisset *et al.* (2016); Wächtler and Santos (2016b); Baillie *et al.* (2016). 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 . To preserve the high density of the coherent component, our measurements focus on short timescales with ms and ms. After preparing the system at the desired , we abruptly switch off the ODT, let the gas expand for a variable , 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 (), we compute the second moments along . Figure 5 (a) exemplifies the ToF evolution of at , , and . When entering the 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 (i. e. ). From , we extract the value of the expansion velocity, , by fitting the data to . Figure 5 (b) shows in an range from to . When the system approach the droplet regime with decreasing scattering length (), undergoes a strong reduction and drops to a minimum equal to at about (). For further lowering of , starts to increase again. A similar behavior is observed for . We note that only the high-density component reveals this intriguing dependency on , 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 sup (). We record the evolution of with and extract the corresponding expansion velocities from the asymptotic behavior of .
Our simulations show a slowing down with a minimum of at (), followed by an increase at lower ^{6}

The expansion behavior can be qualitatively well understood considering the so-called released, or internal energy, .
This is the energy of the system when substracting the energy related with the confinement Pitaevskii and Stringari (2003).
In the MF scenario, as long as the ground state is stable. The BEC expands ballistically and decreases for decreasing and .
In the unstable regime, the expansion velocity depends crucially on the value of 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 is expected to be suppressed. Assuming a fixed (*i.e.* no 3B losses), one can show that
decreases with decreasing and eventually reaches for , marking the onset of the self-bound (SB) solution (e.g. for ) sup ().
However, in stark contrast to the MF case, increases with decreasing in the droplet regime. We note that is then shifted to lower values when 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 for decreasing at a fixed , and the increase of for decreasing in the droplet regime. In the crossover regime, the system smoothly evolves towards a fully self-bound state () 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, R. van Bijnen, R. N. Bisset and B. Blakie, E. Demler, R. Grimm, T. Pfau and the Stuttgart Dy team for multiple useful discussions. We thank G. Faraoni for her support in the final stage of the experiment. The Innsbruck group thanks the European Research Council within the ERC Consolidator Grant RARE no. 681432. LC is supported within the Marie Curie Project DIPPHASE no. 706809 of the European Commission. The Hannover group thanks the support of the DFG (RTG 1729). Both groups acknowledge the support of the DFG/FWF (FOR 2247).

*

## Appendix A Supplementary material

### a.1 Bose-Einstein condensation of

We prepare an ultracold gas of the isotope following a similar trapping and cooling scheme as the one employed for Er Aikawa *et al.* (2012); Frisch *et al.* (2012). We load a crossed-ODT from a narrow-line MOT of atoms at about .
At the end of the MOT sequence, the atoms are automatically spin-polarized in their lowest Zeeman sub-level Frisch *et al.* (2012).
The dipole orientation follows the one of the external applied magnetic field, . In our experiment, the latter is controlled by independent tuning of the components along the experimental coordinate system , 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 nm. One beam propagates horizontally along the -axis, and the other propagates vertically and nearly collinear to the -axis. The -beam has a maximum power of 10 W and an elliptical profile defined by its waists of . The -beam has a maximum power of 27 W, a vertical waist , and a tunable horizontal waist, . The latter can be conveniently tuned from to 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 ( 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 atoms by means of evaporative cooling in the crossed ODT at G (). Typically, we first rapidly (in 600 ms) reduce the power and aspect ratio of the -beam from 24 W to 4 W and to , respectively. We further decrease the power of the -beam from 4 W to 0.3 W in s in an exponential manner and then exponentially increase the aspect ratio from 1.6 to 8 in 2.5 s. The final trap frequencies are typically of Hz. We finally obtain BECs of typically atoms with more than 80% condensed fraction and a temperature nK. We typically measure and the condensed fraction from a bimodal fit of the 2D column density distribution measured along -imaging with ms. is extracted from the evolution of the thermal size of the bimodal fit with varying from to ms.

### a.2 Experimental setup and axes

In our setup we define the orthonormal -coordinate system in the following way: the vertical axis is aligned with gravity and the 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 -coordinate system with and rotated by a small angle as compared to . The magnetic field components , each created by each pair of coils, can be controlled independently. We estimate to be about 15 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 causes a small reduction of the mean DDI energy and corresponding small shifts of the expected characteristics compared to the one predicted for : the MF collapse threshold should appear at a lower and, for a given , and are shifted respectively down and up. We have experimentally evaluated the shift of deep in the stable BEC regime (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 compared to and compared to in the -plane. The -imaging axis is tilted by , mainly in the -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 because of a mixed projection of and , the two first velocity components in the -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.

### a.3 Precision measurements of the -to- conversion in a three-dimensional optical lattice

Precise determination of the -to- 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 should thus rely on experimental investigations.
We perform lattice modulation spectroscopy in a three-dimensional optical lattice. From the measurements of the energy gap in the Mott insulator state we extract .
Our lattice experiments focuses in the region of low -field [ G] and are based on a lattice setup and procedure similar to the one described in Ref. Baier *et al.* (2016).
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 , given in unit of the respective recoil energies kHz. At these lattice dephts, the gas is in the Mott insulator state. We then vary to the desired value by rapidly changing in 2 ms, either just before or just after loading the lattice. We use the latter option for the smallest values at which 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 at a variable frequency for 90 ms with a peak-to-peak amplitude of about . Finally, we ramp down the lattice depths to zero in 150 ms, and measure the recovered condensed fraction as a function of from -imaging ToF picture. For the smallest values considered, we also ramp back to 2 G in 2 ms before switching off the lattice-beams in order to again minimize 3B loss effects.

When varying , we observe a resonant depletion of the condensate due to particle-hole excitations. The resonance condition in the Mott-insulator regime is given by

(2) |

Here , and are respectively the on-site contact interaction, the on-site dipolar interaction and the nearest-neighbor dipolar interaction along in the corresponding extended Bose-Hubbard model. and can be accurately predicted from the knowledge of the lattice depths and dipole orientation and in our typical experimental condition, they are equal to Hz and Hz respectively. By subtracting the theoretical dipolar contributions to the measured frequency, we extract , which is directly proportional to the scattering length . A precise mapping of in the ultralow -field region is then obtained by repeating the above measurements at various values; see Fig. 1.

In the region [] G, our lattice spectroscopy reveals the presence of two FRs, one located at about zero 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 G and G.

The scattering length of Er can be parametrized by the following simple expression Chin *et al.* (2010)

(3) |

in which the specific positions () and widths () 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 mG, mG, G, mG. The -dependent background scattering length accounts for the overlapping resonances and reads as with . We also account for the small effect of the two next resonances, whose positions and widths 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. 3.

### a.4 Ramps in scattering length

Our measurements rely on controlled variations of the scattering length . In our experiments, we either adiabatically change using ms or we quench it using ms. The adiabatic condition for reads as

(4) |

As shown in Fig. 6, we use two different types of time variations of and thus of : (i) a simple linear ramp in and (ii) we design a specific variation in order to minimize the adiabaticity parameter . The resulting shows an exponential-type variation with . The adiabaticity condition of Eq. (4) is more stringent for lower . For ramping down to , we find that (i) verifies Eq. (4) for ms and (ii) for ms. Data from Figs. 1, 3 and 4 (b-d) (resp. Figs. 2, 4 (a) and 5) use ramp (i) (resp. (ii)).

For our theoretical description, we use a linear change of , similar to case (ii).

### a.5 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 . We first prepare a non-degenerate thermal sample of Er atoms at nK in an harmonic trap. We then record the decay of the atom number, , as a function of in a range from 0 to 1 s. is obtained from a Gaussian fit to the measured density distribution.

To fit the time evolution of , we use the integrated 3B rate equation, which reads as . Here, is the mean square density on the cloud. To describe the scaling of with , we use its prediction for an ideal gas at thermal equilibrium at whose state occupancies follow Boltzmann law and take into account the anti-evaporation effect Weber *et al.* (2003). Then is extracted from a fit of along:

(5) |

where is the atom number at ms and is defined via the relation

(6) |

with the Boltzmann constant and .

We account for the -dependence of near a FR by repeating the measurement at different . We check that the measured does not depend on the orientation and measure its B-dependency using and for varying in 0.1 to 1.9 G. Note that, due to the bosonization effect, the in a quantum degenerate bosonic gas is equal to . Figure 7 shows in a quantum degenerate bosonic gas of as a function of using the measured -to- conversion. Despite the existence of many coupled molecular potential in Er, we measure a comparable low as compared to the typical values reported in alkali atoms. varies between ms at and ms at .

### a.6 Time-of-Flight measurements

For our ToF measurements, we abruptly extinguish the ODT in about . 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 , remains unchanged and the dynamics occur at the original and dipole orientation. At time before the image is taken, 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 of the imaging field is chosen constant and equal to 0.31 G for all considered. is set to be 12 ms for -imaging and 15 ms for -imaging where the change in is more drastic for the typical dipole orientation () used in this experiment.

We note that our resolution limit for both - and -imaging is estimated to be . Moreover the effective pixel sizes are set to and in our setup. These limit the size of the structure we are able to observe as well as the minimal we can use, which is typically ms.

### a.7 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 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 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.

### a.8 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 of the BEC, along with smaller amplitude oscillations of the radial size in phase opposition. We extract its frequency by studying the larger amplitude oscillation of . For this, we probe the ToF density distribution of the gas with -imaging after a ToF of 30ms. We focus on the sizes of the central, high-density component of the cloud, which we study as a function of for different . We note that the precise shape of the column density profile is expected to change as a function of and this in a different way for the two axis and 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 . Typically, we use a bimodal MF TF plus Gauss fit for . For 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 with at fixed and thus enable to determine .

To fit , we use a damped-sine function of . Typically we fit the evolution of for varying from 0 to few hundreds of ms, depending on the damping rate observed. The upper value of 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 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.

### a.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 . 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 : (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 writes Pitaevskii and Stringari (2003).(ii) 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 considered its width is typically at least 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.

The quality of the bimodal fit is evaluated by the norm of the residue . For both distributions, it satisfies .

### a.10 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 with . There our aim was to extract the mean density and we did not expand on the mere description of the .

We note that is well accounted for by a double exponential decay evolution of respective amplitude and . is the initial atom number. Each decay corresponds to a different time constant, respectively and , and starts after a different delay time, respectively and . We fix . For ms is approximately constant and equal to ms. The absence of evolution for indicates that the magnetic field may not have reached the target value in the first ms. For ms, can be set to 0. , and (and ) depend on , typically decreasing with it. The evolution is illustrated in Fig. 8 for ms.

### a.11 Gaussian Ansatz including

A good qualitative (and to a large extent quantitative) insight in the physics of dipolar condensates in the presence of LHY stabilization may be gained from a simplified Gaussian ansatz of the form

(7) |

where the variational parameters are the number of atoms , the global phase , the widths , and the phase curvatures . The Lagrangian density reads:

(8) |

We insert ansatz (7) into Eq. (8), obtaining the Lagrangian :

(9) | |||||

with

(10) |

The variational parameters are determined from the corresponding Euler-Lagrange equations Filho *et al.* (2001):

(11) |

where , and . Introducing the dimensionless units , , , with , and after some algebra, we obtain a close set of equations for the Gaussian widths and the number of atoms:

(12) |

(13) |

with , and

(14) | |||||

with and .

Due to its simplicity, Eqs. (12) and (13) 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.

### a.12 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, . In absence of losses, we may evaluate by means of the Gaussian Ansatz for the ground-state of a trapped BEC with scattering length and fixed . Figure 9 shows the results for for different