Transient supersolid properties in an array of dipolar quantum droplets

TL;DR

This paper explores the emergence of supersolid properties in dipolar Bose-Einstein condensates, demonstrating phase coherence and density modulation through theoretical analysis and experimental verification.

Contribution

It provides a combined theoretical and experimental study revealing the phase diagram and supersolid behavior in dipolar quantum droplets.

Findings

Identification of three distinct phases: BEC, incoherent droplets, and coherent droplets.

Experimental observation of in-situ density modulation.

Confirmation of phase coherence via matter wave interference.

Abstract

We study theoretically and experimentally the emergence of supersolid properties in a dipolar Bose-Einstein condensate. The theory reveals a ground state phase diagram with three distinct regimes - a regular Bose-Einstein condensate, incoherent and coherent arrays of quantum droplets. In the latter the droplets are connected by a finite superfluid density, which leads - in addition to the periodic density modulation - to a robust phase coherence throughout the whole system. We further theoretically demonstrate that we are able to dynamically approach the ground state in our experiment and that its lifetime is only limited by three-body losses. Experimentally we probe and confirm the signatures of the phase diagram by observing the in-situ density modulation as well as the phase coherence using matter wave interference.

Figures

Equations

$$SW=\sum _{|k|=0.7\text{upmu}{\text{m}}^{-1}}^{2.1\text{upmu}{\text{m}}^{-1}}S(k)/\sum _{|k|<0.7\text{upmu}{\text{m}}^{-1}}S(k)$$

$$I(z)\propto e^{-\frac{2z^2}{Z_0^2}}\left(1+A_1\cos \left({\Phi }_1+\frac{2\pi z}{D}\right)\right).$$

$$I(z)\propto e^{-\frac{2z^2}{Z_0^2}}\left(1+\sum _{n=1}^{N-1}{A}_n\cos \left({\Phi }_n+\frac{2\pi nz}{D}\right)\right),$$

$$\frac{2}{N}\sum _{j=n+1}^N{e}^{i({\varphi }_j-{\varphi }_{j-n})}.$$

$$I(z)\propto \frac{1}{N}\frac{{\sin }^2(N\pi z/D)}{{\sin }^2(\pi z/D)}.$$

$$\frac{1}{(N-1){⟨p⟩}^2}\sum _{i=1}^N{\left|p_i-⟨p⟩\right|}^2$$

Full text

Whether a material is solid, liquid or gaseous in classical physics depends on the strength of the interactions with respect to the motional energy of the particles. In analogy to this behavior the interplay between quantum fluctuations and interparticle interactions also leads to new phases of matter in quantum mechanics. One example of such quantum phases is the supersolid Boninsegni and Prokof’ev (2012); Leggett (1970); Andreev and Lifshitz (1969); Thouless (1969); Chester (1970); Penrose and Onsager (1956) featuring the periodic density modulation of a solid together with the dissipationless flow of a superfluid. While these properties are normally thought of as mutually exclusive, it was shown that they can actually coexist Leggett (1970). More formally speaking a supersolid features both on- and off-diagonal long-range order in its density matrix Boninsegni and Prokof’ev (2012). Originally, supersolidity was mainly discussed in the context of 4He, for which it remains elusive in the experiments Chan et al. (2013). The concept of supersolidity has since been generalized to other superfluid systems and supersolid properties have been observed in ultracold atomic systems for spin-orbit coupled Bose-Einstein condensates (BECs) Li et al. (2017) as well as BECs symmetrically coupled to two crossed optical cavities Léonard et al. (2017a, b). In these systems the periodicity of the modulation is induced by the underlying periodic optical potentials. In contrast, there are physical systems where the self-organized structure formation is induced by the intrinsic interactions and therefore phonon modes of the periodic modulation are allowed like in classical solids.

One promising system of this type are dipolar quantum gases Baranov et al. (2012); Lu et al. (2015), featuring both short-range contact interactions as well as long-range dipole-dipole interactions. These dipolar systems feature a rotonic dispersion relation similar to 4He Santos et al. (2003) which in addition, is fully tunable by changing the contact interaction strength as well as the external confinement along the dipoles. This dispersion relation has been studied experimentally Chomaz et al. (2018); Petter et al. (2018) and has led to the discovery of 2D arrays of quantum droplets Kadau et al. (2016); Ferrier-Barbut et al. (2016a); Petrov (2015). However, it was shown that the 2D arrays in these early experiments were incoherent, excited states of the system. For the considered geometries the ground state was always a single droplet Wächtler and Santos (2016); Bisset et al. (2016); Chomaz et al. (2016); Ferrier-Barbut et al. (2018). In contrast to this, in Wenzel et al. (2017) we pointed out that the ground state in strongly-confined 2D geometries is made up of droplet arrays, but experimentally observed that these arrays rapidly loose their relative coherence during their dynamical formation process. So while each droplet is coherent by itself, there is no global phase coherence between different droplets. In this work we always refer to this global coherence of the system. Furthermore, increasing the overlap of the droplet wavefunctions through an increase of the remaining weakly confining trapping direction was proposed as a way to establish more robust phase coherence in experiments.

A recent theoretical study Roccuzzo and Ancilotto (2018) examined a similar elongated trapping configuration with periodic boundaries along one axis perpendicular to the magnetic field. In this work it was shown that close to the softening of the roton mode droplets form, which are immersed in the dilute superfluid background of a BEC. This special case of coherent droplets forms only in a very narrow range of interaction strengths, while for smaller contact interaction strengths the background vanishes, leading to isolated droplets. The coherent droplets were shown to exhibit two distinct excitation modes, a phonon and a phase mode, which are hallmarks of supersolidity.

At the same time it was observed experimentally that phase coherent droplets can exist for a narrow range of contact interaction strengths Tanzi et al. (2018). However, in this experiment only the phase information after time-of-flight expansion was accessible 111During the writing of this work we became aware of another related time-of-flight investigation with Er atoms. Private Communication, F. Ferlaino and a detailed theoretical explanation of the observations is still lacking.

Here we present a comprehensive study of the supersolid properties of a trapped dipolar gas. First, we show theoretically that within the framework of the extended Gross-Pitaevskii equation (eGPE) a phase-coherent and density-modulated state can be reached dynamically for our elongated trap geometry. This state is found to be very close to the ground state of the system. Second, we experimentally realize such a state and observe its properties, both in-situ and in time-of-flight. With these two complementary observation techniques we map out the signatures of the theoretical phase diagram that clearly reveals both a coherent and an incoherent density-modulated regime.

We start by looking at the theoretical ground state of the dipolar system within the framework of the eGPE. To this end, we numerically solve the eGPE using imaginary time propagation Wenzel et al. (2017) to calculate the ground states. These simulations are performed for $2-5\times 10^{4}$ 162Dy atoms in a harmonic potential with trapping frequencies $\omega~{}=~{}2\pi\,$(18.5, 53, 81) Hz similar to Tanzi et al. (2018). Depending on the contact interaction strength this yields three distinct regimes, which are summarized on the right side of Fig. 1a for an atom number of $3.5\,\times\,10^{4}$. While for high scattering lengths of $a_{\text{s}}\geq 95\,a_{\text{0}}$ the ground state is a regular BEC, at low scattering lengths $a_{\text{s}}\lesssim 90\,a_{\text{0}}$ we recover the thoroughly studied regime of isolated droplets Kadau et al. (2016); Ferrier-Barbut et al. (2016a, b); Wenzel et al. (2017). The most interesting regime, e.g. for $a_{\text{s}}=94\,a_{\text{0}}$, is found in between these two limits and consists of droplets that are immersed in a residual BEC background. This background acts as a link between the droplets and therefore establishes phase coherence.

With these observations in mind we quantify the transition in terms of the density link between the individual droplets, in particular between the central droplet and its nearest neighbors. To do this we analyze cuts through the center of the simulated densities $n=|\psi|^{2}$ and calculate the ratio between the first minimum and the central maximum. This simple measure characterizes the density overlap between neighboring droplets and is shown in Fig. 1b. For a BEC there exists no density modulation and therefore we set this ratio to 100%. As soon as the density modulation emerges below $a_{\text{s}}\approx 94.5\,a_{\text{0}}$ this ratio is well-defined and shows a distinct jump followed by a steady decrease for lower scattering lengths. Interpreting the overlap as an order parameter this discontinuity is an indication of a first-order phase transition, further evidenced by the observation of hysteresis in our numerical simulations Sup and by the experimental signatures in Tanzi et al. (2018). Moreover, this behavior is reminiscent of the decrease of the superfluid fraction across the supersolid phase transition that was observed in previous works Roccuzzo and Ancilotto (2018); Macrì et al. (2013). However, note that while the ratio shown in Fig. 1b is a measure of the overlap of the droplets, it does not directly correspond to the actual superfluid fraction. The superfluid fraction of the system was shown to notably exceed the values of this overlap ratio Roccuzzo and Ancilotto (2018).

We extend this study to different atom numbers and show it as a phase diagram in Fig. 1c. We observe a clear phase boundary where the density modulated state becomes lower in energy than a regular BEC. This phase boundary shifts to higher scattering length with increasing atom number. While overall being located close to the roton instability our observed scaling for the simulated ground states appears different from the approximate roton scaling for a trapped gas reported in Tanzi et al. (2018), in particular at higher interaction strengths.

To study how one can dynamically establish a phase coherent state in an experiment, we perform time dependent simulations starting from the BEC ground state at $a_{\text{s}}=110\,a_{\text{0}}$ and then linearly ramp the scattering length to its final value within 30 ms. A snapshot of the system after some evolution time is shown schematically in Fig. 1a, where the dynamically reached states are compared to the calculated ground states. This comparison shows that for $a_{\text{s}}=94\,a_{\text{0}}$ the dynamically calculated wavefunction is phase-coherent and very close to the calculated ground state. In contrast to this for $a_{\text{s}}=90\,a_{\text{0}}$, the phase of the individual droplets is different and also their number does not match the ground state prediction Sup ; Wenzel et al. (2017). Using shorter ramp times in the dynamical simulations, we produce states with a different number of droplets also for $a_{\text{s}}=94\,a_{\text{0}}$. Only for ramp times $\gtrsim$ 20 ms we produce states that are close to the ground state of the system.

For a more quantitative comparison, we compare the overlap of the dynamically simulated wavefunction in real time with the ground state we get from the imaginary time evolution. We observe that incoherent droplets form soon after the end of the interaction ramp, while coherent droplets emerge slowly over several tens of milliseconds. This longer formation time for the coherent state is in agreement with a lower energy difference between this state and the initial BEC. In Fig. 2a we show the difference of the densities in the dynamical case $n(t)$ and the ground state $n_{\text{GS}}$ and Fig. 2b depicts the evolution of the phase difference between the central droplet and its nearest neighbor, which is an indicator for their phase coherence. In the density difference one can see that for the case of isolated droplets the difference is significant and actually increases with time, while for the coherent droplets the difference approaches zero after the time required for the droplets to form. The oscillations visible in the density difference after the formation correspond to a breathing mode along the droplet array. In the phase difference we observe that for the isolated droplets the phase difference increases linearly after the formation. We attribute this to the different chemical potentials of the two droplets. In contrast to this in the coherent regime we observe that the phase difference remains significantly smaller with very little variation over time. Our simulations thus reveal the existence of a state that is both density-modulated and phase coherent and can be reached dynamically by an interaction ramp into a narrow range of interaction strengths. Given the superfluid fraction and excitation spectrum that were calculated in Roccuzzo and Ancilotto (2018) we conclude that this state can be identified as a dipolar supersolid.

As a next step we include realistic three-body losses in the simulations with our experimentally measured loss coefficient $L_{3}=1.5\,\times\,10^{-40}\,\text{m}^{6}/\text{s}$ for $a_{\text{s}}=94\,a_{\text{0}}$ Boe , as well as, for reference with a lower loss rate of $L_{3}/4$. When loss is included a comparison of the densities is challenging because the ground state continuously changes with atom number. Therefore we restrict ourselves to the phase coherence shown in Fig. 2b and again observe only a small phase difference starting to form during the droplet formation process. However, this phase difference is rapidly stabilized and subsequently slowly decreases throughout the remaining lifetime of the state. An alternative way of characterizing the difference of the two wavefunctions is the fidelity, defined as $\mathcal{F}=|\left<\Psi_{\text{GS}}|\Psi(t)\right>|^{2}$. Using this we get a numerical value of $\mathcal{F}\approx 90\%$ after our experimental equilibration time of 15 ms. This shows that we expect to dynamically create a state with transient supersolid properties very close to the actual ground state even in the presence of three-body losses.

In order to investigate the formation of a phase-coherent droplet state experimentally, we prepare a quasi-pure dipolar BEC with approximately $4.5\times 10^{4}$ 162Dy atoms at a temperature below 20 nK in a tubular trap with trap frequencies $\omega=2\pi\,$(19(1), 53(1), 87(1)) Hz and $\textbf{B}\,||\,\hat{z}$, similar to Tanzi et al. (2018). We compensate the gravitational force on the atoms by ramping up a magnetic field gradient, allowing for long times of flight to probe the system. Subsequently we change the scattering length from $\sim$140 $a_{\text{0}}$ to $\sim$110 $a_{\text{0}}$ by ramping the magnetic field in 80 ms closer to a double Feshbach resonances of 162Dy located around 5.1 G Baumann et al. (2014); Boe . In order to reach the droplet regime we subsequently ramp the magnetic field again linearly within 30 ms to the final scattering length in the range between 89 $a_{\text{0}}$ and 98 $a_{\text{0}}$ 222Due to the uncertainty in the measurements of the background scattering length all experimentally quoted scattering lengths used throughout this work exhibit an uncertainty on the order of 15%. For more information on the determination of the scattering length see the Supplemental Material. We then hold the atoms for 15 ms at this field in order to equilibrate. Finally, we probe the resulting state either in-situ using far detuned phase-contrast imaging, or after time-of-flight using absorption imaging. As we are in or close to a regime where droplets are self-bound we boost the time-of-flight expansion velocity by ramping up the scattering length to $\sim$140 $a_{0}$ within 100 $\upmu$s just before the release of the atoms from the trapping potential. This has the additional advantage that density rapidly decreases, minimizing interaction effects in the expansion. Moreover, it acts as a zoom greatly increasing the extend of the interference patterns, and thus, giving access to subtle features beyond simple nearest-neighbor phase coherence. As exemplified in Fig. 3, our experiment thereby reveals both the in-situ density modulation as well as the interference pattern of multiple matter waves emerging after time of flight.

First we analyze in more detail how the in-situ density distribution changes for different atom numbers and final scattering lengths. Above a certain atom number threshold we observe the appearance of multiple droplets aligning along the weak axis of the trap. Experimentally this density modulation takes about 5-10 ms to develop. We attribute this faster formation time compared to the numerical simulations to fluctuations due to residual excitations and finite temperature, which can seed the underlying instabilities driving the phase transition. Due to our finite imaging resolution of 1 $\upmu$m, the smaller size of the droplets and the imaging aberrations arising thereupon, we cannot reliably extract the number of droplets or the overlap between droplets and the background BEC directly from our images. This resolution limited imaging also leads to a larger uncertainty in the extracted atom number compared to the data obtained in time-of-flight. As in our previous work Kadau et al. (2016), we therefore use the absolute value of the Fourier transform of the integrated in-situ density to identify images with a density modulation. To this end, we compare the spectral weight of finite momentum contributions to the weight of the central peak in momentum space Sup . This ratio is plotted as a function of scattering length and atom number in Fig. 5a. Typically, every coordinate is an average of a few experimental runs with 80 realizations in total for every scattering length and the different atom numbers have been realized by binning our experimental data. We observe that a density modulation appears above a certain atom number threshold and compare this to the boundary from the theoretical phase diagram. With this we can clearly map out a region showing a density-modulated state in the experiment. For small scattering lengths and high atom number we see the modulation amplitude decreasing again, which is caused by a washing out of the mean distribution due to fluctuations in the number of droplets and therefore also their spatial separation.

Next we study the phase coherence of the realized droplet states via interference after 30 ms time-of-flight. For our parameters we typically realize an array of several droplets. In contrast to the well known interference of two BECs Davis et al. (1995) our situation thus leads to a more complex interference pattern with multiple frequencies Hadzibabic et al. (2004); Sup . For the evaluation of our data we therefore again turn to the absolute value of the Fourier transform of the integrated image after time-of-flight. Note, that this Fourier transform of the time-of-flight density, rather than the wavefunction, does not yield again the initial in-situ wavefunction but rather provides information about the individual frequency components of the interference pattern. Examples of this for scattering lengths of $a_{\text{s}}=92.5\,a_{\text{0}}$ and $a_{\text{s}}=89\,a_{\text{0}}$ are shown in Fig. 4a and b, respectively. The gray lines are single-shot realizations, while the blue and red lines show the mean of all available realizations for a given atom number and scattering length, with more than 300 realizations for every scattering length in total.

We observe a clear difference between the two scattering lengths, caused by two distinct effects. The data for $a_{\text{s}}=92.5\,a_{\text{0}}$ can be identified with coherent droplets and show stable side peaks at the same position in every realization. The data for $a_{\text{s}}=89\,a_{\text{0}}$ correspond to incoherent droplets and shows strong fluctuations from realization to realization. The change in the side peak positions for the incoherent droplets is caused by a varying number of droplets and therefore a different initial separation from shot to shot, while the stable position of the side peaks for the coherent droplets means that the initial state is very reproducible. In particular, as in the usual double-slit interference, the height of the Fourier peaks are a measure of the interference contrast and their shot-to-shot variation thus encodes the phase coherence of the system. In our case, the individual Fourier peaks characterize nearest-neighbor, next-nearest neighbor and even higher-order coherences Sup . If the initial droplets are phase-coherent the side peaks are expected to always exhibit the same height and therefore the corresponding variance should be low. On the other hand, if the initial droplets are not phase-coherent the peak height should fluctuate from realization to realization and therefore we should observe an increase in the variance of the peak height. Taken together these two effects – fluctuating droplet number and, hence, fringe spacing, as well as incoherent phases between independent droplets – wash out the mean distribution for $a_{\text{s}}=89\,a_{\text{0}}$, while the mean distribution for $a_{\text{s}}=92.5\,a_{\text{0}}$ shows clear side peaks. This already provides clear evidence for the phase coherence between the droplets in the latter case, with the first side peak (labeled as 1 in Fig. 4a) corresponding to nearest neighbor coherence, the second peak (labeled as 2) corresponding to next-nearest neighbor coherence. For atom numbers high enough to yield more than three droplets we can even observe a small signal of next-next-nearest neighbor interference (labeled as 3).

To quantify the phase coherence we therefore calculate the variance of the height of the Fourier transform side peaks (1 and 2 in Fig. 4a) normalized to the peak height Sup . We show these ratios in Fig. 5b and c for different atom numbers and final scattering lengths. The results are in very good agreement with our in-situ results in Fig. 5a, as well as the theoretical phase diagram shown in Fig. 1c. Again we find a sharp phase boundary where the multiple droplet becomes energetically favorable compared to the BEC state. Below this boundary in the respective atom number, we observe no interference and therefore no signal in the Fourier transform. Above this threshold we always see interference, however only for a small range of contact interactions the coherence between droplets is present. The observed boundaries in all these plots are in agreement with the simulated ground state phase diagram in Fig. 1c. Combining the in-situ with the interference results reveals the signatures of the theoretical phase diagram and we see that for a small range of the contact interaction strength there exists a phase of the system showing both a density modulation as well as phase coherence and therefore the hallmark properties of a supersolid state of matter.

In conclusion we have shown theoretically and experimentally that for a narrow range of interaction strengths our dipolar quantum gas of 162Dy atoms exhibits a state that is both density-modulated and phase-coherent. Together with the dynamical study of the phase coherence in Tanzi et al. (2018), this observation is the first step towards the realization and identification of a dipolar supersolid, where in contrast to previous works Léonard et al. (2017a, b); Li et al. (2017) the self-organized density modulation is induced by the intrinsic interactions. In order to finally prove the supersolid character of the observed state beyond the phase coherence demonstrated in this work, an experimental proof of phase rigidity, and hence genuine superfluidity, is required. As a next step we therefore plan to investigate the two types of collective excitations, the phonon and phase modes. Another important aspect we plan to study is to extend the lifetime of the observed states which is currently limited due to three-body losses to approximately 20 ms. This could be accomplished by identifying a region with lower losses in the rich Feshbach spectrum of 162Dy Baumann et al. (2014); Frisch et al. (2014); Maier et al. (2015a, b).

Appendix A Experiment

The complete experimental sequence is described in more detail e.g. in Kadau et al. (2016). In short we create Bose-Einstein condensates of 162Dy in a crossed optical dipole trap (cODT) made up of two laser beams with $\lambda_{\text{cODT}}=1064\,$nm. After the creation of the BEC we change the trap to $\omega=2\pi\,$(19(1), 53(1), 87(1)) Hz within 20 ms.

Imaging is performed along the magnetic field axis ($\hat{z}$ axis) using a microscope objective that allows us to reach a resolution of about $1\,\upmu$m. We can image the atoms either with phase-contrast imaging in-situ or with resonant absorption imaging in time of flight. Both techniques result in similar atom numbers with an overall uncertainty of $<10\%$. To keep the atomic cloud in the focus of our objective we apply a magnetic field gradient to compensate gravitational forces. This applied gradient leads to a shift of the magnetic field by $-428\,$mG at the position of the atoms, which we compensate by ramping up the amplitude of the magnetic offset field at the same time as the gradient.

We observe the phase-coherent droplets (point a in Fig. 5b of the main text) at a magnetic field of 5.266(1) mG, calibrated using RF loss spectroscopy. The whole scattering length range shown in Fig. 5 corresponds to a magnetic field amplitude range of 5.261(1) to 5.271(1) mG.

Appendix B Feshbach resonances

For the measurements in this work we use a specific double Feshbach resonance of 162Dy from the rich spectrum of resonances Baumann et al. (2014); Frisch et al. (2014); Maier et al. (2015a, b) in order to control the short-range contact interaction. The resonances are located at fields of $B_{\text{1}}=5.126(1)\,\text{G}$ and $B_{\text{2}}=5.209(1)\,\text{G}$ with widths of $\Delta B_{\text{1}}=35(1)\,\text{mG}$ and $\Delta B_{\text{2}}=12(1)\,\text{mG}$, respectively. Details about the measurement and calibration of these resonances will be presented in future work Boe . Additionally we also include a broader resonance at $B_{\text{3}}=21.95(5)\,\text{G}$ with a width of $\Delta B_{\text{3}}=2.4(8)\,\text{G}$ Lucioni et al. (2018) into our considerations, since this resonance still has a small effect on the scattering length in the magnetic field range considered in this work. This allows us to calculate the scattering length $a_{\text{s}}$ from the magnetic field we measure, with only the background scattering length as a free parameter. 162Dy features a rather high background scattering length $a_{\text{bg}}=140(20)\,a_{\text{0}}$ Tang et al. (2018, 2015, 2016), which means that away from the resonances the sample is contact-dominated, while closer to the zero-crossing of the resonance we can create a dipolar-dominated sample. In the range of scattering length studied in this work our typical field stability leads to an uncertainty of $\sim 1\,a_{0}$. However, due to the uncertainty in the measurements of the background scattering length, all scattering lengths used throughout this work exhibit an uncertainty on the order of 15%.

Another important parameter controlled by the used Feshbach resonances is the three-body loss coefficient $L_{3}$. We measure a three-body loss coefficient of $L_{3}=1.33\times 10^{-41}\,\text{m}^{6}/\text{s}$ away from the Feshbach resonances for a BEC of 162Dy. Closer to the resonances the loss coefficient increases rapidly Boe leading to increased losses and therefore shorter lifetimes of the observed coherent droplets. For the scattering length where we observe the coherent droplets the losses are already enhanced by about a factor of 9 compared to the value away from the resonance. By using a different Feshbach resonance with lower loss rates from the complicated spectrum of resonances one could therefore greatly increase the lifetime of the observed states.

Appendix C In-situ evaluation

As explained in the main text we observe the appearance of an array of multiple droplets above a certain atom number threshold. The small theoretical size of the droplets compared to our imaging resolution leads to imaging aberrations. These aberrations mean that directly fitting the data, as well as directly counting the droplets is not reliable. We therefore quantify the density modulation by analyzing its contribution to the absolute value of the Fourier transform of the integrated images.

To show the general principle of this analysis we first demonstrate it on the simulated ground states that we used in Fig. 1b of the main text to show the phase transition from BEC to coherent droplets to isolated droplets. For a BEC the absolute value of the Fourier transform is peaked around zero with almost no contribution at higher momentum. Crossing the boundary into the droplet phase we observe a pronounced local maximum at a finite momentum. With decreasing scattering length we observe higher order peaks appearing that get more and more pronounced. To quantify the transition we determine the spectral weight $SW$ at finite momentum instead of the ratio of the overlap used in the main text. For the simulated ground states we therefore calculate the weight of the first order momentum peak and compare it to the central peak, by summing up the corresponding signal $S(k)$ in the range indicated by the gray lines in Fig. 6.

[TABLE]

This spectral weight is shown in Fig. 6 across the phase transition from a BEC to quantum droplets. We observe a distinct jump in the spectral weight at the phase transition. Interpreting this as an order parameter this discontinuity, similar to the measure of the overlap used in the main text, is an indication of a first-order phase transition.

For Fig. 5a of the main text we used this analysis in order to map out the phase boundary between BEC and the density modulated arrays of quantum droplets. In Fig. 7a and b we show the absolute value of the Fourier transform of the integrated in-situ images exemplary for the case of coherent droplets ($a_{\text{s}}=92.5\,a_{\text{0}}$) and incoherent droplets ($a_{\text{s}}=89\,a_{\text{0}}$), respectively. Differently than for the simulated ground states we do not observe any higher order peaks in our experimental data, due to the finite imaging resolution. In order to calculate the spectral weight for the experimental data we sum up the signal with $|k|<0.5\upmu\text{m}^{-1}$ for the main peak and with $0.5\upmu\text{m}^{-1}<|k|<2\upmu\text{m}^{-1}$ for the finite momentum contribution. The boundary between the two ranges is indicated by the dashed gray lines in Fig. 7a and b. For large $|k|$ we also put a boundary due to the imaging aberrations we observe. These imaging aberration are more pronounced for the isolated droplets compared to the coherent droplets.

Appendix D Interference of an array of phase independent matter waves

The interference of two waves, whether they are light waves after a double slit or matter waves from two separate sources Andrews et al. (1997); Hadzibabic et al. (2004), is approximately described by the well-known interference pattern

[TABLE]

In this $A_{1}$ is the contrast of the interference pattern and $\Phi_{1}$ the relative phase of the two waves. In the case of matter waves initially confined in a harmonic trap with trapping frequency $\omega_{z}$, the interference fringe spacing given by $D=ht/(md)$, with $d$ the separation of the initial traps, $t$ the time-of-flight, $Z_{0}=\hbar t/(ml)$ the size of the matter waves in time-of-flight, $m$ the atomic mass and $l=\sqrt{\hbar/(2m\omega_{z})}$ the width of the initial density distribution. Here, we have assumed long times of flight $t$.

In comparison to this, a more complex interference pattern characterized by multiple frequencies arises from the interference of light waves passing through a multi-slit arrangement or for multiple interfering matter waves. Here, our aim is to characterize the phase coherence between these individual matter waves. Following Hadzibabic et al. (2004) we can write the interference pattern $I(z)$ for an array of $N$ independent matter waves according to

[TABLE]

where the amplitudes $A_{n}$ and the phases $\Phi_{n}$ are given by the modulus and the argument of a sum of complex numbers

[TABLE]

In this $\phi_{j}$ is the phase of the single matter wave at position $j$. The contrast amplitudes $A_{n}$ and the phases $\Phi_{n}$ of the interference pattern are thus defined by the phase differences of all matter wave pairs initially separated by the distance $nd$. Here the index $n$ runs over all the possible pairs of interfering waves, starting from the nearest neighbors with $n=2$ and ending with the interference of the two outermost waves with $n=N-1$. This is schematically indicated in the inset of Fig. 4a of the main text.

The result in Eq. 3 can be obtained from the general expansion of matter waves, assuming that they all contain the same average number of particles and that the atoms do not interact during the expansion.

For only two interfering matter waves the interference pattern simplifies to Eq. 2, while for infinitely many interfering waves it corresponds to the interference pattern of a grating given by

[TABLE]

Using Eq. 3 we can generate 300 idealized interference patterns for four phase coherent matter waves (Fig. 8a), or for four matter waves with randomly distributed phases (Fig. 8b). The blue solid line corresponds to the mean interference pattern in the coherent case ($\Phi_{n}=0$) and the red solid line to the mean interference pattern in the incoherent case (random distributed $\Phi_{n}$). Examples of single shot realizations are shown in gray. We can see that the interference pattern is perfectly reproducible in the case of coherent droplets, while there are strong fluctuations in the incoherent case, which lead to a washed out mean distribution even if the single-shots show a large contrast. This means that the appearance of high contrast interference pattern alone is not enough to proof phase coherence.

To analyze this in more detail we can use the same evaluation with the Fourier transform as for the experimental data. For Eq. 3 this leads to the spectra shown in Fig. 8c and d, where the peaks of the absolute value of the Fourier transform are located at multiples of $d$. The height of these side peaks are proportional to $A_{n}$. Again the blue and red line correspond to the respective mean distribution, while the gray lines are exemplary single realizations. As expected from Eq. 4 the single realizations for the phase coherent waves show stable side peaks with very low variance in the peak heights for different realizations, while the peak heights for the incoherent matter waves show strong fluctuations. The Fourier transform - and specifically the variance of the peak heights - thus allows us to individually quantify the coherences of order $n$.

Appendix E Revealing the experimental signatures of the phase diagram

In Fig. 5b and c of the main text we show the variance of the height of the side peaks as a measure of the phase coherence between neighboring droplets, as it was explained in the previous section and in Fig. 9. In order to get a result independent of the peak height, we normalize it to the corresponding mean peak height

[TABLE]

with $N$ the number of realizations, $p_{i}$ the peak height of the $i$-th realization and $\left<p\right>$ the mean peak height. In this evaluation we subtract the background signal which we observe in the Fourier transforms even in the absence of interference from the mean peak height.

Instead of looking at the variance of the heights of the side peaks in the Fourier transforms (see Fig. 5b and c of the main text), we can also directly study the peak height as an indicator for the phase coherence. To make the analysis more robust, rather than directly identifying a peak height in the noisy experimental data, we calculate the weight of the mean distribution around the side peaks and compare this to the weight of the main peak. We show this ratio in Fig. 9a-c for different atom numbers and final scattering lengths, up to the signal resulting from the next-next-nearest neighbor interference. The black dashed line is again the guide to the eye from Fig. 1c of the main text. These results are in very good agreement with the results obtained through the analysis of the variance in the peak heights that are shown in Fig. 5b-c of the main text, as well as the theoretical phase diagram shown in Fig. 1c of the main text.

Appendix F Hysteresis in the numerical simulations

The observed discontinuity in the measure of the overlap of neighboring droplets in Fig. 1b of the main text is an indication of a first-order phase transition. This first-order nature of the transition is further evidenced by the observation of hysteresis in Tanzi et al. (2018); Kadau et al. (2016).

To further study the nature of the transition, we numerically solve the eGPE using imaginary time propagation Wenzel et al. (2017) to calculate the ground states, once with a regular BEC as an initial wavefunction (labelled ”down” in Fig. 10) and once with a density modulated, multiple droplet initial wavefunction (labelled ”up”). These simulations are performed for $3.5\times 10^{4}$ 162Dy atoms in a harmonic potential with trapping frequencies $\omega~{}=~{}2\pi\,$(18.5, 53, 81) Hz similar to Tanzi et al. (2018). Again we recover the three distinct regimes of a regular BEC and the coherent and incoherent quantum droplets depending on the contact interaction strength. However, we observe that the boundary between BEC and density modulated droplets is shifted depending on the direction from which we approach the transition. This hysteretic behavior is shown in Fig. 10, where we have used the peak density as a simple probe to identify the transition point. Note that while the peak density can be used to identify the transition between BEC and quantum droplet regime, it can not distinguish between incoherent and coherent droplets. This is because the peak density is significantly higher for the quantum droplets compared to the BEC, but it is similar for the two different droplet regimes.