Long-lived and transient supersolid behaviors in dipolar quantum gases

TL;DR

This study combines theory and experiments to demonstrate long-lived supersolid states in dipolar quantum gases of erbium and dysprosium, revealing differences in their stability and lifetime due to atom loss rates.

Contribution

It provides the first detailed experimental observation of long-lived supersolid behavior in dipolar quantum gases and identifies the parameter regimes for different phases.

Findings

Supersolid states coexist with density modulation and phase coherence.

Dysprosium exhibits longer supersolid lifetime than erbium.

Long-lived supersolids achieved via evaporative cooling from thermal gases.

Abstract

By combining theory and experiments, we demonstrate that dipolar quantum gases of both $^{166}$Er and $^{164}$Dy support a state with supersolid properties, where a spontaneous density modulation and a global phase coherence coexist. This paradoxical state occurs in a well defined parameter range, separating the phases of a regular Bose-Einstein condensate and of an insulating droplet array, and is rooted in the roton mode softening, on the one side, and in the stabilization driven by quantum fluctuations, on the other side. Here, we identify the parameter regime for each of the three phases. In the experiment, we rely on a detailed analysis of the interference patterns resulting from the free expansion of the gas, quantifying both its density modulation and its global phase coherence. Reaching the phases via a slow interaction tuning, starting from a stable condensate, we observe that $^{166}$Er and $^{164}$Dy exhibit a striking difference in the lifetime of the supersolid properties, due to the different atom loss rates in the two systems. Indeed, while in $^{166}$Er the supersolid behavior only survives a few tens of milliseconds, we observe coherent density modulations for more than $150\,$ms in $^{164}$Dy. Building on this long lifetime, we demonstrate an alternative path to reach the supersolid regime, relying solely on evaporative cooling starting from a thermal gas.

Tables

Table 1. Table 1 : Extracted 1/10-lifetime of Er 166 superscript Er 166 {}^{166}{\rm Er} atom number decay for t h ≥ 0 subscript 𝑡 h 0 t_{\rm h}\geq 0 and remaining atom number at long holding time for data in Fig.  AS6 .

$a_s\left(a_0\right)$	$t_N$ (ms)	$N_r\left({10}^4\right)$	$t_{\Phi }$ (ms)
53.3(2)	32(5)	1.03(5)	-
54.0(2)	51(9)	1.29(11)	25(6)
54.2(2)	46(12)	1.7(2)	32(9)

Table 2. Table 2 : Extracted 1/10-lifetime of Dy 164 superscript Dy 164 {}^{164}{\rm Dy} atom number decay and A Φ subscript 𝐴 Φ A_{\Phi} decay for data in Fig.  AS7 .

$B$ (G)	$t_N$ (ms)	$t_{\Phi }$ (ms)
1.8	300(12)	12(5)
2.04	728(34)	152(13)
2.1	926(36)	133(25)

Equations

$$S_i$$

$$=$$

$$J_i$$

Full text

I Introduction

Supersolidity is a paradoxical quantum phase of matter where both crystalline and superfluid order coexist Andreev and Lifshitz (1969); Chester (1970); Leggett (1970). Such a counter-intuitive phase, featuring rather antithetic properties, has been originally considered for quantum crystals with mobile bosonic vacancies, the latter being responsible for the superfluid order. Solid 4He has been long considered a prime system to observe such a phenomenon Kirzhnits and Nepomnyashchii (1971); Schneider and Enz (1971). However, after decades of theoretical and experimental efforts, an unambiguous proof of supersolidity in solid 4He is still missing Balibar (2010); Boninsegni and Prokof’ev (2012).

In search of more favorable and controllable systems, ultracold atoms emerged as a very promising candidate, thanks to their highly tunable interactions. Theoretical works point to the existence of a supersolid ground state in different cold-atom settings, including dipolar Lu et al. (2015a) and Rydberg particles Henkel et al. (2010); Cinti et al. (2010), cold atoms with a soft-core potential Boninsegni (2012), or lattice-confined systems Boninsegni and Prokof’ev (2012). Breakthrough experiments with Bose–Einstein condensates (BECs) coupled to light have recently demonstrated a state with supersolid properties Léonard et al. (2017); Li et al. (2017). While in these systems indeed two continuous symmetries are broken, the crystal periodicity is set by the laser wavelength, making the supersolid incompressible.

Another key notion concerns the close relation between a possible transition to a supersolid ground state and the existence of a local energy minimum at large momentum in the excitation spectrum of a non-modulated superfluid, known as roton mode Landau (1941). Since excitations corresponding to a periodic density modulation at the roton wavelength are energetically favored, the existence of this mode indicates the system’s tendency to crystallize  Nozières (2004) and it is predicted to favor a transition to a supersolid ground state Kirzhnits and Nepomnyashchii (1971); Schneider and Enz (1971); Henkel et al. (2010).

Remarkably, BECs of highly magnetic atoms, in which the particles interact through the long-range and anisotropic dipole-dipole interaction (DDI), appear to gather several key ingredients for realizing a supersolid phase. First, as predicted more than fifteen years ago O’Dell et al. (2003); Santos et al. (2003) and recently demonstrated in experiments Chomaz et al. (2018); Petter et al. (2018), the partial attraction in momentum space due to the DDI gives rise to a roton minimum. The corresponding excitation energy, i. e. the roton gap, can be tuned in the experiments down to vanishing values. Here, the excitation spectrum softens at the roton momentum and the system becomes unstable. Second, there is a non-trivial interplay between the trap geometry and the phase diagram of a dipolar BEC. For instance, our recent observations have pointed out the advantage of axially-elongated trap geometries (i. e. cigar-shaped) compared to the typically considered cylindrically-symmetric ones (i. e.  pancake-shaped) in enhancing the visibility of the roton excitation in experiments. Last but not least, while the concept of a fully softened mode is typically related to instabilities and disruption of a coherent quantum phase, groundbreaking works in the quantum-gas community have demonstrated that quantum fluctuations can play a crucial role in stabilizing a dipolar BEC Kadau et al. (2016); Ferrier-Barbut et al. (2016); Chomaz et al. (2016); Wächtler and Santos (2016a, b); Schmitt et al. (2016); Ferrier-Barbut et al. (2018a). Such a stabilization mechanism enables the existence, beyond the mean-field instability, of a variety of stable ground states, from a single macro-droplet Bisset et al. (2016); Chomaz et al. (2016); Wächtler and Santos (2016b) to striped phases Wenzel et al. (2017), and droplet crystals Baillie and Blakie (2018); see also related works Petrov (2015a); Cabrera et al. (2018); Semeghini et al. (2018); Cheiney et al. (2018). For multi-droplet ground states, efforts have been devoted to understand if a phase coherence among ground-state droplets could be established Wenzel et al. (2017); Baillie and Blakie (2018). However, previous experiments with ${}^{164}{\rm Dy}$ have shown the absence of phase coherence across the droplets Wenzel et al. (2017), probably due to the limited atom numbers.

Droplet ground-states, quantum stabilization, and dipolar rotons have raised a huge excitement with very recent advancements adding key pieces of information to the supersolid scenario. The quench experiments in an ${}^{166}{\rm Er}$ BEC at the roton instability have revealed out-of-equilibrium modulated states with an early-time phase coherence over a timescale shorter than a quarter of the oscillation period along the weak-trap axis Chomaz et al. (2018). In the same work, it has been suggested that the roton softening combined with the quantum stabilization mechanism may open a promising route towards a supersolid ground state. A first confirmation came from a recent theoretical work Roccuzzo and Ancilotto (2018), considering an Er BEC in an infinite elongated trap with periodic boundary conditions and tight transverse confinement. The supersolid phase appears to exist within a narrow region in interaction strength, separating a roton excitation with a vanishing energy and an incoherent assembly of insulating droplets. Almost simultaneously, experiments with ${}^{162}{\rm Dy}$ BECs in a shallow elongated trap, performing a slow tuning of the contact interaction, reported on the production of stripe states with phase coherence persisting up to half of the weak trapping period Tanzi et al. (2018). More recently, such observations have been confirmed in another ${}^{162}{\rm Dy}$ experiment Böttcher et al. (2019). Here, theoretical calculations showed the existence of a phase-coherent droplet ground-state, linking the experimental findings to the realization of a state with supersolid properties. The results on ${}^{162}{\rm Dy}$ show however transient supersolid properties whose lifetime is limited by fast inelastic losses caused by three-body collisions Tanzi et al. (2018); Böttcher et al. (2019). These realizations raise the crucial question of whether a long-lived or stationary supersolid state can be created despite the usually non-negligble atom losses and the crossing of a discontinuous phase transition, which inherently creates excitations in the system.

In this work, we study both experimentally and theoretically the phase diagram of degenerate gases of highly magnetic atoms beyond the roton softening. Our investigations are carried out using two different experimental setups producing BECs of 166Er  Aikawa et al. (2012); Chomaz et al. (2016) and of 164Dy  Trautmann et al. (2018), and rely on a fine tuning of the contact-interaction strength in both systems. In the regime of interest, these two atomic species have different contact-interaction scattering lengths, $a_{\rm s}$, whose precise dependence on the magnetic field is known only for Er Baier et al. (2016); Chomaz et al. (2016, 2018), and different three-body-loss rate coefficients. Moreover, Er and Dy possess different magnetic moments, $\mu$, and masses, $m$, yielding the dipolar lengths, $a_{\rm dd}=\mu_{0}\mu^{2}m/12\pi\hbar^{2}$, of $65.5\,a_{0}$ and $131\,a_{0}$, respectively. Here $\mu_{0}$ is the vacuum permeability, $\hbar=h/2\pi$ the reduced Planck constant, and $a_{0}$ the Bohr radius. For both systems, we find states showing hallmarks of supersolidity, namely the coexistence of density modulation and global phase coherence. For such states, we quantify the extent of the $a_{\rm s}$-parameter range for their existence and study their lifetime. For ${}^{166}{\rm Er}$, we find results very similar to the one recently reported for ${}^{162}{\rm Dy}$ Tanzi et al. (2018); Böttcher et al. (2019), both systems being limited by strong three-body losses, which destroy the supersolid properties in about half of a trap period. However, for ${}^{164}{\rm Dy}$, we have identified an advantageous magnetic-field region where losses are very low and large BECs can be created. In this condition, we observe that the supersolid properties persist over a remarkably long time, well exceeding the trap period. Based on such a high stability, we finally demonstrate a novel route to reach the supersolid state, based on evaporative cooling from a thermal gas.

II Theoretical description

As a first step in our study of the supersolid phase in dipolar BECs, we compute the ground-state phase diagram for both ${}^{166}{\rm Er}$ and ${}^{164}{\rm Dy}$ quantum gases. The gases are confined in a cigar-shaped harmonic trap, as illustrated in Fig. 1(a). Our theory is based on numerical calculations of the extended Gross-Pitaevskii equation (eGPE) sup , which includes our anisotropic trapping potential, the short-range contact and long-range dipolar interactions at a mean-field level, as well as the first-order beyond-mean-field correction in the form of a Lee-Huang-Yang (LHY) term Wächtler and Santos (2016a, b); Bisset et al. (2016); Chomaz et al. (2016, 2018). We note that, while both the exact strength of the LHY term and its dependence on the gas characteristics are under debate Schmitt et al. (2016); Chomaz et al. (2018); Ferrier-Barbut et al. (2018b); Cabrera et al. (2018); Petter et al. (2018), the importance of such a term, scaling with a higher power in density, is essential for stabilizing states beyond the mean-field instability Schmitt et al. (2016); Chomaz et al. (2018); Ferrier-Barbut et al. (2018b); see also Gammal et al. (2000); Bulgac (2002); Lu et al. (2015a); Petrov (2015b).

Our theoretical results are summarized in Fig. 1. By varying the condensed-atom number, $N$, and $a_{\rm s}$, the phase diagram shows three very distinct phases. To illustrate them, we first describe the evolution of the integrated in-situ density profile $n(y)$ with fixed $N$ for varying $a_{\rm s}$, Fig. 1(b). The first phase, appearing at large $a_{\rm s}$, resembles a regular dilute BEC. It corresponds to a non-modulated density profile of low peak density and large axial size, $\sigma_{y}$, exceeding several times the corresponding harmonic oscillator length ($\ell_{y}=\sqrt{\hbar/m\omega_{y}}$), see Fig. 1(e) and the region denoted BEC in (f) and (g). The second phase appears when decreasing $a_{\rm s}$ down to a certain critical value, $a_{\rm s}^{*}$. Here, the system undergoes an abrupt transition to a periodic density-modulated ground state, consisting of an array of overlapping narrow droplets, each of high peak density. Because the droplets are coupled to each other via a density overlap, later quantified in terms of the link strength $S$, particles can tunnel from one droplet to a neighboring one, establishing a global phase coherence across the cloud; see Fig. 1(d). Such a phase, in which periodic density modulation and phase coherence coexist, is identified as the supersolid one Cinti et al. (2010); Roccuzzo and Ancilotto (2018); SSP region in (f) and (g). When further decreasing $a_{\rm s}$, we observe a fast reduction of the density overlap, which eventually vanishes; see Fig. 1(c). Here, the droplets become fully separated. Under realistic experimental conditions, it is expected that the phase relation between such droplets cannot be maintained; see later discussion. We identify this third phase as the one of an insulating droplet array Bisset et al. (2016); Wächtler and Santos (2016); Wenzel et al. (2017); ID region in (f) and (g). The number of droplets in the array decreases with lowering either $a_{\rm s}$ (see (b)) or $N$, eventually resulting in a single droplet of high peak density, as in Refs. Wächtler and Santos (2016b); Bisset et al. (2016); see dark blue region in (f). The existence of these three phases (BEC-SSP-ID) is consistent with recent calculations considering an infinitely elongated Er BEC Roccuzzo and Ancilotto (2018) and a cigar-shaped ${}^{162}{\rm Dy}$ BEC Böttcher et al. (2019), illustrating the generality of this behavior in dipolar gases.

To study the supersolid character of the density-modulated phases, we compute the average of the wavefunction overlap between neighboring droplets, $S$. As an ansatz to extract $S$, we use a Gaussian function to describe the wavefunction of each individual droplet. This is found to be an appropriate description from an analysis of the density profiles of Fig. 1(b-d); see also Wenzel et al. (2018). For two droplets at a distance $d$ and of identical Gaussian widths, $\sigma_{y}$, along the array direction, $S$ is simply $S=\exp(-d^{2}/4\sigma_{y}^{2})$. Here, we generalize the computation of the wavefunction overlap to account for the difference in widths and amplitudes among neighboring droplets. This analysis allows to distinguish between the two types of modulated ground states, SSP and ID in Fig. 1(f–g). Within the Josephson-junction picture Josephson (1962); Javanainen (1986, 1986); Raghavan et al. (1999), the tunneling rate of atoms between neighboring droplets depends on the wavefunction overlap, and an estimate for the single-particle tunneling rate can be derived within the Gaussian approximation Wenzel et al. (2018); see also sup . The ID phase corresponds to vanishingly small values of $S$, yielding tunneling times extremely long compared to any other relevant time scale. In contrast, the supersolid phase is identified by a substantial value of $S$, with a correspondingly short tunneling time.

As shown in Fig. 1(f–g), a comparative analysis of the phase diagram for ${}^{166}{\rm Er}$ and ${}^{164}{\rm Dy}$ reveals similarities between the two species (see also Ref. Böttcher et al. (2019)). A supersolid phase is found for sufficiently high $N$, in a narrow region of $a_{\rm s}$, upper-bounded by the critical value $a_{\rm s}^{*}(N)$. For intermediate $N$, $a_{\rm s}^{*}$ increases with increasing $N$. We note that, for low $N$, the non-modulated BEC evolves directly into a single droplet state for decreasing $a_{\rm s}$ foo (a). In this case, no supersolid phase is found in between, see also Refs. Wächtler and Santos (2016b); Bisset et al. (2016). Despite the general similarities, we see that the supersolid phase for ${}^{164}{\rm Dy}$ appears for lower atom number than for Er and has a larger extension in $a_{\rm s}$. We note that, at large $N$ and for decreasing $a_{\rm s}$, Dy exhibits ground states with a density modulation appearing first in the wings, which then progresses inwards until a substantial modulation over the whole cloud is established foo (b); see inset (g). In this regime, we also observe that $a_{\rm s}^{*}$ decreases with increasing $N$. These type of states have not been previously reported and, although challenging to access in experiments because of the large $N$, they deserve further theoretical investigations.

III Experimental sequence for 166Erbium and 164Dysprosium

To experimentally access the above-discussed physics, we produce dipolar BECs of either ${}^{166}{\rm Er}$ or ${}^{164}{\rm Dy}$ atoms. These two systems are created in different setups and below we summarize the main experimental steps; see also sup .

Erbium – We prepare a stable ${}^{166}{\rm Er}$ BEC following the scheme of Ref. Chomaz et al. (2018). At the end of the preparation, the Er BEC contains about $N=8\times 10^{4}$ atoms at $a_{\rm s}=64.5\,a_{0}$. The sample is confined in a cigar-shaped optical dipole trap with harmonic frequencies $\omega_{x,y,z}=2\pi\times(145,31.5,151)\,$Hz. A homogeneous magnetic field, $B$, polarizes the sample along $z$ and controls the value of $a_{\rm s}$ via a magnetic Feshbach resonance (FR) Chomaz et al. (2016, 2018); sup . Our measurements start by linearly ramping down $a_{\rm s}$ within $20\,$ms and waiting additional $15\,$ms so that $a_{\rm s}$ reaches its target value sup . We note that ramping times between $20$ and $60\,$ms have been tested in the experiment and we do not record a significant difference in the system’s behavior. After the 15 ms stabilization time, we then hold the sample for a variable time $t_{\rm h}$ before switching off the trap. Finally, we let the cloud expand for 30 ms and perform absorption imaging along the $z$ (vertical) direction, from which we extract the density distribution of the cloud in momentum space, $n(k_{x},k_{y})$.

Dysprosium – The experimental procedure to create a ${}^{164}{\rm Dy}$ BEC follows the one described in Ref. Trautmann et al. (2018); see also sup . Similarly to Er, the Dy BEC is also confined in a cigar-shaped optical dipole trap and a homogeneous magnetic field ${B}$ sets the quantization axis along $z$ and the value of $a_{\rm s}$. For Dy, we will discuss our results in terms of magnetic field, $B$, since the $a_{\rm s}$-to-$B$ conversion is not well known in the magnetic-field range considered sup ; Tang et al. (2015); Schmitt et al. (2016); Ferrier-Barbut et al. (2018b). In a first set of measurements, we first produce a stable BEC of about $N=3.5\times 10^{4}$ condensed atoms at a magnetic field of $B=2.5\,$G and then probe the phase diagram by tuning $a_{\rm s}$. Here, before ramping the magnetic field to access the interesting $a_{\rm s}$ regions, we slowly increase the power of the trapping beams within $200\,$ms. The final trap frequencies are $\omega_{x,y,z}=2\pi\times(300,16,222)\,$Hz. After preparing a stable BEC, we ramp $B$ to the desired value within $20\,$ms and hold the sample for $t_{\rm h}$ sup . In a second set of measurements, we study a completely different approach to reach the supersolid state. As discussed later, here we first prepare a thermal sample at a ${B}$ value where supersolid properties are observed and then further cool the sample until a transition to a coherent droplet-array state is reached. In both cases, at the end of the experimental sequence, we perform absorption imaging after typically $27\,$ms of time-of-flight (TOF) expansion. The imaging beam propagates horizontally under an angle $\alpha$ of $\approx 45^{\rm o}$ with respect to the weak axis of the trap ($y$). From the TOF images, we thus extract $n(k_{Y},k_{z})$ with $k_{Y}=\cos(\alpha)k_{y}+\sin(\alpha)k_{x}$.

A special property of ${}^{164}{\rm Dy}$ is that its background scattering length is smaller than $a_{\rm dd}$. This allows to enter the supersolid regime without the need of setting $B$ close to a FR, as done for ${}^{166}{\rm Er}$ and ${}^{162}{\rm Dy}$, which typically causes severe atom losses due to increased three-body loss coefficients. In contrast, in the case of ${}^{164}{\rm Dy}$, the supersolid regime is reached by ramping $B$ away from the FR pole used to produce the stable BEC via evaporative cooling, as the $a_{\rm s}$-range of Fig. 1(g) lies close to the background $a_{\rm s}$ reported in Ref. Tang et al. (2015); see also sup . At the background level, three-body loss coefficients below $1.3\times 10^{-41}\text{m}^{6}\text{s}^{-1}$ have been reported for ${}^{164}{\rm Dy}$ Schmitt et al. (2016).

IV Density Modulation and phase coherence

The coexistence of density modulation and phase coherence is the key feature that characterizes the supersolid phase and allows to discriminate it from the BEC and ID cases. To experimentally probe this aspect in our dipolar quantum gases, we record their density distribution after a TOF expansion for various values of $a_{\rm s}$ across the phase diagram. As for a BEC in a weak optical lattice Greiner et al. (2002) or for an array of BECs Greiner et al. (2001); Paredes et al. (2004); Hadzibabic et al. (2004), the appearance of interference patterns in the TOF images is associated with a density modulation of the in-situ atomic distribution. Moreover, the shot-to-shot reproducibility of the patterns (in amplitude and position) and the persistence of fringes in averaged pictures, obtained from many repeated images taken under the same experimental conditions, reveals the presence of phase coherence across the sample Hadzibabic et al. (2004).

Figure 2 exemplifies snapshots of the TOF distributions for Er, measured at three different $a_{\rm s}$ values; see (a-c). Even if very close in scattering length, the recorded $n(k_{x},k_{y})$ shows a dramatic change in behavior. For $a_{\rm s}=54.7(2)\,a_{0}$, we observe a non-modulated distribution with a density profile characteristic of a dilute BEC. When lowering $a_{\rm s}$ to $53.8(2)\,a_{0}$, we observe the appearance of an interference pattern in the density distribution, consisting of a high central peak and two almost symmetric low-density side peaks Not (a). Remarkably, the observed pattern is very reproducible with a high shot-to-shot stability, as shown in the repeated single snapshots and in the average image (b and e). This behavior indicates a coexistence of density modulation and global phase coherence in the in-situ state, as expected in the supersolid phase. This observation is consistent with our previous quench experiments Chomaz et al. (2018) and with the recent ${}^{162}{\rm Dy}$ experiments Tanzi et al. (2018); Böttcher et al. (2019). When further lowering $a_{\rm s}$ to $53.3(2)\,a_{0}$, complicated patterns develop with fringes varying from shot-to-shot in number, position, and amplitude, signalizing the persistence of in-situ density modulation. However, the interference pattern is completely washed out in the averaged density profiles (f), pointing to the absence of a global phase coherence. We identify this behavior as the one of ID states.

Toy Model – To get an intuitive understanding of the interplay between density modulation and phase coherence and to estimate the role of the different sources of fluctuations in our experiment, we here develop a simple toy model, which is inspired by Ref. Hadzibabic et al. (2004); see also Ref. sup . In our model, the initial state is an array of $N_{D}$ droplets containing in total $N$ atoms. Each droplet is described by a one-dimensional Gaussian wavefunction, $\psi_{i}(y)$, of amplitude $\alpha_{i}$, phase $\phi_{i}$, width $\sigma_{i}$, and center $y_{i}$. To account for fluctuations in the experiments, we allow $\alpha_{i}$, $d_{i}=y_{i}-y_{i-1}$, and $\sigma_{i}$ to vary by 10% around their expectation values. The spread of the phases ${\phi_{i}}$ among the droplets is treated specially as it controls the global phase coherence of the array. By fixing $\phi_{i}=0$ for each droplet or by setting a random distribution of $\phi_{i}$, we range from full phase coherence to the incoherent cases. Therefore, the degree of phase incoherence can be varied by changing the standard deviation of the distribution of $\phi_{i}$.

To mimic our experiment, we compute the free evolution of each individual $\psi_{i}$ over $30\,$ms, and then compute the axial distribution $n(y,t)=|\sum_{i}\psi_{i}(y,t)|^{2}$, from which we extract the momentum distribution $n(k_{y})$, also accounting for the finite imaging resolution sup . For each computation run, we randomly draw $N_{D}$ values for $\phi_{i}$, as well as of $\sigma_{i}$, $d_{i}$ and $\alpha_{i}$ and extract $n(k_{y})$. We then collect a set of $n(k_{y})$ by drawing these values multiple times using the same statistical parameters and compute the expectation value, $\langle n(k_{y})\rangle$; see Fig. 2(j-l). The angled brackets denote the ensemble average.

The results of our toy model show large similarity with the observed behavior in the experiment. In particular, while for each single realization one can clearly distinguish multi-peak structures regardless of the degree of phase coherence between the droplets, the visibility of the interference pattern in the averaged $n(k_{y})$ survives only if the standard deviation of the phase fluctuations between droplets is small (roughly, below $0.3\pi$). In the incoherent case, we note that the shape of the patterns strongly vary from shot to shot. Interestingly, the toy model also shows that the visibility of the coherent peaks in the average images is robust against the typical shot-to-shot fluctuations in droplet size, amplitude and distance that occur in the experiments; see (j,k).

Probing density modulation and phase coherence – To separate and quantify the information on the in-situ density modulation and its phase coherence, we analyse the measured interference patterns in Fourier space  Takeda et al. (1982); Kohstall et al. (2011); Chomaz et al. (2015); Böttcher et al. (2019). Here, we extract two distinct averaged density profiles, $n_{\mathcal{M}}$ and $n_{\Phi}$. Their structures at finite $y$-spatial frequency (i. e. in Fourier space) quantify the two above-mentioned properties.

More precisely, we perform a Fourier transform (FT) of the integrated momentum distributions, $n(k_{y})$, denoted $\mathcal{F}[n](y)$. Generally speaking, modulations in $n(k_{y})$ induce peaks at finite spatial frequency, $y=y^{*}$, in the FT norm, $\left|\mathcal{F}[n](y)\right|$; see Fig. 2(g-i) and (m-o). Following the above discussion (see also Refs. Hadzibabic et al. (2004); Hofferberth et al. (2007)), such peaks in an individual realization hence reveal a density modulation of the corresponding in-situ state, with a wavelength roughly equal to $y^{*}$. Consequently, we consider the average of the FT norm of the individual images, $n_{\mathcal{M}}(y)=\langle|\mathcal{F}[n](y)|\rangle$ as the first profile of interest. The peaks of $n_{\mathcal{M}}$ at finite $y$ then indicate the mere existence of an in-situ density modulation of roughly constant spacing within the different realizations. As second profile of interest, we use the FT norm of the average profile $\langle n(k_{y})\rangle$, $n_{\Phi}(y)=|\mathcal{F}[\langle n\rangle](y)|$. Connecting to our previous discussion, the peaks of $n_{\Phi}$ at finite $y$ point to the persistence of a modulation in the average $\langle n(k_{y})\rangle$, which we identified as a hallmark for a global phase coherence within the density modulated state. In particular we point out that a perfect phase coherence, implying identical interference patterns in all the individual realizations, yields $n_{\mathcal{M}}=n_{\Phi}$ and, thus, identical peaks at finite $y$ in both profiles. We note that, by linearity, $n_{\Phi}$, also matches the norm of the average of the full FT of the individual images, i. e. $n_{\Phi}(y)=|\langle\mathcal{F}[n](y)\rangle|$; see also Ref. sup .

Figure 2 (g-i) and (m-o) demonstrates the significance of our FT analysis scheme by applying it to the momentum distributions from the experiment (d-f) and the ones from the toy model (j-l), respectively. As expected, for the BEC case, both $n_{\mathcal{M}}$ and $n_{\Phi}$ show a single peak at zero spatial frequency, $y=0$, characterizing the absence of density modulation (g). In the case of phase-coherent droplets (e), we observe that $n_{\mathcal{M}}$ and $n_{\Phi}$ are superimposed and both show two symmetric side peaks at finite $y$, in addition to a dominant peak at $y=0$; see (h). In the incoherent droplet case, we find that, while $n_{\mathcal{M}}$ still shows side peaks at finite $y$, the ones in $n_{\Phi}$ wash out from the averaging, (f, i, l, o). For both coherent and incoherent droplet arrays, the toy-model results show behaviors matching the above description, providing a further justification of our FT-analysis scheme; see (j-o). Our toy model additionally proves two interesting features. First, it shows that the equality $n_{\mathcal{M}}=n_{\Phi}$, revealing the global phase coherence of a density modulated state, is remarkably robust to noise in the structure of the droplet arrays; see (j, m). Second, our toy model however shows that phase fluctuations across the droplet array on the order of 0.2$\pi$ standard deviation are already sufficient to make $n_{\Phi}$ and $n_{\mathcal{M}}$ to deviate from each other; see (k, n). The incoherent behavior is also associated with strong variations in the side peak amplitude of the individual realizations of $|\mathcal{F}[n]|$, connecting, e. g. to the observations of Refs. Böttcher et al. (2019).

Finally, to quantify the density modulation and the phase coherence, we fit a three-Gaussian function to both $n_{\mathcal{M}}(y)$ and $n_{\Phi}(y)$ and extract the amplitudes of the finite-spatial-frequency peaks, $A_{\mathcal{M}}$ and $A_{\Phi}$, for both distributions, respectively. Note that for a BEC, which is a phase coherent state, $A_{\Phi}$ will be zero since it probes only finite-spatial-frequency peaks; see Fig. 2 (g-i,m-o).

V Characterization of the Supersolid State

We are now in the position to study two key aspects, namely (i) the evolution of the density modulation and phase coherence across the BEC-supersolid-ID phases, and (ii) the lifetime of the coherent density modulated state in the supersolid regime.

Evolution of the supersolid properties across the phase diagram – The first type of investigations is conducted with ${}^{166}{\rm Er}$ since for this species the scattering length and its dependence on the magnetic field has been precisely characterized Chomaz et al. (2016, 2018). After preparing the sample, we ramp $a_{\rm s}$ to the desired value and study the density patterns as well as their phase coherence by probing the amplitudes $A_{\mathcal{M}}$ and $A_{\Phi}$ as a function of $a_{\rm s}$ after $t_{\rm h}=5\,$ms. As shown in Fig. 3(a), in the BEC region (i. e. for large $a_{\rm s}$), we observe that both $A_{\mathcal{M}}$ and $A_{\Phi}$ are almost zero, evidencing the expected absence of a density modulation in the system. As soon as $a_{\rm s}$ reaches a critical value $a_{\rm s}^{*}$, the system’s behavior dramatically changes with a sharp and simultaneous increase of both $A_{\mathcal{M}}$ and $A_{\Phi}$. While the strength of $A_{\mathcal{M}}$ and $A_{\Phi}$ varies with decreasing $a_{\rm s}$ – first increasing then decreasing – we observe that their ratio $A_{\Phi}/A_{\mathcal{M}}$ remains constant and close to unity over a narrow $a_{\rm s}$-range below $a_{\rm s}^{*}$ of $\gtrsim 1\,a_{0}$ width; see inset. This behavior pinpoints the coexistence in the system of phase coherence and density modulation, as predicted to occur in the supersolid regime. For $(a_{\rm s}-a_{\rm s}^{*})<-1\,a_{0}$, we observe that the two amplitudes depart from each other. Here, while the density modulation still survives with $A_{\mathcal{M}}$ saturating to a lower finite value, the global phase coherence is lost with $A_{{\Phi}}/A_{\mathcal{M}}<1$, as expected in the insulating droplet phase.

To get a deeper insight on how our observations compare to the phase-diagram predictions (see Fig. 1), we study the link strength $S$ as a function of $a_{\rm s}$; see Fig. 3(b). Since $S$ quantifies the density overlap between neighboring droplets and is related to the tunneling rate of atoms across the droplet array, it thus provides information on the ability of the system to establish or maintain a global phase coherence. In this plot, we set $S=0$ in the case where no modulation is found in the ground state. At the BEC-to-supersolid transition, i. e. at $a_{\rm s}=a_{\rm s}^{*}$, a density modulation abruptly appears in the system’s ground-state with $S$ taking a finite value. Here, $S$ is maximal, corresponding to a density modulation of minimal amplitude. Below the transition, we observe a progressive decrease of $S$ with lowering $a_{\rm s}$, pointing to the gradual reduction of the tunneling rate in the droplet arrays. Close to the transition, we estimate a large tunneling compared to all other relevant timescales. However, we expect this rate to become vanishingly small, on the sub-Hertz level sup , when decreasing $a_{\rm s}$ $1$–$2\,a_{0}$ below $a_{\rm s}^{*}$. Our observation also hints to the smooth character of the transition from a supersolid to an ID phase, suggesting a crossover-type behavior.

The general trend of $S$, including the extension in $a_{\rm s}$ where it takes non-vanishing values, is similar to the $a_{\rm s}$-behavior of $A_{\mathcal{M}}$ and $A_{\Phi}$ observed in the experiments  foo (c). We observe in the experiments that the $a_{\rm s}$ dependence at the BEC-to-supersolid transition appears sharper than at the supersolid-to-ID interface, potentially suggesting a different nature of the two transitions. However, more investigations are needed since atom losses, finite-temperature and finite-size effects can affect, and in particular smoothen, the observed behavior Fisher and Ferdinand (1967); Imry and Bergman (1971); Imry (1980). Moreover, dynamical effects, induced by e. g.  excitations created at the crossing of the phase transitions or atoms losses during the time evolution, can also play a substantial role in the experimental observations, complicating a direct comparison with the ground state calculations. The time dynamics as well as a different scheme to achieve a state with supersolid properties will be the focus of the remainder of the manuscript.

Lifetime of the supersolid properties – Having identified the $a_{\rm s}$ range in which our dipolar quantum gas exhibits supersolid properties, the next central question concerns the stability and lifetime of such a fascinating state. Recent experiments on ${}^{162}{\rm Dy}$ have shown the transient character of the supersolid properties, whose lifetime is limited by three-body losses Tanzi et al. (2018); Böttcher et al. (2019). In these experiments, the phase coherence is found to survive up to 20 ms after the density modulation has formed. This time corresponds to about half of the weak-trap period. Stability is a key issue in the supersolid regime, especially since the tuning of $a_{\rm s}$, used to enter this regime, has a twofold consequence on the inelastic loss rate. First, it gives rise to an increase in the peak density (see Fig. 1b-d) and, second, it may lead to an enhancement of the three-body loss coefficient.

We address this question by conducting comparative studies on ${}^{166}{\rm Er}$ and ${}^{164}{\rm Dy}$ gases. These two species allow us to tackle two substantially different scattering scenarios. Indeed, the background value of $a_{\rm s}$ for ${}^{166}{\rm Er}$ (as well as for ${}^{162}{\rm Dy}$) is larger than $a_{\rm dd}$. Thus, reaching the supersolid regime, which occurs at $a_{\rm dd}/a_{\rm s}\approx 1.2-1.4$ in our geometry, requires to tune ${B}$ close to the pole of a FR. This tuning also causes an increase of the three-body loss rate. In contrast, ${}^{164}{\rm Dy}$ realizes the opposite case with the background scattering length smaller than $a_{\rm dd}$. This feature brings the important advantage of requiring tuning $B$ away from the FR pole to reach the supersolid regime. As we will describe below, this important difference in scattering properties leads to a strikingly longer lifetime of the ${}^{164}{\rm Dy}$ supersolid properties with respect to 166Er and to the recently observed behavior in 162Dy  Tanzi et al. (2018); Böttcher et al. (2019).

The measurements proceed as follows. For both ${}^{166}{\rm Er}$ and ${}^{164}{\rm Dy}$, we first prepare the quantum gas in the stable BEC regime and then ramp $a_{\rm s}$ to a fixed value in the supersolid regime for which the system exhibits a state of coherent droplets (i. e. $A_{\Phi}/A_{\mathcal{M}}\approx 1$); see previous discussion. Finally, we record the TOF images after a variable $t_{\rm h}$ and we extract the time evolution of both $A_{\Phi}$ and $A_{\mathcal{M}}$. The study of these two amplitudes will allow us to answer the question whether the droplet structure – i. e. the density modulation in space – persists in time whereas the coherence among droplets is lost ($A_{\mathcal{M}}>A_{\Phi}\rightarrow 0$) or if the density structures themselves vanish in time ($A_{\mathcal{M}}\approx A_{\Phi}\rightarrow 0$).

As shown in Fig. 4, for both species, we observe that $A_{\Phi}$ and $A_{\mathcal{M}}$ decay almost synchronously with a remarkably longer lifetime for ${}^{164}{\rm Dy}$ (b) than ${}^{166}{\rm Er}$ (a). Interestingly, $A_{\Phi}$ and $A_{\mathcal{M}}$ remain approximately equal during the whole time dynamics; see inset (a-b). This behavior indicates that it is the strength of the density modulation itself and not the phase coherence among droplets that decays over time. Similar results have been found theoretically in Ref. foo (d). We connect this decay mainly to three-body losses, especially detrimental for ${}^{166}{\rm Er}$, and possible excitations created while crossing the BEC-to-supersolid phase transition sup . For comparison, the inset shows also the behavior in the ID regime for ${}^{166}{\rm Er}$, where $A_{\Phi}/A_{\mathcal{M}}<1$ already at short $t_{\rm h}$ and remains so during the time evolution sup .

To get a quantitative estimate of the survival time of the phase coherent and density modulated state, we fit a simple exponential function to $A_{\Phi}$ and extract $t_{\Phi}$, defined as the $1/10$ lifetime; see Fig. 4. For ${}^{166}{\rm Er}$, we extract $t_{\Phi}=38(6)\,$ms. For $t_{\rm h}>t_{\Phi}$, the interference patterns become undetectable in our experiment and we recover a signal similar to the one of a non-modulated BEC state (as in Fig. 2(a,d)). These results are consistent with recent observations of transient supersolid properties in ${}^{162}{\rm Dy}$ Tanzi et al. (2018). For ${}^{164}{\rm Dy}$, we observe that the coherent density-modulated state is remarkably long-lived. Here, we find $t_{\Phi}=152(13)\,$ms.

The striking difference in the lifetime and robustness of the supersolid properties between ${}^{166}{\rm Er}$ and ${}^{164}{\rm Dy}$ becomes even more visible when studying $t_{\Phi}$ as a function of $a_{\rm s}$ (${B}$ for Dy). As shown in Fig. 5, $t_{\Phi}$ for Er remains comparatively low in the investigated supersolid regime and slightly varies between 20 and 40 ms. Similarly to the recent studies with ${}^{162}{\rm Dy}$, this finding reveals the transient character of the state and opens the question of whether a stationary supersolid state can be reached with these species. On the contrary, for ${}^{164}{\rm Dy}$ we observe that $t_{\Phi}$ first increases with ${B}$ in the range from $1.8\,$G to about $1.98\,$G. Then, for ${B}>1.98\,$G, $t_{\Phi}$ acquires a remarkably large and almost constant value of about $150\,$ms over a wide ${B}$-range. This shows the long-lived character of the supersolid properties in our ${}^{164}{\rm Dy}$ quantum gas. We note that over the investigated range, $a_{\rm s}$ is expected to monotonously increase with ${B}$ sup . Such a large value of $t_{\Phi}$ exceeds not only the estimated tunneling time across neighbouring droplets but also the weak-axis trap period, which together set the typical timescale to achieve global equilibrium and to study collective excitations.

VI Creation of states with supersolid properties by evaporative cooling

The long-lived supersolid properties in ${}^{164}{\rm Dy}$ motivate us to explore an alternative route to cross the supersolid phase transition, namely by evaporative cooling instead of interaction tuning. For this set of experiments, we have modified the waists of our trapping beams in order to achieve quantum degeneracy in tighter traps with respect to the one used for condensation in the previous set of measurements. In this way, the interference peaks in the supersolid region are already visible without the need to apply a further compression of the trap since the side-to-central-peak distance in the momentum distribution scales roughly as $1/\ell_{z}$ Chomaz et al. (2018). Forced evaporative cooling is performed by reducing the power of the trapping beams piecewise-linearly in subsequent evaporation steps until a final trap with frequencies $2\pi\times(225,37,134)\,$Hz is achieved. During the whole evaporation process, which has an overall duration of about 3 s, the magnetic field is kept either at ${B}=2.43\,$G, where we observe long-lived interference patterns, or at ${B}=2.55\,$G, where we produce a stable non-modulated BEC. We note that these two $B$ values are very close without any FR lying in between sup .

Figure 6 shows the phase transition from a thermal cloud to a final state with supersolid properties by evaporative cooling. In particular, we study the phase transition by varying the duration of the last evaporation ramp, while maintaining the initial and final trap-beam power fixed. This procedure effectively changes the atom number and temperature in the final trap while keeping the trap parameters unchanged, which is important to not alter the final ground-state phase diagram of the system. At the end of the evaporation, we let the system equilibrate and thermalize for $t_{\rm h}=100\,$ms, after which we switch off the trap, let the atoms expand for $26.5\,$ms, and finally perform absorption imaging. We record the TOF images for different ramp durations, i. e. for different thermalization times. For a short ramp, too many atoms are lost such that the critical atom number for condensation is not reached, and the atomic distribution remains thermal; see Fig. 6(a).

By increasing the ramp time, the evaporative cooling becomes more efficient and we observe the appearance of a bimodal density profile with a narrow and dense peak at the center, which we identify as a regular BEC; see Fig. 6(b). By further cooling, the BEC fraction increases and the characteristic pattern of the supersolid state emerges; see Fig. 6(c-d). The observed evaporation process shows a strikingly different behavior in comparison with the corresponding situation at $B=2.55\,$G, where the usual thermal-to-BEC phase transition is observed; see Fig. 6(i-l).

We finally probe the lifetime of the supersolid properties by extracting the time evolution of both the amplitudes $A_{\Phi}$ and $A_{\mathcal{M}}$, as previously discussed. We use the same experimental sequence as the one in Fig. 6(d) – i. e. $300\,$ms duration of the last evaporation ramp and $100\,$ms of equilibration time – and subsequently hold the sample in the trap for a variable $t_{\rm h}$. As shown in Fig. 7(a), we observe a very long lifetime with both amplitudes staying large and almost constant over more than 200 ms. At longer holding time, we observe a slow decay of $A_{\Phi}$ and $A_{\mathcal{M}}$, following the one of the atom number. Moreover, during the dynamics, the ratio $A_{\Phi}/A_{\mathcal{M}}$ stays constant. The long lifetime of the phase-coherent density modulation is also directly visible in the persistence of the interference patterns in the averaged momentum density profiles (similar to Fig. 2(e)), both at intermediate and long times; see Fig.7(b) and (c), respectively. For even longer $t_{\rm h}$, we can not resolve anymore interference patterns in the TOF images. Here, we recover a signal consistent with a regular BEC of low $N$.

Achieving the coherent-droplet phase via evaporative cooling is a very powerful alternative path to supersolidity. We speculate that, for instance, excitations, which might be important when crossing the phase transitions by interaction tuning, may be small or removed by evaporation when reaching this state kinematically. Other interesting questions, open to future investigations, are the nature of the phase transition, the critical atom number, and the role of non-condensed atoms.

VII Conclusions

For both ${}^{166}{\rm Er}$ and ${}^{164}{\rm Dy}$ dipolar quantum gases, we have identified and studied states showing hallmarks of supersolidity, namely global phase coherence and spontaneous density modulations. These states exist in a narrow scattering-length region, lying between a regular BEC phase and a phase of an insulating droplet array. While for ${}^{166}{\rm Er}$, similarly to the recently-reported ${}^{162}{\rm Dy}$ case Tanzi et al. (2018); Böttcher et al. (2019), the observed supersolid properties fade out over a comparatively short time because of atom losses, we find that ${}^{164}{\rm Dy}$ exhibits remarkably long-lived supersolid properties. Moreover, we are able to directly create stationary states with supersolid properties by evaporative cooling, demonstrating a powerful alternative approach to interaction tuning on a BEC. This novel technique provides prospects of creating states with supersolid properties while avoiding additional excitations and dynamics. The ability to produce long-lived supersolid states paves the way for future investigations on quantum fluctuations and many-body correlations as well as of collective excitations in such an intriguing many-body quantum state. A central goal of these future investigations lies in proving the superfluid character of this phase, beyond its global phase coherence Boninsegni and Prokof’ev (2012); Scarola et al. (2006); Lu et al. (2015b); Cinti and Boninsegni (2017); Roccuzzo and Ancilotto (2018).

Note added– During the course of the present work, we became aware of a related work, reporting a theoretical study of the ground-state phase diagram based on Monte-Carlo calculations Kora and Boninsegni (2019).

VIII Acknowledgments

Acknowledgements.

We thank R. Bisset, B. Blakie, M. Boninsegni, G. Modugno, T. Pfau, and in particular L. Santos for many stimulating discussions. Part of the computational results presented have been achieved using the HPC infrastructure LEO of the University of Innsbruck. We acknowledge support by the Austrian Science Fund FWF through the DFG/FWF Forschergruppe (FOR 2247/PI2790), by the ERC Consolidator Grant (RARE, no. 681432), and by a NFRI Grant (MIRARE, no. ÖAW0600) from the Austrian Academy of Science. G. D. and M. S. acknowledge support by the Austrian Science Fund FWF within the DK-ALM (no. W1259-N27).

• Correspondence and requests for materials should be addressed to Francesca.Ferlaino@uibk.ac.at.

Appendix A Supplemental Material

Ground state calculations

We perform numerical calculations of the ground state following the procedure detailed in the supplementary information of Ref. Chomaz et al. (2018). The calculations are based on the conjugate-gradients technique to minimize the energy functional of an eGPE Ronen et al. (2006). In particular, the eGPE accounts for the effect of quantum fluctuations, by including the LHY term $\Delta\mu[n]=32g(na_{\rm s})^{3/2}(1+3\epsilon_{\rm dd}^{2}/2)/3\sqrt{\pi}$ in the system’s Hamiltonian (here $g=4\pi\hbar^{2}a_{\rm s}/m$ and $n=|\psi|^{2}$ is the spatial density of the macroscopic state $\psi$). $\Delta\mu[n]$ has been obtained under a local density approximation in Refs.  Lima and Pelster (2011, 2012). The relevance of the LHY correction has been demonstrated in various studies of dipolar Bose gases close to the mean-field instability Wächtler and Santos (2016a, b); Bisset et al. (2016); Chomaz et al. (2016); Schmitt et al. (2016); Chomaz et al. (2018) as it brings an additional repulsive potential, stabilizing the gas against mean-field collapse at large density. We note that the exact functional form of the potential, originating from beyond mean-field effects, has been questioned by several experimental results in finite-size trapped systems Schmitt et al. (2016); Chomaz et al. (2018); Ferrier-Barbut et al. (2018b); Cabrera et al. (2018), calling for further theory developments Cikojević et al. (2018).

Our numerical calculations provide us with the three-dimensional ground-state wavefunctions $\psi(\boldsymbol{r})$. From this, we compute the axial in-situ density profile along the trap’s weak axis, $n(y)=\int|\psi(\boldsymbol{r})|^{2}dxdz$ and find density profiles, corresponding to the BEC, the supersolid or the ID phase, that we plot in Fig. 1. From the density profiles that exhibit a density modulation, we evaluate $S$ by performing Gaussian fits to each droplet, i. e. to $n(y)$ with $y$ ranging between two neighboring local density minima. From these Gaussian fits, we evaluate the sets of centers $\{y^{(0)}_{i}\}_{i}$ and widths $\{\sigma_{i}\}_{i}$ corresponding to the macroscopic Gaussian wavefunctions $\{\psi_{i}\}_{i}$ associated to the individual droplets in the array. We then approximate the droplet wavefunction via $\psi_{i}(y)\approx\sqrt{n(y\approx y^{(0)}_{i})}=\alpha_{i}\exp\left(-(y-y^{(0)}_{i})^{2}/2\sigma_{i}^{2}\right)$ with $\alpha_{i}$ a normalization coefficient such that $\int|\psi_{i}(y)|^{2}dy=1$. We then evaluate the wavefunction overlap $S_{i}$ between the neighboring droplets $i-1$ and $i$ via:

[TABLE]

The latter equation is obtained via an analytical evaluation of the Gaussian integral. The characteristic link strength defined in the paper is then computed by averaging $S_{i}$ over all droplet links in the array: $S=\langle S_{i}\rangle_{i}$. In our calculation, we only consider as droplets all density peaks of at least 10 % of the global density maximum.

Link Strength and estimate of tunneling rate

Generally speaking, the wavefunction overlap between neighboring droplets relates to a tunneling term, which sets a particle exchange term between two neighboring droplets Josephson (1962); Barone and Paternò (1982); Javanainen (1986); Raghavan et al. (1999). Following the work of Ref. Wenzel et al. (2018), we perform a first estimate of the tunneling coefficient by simply considering the single-particle part of the Hamiltonian and evaluate it between two neighboring droplets. We note that, in our particular setting where the density modulation is not externally imposed but arises from the mere interparticle interactions, the inter-droplet interaction may also play a crucial role. To perform a more precise estimation of the tunneling between droplets, one would certainly need to properly account for this effect. Here, we stress that our approach simply gives a rough idea of the magnitude of tunneling while it does not aim to be a quantitative description of it. This consideration calls for further studies making a systematic analysis of the full Hamiltonian and of the full phase diagram within the Josephson junction formalism and beyond.

Generalizing the description of Ref. Wenzel et al. (2018) to neighboring droplets of different sizes and amplitudes, which are described by a three-dimensional wavefunction $\psi_{i}(\boldsymbol{r})$ approximated to a three-dimensional Gaussian of widths $(\sigma_{i,x},\sigma_{i,y},\sigma_{i,z})$ with $\sigma_{i,y}=\sigma_{i}$, our estimate writes:

[TABLE]

where $\ell_{x,y,z}=\sqrt{\hbar/m\omega_{x,y,z}}$ are the harmonic oscillator lengths.

In general, the tunnelling coefficients set two typical rates relevant for equilibration processes. The first one is the bare single-particle tunneling rate, which is equal to $J_{i}/h$, while the second accounts for the bosonic enhancement from the occupation of the droplet modes and writes $\tilde{t}_{i}=\sqrt{N_{i}N_{i-1}}|J_{i}|/h$ where $N_{i}$ is the number of atoms in droplet $i$. In our analysis, we then define the average rates over the droplet arrays as characteristic rates $J/h=\langle J_{i}\rangle_{i}/h$, and $\tilde{t}=\langle\tilde{t}_{i}\rangle_{i}$; see e.g. Hadzibabic et al. (2004). While the ground state evolves from a BEC to a supersolid to an ID, the relevant timescale for achieving (global) equilibrium crosses from being set by the trap frequencies to the above-mentioned tunneling rates.

Using our approximate model, we here give a first estimate of the rates $J/h$ and $\tilde{t}$ as a function of $a_{\rm s}$, for the parameters of Fig. 1(b-d) of the main text (i.e. Er quantum gas with $N=5\times 10^{4}$ atoms). Here we find that, for $a_{\rm s}=a_{\rm s}^{*}$, $J/h\sim 400\,$Hz and $\tilde{t}\sim 10\,$MHz while for $a_{\rm s}=a_{\rm s}^{*}-2.5\,a_{0}$, $J/h\sim 10^{-7}\,$Hz and $\tilde{t}\sim 10^{-3}\,$Hz.

Toy model for the interference pattern

As described in the main text we use a simple toy model, adapted from Ref. Hadzibabic et al. (2004), to identify the main features of the TOF interference patterns obtained from an insitu density-modulated state. As a quick reminder, our model considers a one-dimensional array of $N_{D}$ Gaussian droplets, described by a single classical field, $\psi_{i}$, thus neglecting quantum and thermal fluctuations. We compute the TOF density distribution from the free-expansion of the individual $\psi_{i}$ during a time $t$ via $n(y,t)=|\sum_{i}\psi_{i}(y,t)|^{2}$. In our calculations, we also account for the finite imaging resolution by convolving the resulting $n(y,t)$ with a gaussian function of width $\sigma_{\rm im}$. Here we allow the characteristics of the individual $\psi_{i}$ to fluctuate. In this aim, we introduce noise on the corresponding parameter with a normal distribution around its expectation value and with a variable standard deviation (only $\phi_{i}$ can also have a uniform distribution). We then perform a Monte-Carlo study and perform ensemble averages, similar to our experimental analysis procedure. We note that, in this simple implementation, the noise on the different parameters – droplet amplitudes, widths and distances – are uncorrelated.

In the main text, we present results for a single set of parameters, namely $N_{D}=4$, $d\equiv\langle d_{i}\rangle_{i}=2.8\,\mu{\rm m}$ (mean droplet distance), $\sigma_{y}\equiv\langle\sigma_{i}\rangle_{i}=0.56\,\mu{\rm m}$ (mean droplet size), $t=30\,$ms, and $\sigma_{\rm im}=3\,\mu{\rm m}$, typical for our experimental Er setting and the corresponding theory expectations in the supersolid regime. $\langle\cdot\rangle_{i}$ denotes the average over the droplets. In this section, we have a deeper look at the impact of the different parameters on both the TOF signal and our FT analysis. We study both the fully phase coherent and fully incoherent case, and the unchanged parameters are set as in Fig. 2(j,m) and (l,o).

In Fig. AS1, we first exemplify the TOF and FT profiles for a varying number of droplets, between 2 and 8, which cover the range of relevant $N_{D}$ over the phase diagram of Fig. 1. The results remain remarkably similar to the realization of Fig. 2 with only slight quantitative changes. The main difference lies in the individual interference patterns obtained in the phase incoherent case. With increasing $N_{D}$, those profiles become more complex and made of a larger number of peaks (see (c,g,k)). Yet, in this incoherent case, a similar (non-modulated) profile is recovered in the averaged $n(k_{y})$ for all $N_{D}$. Additionally, we note that for the coherent case with $N_{D}=8$, the side peaks in the FT analysis (see (j)) become less visible. By performing additional tests, we attribute this behavior to the limited TOF duration, $t$, used in our experiment yielding a typical length scale, $\sqrt{\hbar t/m}$ ($=3.4\mu{\rm m}$), which becomes small compared to the system size ($\approx(N_{D}-1)d+\sigma_{y}$) for large $N_{D}$. This intermediate regime in the TOF expansion leads to more complex features, including smaller-sized motifs, in the interference patterns. Finally, when accounting for our imaging resolution, it yields a broadening of the structure observed in the TOF images and less visible peaks in the FT (see (i)). We note that our experiments, because of limited $N$ and additional losses, should rather lie in the regime $2\leq N_{D}\leq 5$; see Fig. 1(b).

We then investigate the evolution of the interference patterns and their FT analysis for a varying mean droplet size, $\sigma_{y}$, while keeping their mean distance, $d$, fixed. This study is particularly relevant recalling that, within the Josephson junction formalism (see main text and corresponding section of this Supplemental Material), the key parameter controlling the tunneling rate between the droplets is set by the ratio $\sigma_{y}/d$, and the link strength parameter that we use to characterize the supersolid regime scales roughly as $\exp(-(d/2\sigma_{y})^{2})$. Thus, in our experiment, $\sigma_{y}/d$ is intrinsically expected to decrease with the scattering length (see Fig. 3). Performing a direct estimate of the average droplet link from the initial state of our toy model, we find $S=0.004$ for the calculations of Fig. 2(j-o), lying in an expected supersolid regime yet rather close to the supersolid-to-ID transition. Figure AS2 investigates the effect of smaller and larger values of $\sigma_{y}/d$ (and consequently of $S$) on the TOF and FT profiles while independently assuming phase coherence or incoherence. Qualitatively, the features remain similar as in Fig. 2(j-o). In the coherent case, side peaks are visible in the individual as well as in the mean $n(k_{y})$ (see (a,e,i)) and yield side peaks in the FT profiles, with $n_{\mathcal{M}}\approx n_{\mathcal{\Phi}}$ (see (b,f,j)). Increasing (decreasing) $\sigma/d$ mainly results in a stronger (weaker) signal both in the interference pattern and their FT analysis. Within our toy model, we find that, already for $\sigma/d=0.25$, the signal nearly vanishes; see (i,j). Even if, given the approximations used in our toy model, this exact value may not fully hold for our experimental conditions, we expect a similar trend. It is interesting to keep in mind that this effect may limit our capacity of detecting an underlying supersolid state via matter-wave interference in experiments. In the incoherent case, the effect of decreasing $\sigma_{y}/d$ mainly results in a broader shape of the mean density profile, while it remains non-modulated; see (c,g,k). In the FT analysis $n_{\Phi}$ remains structure-less independently of $\sigma_{y}/d$ while the structures in $n_{\mathcal{M}}$ becomes sharper with decreasing $\sigma_{y}/d$, as in the coherent case; see (d,h,l).

Finally, we investigate how a possible shot-to-shot noise on the position of the central interference peak could affect our observables of the density modulation and phase coherence. In the experiments, such fluctuations may occur, for instance, because of beam-pointing fluctuations or excitations of the gas. Although we compensate for such effects by recentering the individual images (see Imaging Analysis section), residual effects may remain, in particular due to center misestimation in the mere presence of the interference patterns of interest. To investigate this aspect, we repeat our toy model calculations now including noise in the global droplet array position and using a standard deviation of $2\,\mu{\rm m}$ for two values of $\sigma_{y}/d$; see Fig.AS3. Again, qualitatively the observed features remains similar to our prediction in the main text. The main effect lies in the appearance of a small discrepancy in the coherent case between $n_{\Phi}$ and $n_{\mathcal{M}}$, while the structure in the incoherent case remains similar. As the center misestimation should be the most severe in the latter case (due to the variability of the interference patterns observed here), our test shows the robustness of our analysis procedure against this issue.

Imaging Analysis: ${}^{164}{\rm Dy}$ and ${}^{166}{\rm Er}$

The density distributions in momentum space are extracted from the TOF images using the free-expansion expectation. In the Dy case, the thermal component is subtracted from the individual distribution by cutting out the central region of the cloud and performing an isotropic Gaussian fit on the outer region. This subtraction is beneficial because of the large thermal fraction. In the ${}^{166}{\rm Er}$ case, such a subtraction is on the contrary complicated because of the weak thermal component and this pre-treatment may lead to improper estimation of $A_{\mathcal{M}}$ and $A_{\Phi}$ in the later analysis. The obtained momentum density distributions are then recentered and integrated numerically along $k_{z}$($k_{x}$) between $[-2.0,+2.0]\,\mu$m*-1* ($[-1.28,+1.28]\,\mu$m*-1*) to obtain $n(k_{Y})$ ($n(k_{y})$) for ${}^{164}{\rm Dy}$ (${}^{166}{\rm Er}$). The recentering procedure uses the result a single Gauss fit to the TOF images. The fit is performed after convoluting each image with a Gaussian function of width $0.5\,\mu{\rm m}$ whose purpose is to reduce the impact of the interference pattern on the center estimation foo (e).

In order to characterise the system’s state, we use the Fourier transform, $\mathcal{F}[n](y)$ of the single density profile, $n(k_{y})$. We then compute two average profiles, $n_{\mathcal{M}}$ and $n_{\Phi}$, relying on ensemble average over all measurements under the same experimental conditions; see below for a detailed discussion on $n_{\mathcal{M}}$ and $n_{\Phi}$. In all the measurements reported in this work we use averages over typically 15 to 100 realizations.

To quantify both the existence of a density modulation and global phase coherence on top of this modulation, we fit both $n_{\mathcal{M}}(y)$ and $n_{\Phi}(y)$ with a triple-Gaussian function, where one Gaussian accounts for the central peak and the other Gaussians are accounting for the symmetric side peaks. The amplitudes of the latter give $A_{\mathcal{M}}$ and $A_{\Phi}$, respectively. The distance between the side peaks and the central one is allowed to vary between $[2.5,2.7]\,\mu$m ($[2.3,2.5]\,\mu$m) in the case of ${}^{164}{\rm Dy}$ (${}^{166}{\rm Er}$).

Details on the Fourier analysis

In our analysis we rely on two averaged profiles, named $n_{\mathcal{M}}$ or $n_{\Phi}$, to quantify both the density modulation and its phase coherence. Here we detail the meaning of the average performed.

The Fourier transform (FT) of the integrated momentum distributions, $n(k_{y})$, which reads $\mathcal{F}[n](y)=\left|\mathcal{F}[n](y)\right|\,\exp(i\arg\left(\mathcal{F}[n](y)\right))$ sets the ground for our analysis. As stated in the main text, an in-situ density modulation of wavelength $y^{*}$ yields patterns in $n(k_{y})$ and consequently induce peaks at $y\approx y^{*}$, in the FT norm, $\left|\mathcal{F}[n](y)\right|$, see Fig. 2(g-i) and (m-o). Spatial variations of the phase relation within the above-mentioned density modulation translate into phase shifts of the interference patterns, which are stored in the FT argument at $y\approx y^{*}$, $\arg\left(\mathcal{F}[n](y^{*})\right)$; see also Ref. Hadzibabic et al. (2004); Hofferberth et al. (2007).

The first average that we use is $n_{\mathcal{M}}(y)=\langle|\mathcal{F}[n](y)|\rangle$, i. e. the average of the FT norm of the individual images. As the phase information contained in $\arg\left(\mathcal{F}[n](y)\right)$ is discarded from $n_{\mathcal{M}}$ when taking the norm, the peaks in $n_{\mathcal{M}}$ probe the mere existence of an insitu density modulation of roughly constant spacing within the different realizations. The second average of interest is $n_{\Phi}(y)=|\langle\mathcal{F}[n](y)\rangle|$, i. e. the average of the full FT of the individual images. In contrast to $n_{\mathcal{M}}$, $n_{\Phi}$ keeps the phase information of the individual realizations contained in $\arg\left(\mathcal{F}[n](y^{*})\right)$. Consequently, peaks in $n_{\Phi}$ indicate that the phase relation is maintained over the density modulation, in a similar way for all realizations. Their presence thus provides information on the global phase coherence of a density-modulated state.

Experimental sequence: ${}^{164}{\rm Dy}$ and ${}^{166}{\rm Er}$

166**Erbium - The BEC of ${}^{166}{\rm Er}$ is prepared similarly to Refs. Aikawa et al. (2012); Chomaz et al. (2016, 2018); Petter et al. (2018). We start from a magneto-optical trap with $2.4\times 10^{7}$ 166Er  atoms at a temperature of $10\mu{\rm K}$, spin-polarized in the lowest Zeeman sub-level. In a next step we load about $3\times 10^{6}$ atoms into a crossed optical dipole trap (ODT) operated at 1064 nm. We evaporatively cool the atomic cloud by reducing the power and then increasing the ellipticity of one of the ODT beams. During the whole evaporation a constant magnetic field of $B=1.9\,$G ($a_{\rm s}=80\,a_{0}$) along $z$ is applied. We typically achieve BEC with $1.4\times 10^{5}$ atoms and a condensed fraction of $70\%$. In a next step the ODT is reshaped in 300 ms into the final trapping frequencies $\omega_{x,y,z}=2\pi\times(145,31.5,151)\,$Hz. Consecutively, we ramp $B$ linearly to $0.62\,$G ($64.5\,a_{0}$) in 50 ms and obtain a BEC with $8.5\times 10^{4}$ atoms, which are surrounded by $3.5\times 10^{4}$ thermal atoms. This point marks the start of the ramp to the final $a_{s}$.

164**Dysprosium - For the production of a ${}^{164}\text{Dy}$ BEC we closely follow the scheme presented in Trautmann et al. (2018). Starting from a $3\,$s loading phase of our 5-beam MOT in open-top configuration Ilzhöfer et al. (2018), we overlap a $1064\,$nm single-beam dipole trap with a $\nicefrac{{1}}{{e^{2}}}$-waist of about $22\,\mu$m, for $120\,$ms. Eventually, we transfer typically $8\times 10^{6}$ atoms utilizing a time averaging potential technique to increase the spatial overlap with the MOT. After an initial $1.1\,$s evaporative cooling phase by lowering the power of the beam, we add a vertically propagating beam, derived from the same laser, with a $\nicefrac{{1}}{{e^{2}}}$-waist of about $130\,\mu$m to form a crossed optical dipole trap for additional confinement. Subsequently, we proceed forced evaporative cooling to reach quantum degeneracy by nearly exponentially decreasing the laser powers in the two dipole-trap beams over $3.6\,$s. We achieve BECs of ${}^{164}\text{Dy}$ with typically $10^{5}$ atoms and condensate fractions of about $40\%$. During the entire evaporation sequence the magnetic field is kept constant at $2.5\,$G pointing along the vertical ($z$-) axis.

To be able to condense directly into the supersolid, we modify the dipole trap to condense at a stronger confinement of $\omega_{x,y,z}=2\pi\times\left(225,37,134\right)\,$Hz. After a total evaporative cooling duration of $3.1\,$s, we achieve Bose-Einstein condensation at $2.55\,$G and reach a state with supersolid properties at $2.43\,$G, keeping the magnetic field constant throughout the entire evaporation sequence for both cases.

Time of flight and imaging for ${}^{166}{\rm Er}$ and ${}^{164}{\rm Dy}$ - In order to probe the momentum distribution of the Dy (Er) gases, we switch off the confining laser beams and let the atoms expand freely for $18\,$ms (15 ms), while keeping the magnetic field constant. Consecutively the amplitude of $B$ is increased to a fixed amplitude of 5.4 G (0.6 G). In the case of ${}^{164}{\rm Dy}$, the magnetic field orientation is rotated in order to point along the imaging axis. This ensures constant imaging conditions for different $a_{\rm s}$. After an additional $9\,$ms (15 ms) we perform a standard absorption imaging.

Tuning the scattering length in ${}^{166}{\rm Er}$ and ${}^{164}{\rm Dy}$

166**Erbium - All measurements start with a BEC at $64.5\,a_{0}$. In order to probe the BEC-supersolid-ID region, we linearly ramp $a_{s}$ to its target value in $t_{\rm r}=20\,$ms by performing a corresponding ramp in $B$. Due to a finite time delay of the magnetic field in our experimental setup and the highly precise values of $a_{s}$ needed for the experiment, we let the magnetic field stabilize for another 15 ms before $t_{h}=0\,$ starts. By this, we ensure that the residual lowering of $a_{s}$ during the entire hold time is $\lesssim 0.3\,a_{0}$. In the main text, we always give the $a_{s}$ at $t_{\rm h}=0\,$. Furthermore, we estimate our magnetic field uncertainty to be $\pm 2.5\,$mG, leading to a $\pm 0.2\,a_{0}$ uncertainty of $a_{s}$ in our experiments.

To choose the best ramping scheme, we have performed experiments varying $t_{\rm r}$ from 0.5 ms to 60 ms, ramping to a fixed $a_{\rm s}$ lying in the supersolid regime, and holding for $t_{\rm h}=5\,$ms after a fixed 15 ms waiting time. We record the evolution of $A_{\Phi}$ as a function of $t_{\rm r}$; see Fig. AS4. When increasing $t_{\rm r}$, we first observe that $A_{\Phi}$ increases, up to $t_{\rm r}=20$ ms, and then $A_{\Phi}$ gradually decreases. The initial increase can be due to diabatic effects and larger excitation when fast-crossing the phase transition. On the other hand, the slow decrease at longer $t_{\rm r}$ can be explained by larger atom loss during the ramp. We then choose $t_{\rm r}=20\,$ms as an optimum value where a supersolid behavior develops and maintains itself over a significant time while the losses are minimal.

164**Dysprosium - As the value of the background scattering, $a_{\rm bg}$ length for 164Dy is still under debate Tang et al. (2015); Schmitt et al. (2016); Ferrier-Barbut et al. (2018b), we discuss the experimental settings in terms of magnetic field. Yet, to gain a better understanding of the tunability of $a_{\rm s}$ in our experiment, we first perform a Feshbach spectroscopy scan on a BEC at $T=60\,$nK. After evaporative cooling at $B=2.5\,$G, we jump to $B$ varying from $1\,$G to $7.5\,$G and we hold the sample for $100\,$ms. Finally, we switch off the trap, let the cloud expand for $26\,ms$ and record the total atom number as a function of $B$. We then fit the observed loss features with a gaussian fit to obtain the position $B_{0,i}$ and width $\Delta B_{i}$ of the FRs, numbered ${i}$. We finally use the standard Feshbach resonance formula to estimate the $a_{\rm s}$-to-$B$ dependence via $a_{\rm s}(B)=a_{\rm bg}\prod_{i}\left(1-\Delta B_{i}/(B-B_{0,i})\right)$. Here we account for 8 FRs located between 1.2 G and 7.2 G. Depending on the background scattering length $a_{\rm bg}$, the overall magnitude of $a_{\rm s}(B)$ changes. We can get an estimate of $a_{\rm bg}$ from literature. In Fig. AS5, we use the value of $a_{\rm s}$ from Ref. Tang et al. (2015) obtained at $1.58$ G close to the $B$-region investigated in our experiment, $a_{\rm s}=92(8)\,a_{0}$. By reverting the $a_{\rm s}(B)$ formula, we set $a_{\rm bg}=87(8)\,a_{0}$. For the measurements of Figs. 4-5, we ramp $B$ linearly from 2.5 G in 20 ms to a final value ranging from 1.8 to 2.1 G, for which we estimate $a_{\rm s}$ ranging from $97(9)\,a_{0}$ to $105(10)\,a_{0}$. We calibrate our magnetic field using RF spectroscopy, with a stability of about 2 mG. In the Dy case, we do not apply an additional stabilization time. This is justified because of the more mellow $a_{\rm s}$-to-$B$ dependence in the $B$-range of interest as well as of the wider $a_{\rm s}$-range of the superoslid regime (see Fig. 1) compared to the Er case. For the measurements of Figs. 6–7, we use two $B$-values, namely 2.43 G and 2.55 G, at which we perform the evaporative cooling scheme. Here we estimate $a_{\rm s}=109(10)\,a_{0}$ and $a_{\rm s}=134(12)\,a_{0}$, respectively.

Atom losses in ${}^{166}{\rm Er}$ and ${}^{164}{\rm Dy}$

As pointed out in the main text, in the time evolution of the quantum gases in both the supersolid and the ID regime, inelastic atom losses play a crucial role. The atom losses are increased in the above mentioned regime as (i) higher densities are required so that a stabilization under quantum fluctuation effects becomes relevant and (ii) the magnetic field may need to be tune close to a FR pole to access the relevant regime of interaction parameters. (i) is at play for all magnetic species but more significant for 166Er due to the smaller value of $a_{\rm dd}$. (ii) is relevant for both ${}^{166}{\rm Er}$ and ${}^{162}{\rm Dy}$ but conveniently avoided for ${}^{164}{\rm Dy}$ thanks to the special short-range properties of this isotope.

To quantify the role of these losses, we report here the evolution of the number of condensed atoms, $N$, as a function of the hold time in parallel to the phase coherent character of the density modulation observed. We count $N$ by fitting the thermal fraction of each individual image with a two-dimensional Gaussian function. To ensure that only the thermal atoms are fitted, we mask out the central region of the cloud associated with the quantum gas. Afterwards we subtract this fit from the image and perform a numerical integration of the resulting image (so called pixel count) to obtain $N$.

166**Erbium - In the Er case, a 15 ms stabilization time is added to ensure that $a_{\rm s}$ is reached up to $0.3\,a_{0}$. During this time, i. e. for $t_{\rm h}<0$, we suspect that the time-evolution of the cloud properties is mainly dictated by the mere evolution of the scattering length. Therefore, in the main text, we report on the time evolution for $t_{\rm h}\geq 0$. We note that because of the narrow $a_{\rm s}$-range for the supersolid regime, the long stabilization time for $a_{\rm s}$ is crucial. However, because of the significant role of the atom losses in our system, in particular for ${}^{166}{\rm Er}$, the early evolution of $N$ and the cloud’s properties are intimately connected. Therefore, the early time evolution at $t_{h}<0\,$ is certainly of high importance for our observations at $t_{h}\geq\,0$.

To fully report on this behavior, we show the evolution of $N$ and $A_{\Phi}$ during both the stabilization and the holding time in Fig. AS6 for three different $a_{\rm s}$ values – either in the ID (a, b) or supersolid regime (c-f). The time evolution shows significant atom loss, prominent already during the stabilization time, and levels off towards a remaining atom number at longer holding times in which we recover small BECs. Simultaneously, in each case reported here, we observe that during the stabilization time $A_{\Phi}$ increases and a coherent density modulated state grows. This density modulation starts to appear at a typical atom number of $N\gtrsim 6\times 10^{4}$ and consecutively decays. For the lower $a_{\rm s}=53.3(2)\,a_{0}$ case, we observe that the coherent state does not survive the $a_{\rm s}$ stabilization time, and decays faster than the atoms loss; see Fig. AS6 (a, b). This behavior corresponds to the ID case discussed in the main text. The central point of the present work is to identify a parameter range where the coherence of the density modulated state survives for $t_{\rm h}>0\,$ and its decay time scale is similar to the one of the atom loss. In order to quantify a timescale for the atom number decay, we fit an exponential decay to $t_{\rm h}\geq 0\,$ms. Here we allow an offset $N_{r}$ of the fit, accounting for the BEC recovered at long holding times. In Table 1, we report on the typical 1/10-decay times of the atom number, which are up to 50 ms. These values are of the order as the extracted $t_{\Phi}$, see Table 1 and Fig. 5 of the main text. This reveals that in ${}^{166}{\rm Er}$ the extracted lifetime of the coherent density modulated states are mainly limited by atom loss.

Furthermore we note that the extracted $N_{r}$ values for the recovered BECs are smaller than $2\times 10^{4}$, which is consistent with the BEC region found in the phase diagram of Fig. 1(f).

164**Dysprosium - Differently from the 166Er case, for ${}^{164}{\rm Dy}$, we operate in a magnetic-field range in which the three-body collision coefficients are small and only moderate atom losses occur. This enables the observation of an unprecendented long-lived supersolid behavior. To understand the effects limiting the supersolid lifetime, we study the lifetime of the condensed-atom number for different $B$. We perform this detailed study for the data of Fig. 5 of the main text, which are obtained after preparing a stable BEC and then ramping $B$ to the target value. Fig. AS7 shows the parallel evolution of $N$ and $A_{\Phi}$ for three different magnetic field values $1.8\,$G, $2.04\,$G and $2.1\,$G. Here we observe that, for all $B$ values, $A_{\Phi}$ seems to decay faster than the atom number. This suggests that the lifetime of the density-modulated state in our ${}^{164}{\rm Dy}$ experiment is not limited by atom losses. To confirm this observation, we extract the 1/10 lifetimes of both $N$ and $A_{\Phi}$; see Table 2. The values confirm our observation and shows an atom number lifetime larger than $t_{\Phi}$ at least by a factor of $\approx 5$. In addition, we find that the ratio $t_{N}/t_{\Phi}$ varies, indicating that atom losses are not the only mechanism limiting the lifetime of the supersolid properties in Dy.