Observation of the Roton Mode in a Dipolar Quantum Gas

Observation of the Roton Mode in a Dipolar Quantum Gas


The concept of a roton, a special kind of elementary excitation, forming a minimum of energy at finite momentum, has been essential to understand the properties of superfluid He. In quantum liquids, rotons arise from strong interparticle interactions, whose microscopic description remains debated. In the realm of highly-controllable quantum gases, a roton mode has been predicted to emerge due to dipolar interparticle interactions despite of their weakly-interacting character. Yet it has remained elusive to observations. Here we report momentum-distribution measurements in dipolar quantum gases of highly-magnetic erbium atoms, revealing the existence of the long-sought roton. We observe the appearance of peculiar peaks at well-defined momentum matching the inverse of the tight confinement length as expected for dipolar rotons. Our combined theoretical and experimental work demonstrates unambiguously the roton softening of the excitation spectrum and provides a further step in the quest towards supersolidity.


Quantum properties of matter continuously challenge our intuition, especially when many-body effects emerge at a macroscopic scale. In this regard, the phenomenon of superfluidity is a paradigmatic case, which continues to reveal fascinating facets since its discovery in the late 1930s (1); (2). A major breakthrough in understanding superfluidity thrived on the concept of quasiparticles, introduced by Landau in 1941 (3). Quasiparticles are elementary excitations of momentum , whose energies define the dispersion (energy-momentum) relation .

To explain the special thermodynamic properties of superfluid , Landau postulated the existence of two types of low-energy quasiparticles: phonons, referring to low- acoustic waves, and rotons, gapped excitations at finite initially interpreted as elementary vortices. These two types of excitations were later unified in a unique dispersion relation (4), which continuously evolves from linear at low (phonons) to parabolic-like with a minimum (roton) at a finite . Neutron scattering experiments confirmed Landau‘s remarkable intuition (5). In liquid He, scales as the inverse of the interatomic distance. This manifests a tendency of the system to establish a local order, which is driven by the strong correlations among the atoms. Yet, the same strongly-correlated nature of helium handicaps a microscopic understanding from first principles of the roton properties (1); (6).

In the realm of low-temperature quantum physics, ultra-cold quantum gases realise the other extreme limit for which the interparticle interactions - and correlations - are weak, meaning that classically their range of action is much smaller than the mean interparticle distance (7); (2). Because of this diluteness, roton excitations are typically absent in ordinary quantum gases, i. e. in Bose-Einstein condensates (BECs) with contact interactions (2). However, the degree of tunability in BECs is remarkable and a roton-like softening has been induced in hybrid systems via cavity-mediated interactions (8), and in spin-orbit-coupled BECs (9) and quantum gases in shaken optical lattices (10) by engineering the single-particle dispersion relation.

About 15 years ago, the existence of a roton minimum was theoretically predicted in BECs with dipole-dipole interactions (DDI) (11); (12). The DDI is long-range and anisotropic; in particular it can change sign depending on the dipole configuration, being attractive for head-to-tail dipoles and repulsive for side-by-side ones (Fig. 1a). Despite the weakly-interacting character of the gas, the roton minimum in dipolar BECs (dBECs) is genuinely interaction-induced as in superfluid He. However, in contrast to helium, the dispersion minimum originates from the peculiar -dependence of the DDI rather than from strong correlations (Fig. 1b).

The realisation of a roton minimum in dBECs would allow for an unprecedented degree of control and microscopic understanding of the roton properties unavailable in helium physics. This prospect has triggered remarkable theoretical works, devoted to unveil the special properties of the dipolar roton and to link them to experimental observables (13); (14); (15); (16); (17); (18); (19); (20); (21); (22); (23); (24); (25); (26); (27). Despite the maturity achieved in the theoretical understanding, the observation of dipolar roton modes has remained so far an elusive goal. For a long time, the only dBEC available in experiments consisted of chromium atoms (28); (29), for which the achievable dipolar character is hardly sufficient to support a roton mode. With the advent of the more magnetic lanthanide atoms (30); (31), a broader range of dipolar parameters became available, providing new prospects for the observation of rotons.

.1 Roton mode in dBECs and its signature

To elucidate the mechanism of rotonization in dBECs, it is instructive to first review the case of cylindrically symmetric pancake traps with the dipoles aligned along the symmetry axis (Fig. 1e1) (11); (13); (15); (16); (17); (18); (19); (20); (21); (22); (23); (24); (25); (26); (27). In this quasi-two-dimensional (q2D) geometry, we consider elementary excitations of planar momentum , corresponding to in-plane density modulations of wavelength (Fig. 1c-d). In a mean-field picture, the energy cost to create such an excitation, , arises from both kinetic and interparticle interaction contributions. The latter includes the contact interaction and the DDI, and the DDI changes sign with because of the much stronger confinement along . At low , the repulsive nature of the in-plane DDI prevails, stiffening the dispersion relation (Fig. 1c). In contrast, for , being the transverse confinement length, the three-dimensional (3D) nature of the excitation is reestablished and the attractive part of the DDI dominates, softening (Fig. 1d). This softening is counterbalanced by the contributions of the repulsive interactions, namely the contact interaction at positive scattering length , and of the kinetic energy, increasing with . For very large , the kinetic energy cost prevails and eventually bends up into a single-particle spectrum. For strong-enough DDI, the dipole-induced softening dominates at intermediate momenta, and develops a minimum at  (11) (Fig. 1b).

Figure 1: Roton mode in a dBEC. a, dipole orientations with repulsive (blue) and attractive (red) DDI. b, illustration of the dispersion relation with real (solid lines) and imaginary (norm of the dotted line) parts of a dBEC in constrained geometries (see e1-f1) and varying (dashed arrow) (11). c-d, spatial density modulations associated to small- (c) and large- (d) excitations. The shape along and the color gradient illustrate the density modulation. The color code, following that of a, pictures the dominant DDI contribution to . In b, a roton emerges at when the attractive DDI prevails. In q2D geometries (e1), the roton in-plane momentum distribution shows a ring (e2). f1 illustrates the advantageous q1D geometry, where the roton momentum distribution is double-peaked (f2). e2 and f2 have the same total population.

Ultracold atoms allow to regulate the impact of the DDI in by tuning the value of through a Feshbach resonance (FR) (32). This provides a powerful control knob for the roton physics that is absent in the He case. The figure of merit is the parameter , defined as the ratio between the dipolar length, , and . Here is the mass and the magnetic moment of the atoms. Increasing (decreasing ) mitigates the energy cost for large- excitations and the system can develop a roton minimum in . By further increasing , the roton gap, , decreases and eventually vanishes (Fig. 1b). In the latter case, the system undergoes a local instability (18); (19) and develops a short-wavelength density modulation on the scale of , , resembling the case of the expected superfluid-to-solid transition in overpressurized He (33).

When the dBEC spectrum gets deeply rotonized (i. e.  small ), the momentum distribution can be profoundly modified, providing a signature for the effect. In q2D geometries (Fig. 1e1), the excitation softening develops in plane - both in the and directions. This leads to a radial roton corresponding to a set of in-plane momenta of . The privileged population of these modes translates into the emergence of a ring of radius in the momentum distribution (Fig. 1e2).

Figure 2: BdG excitation spectrum. a, Excitation spectrum of the ground state of a BEC with atoms in a trap with Hz and scattering length , obtained by numerically solving the BdG equations. Roton modes appear as isolated modes lying below the main branch of the spectrum, forming a ’roton finger’. The lowest roton mode and an exemplary phonon mode are highlighted, with panels b, c showing the corresponding excited state density modulation (with black dashed line the ground state, red solid line the excited state) and, b1, c1, the momentum distribution from  ms ToF expansion (Methods). b, b1 phonon mode, c, c1 roton mode.

Here we extend the roton physics to a largely unexplored geometry: that of axially elongated dBECs with dipoles oriented orthogonally to the elongation axis (Fig. 1f1). Although the above energy arguments equally apply to this quasi-one-dimensional (q1D) geometry, the axial elongation provides an important difference. Here, the system develops a linear roton mode corresponding to a narrow set of with component . As a consequence, the privileged population of this mode translates into two marked peaks in the momentum distribution (Fig. 1f2). This geometrical focusing of the roton population greatly enhances the visibility of the effect compared to the q2D case (Fig. 1e2-f2).

.2 Rotonized dispersion relation in q1D dBECs

Figure 3: Appearance of roton peaks in a q1D geometry. a, illustration of the in-trap Er BEC geometry, forming a cigar shape with -elongation and dipoles aligned by an external uniform -field along the transverse direction. b-d, observed momentum distributions obtained by averaging 15-to-25 absorption images with , ms and for different : b, , c, and d, . b1-d1 show the corresponding central cuts with (red dots) and their fits with three-Gaussian distributions (blue lines), from which we extract and (see text and Methods). e, evolution of with for and ms (dots). Error bars are deduced from the fit statistical uncertainties on the amplitudes . We fit an empirically chosen linear step function to identify (line).

The existence of a roton minimum in our q1D geometry is well explained by a simplified model inspired from the uniform q2D calculations (11). Here we consider that the condensate is trapped along and , with harmonic frequencies and , but homogeneous, i.e. unconfined, along the axis . The physics of the dBEC is well captured by a non-local Gross-Pitaevskii equation (NLGPE) (2); (34), which contains the transverse confinement, the short-range interactions, and the DDI (Methods). Within the Thomas-Fermi (TF) approximation, the BEC density takes the simple form , with the homogeneous axial density (12) (Methods). The excitation spectrum is obtained from the linearisation of the NLGPE around the ground-state wavefunction . Due to the homogeneity along , the elementary excitations have a well-defined momentum . Proceeding as in Ref. (11), we obtain an analytic form for in the relevant case of 3D modes, i. e.  for  (Methods). For dominant DDI (), indeed rotonizes. In the vicinity of the roton minimum the dispersion acquires a gapped quadratic form similar to that of helium rotons:


for . Here the roton momentum is , with and , while its gap reads , with . Note that the TF radii and can be evaluated as functions of , and  (12) and in particular . Then, close to the instability, , takes the simple form with and a geometrical factor that depends only on  (Methods). In the case , .

To move beyond this simplified model and closer to realistic experimental conditions, we develop a numerical approach based on the 3D NLGPE for a generic anisotropic harmonic confinement of frequencies  (Methods) (14). For a quantitative understanding of the experiments, our numerical treatment also includes quantum fluctuations (i. e. local Lee-Huang-Yang (LHY) corrections) (35); (36); (37); (38); (39); (40); (41), finite temperature effects (42); (37), and three-body losses (43) (Methods). Our numerical platform offers rich possibilities to investigate the physics at play: it allows for real-time evolution of the quantum gas wavefunction, and provides access to the Bogoliubov-de Gennes (BdG) excitation spectrum.

The excitation spectrum is obtained by linearizing the NLGPE around a stationary state and numerically solving the resulting BdG equations (Methods) (44). The calculated spectrum is qualitatively modified compared to that of an homogeneous system (Fig. 2). In order to depict the spectrum as a quasi-dispersion relation even in the presence of an axial confinement, we associate to each elementary excitation an effective momentum  (17). The spectrum is discrete with phonon-like collective modes at low . For higher the spectrum flattens, but eventually bends upwards again due to the dominant kinetic energy. In our experimentally relevant q1D geometries, instead of developing a smooth minimum, roton excitations appear as isolated low-lying modes at intermediate momenta that depart from the overall spectrum (17). These so-called roton fingers are related to confinement of the roton modes in the inhomogeneous BEC of profile  (16). Using the local-density approximation in Eq. (1), the spatially-dependent spectrum has a minimal roton gap at the trap center, translating into an effective confinement of the mode.

The confinement is evident from our BdG calculations, in which the lowest roton mode forms a short-wavelength density modulation localized at the trap center (Fig. 2c). This contrasts with phonon modes for which the modulation is delocalised over the entire condensate (Fig. 2b). Whereas the density modulation (Fig. 2c) and the finite-momentum peaks (Fig. 1f2) are signatures of the same effect, from now on we will focus on the latter aspect since it is the one we probe in the experiments via time-of-flight (ToF) expansion measurements (2). Accordingly, we compute the momentum distribution from the ToF expanded densities associated with the selected modes (Fig. 2b1-c1) and observe that the roton mode indeed presents as two symmetric peaks localised at positions corresponding to .

To enrich the stationary-state picture, we additionally develop real-time evolution simulations of the 3D NLGPE. This enables to fully account for the experimental procedure, as described in the next Section, as well as the effects due to finite temperature and atom losses (Methods). Also the real-time evolution shows a population of the roton mode with the same signatures as in Fig. 2c-c1, both in trap and in ToF. The extracted agrees well, within 10%, with the BdG approach in its validity regime, i. e.  when a stationary dBEC state is found.

.3 Roton peaks in quench experiments

Our experiment relies on the strong dipolar character of Er () and on the ability to fine tune below . The experimental sequence starts with a dBEC in a trap elongated along the axis. The trap aspect ratio, , can span from about 4 to 30, corresponding to ranging from 150 Hz to 800 Hz. We note that and are kept roughly constant to Hz and , respectively (Methods). An external homogeneous magnetic field, , defines the dipole orientation (magnetization) with respect to the trap axes and sets the values of through a magnetic FR, centered close to G. The precise -to- conversion has been measured in Ref. (39). At the end of the preparation stage, we obtain a stable BEC at (, G) with transverse () magnetization (Fig. 3a) (Methods).

To access the roton regime, we suddenly quench to the desired lower value and we shortly hold the atoms in the trap for a time , typically ms. We note that during , exponentially converges to its set value with a -time of ms (Methods). We then release the atoms from the trap, change back to approximately its initial value, and let the cloud expand for ms, after which we perform resonant absorption imaging. The imaging beam propagates transversely, i. e.  nearly collinear with the -axis. Hence the ToF images probe the two-dimensional density distribution of the expanded cloud, . For long enough to ignore the in-situ width of the cloud, maps the momentum distribution of the gas in trap, , assuming negligible interactions during the expansion. Our real-time simulations, accounting for the precise experimental sequence, confirm this assumption (Methods).

First we investigate as a function of for , and ms. For large enough , shows a single narrow peak with an inverted aspect ratio compared to the trapped gas, typical of a stable BEC (2) (Fig. 3b). We define the center of the distribution as the origin of . In contrast, when the system enters the dominantly dipolar regime by decreasing , changes fundamentally. We observe a sudden appearance of two symmetric finite-momentum peaks along the axis, located at and of similar shape (Fig. 3c-d). Beside their remarkable symmetry, these finite-momentum peaks also exhibit a high shot-to-shot reproducibility, as evidenced by their high contrast in the averaged distributions (Fig. 3c-d). The observed side peaks show the privileged population of specific high-momentum modes () in the excitation spectrum of our trapped dBEC.

To quantitatively investigate the peak structures, we fit a sum of three Gaussian distributions to the central cuts of the average (Fig. 3(b1-d1)) (Methods). From the fit we extract , the amplitudes of the side and central peaks , and derive the contrast of the side peaks, . Typical values of are -. The evolution of with reveals an interesting behaviour (Fig. 3e). For large , we observe an essentially zero and flat contrast. The residual arises from the usual thermal background and the imaging noise. With decreasing , first exhibits a sharp increase and then saturates to an upper value of about . We define the onset of the finite-momentum excitation, , as the value of at which starts to rise. For the parameter here considered, we find , corresponding to .

Our ToF observations directly reveal the presence of the long-sought roton mode. The observed structures closely follow the predicted signature of this mode (Fig. 2c1). In addition, the fitted are in very good agreement with the predicted from our theory (Sec. .4). Finally, the high shot-to-shot repeatability of the remarkable peak structure in suggests the persistence of a macroscopic phase coherence in the gas.

.4 Properties of the roton excitation.

Figure 4: Geometrical scaling of the roton momentum. Measured at the onset of a roton population in the quench experiment (dots) as a function of . Error bars show the statistical uncertainty of the fit. Here short ranging from 3 to 6 ms are used (Methods). The dotted and dashed lines show the theory predictions using the experimental parameters, from the analytic q1D model of Eq. (1) and from our real-time simulations reproducing the experimental procedure, respectively. Both theories use the values being just below the predicted roton softening (Methods). The blue solid line shows a linear fit passing by the origin, to the data. The fitted slope is 1.61(4).
Figure 5: Roton characteristic diagram. Summary of the measured roton properties as a function of and . Here short ranging from 3 to 6 ms are used (Methods). The white-to-red colormap shows the reduced roton momentum . In the blue region no roton excitation is observed (e. g.  Fig. 3b), and in the grey region we do not have a full set of data (Methods). The yellow triangles show the onset of the roton softening, , extracted from contrast measurements at ms (e. g.  Fig. 3d) and the solid black line is a linear fit to . b, as a function of for the reference trap , corresponding to the dashed line in a. Error bars show the statistical uncertainty of the fit.

A major fingerprint of the roton mode in a dipolar gas is its characteristic dependence on the trap geometry. Its geometrical nature has been proven in various contexts, from pancake- (11); (13); (15); (16); (17) to our cigar-shaped traps (Sec. .2), highlighting a universal scaling with . To experimentally investigate such geometrical properties in our elongated dBEC, we repeat the quench measurements for different trap parameters (, ) and we extract the corresponding value of (Fig. 4) (Methods). Our data, plotted as a function of , reveal a marked increase of , and a linear fit to the data gives a slope of 1.61(4). This value is remarkably close to the mean calculated from our simple analytic model using the experimental parameters. To firmly tie our observation to the roton excitation, we perform real-time simulations of the NLGPE, accounting for the specific experimental sequences (Methods). The extracted parameter-free curve (dashed line in Fig. 4) is in excellent agreement with the data and also captures the smooth deviation from a linear slope reflecting the experimental details. The comparison between the experimental observations and our two complementary theories unambiguously establishes the rotonic nature of the finite-momentum peaks.

To gain further insights into the roton properties, we additionally investigate the effects of and . From all our measurements, we construct a characteristic diagram, presented in terms of the dimensionless parameters and (Fig. 5a). The former parameter governs the geometric competition between the attractive and repulsive parts of the DDI, and embodies the interplay between dipolar and short-range interactions (Methods). The white-to-red color scale indicates the values of the reduced roton momentum . Each column of the diagram shows the variation of with , whereas each row gives the evolution of with for a given trap configuration. We observe varying at most by along the lines and rows, which is on the order of our experimental precision. Close to the instability, both dependencies are expected to be mild as remains mainly set by its geometrical nature (see Secs. .1,.2 and e. g. Refs.(11); (15)). The smooth dependence of with is particularly relevant (Fig. 5b) as it contrasts the observed softening from a phonon instability in which the minimum in is expected to appear at and then move to larger values when lowering . From our analytic q1D model for the roton softening, we extract a residual increase of of 17% from to in the geometry  Hz, similar to the experiments.

Besides the geometric determination of , and control the roton gap . In the experiment, we observe that the critical for the onset of the roton population moves to larger values for more elongated traps (blue-to-red border and -measurements in Fig. 5). This increase can be explained by the following argument: a macroscopic roton population appears for small enough  (Methods). This condition is realised when the attractive DDI dominates over the total repulsive (contact and repulsive DDI) interaction. As the ratio of the attractive and repulsive contributions of the DDI depends on , so will . A larger favors the repulsive part of the DDI, thus weaker contact interactions are needed for a similar roton gap, explaining a lower . Predictions from our q1D model reproduce well the observed , both its variation direction and the actual values within  (Methods).

.5 Discussion

In conclusion, we report on the observation of the long-sought roton mode in a dipolar BEC of Er atoms. Our investigations take advantage of an axially elongated geometry, amplifying the roton signature in momentum space. Our results, combining experimental studies, numerical and analytic theories, unambiguously establish the roton softening of the excitation spectrum and the universal geometrical scaling of the excitation.

Our work opens fascinating new ground for the study of roton physics in dipolar gases. In the future, it might be interesting to explore out-of-equilibrium dynamics of the roton mode, to spectroscopically access the roton gap, and to directly probe in-trap density modulated states. These density modulations, providing a signature of the roton softening in real space, are a precursor of a supersolid phase, in which a phase-coherent density-modulated ground-state would arise (45). Although the density modulation is expected to be mean-field unstable against local collapses (46), quantum stabilization may prevent collapse as for the case of recently explored quantum droplets (38); (39); (37); (40); (41). Future investigations may be directed in exploring, both in theory and experiment, the possibility of a stable supersolid state in dBECs.

  • Acknoledgements: We thank M. Baranov and E. Demler for fruitful discussions and G. Natale for his support in the final stage of the experiment. This work is dedicated in memory of D. Jin and her inspiring example. The Innsbruck group is supported through a ERC Consolidator Grant (RARE, no. 681432) and a FET Proactive project (RySQ, no. 640378) of the EU H2020. LC is supported within a Marie Curie Project (DipPhase, no. 706809) of the EU H2020. FW and LS thank the DFG (SFB 1227 DQ-mat). All the authors thank the DFG/FWF (FOR 2247). Part of the computational results presented have been achieved using the HPC infrastructure LEO of the University of Innsbruck.

  • Competing Interests: The authors declare that they have no competing financial interests.

  • Correspondence: Correspondence and requests for materials should be addressed to F.F. (email: francesca.ferlaino@uibk.ac.at).

  • Author contributions: FF, LC, DP, GF, MJM, JHB, SB conceived and supervised the experiment and took the experimental data. LC analysed them. RMWvB developed the BdG calculations. FW and LS performed the real time simulations. LS derived the q1D analytical model. LC, FF, RMWvB and LS wrote the paper with the contributions of all the authors.


.6 Production of BECs

We prepare a BEC similarly to Refs. (39); (31). From a narrow-line magneto-optical trap with Er  atoms, automatically spin-polarized in their absolute lowest Zeeman sub-level  (47), at about , we directly load the atomic gas in a crossed optical dipole trap (ODT) with an efficiency of more than . A uniform magnetic field, , is permanently applied along the vertical axis, fixing the dipole orientation, while its value is varied during the experimental sequence, to tune (Method .8). We achieve condensation by means of evaporative cooling in the crossed ODT at G (). During the evaporation procedure, we first change the power and then the ellipticity of one of the ODT beams (Method .7). The final atom number, typically , and condensed fraction, typically , are assessed by fitting the ToF absorption images of the gas to a bimodal distribution, sum of a TF profile and a broad Gaussian background.

.7 Trapping geometries

The BEC is confined in a harmonic trapping potential , characterized by the frequencies . The ODT results from the crossing of two red-detuned laser beams of nm wavelength at their respective focii. One beam propagates along the -axis and is denoted hODTb; the other, called vODTb, propagates nearly collinear to the -axis. The vODTb has a maximum power of 7 W and an elliptical profile with waists of and along and respectively. The hODTb has a maximum power of 24 W, a vertical waist , and a controllable horizontal waist, . The ellipticity can be tuned from to by time averaging the frequency of the first-order deflection of an Acousto-Optic Modulator (48). By adjusting independently and the powers of the vODTb and of the hODTb , we can widely control the geometry of the trap in a dynamic manner. Precisely, is essentially set by the vODTb power, by that of the hODTb. is controlled by both the power and ellipticity of the hODTb, with . The independent control of and yields an especially easy tuning of the relevant trap aspect ratio, , for our q1D geometry.

We use the tuning of the hODTb power and ellipticity to perform evaporative cooling to quantum degeneracy (Method .6). After reaching condensation, we again modify the beam parameters to shape the trap into a favourable q1D configuration for observing the roton physics (). The different trapping geometries probed in Figs. 3-5 are achieved by changing the hODTb power with and the vODTb power set to its maximum so that and are kept roughly constant. The corresponding are calibrated via the excitation and probing of the center-of-mass oscillation of thermal samples and reported in Extended Data Table 1. We note that the final atom number , BEC fraction , and temperature after the shaping procedure depend on the final configurations, as detailed in Extended Data Table 1. is extracted from the evolution with the ToF duration of the size of the background Gaussian in the TF-plus-Gaussian bimodal fit to the corresponding ToF images of gas. The values of and are used for the initial states of our real-time simulations (Method .14).

.8 Quench of the scattering length

To control we use a magnetic Feshbach resonance between atoms in their absolute ground state, which is centered around G. The -to- conversion has been previously precisely measured via lattice spectroscopy, as reported in Ref. (39). Errors on , taking into account statistical uncertainties of the conversion and effects of magnetic field fluctuations (e.g. from stray fields), are of 3-to-5  for the relevant range - in this work. After the BEC preparation and in order to investigate the roton physics via an interaction quench, we suddenly change the magnetic field set value, , twice. First we perform the quench itself and abruptly change from G () to the desired lower value at the beginning of the hold in trap (ms). Second we prepare the ToF expansion and imaging conditions (see Method .9) and abruptly change from the quenched value to G () at the beginning of the ToF expansion (ms). Due to delays in the experimental setup, e.g. coming from eddy currents in our main chamber, the actual value felt by the atoms responds to a change of via (49). By performing pulsed-radio-frequency spectroscopy measurements (pulse duration 100 s) on a BEC after changing (from 0.4 to 0.2 G), we verify this law and calibrate ms. As a consequence is evolving during and . This effect is fully accounted in the experiment and simulations and, we report the roton properties as a function of the effective value of at the end of . We use ranging from to  ms. The lower bound on comes from the time needed for to effectively reach the regime of interest. Here we then consider the initial evolution for which , being the characteristic collision rate. We estimate that ranges typically from 40 to 90 ms in the initial BECs of Extended Data Table 1 at . We note that during the considered , the atom loss remains below for our less confined geometries Hz and can go up in the tightest traps of Hz. We have checked that the roton, if it exists, has developed within this range of . In Figs. 4-5, the parameter is optimized for each value to obtain the largest visibility of the roton peaks (largest ).

.9 Imaging procedure

The in-trap density modulation associated with the roton excitation has a characteristic wavelength of (e.g. Fig. 2c). This value is much smaller than the axial width of the cloud () and below our imaging resolution (). In our experiments, we employ ToF expansion measurements, accessing the momentum distribution of the gas (2), to probe the roton mode population. We let the gas expand freely for ms, which translates the imaging resolution in space into a momentum resolution of . This means that we can resolve the population of modulation modes with wavelength and the roton mode of interest should be well detectable.

In the experiments, we record 2D absorption pictures of the cloud after ToF expansion by means of standard resonant absorption imaging on the atomic transition at 401 nm. The imaging beam propagates nearly vertically, with a remaining angle of compared to the -axis within the -plane. Thus the ToF images essentially probe the spatial density distribution in the plane. When releasing the cloud (ms), we change back to G (Method .8). This change enables constant and optimal imaging conditions with a fixed probing procedure, i. e. a maximal absorption cross-section. In addition, the associated increase of to allows to minimize the time during which the evolution happens in the small- regime where roton physics develops. Thanks to this we are able to effectively probe the short-time evolution of the gas. In this work, we use the simple mapping:


which neglects the initial size of the cloud in the trap and the effect of interparticle interactions during the free expansion. Using real-time simulations (Method .14), we simulate the experimental sequence and are able to compute both the real momentum density from the in-trap wavefunction and the spatial ToF distribution ms after switching off the trap. Using the mapping of Eq. (2) and our experimental parameters, the two calculated distributions are very similar, and, in particular, the two extracted momenta associated with the roton signal agrees within . This confirms that the interparticle interactions play little role during the expansion and justifies the use of Eq. (2).

.10 Fit procedure for the ToF images

For each data point of Figs. 3-5, we record between 12 and 25 ToF images. By fitting a two-dimensional Gaussian distribution to the individual images, we extract their origin and recenter each image. From the recentered images, we compute the averaged , from which we characterise the linear roton developing along . To do so, we extract a one-dimensional profile by averaging the one-dimensional cut of of fixed within : . To quantitatively analyse the observed roton peaks, we fit a sum of three Gaussian distributions to . One Gaussian accounts for the central peak and its centre is imposed to . The two other Gaussians are symmetrically located at , and we impose their size and amplitude to be identical. We focus on the roton side peaks by constraining and . In all the figures, the statistical error of the fit is characterised by its 95 % confidence interval.

In order to analyze the evolution of the contrast of the side peaks (see Fig. 3e and Fig. 5), we perform a second run of the fitting procedure, in which we constrain more strictly the value of . The interval of allowed values is defined for each trapping geometry (Extended Data Table 1) and is set from the results of the first run of the fitting procedure. We use the results of the (, )-configurations where the peaks are clearly visible (see Fig.5) and we set the allowed -range to that covered by the 95% confidence intervals of the first-fitted in these configurations. This constraint enables that, for , the fitting procedure estimates the residual background population on the relevant momentum range for the roton physics.

In the ToF images, we qualitatively observe that the momentum distribution broadens (i. e. larger ) when increases and decreases. This behaviour effectively limits the visibility of the roton peaks, and explains that the optimal for the roton observation (see Method .8) gets shorter for lower . In Fig. 5, the grey region indicates the (,)-configurations where the whole momentum distribution was observed to be too spread out even at the shortest ms, so that the roton peaks could not be clearly detected. This region was then excluded and no full set of measurements is available. We note that, in this region of the diagram, is the smallest, making the blurring effect due to the distribution braodening more drastic, since the roton and BEC peaks are closer together. The non-detectability of roton peaks can also relate to the phase diagram characteristics. The small--small- region, corresponding to the grey area lies deeper in the roton regime (as embodied by the increase of for decreasing , see main text). It is then possible that a larger number of modes are dynamically destabilized and populated, resulting in a blurring of the roton signal.

.11 Generalized Gross–Pitaevskii equation

Our theory is based on an extended version of the NLGPE


governing the evolution of a macroscopically occupied wavefunction , with corresponding atomic density at position and time . The standard dipolar NLGPE includes the kinetic energy, external trap potential and the mean-field effect of the interactions. These correspond to the three first terms of Eq. (.11), where the mean-field interaction potential takes the form of a convolution of with the binary interaction potential


for two particles separated by  (34). The first term corresponds to contact interactions between the particles with strength . The DDI gives rise to the second term, which depends on both the distance and orientation of the vector compared to the polarisation axis ( axis) of the dipoles (angle ). Most properties of dBECs are well captured by this standard NLGPE (mean-field) (2); (34).

Recent experimental and theoretical results, however, have established the importance of accounting for quantum fluctuations in dipolar condensates (38); (39); (37); (40); (41). Their effect can be included in the NLGPE in a mean field treatment through a Lee-Huang-Yang correction to the chemical potential, , which is obtained under a local density approximation (36); (35). The accuracy of this mean field treatment has been established, e.g., in Refs. (37); (40); (41), and has proven succesful in explaining recent experimental results (38); (39). The final nonlinear term in the extended NLGPE accounts for three-body losses, with an experimentally determined loss parameter , which is dependent on and typically of the order , as reported in Ref. (39).

.12 BdG spectrum calculations

Collective excitations of the dBEC are obtained by linearising the NLGPE around a stationary state , which can be obtained by imaginary time propagation (Method .14). We write , where is the chemical potential associated with state , and are spatial modes oscillating in time with characteristic frequency and (44). Inserting this ansatz in the NLGPE, and retaining only terms up to linear order in we obtain the BdG equations


where the operator , acting on a function and evaluated at point , is defined as


The above equations constitute an eigenvalue problem, which we solve numerically using the Arnoldi method to obtain eigenmodes and corresponding excitation energies . The equations presented here are a generalisation of the BdG equations for dipolar systems as derived in Ref. (44), to include the LHY correction for quantum fluctuations. The LHY term generally serves to stabilise the excitation spectrum, causing the roton instability point to shift to lower scattering lengths. The excited states shown in Fig. 2b-c correspond to the density , for particular pairs of corresponding to phonon and roton modes. Even while the amplitude of both modes is taken to be equal in Figs. 2b and 2c, the roton excitation (Fig. 2c) leads to markedly larger local density modulations than the phonon excitation (Fig. 2b). The ToF signatures in Fig. 2 are computed by letting the wave function of the excitation, , expand ballistically for , i.e. neglecting interactions during the expansion. The resulting density is then plotted (Fig. 2b1, c1).

.13 Analytical dispersion relation for an infinite axially elongated geometry

Equation (1) in the main text results from a similar procedure as that used in Ref. (11) for rotons in infinite q2D traps. We consider a dBEC homogeneous along but harmonically confined with frequencies and along and . For sufficiently strong interactions the BEC is in the TF regime on the plane, in which the BEC wavefunction acquires the form , with , where and are the TF radii, and . The calculation of , and is detailed at the end of this section.

Due to the axial homogeneity, the elementary excitations discussed in the previous section have a defined axial momentum , and take the form . We consider the GPE given by Eq. (.11) without LHY correction and three-body losses, and insert the perturbed solution . After linearization we obtain the BdG equations for :




Employing , and for , we obtain the following equation for the function :


where , , . For the last term of Eq. (12) vanishes. In that case, the lowest-energy solution is given by , whose eigen-energy builds, as a function of , the dispersion with


with . In the vicinity of , the effect of the last term in Eq. (12) may be evaluated perturbatively, resulting in the dispersion


with .

This expression for the dispersion presents a roton minimum for at . Expanding Eq. (14) in the vicinity of the roton minimum, , we obtain Eq. (1) of the main text, with . At the instability, , and .

Employing a similar procedure as in Ref. (12) we obtain that the BEC aspect ratio fulfills:


with and


These two equations fully determine the TF solution for given , , and the ratio . By inserting the expressions of and in , we find for :


whereas simplifies into . As a result, at the instability depends only on the transverse confinement aspect ratio , and we obtain the geometrical factor introduced in the main text:


In order to better compare quantitatively with our quench experiments, we evaluate the 3D TF solution, for the axially trapped condensate prior to the quench. We then use the central density and the TF radii and in the evaluation of the roton spectrum above, employing the value of after the quench. This procedure takes the quench of into account by evaluating an initial (i. e.  at ms) instantaneous roton spectrum, within the assumption of an immediate change of . Further refinements might be required to fully describe the dynamical population of the roton mode, related to the roton softening of the time-dependent instantaneous excitation spectrum depicted in Method .15. Yet, as being inspired from the latter concept, our approach allows for a better description of our experiments than purely stationary calculations.

The evaluation of the 3D TF solution is performed by following the procedure of Ref. (50). We introduce the BEC aspect ratios: , and . For a given , these aspect ratios fulfil the equations:


where , , and . Here we have introduced the functions:


with , and . Once are determined one may then evaluate the TF radii:


and the central density