Density Fluctuations across the Superfluid-Supersolid Phase Transition in a Dipolar Quantum Gas

Phase transitions share the universal feature of enhanced fluctuations near the transition point. Here we show that density fluctuations reveal how a Bose-Einstein condensate of dipolar atoms spontaneously breaks its translation symmetry and enters the supersolid state of matter -- a phase that combines superfluidity with crystalline order. We report on the first direct in situ measurement of density fluctuations across the superfluid-supersolid phase transition. This allows us to introduce a general and straightforward way to extract the static structure factor, estimate the spectrum of elementary excitations and image the dominant fluctuation patterns. We observe a strong response in the static structure factor and infer a distinct roton minimum in the dispersion relation. Furthermore, we show that the characteristic fluctuations correspond to elementary excitations such as the roton modes, which have been theoretically predicted to be dominant at the quantum critical point, and that the supersolid state supports both superfluid as well as crystal phonons.

Fluctuations play a central role in quantum manybody systems. They connect the response and correlation of the system to its excitation spectrum, instabilities, phase transitions and thermodynamic properties. A quantity that is fundamental to the theoretical description of fluctuations in many-body systems is the structure factor, which can be formulated as the Fourier transform of the density-density correlation function [1,2]. Superfluid helium is an important example of a quantum manybody state, where the determination of the structure factor was crucial to understand its elementary excitations and therefore improved our understanding of superfluidity [2][3][4][5]. In the case of quantum gases, the structure factor of Bose-Einstein condensates (BECs) and superfluid Fermi gases is often investigated by Bragg spectroscopy [6][7][8]. In contact-interacting BECs this enabled the study of the spectrum and collective modes [9]. In dipolar BECs, it has provided indications of the roton minimum in the dispersion relation [10] analogous to the neutron and X-ray scattering data for helium [2,11]. A different approach is to look at the condensate density directly in situ, which provides access to finite temperature and quantum fluctuations [12][13][14][15][16][17][18], and enables to extract the static structure factor simultaneously at all momenta.
The roton minimum both in helium and dipolar quantum gases is accompanied by a characteristic peak in the static structure factor close to the roton momentum [2,[19][20][21]. Unlike in helium however, the contact interactions in dipolar quantum gases are tunable [22]. This tunability allows for precise control of the dispersion relation and the controllable softening of the roton minimum. The roton modes associated to this minimum manifest as density modulations on top of the ground-state density distribution [23,24]. An instability in the ground state can appear once the roton minimum is sufficiently soft. Since these modes represent precursors to solidification, dipolar BECs have long been proposed as candidates for the elusive supersolid state of matter, which simultaneously combines crystalline order with superfluidity [25].
Recently, a dipolar supersolid state of matter has been realized through a phase transition from a BEC to an array of coherent quantum droplets [26][27][28][29][30][31][32] by precisely tuning the contact interaction strength. Close to the transition point these droplets are immersed in a superfluid background, and by lowering the scattering length further the superfluid fraction decreases towards a regime of isolated droplets. As the superfluid-supersolid phase transition is governed by intrinsic interactions it is of interest to study the fluctuations that emerge across the transition [21,33,34], facilitate structure formation [15,20,[35][36][37] and give rise to the supersolid state. The low-lying collective modes were shown to be particularly interesting regarding the aspect of supersolidity in this system [29][30][31]38]. Those modes are facilitated by a continuous superfluidity across the droplet array, despite the translational symmetry breaking. Schematic of the dispersion relationhω(k) and associated static structure factor S(k) of an elongated and strongly dipolar BEC [21]. A decrease in scattering length as causes a roton minimum to emerge in the excitation spectrum, associated with a characteristic peak in the static structure factor. The roton momentum krot is indicated where the roton minimum drops near zero. (c) For a given scattering length as we observe a large number of in situ densities nj(r), calculate their mean n(r) and the density fluctuations δnj(r) = nj(r) − n(r) as the deviation of the in situ images from their mean. Investigating the mean power spectrum of the fluctuations |δn(k)| 2 for different scattering lengths across the transition allows us to directly observe the static structure factor as the system passes from BEC to supersolid to isolated droplet states. The colored arrow on top indicates the direction towards lower scattering lengths, passing from BEC to supersolid to isolated droplet regimes. The colormap used for the images shows the normalized amplitude of densities, density fluctuations and mean power spectra, respectively.
Here we provide the first direct in situ observation of density fluctuations across the superfluid-supersolid phase transition in a trapped dipolar quantum gas. By analyzing hundreds of in situ images of the atomic cloud around the phase transition point we spatially resolve characteristic fluctuation patterns that arise across the transition. From the observed fluctuations we determine the static structure factor and estimate the spectrum of elementary excitations. We observe a strong peak in the static structure and an associated roton minimum in the dispersion relation. Moreover, we experimentally determine that the dominant fluctuations at the transition point correspond to two degenerate roton modes [29,38] and that the supersolid state supports both superfluid as well as crystal phonons in a narrow range of scattering lengths. Our study combines the fluctuations with the excitation spectrum of a dipolar supersolid and highlights its bipartite nature between the superfluid BEC and the crystalline isolated droplets.
The rotonic dispersion relation of strongly dipolar BECs [10,23,24] is schematically shown in Fig. 1(a). The system becomes more susceptible to density fluctuations as the roton minimum softens. These density fluctuations are associated with a characteristic peak [2,[39][40][41][42] in the static structure factor, illustrated in Fig. 1(b). This can be understood by considering the general Feynman-Bijl formula S(k) =h 2 k 2 /2mε(k) [2,3], connecting the static structure factor S(k) to the excitation spectrum (k) at zero temperature. As the energy of the roton modes drops near zero, the density fluctuations and thus the structure factor increase dramatically. Eventually, the roton minimum has sufficiently softened in order for the system to enter the roton instability. This roton instability triggers the phase transition to a dipolar supersolid and arrays of isolated quantum droplets.
To study the static structure factor experimentally, we prepare a dipolar BEC with typically 40 × 10 3 162 Dy atoms at temperatures of approximately 20 nK in a cigar-shaped trap with trapping frequencies ω/2π = [30 (1), 89(2), 108(2)] Hz and a magnetic field oriented alongŷ (see Methods). The scattering length is tuned via a double Feshbach resonance [43] to final values between 90 a 0 and 105 a 0 by linearly ramping the magnetic field in 30 ms. We wait for 15 ms to allow for the system to equilibrate and then the atoms are imaged using phase-contrast imaging along theẑ-axis with a resolution of ∼ 1 µm. We find either a BEC, a supersolid phase (SSP) or isolated droplets (ID) for large (a s 105 a 0 ), intermediate (a s 98.4 a 0 ) and small (a s 90 a 0 ) scattering lenghts, respectively [29]. We accumulate enough averages for a statistical evaluation of the structure factor by repeating the experiment around 200 times for every scattering length.
We obtain S(k) experimentally by analyzing the in situ images as illustrated in Fig. 1(c). For every scattering length we center the in situ densities n j (r) to their center of mass and normalize them to the mean atom number. With the former step we remove contributions of the dipole center of mass motion [29] and with the latter we correct for shot to shot total atom number fluctuations ( [13]; Methods) that would otherwise give contributions to S(k) near k = 0. From these in situ images we obtain the mean image n(r) and the density fluctuations δn j (r) = n j (r)− n(r) as the deviation of the in situ images from their mean. With the Fourier transform of the density fluctuation δn j (k) = d 3 rδn j (r)e ik·r we obtain the mean power spectrum of the fluctuations |δn(k)| 2 . In homogeneous systems, the static structure factor can be directly written as S(k) = |δn(k)| 2 /N , where N is the atom number [2,15,17]. In practice, the interpretation is less straightforward [18,[44][45][46] since the expectation values of the density n(r) are spatially dependent due to the finite size and the translational symmetry breaking in the supersolid and droplet regime. Nonetheless S(k) gives insight into the strength of fluctuations [13,18,[47][48][49] and is a quantity that can be continuously evaluated from the BEC via the supersolid to the isolated droplet regime. We note that our evaluation is limited to intermediate momenta between k min /2π 0.08 µm −1 and k max /2π 1 µm −1 due to the finite system size and the finite resolution of our imaging system, respectively [15,18]. We extract the static structure factor S(k x , k y , k z = 0) cut along the k z = 0 plane according to the Fourier slice theorem since the atomic densities are integrated along the line of sight during the imaging process ( [50]; Methods). Due to the cigar-shaped trap geometry, the fluctuations predominantly show structure alongk x (see Fig. 1(c)), which allows us to extract a cut of the mean power spectrum at k y = 0 to obtain the 1D structure factor S(k x ).
Using the above described analysis, we obtain S(k x ) across the phase transition as shown in Fig. 2. In the BEC regime at a s 104 a 0 , we find the structure factor to be relatively flat with the exception of a small peak at around k x /2π 0.25 µm −1 . This peak is an indication that far in the BEC regime roton modes can be excited [29] and consequently that the spectrum features modifications from a purely contact interacting quantum gas.
As the scattering length is reduced, the position and amplitude of this characteristic peak are observed to increase continuously towards the phase transition point (a s 98.4 a 0 ). A continuously growing peak ampli- tude of the structure factor signals enhanced fluctuations, consistent with a softening roton minimum towards the transition point. The structure factor reaches its maximum value as a function of the scattering length at the transition point and is located at the roton momentum k rot /2π 0.29 µm −1 . This value is mainly given by the harmonic oscillator length l y along the magnetic field direction [52]. At the transition point, the enhanced fluctuations provide the roton instability and lead to the formation of supersolid quantum droplets, whose spacing d 3 µm smoothly matches the roton wavelength 2π/k rot [27,29]. The observed increase of the roton momentum towards smaller scattering lengths can be understood in a variational approach of elongated dipolar condenstates [53]. Around the transition point, the in situ densities we observe from shot to shot show droplets immersed in an overall BEC background, constituting the supersolid state of matter [26-32, 38, 54]. Here the density fluctuation patterns become most clear and show spatial oscillations (see Fig. 1(c), middle column). These characteristic fluctuations can directly be attributed to the symmetric and antisymmetric roton modes we found in our previous work [38].
After crossing the phase transition (a s < ∼ 98.4 a 0 ) the Experimentally determined excitation energy ω(kx) according to equation (1) assuming a temperature of 20 nK, for scattering lengths above the phase transition point. A clear roton minimum at finite momentum is observed that softens towards the transition point.
peak amplitude of the structure factor decreases and a shoulder develops at smaller momenta. This shoulder increases further for smaller scattering lengths and eventually leads to a double-peak structure as seen in Fig. 2. The origin of this rising double peak can be understood by means of a principle component analysis.
The maximum of the structure factor for different scattering lengths acts as a measure of the density fluctuation strength across the superfluid to supersolid phase transition. It quickly increases from the BEC side when approaching the phase transition indicating a significant enhancement of the characteristic fluctuations close to the phase transition point. We see that the increase from the BEC side towards the phase transition is sharper than the decrease on the droplet side. The magnitude of the structure factor (S max = 260) can mainly be explained by thermal enhancement of the participating low-energy modes.
To estimate the dispersion relation based on the experimentally determined structure factor we use the relation which extends the Feynman-Bijl formula, S(k) =h 2 k 2 /2mε(k), valid at T = 0, to nonzero temperatures T [2,45]. At nonzero temperatures and small excitation energies the contribution of low-lying modes to the structure factor can easily be enhanced by several orders of magnitude. Close to the transition point where the roton gap ∆ rot is small compared to the temperature of the system (h∆ rot /k B T < ∼ 1), equation (1) can be expanded and the peak of the static structure factor scales as S max ∼ T /∆ 2 rot [19]. Note that equation (1) is an excellent description of the structure factor for a weakly interacting superfluid, where the excitation spectrum is dominated by a single mode, and where the influence of the quantum as well as thermal depletion can be ignored. Although we study a finite system, leading to a discrete excitation spectrum, a continuous approximation to the dispersion relation yields a meaningful estimate for the excitation energies (see Methods). We show the resulting spectrum Fig. 3. To do so we assumed a mean temperature of 20 nK, a conservative approximation to include additional minor heating during the preparation (see Methods). In Fig. 3, one can see a small roton minimum already well above the trap frequency of 30 Hz for a large scattering length. The roton minimum softens and moves towards higher momenta k x as the scattering length is lowered and finally reaches its lowest energy at the phase transition point. After crossing the phase transition point, when the crystalline structure has developed, the excitation spectrum should have a band structure due to the translational symmetry breaking. In this case equation (1) is no longer necessarily justified as several modes contribute to the excitation spectrum, and therefore it is no longer straightforward to extract the excitation spectrum from the measured static structure factor.
To gain a better insight into the modes that dominantly contribute to the fluctuations, we use principal component analysis (PCA) [55] on the density fluctuations for all scattering lengths combined. This modelfree statistical analysis is a general method to extract dominant components or to reduce the dimensionality of a dataset. We study the principal components (PCs) across the phase transition since there is a direct correspondence [56] to the dominant collective excitations obtained with the Bogoliubov-de-Gennes (BdG) formalism [38]. This allows us to identify and compare the most dominant PCs with specific BdG modes and study how their weight behaves across the transition, as shown in Fig. 4 and Fig. 5.
The first principal component is structureless and only represents the global atom number fluctuation [56,57]. The subsequent principal components are shown in Fig. 4(a)-(b). These two components are dominant across the phase transition and represent a periodic spatial pattern. Close to the transition point we can identify these characteristic patterns in many single-shot realizations of the density fluctuations, as shown in the central column of Fig. 1(c). We compare the profiles of these PCs to the antisymmetric and symmetric roton modes from the BdG theory at the transition point in Fig. 4(c)-(d) and find them to be in excellent agreement. These two roton modes are developing into the Goldstone [29] and amplitude (Higgs) [38] modes of the supersolid.
The mean absolute value of the weights for these two PCs are shown in Fig. 4(e) as an indication of their strength across the phase transition. Starting from the BEC side, where they are comparable in strength to other modes, these PCs gain rapidly in strength as the phase transition is approached (see Methods). We note that the maximum of the structure factor behaves similar to the weights of these two PCs as a function of scatter- The gray area indicates the supersolid region previously determined [29]. Error bars indicate the standard error of the mean.
ing length. Leading up to the quantum critical point at a s 98.4 a 0 these two modes have identical weights, in accordance to our previous work [38], in which we showed that the two roton modes remain degenerate while softening towards the phase transition.
Further into the isolated droplet regime the weight of the roton PCs decrease and other PCs become more important because further modes are softening. In Fig. 5 we present the three next higher PCs that correspond to the BEC phonon (a) and the antisymmetric (b) and symmetric (c) crystal phonon, respectively. The quadrupole mode of the BEC has a relatively high weight in the BEC regime and vanishes for small scattering lengths towards the isolated droplet regime. The breathing or lowest phonon modes of the droplet crystal show a clear splitting in the Fourier transform (Fig. 5(d)-(e)) explaining the observed double-peak structure in S(k x ) for low scat- tering lengths. This can be understood as the appearance of the band structure, where excitations are split around the edge of the Brillouin zone. These modes only have an appreciable weight for low scattering length after crossing the phase transition ( Fig. 5(g)). In the experiment the excitation of the crystal breathing mode is further enhanced by the preparation process [27,30]. Note that there is a small region close to the phase transition where both types of modes have a non-vanishing weight. This subtle feature shows the coexistence of both BEC and droplet crystal, which is a prerequisite of the supersolid nature of the phase. One can see that the supersolid state supports both types of excitations -the phonon of the superfluid BEC and the crystalline phonons of the droplets.
In conclusion, we reported the first in situ measurement of the density fluctuations across the superfluid to supersolid phase transition in a dipolar quantum gas. We quantified the fluctuation strength across the transition by the static structure factor S(k x ) using a statistical evaluation of in situ images and found a characteristic peak in S(k x ) that strongly increases towards the phase transition point. We showed that this peak is unambiguously dominated by the low-lying modes of the rotonic dispersion relation. The characteristic fluctuations close to the transition point are stronger compared to the BEC or isolated droplet regime. The large amplitude of the measured static structure factor reveals the important role played by temperature at the phase transition, an aspect which has so far been absent in the discussion of the superfluid-supersolid phase transition. Using principal component analysis we spatially resolved the dominant fluctuations and identified them as two roton modes. Furthermore, we showed that the supersolid supports both superfluid and crystal phonons. Our study provides a promising outlook to extract thermodynamic properties [20] and possibly universal access to the condensate fraction [37] of the supersolid state. Exciting avenues for future work includes using fluctuations as a tool for thermometry of the supersolid state [58], understanding the out of equilibrium dynamics that arise when crossing the phase transition [59,60], and exploring the roles of fluctuations in the Kibble-Zurek mechanism [33,34,61].

COMPETING INTERESTS
The authors declare no competing financial interests.

DATA AVAILABILITY
Correspondence and requests for materials should be addressed to T.P. (email: t.pfau@physik.unistuttgart.de).

A. Experimental protocol
The complete experimental procedure has been described in detail in our previous publications [27,29,62]. After preparing a quasipure BEC of 162 Dy with T 10 nK in a crossed optical dipole trap at 1064 nm, we reshape the trap within 20 ms to its final geometry with trap frequencies of ω/2π = [30(1), 89(2), 108(2)] Hz. The magnetic field is oriented along theŷ-axis and is used to tune the contact interaction strength.
We ramp the magnetic field in a two-step ramp closer to the double Feshbach resonance near 5.1 G to tune the scattering length from its initial background value of a bg = 140(20) a 0 [63][64][65] to a final value between 90 a 0 and 105 a 0 . This corresponds to the droplet and BEC regime, respectively. We expect the preparation scheme to induce some additional heating and thus assume a temperature of 20 nK in our later analysis of S(k x ) and ω(k x ). After 15 ms of free evolution the droplets have formed and equilibrated. We finally probe the atomic system using in situ phase-contrast imaging along the verticalẑ-axis, which is done using a microscope objective featuring a numerical aperture of 0.3. We reach an imaging resolution of 1 µm.

B. Analysis method and principle component analysis
In the following we will describe our analysis procedure. We start with a large set of images that contains 200 averages for each scattering length. We center the images with respect to their center-of-mass to get rid of the otherwise dominating dipole modes. Afterwards, we post-select on atom number for every scattering length and only take images with an atom number that lays in an interval of ±30% around the mean atom number at that scattering length. We have confirmed that changing the tolerance in the post selection does not affect the features of the structure factor. In a next step, we normalize each image to its atom numberñ j = n j /N j and calculate the fluctuations δñ j (r) =ñ j (r) − ñ(r) as the deviation of the normalized imageñ j (r) from the mean image ñ(r) . The structure factor is then given by the power spectrum of these fluctuations multiplied by the mean atom numberN where δñ j (k) = F[δñ j ](k) = d 2 r δñ j (r)e ik·r is the Fourier transform of the normalized fluctuations which were obtained from the line integrated images. According to the Fourier-slice theorem, taking the Fourier transform of two-dimensional integrated atomic densities is connected with the Fourier transform of the full Supplementary Figure 1. Average atom number. Average atom number for the scattering length range in the experiment. Crossing the phase transition from BEC to droplet arrays the increasing density leads to larger three-body losses causing lower atom numbers for smaller scattering lengths. Error bars indicate one standard deviation from the mean.
three-dimensional atomic density distribution. For an arbitrary function f (x, y, z) and its projection p(x, y) = dz f (x, y, z) the Fourier-slice theorem reads As a result, the static structure factor S(k x , k y ) extracted from the integrated densities is in fact a slice through the structure factorS(k x , k y , k z ) one would get if one had access to the full three-dimensional density distribution S(k x , k y ) =S(k x , k y , 0). Further it is worth noticing that normalizing each image to its atom number when comparing the structure factor at different scattering lengths has an effect in the small k regime. This part of the structure factor is highly dependent on the total atom number which is changing across the transition due to higher three-body losses towards the droplet regime (Fig. (1)). Finally we take a cut at k y /2π = 0 to determine the one-dimensional static structure factor S(k x ) we present in Fig. 2 of the main text. We confirmed that averaging within the k y -values that correspond to our imaging resolution is below our statistical error obtained by bootstrapping [51,[66][67][68].
Two natural boundaries arise at small-k and at largek because we are probing a finite system. The small-k cut-off simply comes from the finite size of the atomic cloud on the image. This results in a lowest k-value k min /2π 0.08 µm −1 that a possible excitation must have in a system of size L 12 µm. In contrast, the high-k cut-off has its origin in the finite imaging resolution of the microscope objective which leads to k max /2π 1 µm −1 .
We note that Ref. [27] has shown that the dynamical preparation scheme with the scattering length ramp only leads to states close to the actual ground state of the system for the BEC and supersolid regime. In regime of isolated droplets, the preparation process can lead to a different number of droplets from realization to realization. In comparison to the BdG theory, where only excitation on top of the actual ground state are considered, this could lead to increased fluctuations for the isolated droplets in the experiment. For all scattering lengths, we expect the dynamical nature of the sample preparation to slightly modify the observed fluctuations compared to purely thermally populated collective modes.
In order to get a better intuition of the different contributions to the experimental structure factor, we use a model-free statistical analysis, the principle component analysis (PCA) [55]. PCA has a wide range of general applications [55], from image analysis to dimensional reduction of large dataset. For ultracold atomic systems, it turns out that there is a direct correspondence [56] between the principal components (PCs) and the density variation of a mode obtained by the Bogoliubov-de-Gennes (BdG) formalism. This allows us to identify a specific PC with a certain BdG mode as long as the imaging noise is negligible. One of PCA's properties is that the signal (in our case the centered images n j (r)) can be reconstructed exactly using a superposition of all PCs, using their respective weights. However, a small subset of PCs accounts for most of the information contained within the scattering length scan. This becomes essential when PCA is used for dimensional reduction.
In the experiment we combine all images independent of their scattering length to one large dataset. This allows us not only to illustrate the different contributions to the fluctuations and therefore the structure factor but also to track the weight of different PCs over a certain scattering length range [56]. We confirm by treating the BEC and droplet regime separately that the relevant lowlying PCs do not change significantly except for their order (or variance). The first two PCs, namely the roton modes, do not change their shape over the complete scattering length range. By limiting the scattering length range to the BEC or droplet side we find in addition to the roton modes only the quadrupole or breathing modes, respectively. This is in agreement with the results presented in the main text, where we analyzed the complete dataset and observed a vanishing weight of these components on the respective side of the phase transition.

C. Simulation details
In this section, we briefly summarize simulation details not explained in the main text. A system of dipolar atoms that undergoes the transition from a BEC to a supersolid can be described by means of the extended Gross-Pitaevskii equation (eGPE) [54,69,70] where we define H GP := H 0 + g|ψ| 2 + Φ dip [ψ] + g qf |ψ| 3 and ψ is normalized to the atom number N = d 3 r |ψ(r)| 2 . The term H 0 = −h 2 ∇ 2 /2m + V ext contains the kinetic energy and trap confinement V ext (r) = m(ω 2 x x 2 + ω 2 y y 2 + ω 2 z z 2 )/2. The contact interaction strength g = 4πh 2 a s /m is given by the scattering length a s .
The dipolar mean field po- is the dipolar interaction for aligned dipoles. The strength of the dipolar interaction is given by the parameter g dd = 4πh 2 a dd /m characterized by the dipolar length a dd = µ 0 µ 2 m m/(12πh 2 ). Here, µ m is the magnetic moment and ϑ is the angle between r and the magnetic field axis. Furthermore, the quantity g qf = 32/(3 √ π)ga 3/2 Q 5 ( dd ) represents quantum fluctuations within the local density approximation for dipolar systems [71,72], where dd = g dd /g = a dd /a s is the relative dipolar strength. In our simulations, we use the approximation Q 5 ( dd ) = 1 + 3 2 2 dd [70,[72][73][74][75]. The mean-field dipolar potential is effectively calculated using a Fourier transform, where we use a spherical cutoff for the dipolar potential. The cutoff radius is set to the size of the simulations space such that there is no spurious interaction between periodic images [69,76,77].
In order to study the elementary excitations, we use the BdG formalism as described in Ref. [38] and linearly expand the wavefunction ψ(r, t) = ψ 0 (r) + λ[u(r)e −iωt + v * (r)e iωt ]e −iµt/h around the ground state ψ 0 with the chemical potential µ. This ansatz together with equation (4) leads to a system of linear equations that can be expressed in matrix form. For the actual form of the BdG matrix representation we refer the interested reader to the literature [38,52,75]. We numerically solve these equations to obtain the modes u and v corresponding to the lowest excitation energieshω. Due to our finite-sized system, we obtain a discrete spectrum of elementary excitations, rather than a continuous one. In addition there is a systematic shift between the theoretical and experimental transition point of approximately ∆a s 2.6 a 0 [10,31,43,52,78].

D. Temperature-enhanced static structure factor
In order to obtain the dynamic structure factor from the Bogoliubov-spectrum we first define the theoretical density fluctuation corresponding to the mode j as δn j = f * j ψ 0 with f j = u j + v j and the ground state ψ 0 . The static structure factor then reads +n(ω j , T )δ(ω + ω j )), where δn j (k) is the Fourier transform of the fluctuations corresponding to the j'th mode and Supplementary Figure 2. Dynamic structure factor from numerical simulations. Dynamic structure factor for two different scattering lengths before (a) and just after the transition (b). The finite system size results in discrete modes. The red line is calculated via the Feynman-Bijl formula equation (1) at T = 0. The softening roton mode clearly dominates the dynamic structure factor when approaching the phase transition. The colorbar shows the amplitude of the dynamic structure factor. For illustration purposes the dynamic structure factor was convolved along the ω-axis with a Gaussian of width σ = 0.5 Hz.
The static structure factor S(k) = N −1 dω S(k, ω) is then given after integrating along ω [2] which leads to the result similar to equation (1) of the main text. Here we made use of the identity 2n(ω j , T ) + 1 = coth(hω j /2k B T ). Equation (6) clearly indicates that low-lying modes satisfying the conditionhω j k B T can be drastically enhanced. For our situation this means that the roton modes that soften towards the phase transition dominate the structure factor compared to all other modes. Only a finite number of modes are essentially necessary to quantitatively describe experiments at finite temperature.
In Fig. 2 we show the zero-temperature dynamic structure factor for two different scattering lengths, in the BEC phase and just after the phase transition. Due to our finite system size one can clearly see the discrete mode structure along the k x -and ω-axis. The color bar indicates the amplitude of the dynamic structure factor. The red line shows the excitation energy obtained via the Feynman-Bijl formula, equation (1) of the main text at T = 0. It illustrates that the continuous dispersion relation obtained from this equation yields a meaningful estimate of the discrete Bogoliubov spectrum.
E. Temperature dependence of the experimental excitation spectrum In Fig. 3 of the main text we show the dispersion relation determined from the experimental static structure factor. This is done by solving equation (1) numerically. For dilute, weakly-interacting Bose gases with negligible quantum and thermal depletion this is a valid approximation. Although the system is undergoing a phase transition, the gas parameter na 3 s 10 −5 is still small enough to consider it as dilute and weakly-interacting with negligible quantum depletion, which is the case in our situation. The thermal component at a temperature of 20 nK can be estimated to be less than 5 % [2].
In the situation where, as we have seen above, the structure factor is mainly dominated by the contribution from the two degenerate roton modes, it is possible to expand coth x 1/x in equation (1) for small energies [19], yielding S(k, T ) S(k, 0) 2k B T ε(k) =h 2 k 2 k B T mε(k) 2 .
Equation (7) might then also be used in a second step to calculate back the zero-temperature static structure factor. In fact we find due to the validity of this expansion that the expression (k) h 2 k 2 k B T /mS(k, T ) is a good approximation to the numerical solution of equation (1).
The usage of equation (1) requires knowledge of the temperature in the system. As measuring low temperatures with a vanishing thermal fraction is rather challenging, we show in Fig. 3 how the uncertainty of our temperature of 20(5) nK affects the determined excitation spectrum for two scattering lengths.