Observation of Roton Mode Population in a Dipolar Quantum Gas

TL;DR

This paper reports the first experimental observation of a roton mode in a dipolar quantum gas of erbium atoms, confirming theoretical predictions and advancing understanding of excitations in quantum liquids.

Contribution

It provides the first direct evidence of roton modes in dipolar quantum gases through experimental momentum distribution measurements.

Findings

Roton peaks appear at predicted momenta after interaction quench.

Roton momentum scales inversely with confinement length.

Roton softening and imaginary gap are observed over time.

Abstract

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 $^4$He. In quantum liquids, rotons arise from the 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 magnetic dipole-dipole interactions despite of their weakly-interacting character. This prospect has raised considerable interest; yet roton modes in dipolar quantum gases have remained elusive to observations. Here we report experimental and theoretical studies of the momentum distribution in Bose-Einstein condensates of highly-magnetic erbium atoms, revealing the existence of the long-sought roton mode. By quenching the interactions, we observe the roton appearance of peaks at well-defined momentum. The roton momentum follows the predicted geometrical scaling with the inverse of the confinement length along the magnetisation axis. From the growth of the roton population, we probe the roton softening of the excitation spectrum in time and extract the corresponding imaginary roton gap. Our results provide a further step in the quest towards supersolidity in dipolar quantum gases.

Tables

Table 1. Extended Data Table 1: dBEC parameters for the experimental measurements (Figs. 2-4). The typical statistical uncertainties on ν x subscript 𝜈 𝑥 \nu_{x} and ν z subscript 𝜈 𝑧 \nu_{z} are below 1 % percent 1 1\,\% , and can be up to 10 % percent 10 10\,\% for ν y subscript 𝜈 𝑦 \nu_{y} . The experimental repeatability results in 5 5 5 -to- 10 % percent 10 10\,\% shot-to-shot fluctuations of N 𝑁 N , f 𝑓 f and T 𝑇 T .

$\lambda$	${\nu }_x$ ($\mathrm{Hz}$ )	${\nu }_y$ ($\mathrm{Hz}$ )	${\nu }_z$ ($\mathrm{Hz}$ )	$N$ (${10}^4$)	$f$ (%)	$T$ (nK)
$4.3$	$114$	$35$	$149$	$9$	$66$	$45$
$4.9$	$154$	$40$	$195$	$10$	$65$	$50$
$10.2$	$183$	$30$	$306$	$11$	$62$	$104$
$14.3$	$267$	$32$	$456$	$8.6$	$50$	$150$
$21.3$	$357$	$30$	$638$	$8.4$	$36$	$179$
$29.7$	$432$	$26$	$771$	$7$	$20$	$171$

Table 2. Extended Data Table 2: Roton softening threshold for the geometries realized in the experiment (Table 1 ), deduced from the 𝒞 𝒞 \mathcal{C} fits for the experimental onset of the population ( a s ∗ , c superscript subscript 𝑎 s c a_{\rm s}^{*,\rm c} ) and from imaginary time NS ( a s ∗ , st superscript subscript 𝑎 s st a_{\rm s}^{*,\rm st} ). The fit uncertainties on the experimental values are typically of 0.2 a 0 subscript 𝑎 0 a_{0} , except for the last trap it is of 0.8  a 0 subscript 𝑎 0 a_{0} .

$\lambda$	${\nu }_z$ ($\mathrm{Hz}$ )	$a_s^{\ast ,c}$ ($a_0$)	$a_s^{\ast ,\mathrm{st}}$ ($a_0$)
$4.3$	$149$	$51.4$	$50$
$10.2$	$306$	$50.2$	$46$
$14.3$	$456$	$49$	$43$
$21.3$	$638$	$47.7$	$39$
$29.7$	$771$	$45.2$	$30$

Equations

$$ϵ(k_y{)}^2\simeq {\Delta }^2+\frac{2{\hslash }^2{k}_{\mathrm{rot}}^2}{m}\frac{{\hslash }^2}{2m}(k_y-k_{\mathrm{rot}}{)}^2.$$

$$\tilde{n}(k_x,k_y)={\left(\frac{\hslash t_{\mathrm{f}}}{m}\right)}^2{n}_{\mathrm{tof}}\left(\frac{\hslash k_x{t}_{\mathrm{f}}}{m},\frac{\hslash k_y{t}_{\mathrm{f}}}{m},t_{\mathrm{f}}\right),$$

$$A_{\ast }(a_{\mathrm{s}}^{\mathrm{f}},t_{\mathrm{h}})\propto \exp \left(2{\int }_0^{t_{\mathrm{h}}}dt\mathrm{Im}[\bar{\Delta }(a_{\mathrm{s}}^{\mathrm{f}},t)]/\hslash \right),$$

$$A_{\ast }(a_{\mathrm{s}}^{\mathrm{f}},t_{\mathrm{h}})\propto \exp \left(2\Gamma {\int }_0^{t_{\mathrm{h}}}\delta (a_{\mathrm{s}}(t))dt\right).$$

$$ϵf_{-}(\mathbf{\rho })$$

$$ϵf_{+}(\mathbf{\rho })$$

$$H_{\mathrm{kin}}{f}_{\pm }(\mathbf{\rho })$$

$$H_{\mathrm{int}}[f_{-}(\mathbf{\rho })]$$

$$0$$

$${ϵ}_0(k_y{)}^2=E(k_y{)}^2+E_0^2,$$

$$ϵ(k_y{)}^2\simeq {ϵ}_0(k_y{)}^2-2E_{I}E(k_y),$$

$${\chi }^2\left[\frac{(1-{ϵ}_{dd})(1+\chi {)}^2+3{ϵ}_{dd}}{(1+2{ϵ}_{dd})(1+\chi {)}^2-3{ϵ}_{dd}{\chi }^2}\right]={\lambda }_{\perp }^2,$$

$$Z^2=\frac{g{n}_0}{2{\pi }^{2}m{\nu }_z^2}\left[(1+2{ϵ}_{dd})-\frac{3{ϵ}_{dd}{\chi }^2}{(1+\chi {)}^2}\right]$$

$$E_0^2=\frac{h^2{\nu }_z^2}{6}(1+\chi {)}^2\left({\lambda }_{\perp }^2+\frac{1}{1+2\chi }\right)$$

$$\kappa =k_{\mathrm{rot}}{\ell }_z={\left(\frac{2}{3}\right)}^{1/4}\sqrt{1+\chi }{\left({\lambda }_{\perp }^2+\frac{1}{1+2\chi }\right)}^{1/4}.$$

$$n(\mathbf{r},t)=n_0(t)\left[1-{\left(\frac{x}{X(t)}\right)}^2-{\left(\frac{y}{Y(t)}\right)}^2-{\left(\frac{z}{Z(t)}\right)}^2\right],$$

$$\frac{{ϵ}^2(k_y,y,t)}{(\hslash {\omega }_z{)}^2{l}_z^4}$$

$${\tilde{n}}_1(k_y,t)\sim \int dy\left[e^{2{\int }_0^{t}d{t}^{\prime }\mathrm{Im}\left[ϵ(k_y,y,t^{\prime })\right]}-1\right].$$

$$\mathrm{i}\hslash \frac{\partial \psi (\mathbf{r},t)}{\partial t}$$

$$\displaystyle\hskip 56.9055pt+\Delta\mu[n]-{\mathrm{i}}\hbar\frac{L_{3}}{2}n^{2}\Big{)}\psi(\mathbf{r},t)$$

$$\displaystyle\equiv\Big{(}\hat{H}_{\mathrm{GP}}[\psi]-\mathrm{i}\hbar\frac{L_{3}}{2}n^{2}\Big{)}\psi(\mathbf{r},t),$$

$$U(\mathbf{r})=g\left(\delta (\mathbf{r})+\frac{3{ϵ}_{\mathrm{dd}}}{4\pi }\frac{1-3{\cos }^2\theta }{|\mathbf{r}{|}^3}\right),$$

$$\left(\begin{array}[]{cc}\hat{H}_{\mathrm{GP}}[\psi_{0}]+A&-A\\ A&-\hat{H}_{\mathrm{GP}}[\psi_{0}]-A\end{array}\right)\left(\begin{array}[]{c}u\\ v\end{array}\right)=\epsilon\left(\begin{array}[]{c}u\\ v\end{array}\right),$$

$$\begin{gathered} (Af)(\mathbf{r})=\int \mathrm{d}{\mathbf{r}}^{\prime }{\psi }_0({\mathbf{r}}^{\prime })U(\mathbf{r}-{\mathbf{r}}^{\prime })f({\mathbf{r}}^{\prime }){\psi }_0(\mathbf{r}) \\ +\frac{16}{\sqrt{\pi }}g{a}_{\mathrm{s}}^{3/2}\left(1+\frac{3}{2}{ϵ}_{\mathrm{dd}}^2\right)|{\psi }_0(\mathbf{r}){|}^{3}f(\mathbf{r}). \end{gathered}$$

$$S(\mathbf{k},{\omega }^{\prime })=\sum _{\mathbf{j}}{\left|\int \mathrm{d}\mathbf{r}[{\mathbf{u}}^{\ast }(\mathbf{r})+{\mathbf{v}}^{\ast }(\mathbf{r})]{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{r}}{\psi }_0(\mathbf{r})\right|}^2\delta ({\omega }^{\prime }-\omega ),$$

$$n(\mathbf{r},t)=n_0(t)\left[1-{\left(\frac{x}{X(t)}\right)}^2-{\left(\frac{y}{Y(t)}\right)}^2-{\left(\frac{z}{Z(t)}\right)}^2\right],$$

$$\frac{m}{2}\nabla \left[\frac{\ddot{X}}{X}{x}^2+\frac{\ddot{Y}}{y}{y}^2+\frac{\ddot{Z}}{Z}{z}^2\right]=-\nabla \mu (\mathbf{r},t),$$

$$\mu (\mathbf{r},t)=V(\mathbf{r})+g{n}_0\left[D-A{\left(\frac{x}{X}\right)}^2-B{\left(\frac{y}{Y}\right)}^2-C{\left(\frac{z}{Z}\right)}^2\right],$$

$$F_1({\chi }_x,{\chi }_y)$$

$$F_2({\chi }_x,{\chi }_y)$$

$$F_3({\chi }_x,{\chi }_y)$$

$$\frac{1}{(2\pi {)}^2}\frac{{\ddot{b}}_x}{{\nu }_x^2}$$

$$\frac{1}{(2\pi {)}^2}\frac{{\ddot{b}}_y}{{\nu }_y^2}$$

$$\frac{1}{(2\pi {)}^2}\frac{{\ddot{b}}_z}{{\nu }_z^2}$$

Full text

\DeclareCaptionType

extfig

Observation of Roton Mode Population in a Dipolar Quantum Gas

L. Chomaz1, R. M. W. van Bijnen2, D. Petter1, G. Faraoni1,3, S. Baier1, J. H. Becher1, M. J. Mark1,2, F. Wächtler4, L. Santos4, F. Ferlaino1,2,∗

1Institut für Experimentalphysik,Universität Innsbruck, Technikerstraße 25, 6020 Innsbruck, Austria

2Institut für Quantenoptik und Quanteninformation,Österreichische Akademie der Wissenschaften, 6020 Innsbruck, Austria

3Dipartimento di Fisica e Astronomia, Università di Firenze, Via Sansone 1, 50019 Sesto Fiorentino, Italy

4Institut für Theoretische Physik, Leibniz Universität Hannover, Appelstr. 2, 30167 Hannover, Germany

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 4He Landau41 . In quantum liquids, rotons arise from the strong interparticle interactions, whose microscopic description remains debated Griffin:1993 . In the realm of highly-controllable quantum gases, a roton mode has been predicted to emerge due to magnetic dipole-dipole interactions despite of their weakly-interacting character Santos:2003 . This prospect has raised considerable interest Ronen:2007 ; Bohn:2009 ; Parker:2009 ; Martin:2012 ; Blakie:2012 ; Lasinio:2013 ; Wilson:2010 ; Natu:2014 ; Pitaevskii:2016 ; yet roton modes in dipolar quantum gases have remained elusive to observations. Here we report experimental and theoretical studies of the momentum distribution in Bose-Einstein condensates of highly-magnetic erbium atoms, revealing the existence of the long-sought roton mode. By quenching the interactions, we observe the roton appearance of peaks at well-defined momentum. The roton momentum follows the predicted geometrical scaling with the inverse of the confinement length along the magnetisation axis. From the growth of the roton population, we probe the roton softening of the excitation spectrum in time and extract the corresponding imaginary roton gap. Our results provide a further step in the quest towards supersolidity in dipolar quantum gases Boninsegni:2012 .

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 Griffin:1993 ; Pitaevskii:2016 . A major breakthrough in understanding superfluidity thrived on the concept of quasiparticles, introduced by Landau in 1941 Landau41 . Quasiparticles are elementary excitations of momentum $k$, whose energies $\epsilon$ define the dispersion (energy-momentum) relation $\epsilon(k)$.

To explain the special thermodynamic properties of superfluid ${}^{4}\rm{He}$, Landau postulated the existence of two types of low-energy quasiparticles: phonons, referring to low-$k$ acoustic waves, and rotons, gapped excitations at finite $k$ initially interpreted as elementary vortices. The dispersion relation continuously evolves from linear at low $k$ (phonons) to parabolic-like with a minimum (roton) at a finite $k=k_{\rm rot}$. Neutron scattering experiments confirmed Landau‘s remarkable intuition Henshaw61 . In liquid 4He, $k_{\rm rot}$ 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 Griffin:1993 .

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 typically weak, meaning that classically their range of action is much smaller than the mean interparticle distance. Because of this diluteness, roton excitations are absent in ordinary quantum gases, i. e. in Bose-Einstein condensates (BECs) with contact (short-range) interactions Pitaevskii:2016 . However, about 15 years ago, seminal theoretical works predicted the existence of a roton minimum both in BECs with magnetic dipole-dipole interactions (DDIs) Santos:2003 and in BECs irradiated by off-resonant laser light Odell2003 . Following the lines of the latter proposal, a roton softening has been recently observed in BECs coupled to an optical cavity Mottl12 . Here, the excitation arises from the infinite-range photon-mediated interactions and the inverse of the laser wavelength sets the value of $k_{\rm rot}$. Additionally, roton-like softening has been created in spin-orbit-coupled BECs Khamehchi:2014 ; Ji2015 and quantum gases in shaken optical lattices LiChung15 by engineering the single-particle dispersion relation.

Our work focuses on dipolar BECs (dBECs). As in superfluid 4He, the roton spectrum in such systems is a genuine consequence of the underlying interactions among the particles. However, in contrast to helium, the emergence of a minimum at finite momentum does not require strong inter-particle interactions. It instead exists in the weakly-interacting regime and originates from the peculiar anisotropic and long-range character of the DDI in real and momentum space Santos:2003 ; Ronen:2007 ; Bohn:2009 ; Parker:2009 ; Martin:2012 ; Blakie:2012 ; Lasinio:2013 ; Wilson:2010 ; Natu:2014 ; Pitaevskii:2016 . 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 Griesmaier:2005 , for which the achievable dipolar character is hardly sufficient to support a roton mode. With the advent of the more magnetic lanthanide atoms Lu:2011 ; Aikawa:2012 , a broader range of dipolar parameters became available, opening the way to access the regime of dominant DDI. In this regime, novel exciting many-body phenomena have been recently observed, as the formation of droplet states stabilised by quantum fluctuations Kadau:2016 ; Igor:2016 ; Chomaz:2016 , which may become self-bound Schmitt:2016 . Lanthanide dBECs hence open new roads toward the long-sought observation of roton modes.

Prior to this work, dipolar rotons have been mostly connected to pancake-like geometries Santos:2003 ; Ronen:2007 ; Bohn:2009 ; Parker:2009 ; Blakie:2012 ; Lasinio:2013 ; Wilson:2010 ; Natu:2014 ; Pitaevskii:2016 . Here, we extend the study of roton physics to the case of a cigar-like geometry with trap elongation along only one direction ($y$) transverse to the magnetisation axis ($z$) (Fig. 1a). The anisotropic character of the DDI (Fig. 1b, inset) together with the tighter confinement along $z$ is responsible for the rotonization of the excitation spectrum along $y$ (Fig. 1b). To illustrate this phenomenon, we consider an infinite cigar-shaped dBEC and focus on its axial elementary excitations, of momentum $k_{y}$. These excitations correspond in real space to a density modulation along $y$ of wavelength $2\pi/k_{y}$. For low $k_{y}$, the atoms sit mainly side-by-side and the repulsive nature of the DDI prevails, stiffening the phononic part of the dispersion relation (Fig. 1b1). In contrast, for $k_{y}\ell_{z}\gtrsim 1$, $\ell_{z}$ being the characteristic $z$ confinement length, the excitation favours head-to-tail alignments and the DDI contribution to $\epsilon(k_{y})$ eventually changes sign Santos:2003 (Fig. 1b2). The resulting softening of $\epsilon(k_{y})$ is counterbalanced by the contributions of the repulsive contact interaction, and of the kinetic energy, which ultimately dominates at very large $k_{y}$. For strong enough DDI, this competition gives rise to a roton minimum in $\epsilon(k_{y})$, occurring at a momentum $k_{y}=k_{\rm rot}$ set by the geometrical scaling $k_{\rm rot}\sim 1/\ell_{z}$ (see below and e.g. Santos:2003 ; Blakie:2012 ; Lasinio:2013 ).

Similar to the helium case, the roton energy gap, $\Delta=\epsilon(k_{\rm rot})$, depends on the density and on the strength of the interactions. In ultracold gases, both quantities can be controlled. In particular, the scattering length $a_{\rm s}$, setting the strength of the contact interaction, can be tuned using Feshbach resonances Pitaevskii:2016 . As $a_{\rm s}$ is reduced, $\Delta$ decreases, vanishes, and eventually becomes imaginary (Fig. 1b). In the latter case, the system undergoes a roton instability and the population at $k_{y}=0$ is transferred in $\pm k_{\rm rot}$ at an exponential rate Bohn:2009 ; Parker:2009 . The population of the roton mode is then readily visible in the momentum distribution of the gas (Fig. 1c-d). In the extensively-studied pancake geometries, the roton population in $k$-space spreads over a ring of radius $k=k_{\rm rot}$ because of the radial symmetry of the confinement (Fig. 1d). Such a spread can be avoided using a cigar geometry. Here, the roton population focuses in two prominent peaks at $k_{y}=\pm k_{\rm rot}$, enhancing the visibility of the effect (Fig. 1c).

We explore the above-described physics using strongly magnetic ${}^{166}{\rm Er}$ atoms. The experiment starts with a stable dBEC in a cigar-shaped harmonic trap of frequencies $\nu_{x,y,z}$, elongated along the $y$ axis. The trap aspect ratio, $\lambda=\nu_{z}/\nu_{y}$, can be tuned from about 4 to 30, corresponding to $\nu_{z}$ ranging from 150 Hz to 800 Hz, whereas $\nu_{y}$ and $\nu_{z}/\nu_{x}$ are kept constant at about $35\,$Hz and $1.6$, respectively (Methods). An external homogeneous magnetic field, $B$, fixes the dipole orientation (magnetisation) with respect to the trap axes and sets the values of $a_{\rm s}$ through a magnetic Feshbach resonance, centered close to $B=0\,$G. In previous experiments, we precisely calibrated the $B$-to-$a_{\rm s}$ conversion for this resonance Chomaz:2016 . The BEC is prepared at $a_{\rm s}^{\rm i}=61\,a_{0}$ ($B=0.4\,$G) with transverse ($z$) magnetisation. The characteristic dipolar length, defined as $a_{\rm dd}=\mu_{0}\mu_{m}^{2}m/12\pi\hbar^{2}$, is $65.5\,a_{0}$, with $m$ the mass and $\mu_{m}$ the magnetic moment of the atoms, $\hbar=h/2\pi$ is the reduced Planck constant, $\mu_{0}$ the vacuum permeability and $a_{0}$ the Bohr radius.

To excite the roton mode, we quench $a_{\rm s}$ to a desired lower value, $a_{\rm s}^{\rm f}$, and shortly hold the atoms in the trap for a time $t_{\rm h}$. We measure that $a_{\rm s}$ converges to its set value with a characteristic time constant of $1\,{\rm ms}$ during $t_{\rm h}$. We then release the atoms from the trap, change $a_{\rm s}$ back close to its initial value and let the cloud expand for $30\,$ms. We probe the momentum distribution $\tilde{n}(k_{x},k_{y})$ by performing standard resonant absorption imaging on the expanded cloud (Methods). The measurement is then repeated at various values of $a_{\rm s}<a_{\rm dd}$ in a fixed trap geometry. The momentum distribution shows a striking behaviour (Fig. 2). For large enough $a_{\rm s}$, $\tilde{n}(k_{x},k_{y})$ shows a single narrow peak with an inverted aspect ratio compared to the trapped gas, typical of a stable BEC Pitaevskii:2016 (Fig. 2a). We define the center of the distribution as the origin of $k$. In contrast, when further decreasing $a_{\rm s}$ below a critical value $a_{\rm s}^{*}$, we observe a sudden appearance of two symmetric finite-momentum peaks, of similar shape and located at $k_{y}=\pm k_{\rm rot}$ (Fig. 2b-c). By repeating the experiment several times, we observe that the peaks consistently appear at the same locations, and they are visible in the averaged distributions. To quantitatively investigate the peak structures, we fit a sum of three Gaussian distributions to the central cuts of the average $\tilde{n}(k_{x},k_{y})$ (Fig. 2a1-c1). From the fit we extract the central momentum, $k_{\rm rot}$, and the amplitude, $A_{*}$, of the side peaks.

A major fingerprint of the roton mode in dBECs is its geometrical nature, leading to a universal scaling $k_{\rm rot}\sim 1/\ell_{z}$ (see e.g. Santos:2003 ; Ronen:2007 ; Blakie:2012 ; Lasinio:2013 ). In addition, the dependency of $k_{\rm rot}$ on $a_{\rm s}$ close to the instability is expected to be mild as $k_{\rm rot}$ remains mainly set by its geometrical nature Santos:2003 ; Blakie:2012 . We investigate both dependencies in the experiment. In a first set of experiments, we repeat the quench measurements for various trap parameters and extract $k_{\rm rot}$. We clearly observe the expected geometrical scaling. $k_{\rm rot}$ shows a marked increase with $1/\ell_{z}=2\pi\sqrt{m\nu_{z}/h}$, matching well with a linear progression with a slope of 1.61(4) (Fig. 2d). Note that no dependence of $k_{\rm rot}$ on $\nu_{y}$ is observed. In a second set of experiments, we fix the trap geometry and explore the dependence of $k_{\rm rot}$ on $a_{\rm s}$. We observe that, within our experimental uncertainty, $k_{\rm rot}$ stays constant when decreasing $a_{\rm s}$  (Fig. 2e). This behaviour contrasts the one expected for a phonon-driven modulation instability that exhibits a strong $a_{\rm s}$-dependence Nguyen:2017 .

To gain a deeper understanding of the roton excitations in our system and its dynamical population, we develop both an analytical model and full numerical simulations and compare the findings with our experimental data. Our analytical model starts by calculating the roton spectrum of the stationary BEC, generalizing the results of Ref. Santos:2003 to a non-radially symmetric configuration. Since the roton wavelength is much smaller than the extension of our 3D BEC along $y$, this mode can be evaluated using a local density approximation in $y$. Hence, for our model, we consider a dBEC homogeneous along $y$, of axial density $n_{0}$, harmonically confined along $x$ and $z$. To analytically evaluate the roton spectrum, we approximate the BEC wavefunction using the Thomas-Fermi (TF) approximation. For dominant DDI, $\epsilon_{\rm dd}=a_{\rm dd}/a_{\rm s}\geq 1$, we find that $\epsilon(k_{y})$ indeed rotonizes (Methods). In the vicinity of the roton minimum and for $\epsilon_{\rm dd}\sim 1$, the computed dispersion acquires a gapped quadratic form similar to that of helium rotons:

[TABLE]

The roton momentum reads as $k_{\rm rot}=\sqrt{2m}(E_{0}^{2}-\Delta^{2})^{1/4}/\hbar$, with $\Delta=\sqrt{E_{0}^{2}-E_{I}^{2}}$, $E_{I}=2gn_{0}(\epsilon_{dd}-1)/3$, $g=4\pi\hbar^{2}a_{\rm s}/m$, and $E_{0}^{2}=2g\epsilon_{dd}n_{0}\frac{\hbar^{2}}{2m}\left(X^{-2}+Z^{-2}\right)$. The TF radii verify $X^{2},Z^{2}\propto gn_{0}$ so that $E_{0}\propto h\nu_{z}$, with a scale set by $\epsilon_{\rm dd}$ and $\nu_{z}/\nu_{x}$ but independent of $gn_{0}$ Odell2003 . Close to the instability ($\Delta\simeq 0$), the roton momentum thus follows a simple geometrical scaling, $k_{\rm rot}=\kappa/\ell_{z}$ with the geometrical factor $\kappa$ depending on $\nu_{z}/\nu_{x}$ alone. For completeness, we have also performed full 3D numerical calculations of the static Bogoliubov spectrum of a finite trapped dBEC, confirming the existence and scaling of the roton mode (Supplementary Information).

The above stationary description accounts for the existence of the roton mode in the cigar geometry used in experiments, and predicts the scaling of $k_{\rm rot}$ and $\Delta$ with the system parameters. However the quench of $a_{\rm s}$ introduces a dynamic, which is crucial to quantitatively reproduce the experimental observations. The reduction of $a_{\rm s}$ decreases the contact interaction and additionally induces a compression of the cloud. This yields a dynamical modification of the local roton dispersion relation, and the roton may be destabilised during the evolution. This dynamical destabilisation is well accounted for by a self-similar model (SSM) describing the evolution of the cloud shape after the quench Pitaevskii:2016 . In particular, we consider a 3D harmonic confinement, and, starting from the stationary TF profile at $a_{\rm s}^{\rm i}$, we evaluate the evolution of the cloud along the change of $a_{\rm s}$ assuming that the TF shape is maintained but with time-dependent TF radii (Methods). We then estimate the local (along $y$) instantaneous roton spectrum $\epsilon(k_{y},y,t_{\rm h})$ using a local density approximation, i.e. evaluating Eq. (1) with the experimentally calibrated $a_{\rm s}(t_{\rm h})$, and the $n_{0}(y,t_{\rm h})$, $X(y,t_{\rm h})$, $Z(y,t_{\rm h})$ estimated from the 3D profile (Methods). We indeed find that the roton gap decreases and eventually turns imaginary ($\Delta(t_{\rm h})^{2}<0$) (Fig. 3a). When this occurs the population at $k_{y}\approx k_{\rm rot}$ exponentially grows with an instantaneous local rate $2\operatorname{Im}[\epsilon(k_{y},y,t_{\rm h})]/\hbar$, $\operatorname{Im}[.]$ indicating the imaginary part, giving rise to two symmetric side peaks in the axial momentum distribution $\tilde{n}_{1}(k_{y},t_{\rm h})\propto\int\exp\left(2\int_{0}^{t_{\rm h}}\operatorname{Im}[\epsilon(k_{y},y,t)]dt\right)dy$  (Fig. 3b) (Methods). The center of the peaks, $\pm k_{\rm rot}(t_{\rm h})$, also evolves with $t_{\rm h}$ but quickly converges after a few ms to its final value, $k_{\rm rot}$ (Fig. 3c). The measured $k_{\rm rot}$ and the calculated values from our parameter-free theory are in remarkable agreement (Fig. 2d-e and 3c).

Our SSM quantitatively explains the experimental observations and provides us with a physical understanding. For completeness, we perform numerical simulations (NS) of the system dynamics. We calculate the time evolution of the generalized non-local Gross-Pitaevskii equation, which accounts also for quantum fluctuations, three-body loss processes, and finite temperature, not included in the SSM (Supplementary Information). The first two contributions limit the peak density of the atomic cloud and stabilise the dBEC against collapse Ronen:2007 ; Bohn:2009 ; Parker:2009 ; Igor:2016 ; Chomaz:2016 , whereas the latter term thermally seeds the initial roton-mode population. The NS confirm both the SSM results and the experimental observations. Indeed the calculations show that, few ms after the quench, the system develops roton peaks in momentum space (Fig. 3d) and short-wavelength density modulations at the center of the BEC (Fig. 3e), showing the predicted roton confinement Lasinio:2013 . The extracted value of $k_{\rm rot}$ and its geometrical scaling are in very good quantitative agreement with both the experimental data and the SSM calculations (Fig. 2d-e and 3c). Interestingly, the timescales for the emergence of roton peaks in the NS are few ms longer than those observed in the experiment. The origin of this time shift remains an open question, whose answer could e.g. require a refinement of current models of quantum fluctuations. However, despite this delay, the growth rate of the roton population in the NS matches both the experimental observations and the SSM predictions in the early dynamics (Fig. 3b).

In the experiment, we study the time evolution of the roton mode in a fixed geometry $(\nu_{z},\lambda)$. In a first set of measurements, we fix $a_{\rm s}^{\rm f}$ and follow the dynamics by recording the momentum distribution at various $t_{\rm h}$. We observe that $k_{\rm rot}$ does not change significantly while the roton population initially grows, in excellent agreement with the theories (Fig. 3b-c). The growth rate is a particularly relevant quantity as it is directly connected to the imaginary excitation energy in the Bogoliubov description (for the roton, its gap), as revealed by our theory. In a second set of experiments, we thus systematically study the growth rate of the roton population for various $a_{\rm s}^{\rm f}$ (Fig. 4). Our data show that the roton mode begins to populate from a critical value of the scattering length, $a_{\rm s}^{*}$. For $a_{\rm s}^{\rm f}>52a_{0}$, we do not observe the roton peaks at any time. For $a_{\rm s}^{\rm f}\lesssim 52\,a_{0}$, after a time delay, $A_{*}$ undergoes an abrupt increase. At longer time, $A_{*}$ then saturates and eventually slowly decreases while atoms are coincidentally lost (Fig. 4a). By further lowering $a_{\rm s}^{\rm f}$, the roton population exhibits a faster growth rate and shorter time delay.

From the growth rate of the roton population we now extract an overall roton gap $\bar{\Delta}(a_{\rm s})$ (Methods). In brief, the dynamical variation of the gap is determined by the leading time dependence of $a_{\rm s}$ and $A_{*}(a_{\rm s}^{\rm f},t_{\rm h})\propto\exp\left(2\int_{0}^{t_{\rm h}}dt\operatorname{Im}[\bar{\Delta}(a_{\rm s}(t))]/\hbar\right)$. To investigate the scaling $\bar{\Delta}(a_{\rm s})$, we first consider the analytical expression of the elementary local roton gap, $\Delta$ as a function of $a_{\rm s}$ (see Eq.(1)). Developing $\Delta$ in the vicinity of $a_{\rm s}=a_{\rm s}^{*}$ yields $\Delta(a_{\rm s})\propto\left(a_{\rm s}^{*}-a_{\rm s}\right)^{1/2}$. Here we confirm this scaling by considering the generic power-law dependence $\operatorname{Im}[\bar{\Delta}(a_{\rm s})]/\hbar=\Gamma\delta(a_{\rm s})$ with $\delta(a_{\rm s})=\left(\frac{a_{\rm s}^{*}-a_{\rm s}}{a_{0}}\right)^{\beta}\left(a_{\rm s}\leq a_{\rm s}^{*}\right)$. Consequently, $A_{*}$ grows as $\exp\left(2\Gamma\int_{0}^{t_{\rm h}}\delta(a_{\rm s}(t))\,dt\right)\,=\,\exp\left(2\Gamma T_{\rm h}\right)$. By rescaling the time variable as $t_{\rm h}\rightarrow T_{\rm h}=\int_{0}^{t_{\rm h}}\delta(a_{\rm s}(t))dt$, we determine the parameters $a_{\rm s}^{*}$ and $\beta$ such that all the curves of Fig. 4a fall on top of each other (Fig. 4b). The best overlap of the experimental curves is found for $a_{\rm s}^{*}=53.0(4)\,a_{0}$ and $\beta=0.55(8)$. We then perform an exponential fit to all the rescaled data for $T_{\rm h}<3\,$ms and extract $\Gamma=465(83)\,{\rm s}^{-1}$. The same analysis, applied to the calculated population growth from the SSM for different $a_{\rm s}^{\rm f}$, gives $a_{\rm s}^{*}\simeq 52.8\,a_{0}$, $\beta\simeq 0.55$, and $\Gamma\simeq 472\,{\rm s}^{-1}$, which are very close to the experimental values. This time-resolved study allows to readily extract the imaginary roton gap, $\operatorname{Im}[\bar{\Delta}]$, from both the experiment and the SSM, as a function of $a_{\rm s}$ and $t_{\rm h}$ (Fig. 4c). Our observations show the softening of the roton mode at $a_{\rm s}=a_{\rm s}^{*}$ and the expected increase of $\operatorname{Im}[\bar{\Delta}]$ for $a_{\rm s}<a_{\rm s}^{*}$, nicely grasping the dynamics during the growth of the roton population (inset of Fig. 4c). The observed behaviour shows again a remarkable agreement with the theory predictions.

To conclude, our work demonstrates the power of weakly-interacting dipolar quantum gases to access the regime of large-momentum, yet low-energy, excitations dressed by interactions. This newly accessible regime, which is largely unexplored in ultracold gases, raises fundamental questions and opens novel directions. Future key developments are to study the impact of such low-lying excitations on the superfluid behavior of a dBEC Lasinio:2013 ; Blakie:2012 ; Wilson:2010 ; Pitaevskii:2016 , the interplay between phononic and rotonic modes in the stability diagram of the quantum gas Ronen:2007 ; Martin:2012 , and the possible role of roton excitations as triggering mechanism of the recently observed instability leading to the formation of metastable droplet arrays in elongated traps Igor:2016 ; Wenzel:2017 . Of particular interest is the prospect of creating a supersolid and striped ground-states in dBECs Boninsegni:2012 ; Wenzel:2017 . Indeed, the short-wavelength density modulation tied to the roton softening together with the quantum stabilization mechanism can favor the formation of such an exotic phase of matter, in which crystalline order coexist with phase coherence. Contrary to recent experiments Li:2017 ; Leonard:2017 , where the density modulation is imposed by external fields yielding a supersolid-like arrangement of infinite stiffness, a dipolar supersolid would be compressible. Experiments on dBECs provide hence the exciting opportunity to unveil similarities and differences among complementary approaches to supersolidity.

Methods

.1 Trapping geometries.

The BEC is confined in a harmonic trapping potential $V(\mathbf{r})=2m\pi^{2}(\nu_{x}^{2}x^{2}+\nu_{y}^{2}y^{2}+\nu_{z}^{2}z^{2})$, characterized by the frequencies $(\nu_{x},\nu_{y},\nu_{z})$. The ODT results from the crossing of two red-detuned laser beams of $1064\,$nm wavelength at their respective focii. One beam, called vODTb, propagates nearly collinear to the $z$-axis and the other, denoted hODTb propagates along the $y$-axis. By adjusting independently the parameters of the vODTb and hODTb, we can widely and dynamically control the geometry of the trap (see Supplementary Information). In particular, $\nu_{y}$ is essentially set by the vODTb power while $\nu_{z}$ (and $\nu_{x}$) independently by that of the hODTb. This yields an easy tuning of the trap aspect ratio, $\lambda=\nu_{z}/\nu_{y}$, relevant for our cigar-like geometry.

After reaching condensation (see Supplementary Information), we modify the beam parameters to shape the trap into an axially elongated configuration, favourable for observing the roton physics ($\nu_{y}\ll\nu_{x},\nu_{z}$). The trapping geometries probed in the experiments, whose $(\lambda,\nu_{x},\nu_{y},\nu_{z})$ are reported in Supplementary Table 1, are achieved by changing the hODTb power with the vODTb power set to 7W so that $\nu_{y}$ and $\nu_{z}/\nu_{x}$ are kept roughly constant. Only the green triangle in Fig. 2d is obtained in a distinct configuration, with a vODTb power of 2W, leading to $(\nu_{x},\nu_{y},\nu_{z})=(156,17,198)\,$Hz and $\lambda=11.6$. The $(\lambda,\nu_{x},\nu_{y},\nu_{z})$ are experimentally calibrated via exciting and probing the center-of-mass oscillation of thermal samples. We note that the final atom number $N$, BEC fraction $f$, and temperature $T$ after the shaping procedure depend on the final configurations, as detailed in Supplementary Table 1.

.2 Quench of the scattering length $a_{\rm s}$.

To control $a_{\rm s}$ we use a magnetic Feshbach resonance between ${}^{166}{\rm Er}$ atoms in their absolute ground state, which is centered around $B=0\,$G Chin2010fri . The $B$-to-$a_{\rm s}$ conversion has been previously precisely measured via lattice spectroscopy, as reported in Ref. Chomaz:2016 . Errors on $a_{\rm s}$, taking into account statistical uncertainties of the conversion and effects of magnetic field fluctuations (e.g. from stray fields), are of 3-to-5 $a_{0}$ for the relevant $a_{\rm s}$ range $27$-$67$ $a_{0}$ 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, $B_{\rm set}$, twice. First we perform the quench itself and abruptly change $B_{\rm set}$ from $0.4\,$G ($a_{\rm s}^{\rm i}=61\,a_{0}$) to the desired lower value (corresponding to $a_{\rm s}^{\rm f}$) at the beginning of the hold in trap ($t_{\rm h}=0\,$ms). Second we prepare the ToF expansion and imaging conditions (see Method Imaging procedure) and abruptly change $B_{\rm set}$ from the quenched value back to $0.3\,$G ($a_{\rm s}=57\,a_{0}$) at the beginning of the ToF expansion ($t_{\rm f}=0\,$ms). Due to delays in the experimental setup, e.g. coming from eddy currents in our main chamber, the actual $B$ value felt by the atoms responds to a change of $B_{\rm set}$ via $B(t)=B_{\rm set}(t)+\tau dB/dt$ Lahaye:2008 . By performing pulsed-radio-frequency spectroscopy measurements (pulse duration 100 $\mu$s) on a BEC after changing $B_{\rm set}$ (from 0.4 to 0.2 G), we verify this law and extract a time constant $\tau=0.98(5)\,$ms. Then $a_{\rm s}$ is also evolving during $t_{\rm h}$ and $t_{\rm f}$ on a similar timescale. This effect is fully accounted for in the experiments and simulations, and we report the roton properties as a function of the effective value of $a_{\rm s}$ at $t_{\rm h}$. We use $t_{\rm h}$ ranging from $2$ to $7$ ms. The lower bound on $t_{\rm h}$ comes from the time needed for $a_{\rm s}$ to effectively reach the regime of interest. We then consider the initial evolution for which $t_{\rm h}/\nu_{y}\,;\,t_{\rm h}/\tau_{\rm coll}\ll 1$, $1/\tau_{\rm coll}$ being the characteristic collision rate. We estimate that $\tau_{\rm coll}$ ranges typically from 40 to 90 ms in the initial BECs of Supplementary Table 1 at $a_{\rm s}^{\rm i}=61\,a_{0}$. Experimentally we observe that the roton, if it ever develops, has developed within the considered range of $t_{\rm h}$.

.3 Imaging procedure.

In our experiments, we employ ToF expansion measurements, accessing the momentum distribution of the gas Pitaevskii:2016 , to probe the roton mode population. We let the gas expand freely for $t_{\rm f}=30\,$ms, which translates the spatial imaging resolution ($\sim\,3.7\,\mu{\rm m}$) into a momentum resolution of $\sim 0.32\,\mu{\rm m}^{-1}$ while our typical roton momentum is $k_{\rm rot}\gtrsim 2\,\mu{\rm m}^{-1}$. After ToF expansion, we record 2D absorption pictures of the cloud 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 $\sim 15^{\rm o}$ compared to the $z$-axis within the $xz$-plane. Thus the ToF images essentially probe the spatial density distribution $n_{\rm tof}\left(x,y,t_{\rm f}\right)$ in the $xy$ plane. When releasing the cloud ($t_{\rm f}=0\,$ms), we change $B$ back to $B=0.3\,$G (Method Quench of the scattering length $a_{\rm s}$). 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 $a_{\rm s}$ to $57\,a_{0}$ allows to minimize the time during which the evolution happens in the small-$a_{\rm s}$ regime where roton physics develops, such that we effectively probe only the short-time evolution of the gas. In this work, we use the simple mapping:

[TABLE]

which neglects the initial size of the cloud in trap and the effect of interparticle interactions during the free expansion. Using real-time simulations (see Supplementary Information), we simulate the experimental sequence and are able to compute both the real momentum distribution from the in-trap wavefunction and the spatial ToF distribution $30\,$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 agree within $5\%$. This confirms that the interparticle interactions play little role during the expansion and justifies the use of Eq. (2).

.4 Fit procedure for the ToF images.

For each data point of Figs. 2-4, we record between 12 and 25 ToF images. By fitting a two-dimensional Gaussian distribution to the individual images, we extract their origin $(k_{x},k_{y})=(0,0)$ and recenter each image. From the recentered images, we compute the averaged $\tilde{n}(k_{x},k_{y})$, from which we characterise the linear roton developing along $k_{y}$. To do so, we extract a one-dimensional profile $\tilde{n}_{1}(k_{y})$ by averaging $\tilde{n}(k_{x},k_{y})$ over $k_{x}$ within $|k_{x}|\leq k_{\rm m}=3.5\mu{\rm m}^{-1}$ : $\tilde{n}_{1}(k_{y})=\int_{-k_{\rm m}}^{k_{\rm m}}\tilde{n}(k_{x},k_{y})dk_{x}/\int_{-k_{\rm m}}^{k_{\rm m}}dk_{x}$. To quantitatively analyse the observed roton peaks, we fit a sum of three Gaussian distributions to $\tilde{n}_{1}(k_{y})$. One Gaussian accounts for the central peak and its centre is constrained to $k_{0}\sim 0$. Its amplitude is denoted $A_{0}$. The two other Gaussians are symmetrically located at $k_{0}\pm k_{y}^{*}$, and we constrain their sizes and amplitudes to be identical, respectively equal to $\sigma^{*}$ and $A_{*}$ . We define the roton population contrast via $\mathcal{C}=A_{*}/A_{0}$. We focus on the roton side peaks by constraining $k_{y}^{*}>0.5\,\mu{\rm m}^{-1}$ and $\sigma^{*}<3\,\mu{\rm m}^{-1}$ (peak at finite momentum of moderate extension compared to the overall distribution).

In order to analyse the onset and evolution of the roton population (see Fig.4), we perform a second run of the fitting procedure, in which we constrain the value of $k_{y}^{*}$ more strictly. The interval of allowed values is defined for each trapping geometry based on the results of the first run of the fitting procedure. We use the results of the ($a_{\rm s}^{\rm f}$, $t_{\rm h}$)-configurations where the peaks are clearly visible and we set the allowed $k_{y}^{*}$-range to that covered by the 95% confidence intervals of the first-fitted $k_{y}^{*}$ in these configurations. This constraint enables that the fitting procedure estimates the residual background population on the relevant momentum range for the roton physics even for $a_{\rm s}>a_{\rm s}^{*}$ (see e.g. Fig.4a).

.5 Time-rescaling analysis and roton gap estimate.

In Fig. 4a-b, we systematically analyse the time evolution of the roton population for various $a_{\rm s}^{\rm f}$ and link it to the roton spectrum in a quench picture. The roton population is embodied by the amplitude $A_{*}$ of the three-Gaussian fit (Method Fit procedure for the ToF images), which measures the density $\tilde{n}_{1}(k_{\rm rot})$. $A_{*}$ is observed to initially increase if $a_{\rm s}^{\rm f}<a_{\rm s}^{*}$, its growth rate increases for decreasing $a_{\rm s}^{\rm f}$ and $a_{\rm s}^{*}$ depends on the trap geometry.

If dynamically unstable, the $k_{\rm rot}$ component is expected to grow exponentially with an instantaneous rate $2\operatorname{Im}[\epsilon(k_{\rm rot},t_{\rm h})]/\hbar$. We hence expect an initial growth of $A_{*}$ of the form:

[TABLE]

$\bar{\Delta}(a_{\rm s}^{\rm f},t_{\rm h})$ being the instantaneous value of the overall roton gap, corresponding to an average over the cloud, after quenching $a_{\rm s}$ to $a_{\rm s}^{\rm f}$. Our results show that the most relevant effect of the quench on the roton spectrum is given by the reduction of $a_{\rm s}$ itself. Hence the time dependence of $\bar{\Delta}$ is determined by $a_{\rm s}(t)$ and, by monitoring $A_{*}(a_{\rm s}^{\rm f},t_{\rm h})$, one can readily extract the scattering-length dependent gap $\bar{\Delta}(a_{\rm s})$, $\bar{\Delta}(a_{\rm s}^{\rm f},t)=\bar{\Delta}(a_{\rm s}(t))$. To investigate the scaling $\bar{\Delta}(a_{\rm s})$, we first consider the simpler case of the static and local roton gap $\Delta$. From Eq. (1) of the main text characterising this elementary gap, assuming $n_{0}$, $X$ and $Z$ fixed, and using the fact that $\Delta(a_{\rm s}=a_{\rm s}^{*})=0$, one easily obtain that for $a_{\rm s}$ in the vicinity of $a_{\rm s}^{*}$, $\Delta(a_{\rm s})\propto(a_{\rm s}^{*}-a_{\rm s})^{1/2}$. Here we verify this scaling for the global and dynamical quantity $\bar{\Delta}(a_{\rm s})$, considering the generic power-law dependence $\operatorname{Im}[\Delta(a_{\rm s})]/\hbar=\Gamma\delta(a_{\rm s})$ with $\delta(a_{\rm s})=\left(\frac{a_{\rm s}-a_{\rm s}^{*}}{a_{0}}\right)^{\beta}\left(a_{\rm s}(t)\leq a_{\rm s}^{*}\right)$, in which the parameters $\Gamma$, $a_{\rm s}^{*}$ and $\beta$ mainly depend on the trap geometry.

Our full set of data $A_{*}(a_{\rm s}^{\rm f},t_{\rm h})$ for the geometry $(\nu_{z},\lambda)=(149\,{\rm Hz},4.3)$ enables us to assess this scaling. Indeed Eq. (3) then reads:

[TABLE]

This defines a time rescaling $t_{\rm h}\rightarrow T_{\rm h}=\int_{0}^{t_{\rm h}}\delta(a_{\rm s}(t))dt$ along which all our experimental data $A_{*}(a_{\rm s},t_{\rm h})$ should collapse in a unique curve, marked by an initial exponential growth of rate $2\Gamma$. The relevant values of $a_{\rm s}^{*}$, $\beta$ and $\Gamma$ are the one that result in the best overlap of the data for the initial growth of $A_{*}$ (minimal spread in $T_{\rm h}$).

In order to determine $a_{\rm s}^{*}$ and $\beta$, we then plot our full dataset as a function of $T_{\rm h}$ for various trial values of these parameters and evaluate the relevance of the trial couple $(a_{\rm s}^{*},\beta)$. Precisely, we assess the dispersion in $T_{\rm h}$ of the full dataset for few fixed $A_{*}=A_{*}^{i}$. We use a panel of 10 values $A_{*}^{i}$ within [80,200]$\mu{\rm m}^{2}$, corresponding to the range of the initial growth of $A_{*}$ for the geometry of Fig. 4a-b. We interpolate the experimental data $A_{*}(T_{\rm h})$ for each $a_{\rm s}^{\rm f}$ using piecewise cubic polynomial interpolation, extract the corresponding set of $T_{\rm h}^{i}$ at which $A_{*}(T_{\rm h}^{i})=A_{*}^{i}$. For each $i$, we evaluate the spread of the $T_{\rm h}^{i}(a_{\rm s}^{\rm f})$ by two complementary quantities: (a) the square-root of their variance and (b) the discrepancy between the maximal and minimal value of the set. We finally estimate the accuracy of a unified dependence $A_{*}(T_{\rm h})$ for a given ($a_{\rm s}^{*},\beta$) via the geometrical average of (a) and (b) for all the ten $A_{*}^{i}$ values. We fit the relevant $a_{\rm s}^{*}$ and $\beta$ by minimizing this averaged spread. For the experimental data of Fig. 4a-b, this results in $a_{\rm s}^{*}=53.0(4)a_{0}$ and $\beta=0.55(8)$.

Using these values of $a_{\rm s}^{*}$ and $\beta$, we observe that the initial growth of $A_{*}$ extends for $T_{\rm h}\leq 3\,$ms (before saturating and decreasing for longer $T_{\rm h}$). We estimate $\Gamma$ by performing an exponential fit on the full set of data $A_{*}(T_{\rm h})$ with $T_{\rm h}\leq 3\,$ms. This gives $\Gamma=465(83)\,{\rm s}^{-1}$. Note that Fig. 4a-b only show 5 values of $a_{\rm s}^{\rm f}$ while the reported analysis considers all available experimental data, i.e. 11 values between $31\,a_{0}$ and $52\,a_{0}$.

From the formula with the estimated $a_{\rm s}^{*}$, $\beta$ and $\Gamma$, we can then compute the global roton gap $\bar{\Delta}$. The extracted value is only meaningful for $\bar{\Delta}^{2}<0$. For the experiments, we also restrict its relevance to $T_{\rm h}\leq 3\,$ms so that, e.g., in the inset of Fig. 4c, $\operatorname{Im}[\bar{\Delta}]$ is shown up to $t_{\rm h}=4\,$ms after which the roton population is observed to deviate from the exponential growth picture (see Figs. 3b and 4a). Note that the $\operatorname{Im}[\bar{\Delta}]$ estimates at the different $t_{\rm h}$ are not independent.

In Fig. 2, we report on the roton momentum as a function of the system characteristics ($a_{\rm s}$, $\ell_{z}$). Here we estimate $k_{\rm rot}$ from a $t_{\rm h}$ value that is individually optimised for each $a_{\rm s}$ and $\ell_{z}$ investigated (largest visibility). We point out that the selected $t_{\rm h}$ correspond, in the $A_{*}$ dynamics, to the late stage of the exponential growth, close to the maximum (i.e. $T_{\rm h}\sim 3\,$ms). Here, the atom loss remains at a few percents level.

In Fig. 2d, we show the value of $k_{\rm rot}$ at the onset of the population of the roton mode, i.e. at $a_{\rm s}=a_{\rm s}^{*}$. We estimate $a_{\rm s}^{*}$ for each trap geometry by analysing the evolution of the roton population with $a_{\rm s}$. We employ a simplified approach with respect to the one in Fig. 4 (see Supplementary Information), which we estimate to lead to a maximum underestimate for $a_{\rm s}^{*}$ of about $1.5\,a_{0}$, lying within our experimental uncertainty on $a_{\rm s}$; see Supplementary Table 2.

.6 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. Santos:2003 for rotons in infinite q2D traps. We consider a dBEC homogeneous along $y$ but harmonically confined with frequencies $\nu_{x}$ and $\nu_{z}$ along $x$ and $z$. For sufficiently strong interactions the BEC is in the TF regime on the $xz$ plane, in which the BEC wavefunction acquires the form $\psi_{0}(\boldsymbol{\rho})=\sqrt{n(\boldsymbol{\rho})}$, with $n(\boldsymbol{\rho})=n_{0}\left(1-(x/X)^{2}-(z/Z)^{2}\right)$, where $X$ and $Z$ are the TF radii, and $\boldsymbol{\rho}=(x,z)$. The calculation of $n_{0}$, $X$ and $Z$ is detailed at the end of this section.

Due to the axial homogeneity, the elementary excitations of the Bogoliubov-de Gennes spectrum have a defined axial momentum $k_{y}$, and take the form $\delta\psi(\boldsymbol{r},t)=u({\boldsymbol{\rho}})e^{ik_{y}y-i\epsilon t/\hbar}-v({\boldsymbol{\rho}})e^{-ik_{y}y+i\epsilon t/\hbar}$, where $u,v$ denote the amplitudes of the spatial modes oscillating in time with characteristic frequency $\epsilon/\hbar$ (see Supplementary Information for more details). We consider the standard GPE without LHY correction and three-body losses, and insert the perturbed solution $\psi(\boldsymbol{r},t)=\left(\psi_{0}(\boldsymbol{\rho})+\eta\delta\psi(\boldsymbol{r},t)\right)e^{-i\mu t/\hbar}$ , where $\mu$ is the chemical potential associated to $\psi_{0}$ and $\eta\ll 1$. After linearization we obtain the BdG equations for $f_{\pm}(\boldsymbol{\rho})=u(\boldsymbol{\rho})\pm v(\boldsymbol{\rho})$:

[TABLE]

where

[TABLE]

where $U(\mathbf{r})$ is the binary interaction potential for two particles separated by $\mathbf{r}$. $U(\mathbf{r})=g\left(\delta(\mathbf{r})+\frac{3\epsilon_{\rm dd}}{4\pi}\frac{1-3\cos^{2}\theta}{|\mathbf{r}|^{3}}\right)$ includes both contact and dipolar interactions (see Supplementary information).

Employing $f_{+}(\boldsymbol{\rho})=W(\boldsymbol{\rho})\psi_{0}(\boldsymbol{\rho})$, and for $k_{y}\gg 1/X,1/Z$, we obtain the following equation for the function $W(\boldsymbol{\rho})$:

[TABLE]

where ${\tilde{x}}=x/X$, ${\tilde{z}}=z/Z$, $E(k_{y})=\hbar^{2}k_{y}^{2}/2m$. For $\epsilon_{dd}=1$ the last term of Eq. (9) vanishes. In that case, the lowest-energy solution is given by $W=1$, whose eigen-energy builds, as a function of $k_{y}$, the dispersion $\epsilon_{0}(k_{y})$ with

[TABLE]

with $E_{0}^{2}=2g\epsilon_{dd}n_{0}\frac{\hbar^{2}}{2m}\left(\frac{1}{X^{2}}+\frac{1}{Z^{2}}\right)$. In the vicinity of $\epsilon_{\rm dd}=1$, the effect of the last term in Eq. (9) may be evaluated perturbatively, resulting in the dispersion

[TABLE]

with $E_{I}=\frac{2}{3}gn_{0}(\epsilon_{dd}-1)$.

This expression for the dispersion presents a roton minimum for $\epsilon_{dd}>1$ at $k_{\rm rot}=\frac{1}{\hbar}\sqrt{2mE_{I}}$. Expanding Eq. (11) in the vicinity of the roton minimum, $\epsilon(k_{y})^{2}\simeq\epsilon(k_{\rm rot})^{2}+\frac{1}{2}\left[\frac{d^{2}\epsilon^{2}(k_{y})}{dk_{y}^{2}}\right]_{k_{y}=k_{\rm rot}}$, we obtain Eq. (1) of the main text, with $\Delta=\epsilon(k_{\rm rot})=\sqrt{E_{0}^{2}-E_{I}^{2}}$. At the instability, $\Delta=0$, and $k_{\rm rot}=\frac{1}{\hbar}\sqrt{2mE_{0}}$.

Employing a similar procedure as in Ref. Odell2003 we obtain that the BEC aspect ratio $\chi=Z/X$ fulfills:

[TABLE]

with $\lambda_{\perp}=\nu_{x}/\nu_{z}$ and

[TABLE]

These two equations fully determine the TF solution for given $\epsilon_{dd}$, $gn_{0}$, and $\lambda_{\perp}$. By inserting the expressions of $X^{2}$ and $Z^{2}$ in $E_{0}$, we find for $\epsilon_{\rm dd}\simeq 1$:

[TABLE]

whereas $\chi$ simplifies into $\chi=\lambda_{\perp}(1+\sqrt{1+1/\lambda_{\perp}})$. As a result, at the instability $k_{\rm rot}\ell_{z}$ depends only on the transverse confinement aspect ratio $\lambda_{\perp}$, giving the geometrical factor $\kappa$:

[TABLE]

.7 Self similar model and instantaneous roton spectrum.

The dynamics of the condensate during and after the quench is crucial to understand the resulting momentum peaks and their direct relation to the growth of roton excitations. Due to its large characteristic momentum, the roton spectrum may be evaluated at any time using local density approximation. On the other side, the evolution of the local density is determined by the global dynamics of the condensate induced by the quench. Hence, interestingly, the analysis of the effect of the quench on the roton spectrum may be performed by combining a self-similar theory of the global dynamics Giovanazzi:2006 , describing the evolution at low $k_{y}$, and the model of the roton spectrum developed in the Method Analytical dispersion relation for an infinite axially elongated geometry, accounting for the high-$k_{y}$ region of interest. We note that this analysis approximates the real time evolution provided by the standard NLGPE Baranov:2008 ; Pitaevskii:2016 ; Blakie:2008 . This treatment is valid for moderate-enough quenches so that each of the descriptions applies (see below for an in-depth discussion) and as long as the high-$k$ modes population only minimally impact the BEC dynamics, that is for short-enough time scales.

We assume that the condensate preserves its TF shape during the evolution:

[TABLE]

where $X(t)=b_{x}(t)X_{0}$, $Y(t)=b_{y}(t)Y_{0}$, and $Z(t)=b_{z}(t)Z_{0}$ are the re-scaled TF radii, with $b_{x,y,z}(t)$ the scaling coefficients ($b_{x,y,z}(0)=1$). Their evolution after the quench can be deduce from solving the hydrodynamics equations (see Supplementary information). Prior to the quench of $a_{\rm s}$, the stationary TF solution is obtained from the self-consistent equations $\chi_{x}(0)=\frac{\omega_{x}}{\omega_{z}}\sqrt{C(0)/A(0)}$ and $\chi_{y}(0)=\frac{\omega_{y}}{\omega_{z}}\sqrt{C(0)/B(0)}$. Solving these equations and using normalisation provides $X(0)$, $Y(0)$, $Z(0)$, and $n(0)$ for known numbers of atoms $N$ and trap frequencies $\omega_{x,y,z}$. We use this stationary solution as the initial condition at the start of the quench ($t=0$), and solve the system of differential equations resulting from the hydrodynamics evolution to obtain the scaling coefficients. An example for the relevant parameters of Fig. 3 is shown in Supplementary Fig. 3.

For a given position $y$ we then evaluate the local TF profile: n({\mathbf{r}},t)=n_{0}(y,t)\allowbreak\big{[}1-(x/X(y,t))^{2}-(z/Z(y,t))^{2}\big{]}, where $n_{0}(y,t)=n_{0}(t)\left[1-(y/Y(t))^{2}\right]$, $X(y,t)=X(t)\left[1-(y/Y(t))^{2}\right]^{1/2}$, and $Z(y,t)=Z(t)\left[1-(y/Y(t))^{2}\right]^{1/2}$, with $n_{0}(t)=n_{0}(0)/b_{x}b_{y}b_{z}$. We may then employ a local density approximation and evaluate the local (in $y$) instantaneous roton spectrum using the results of the Method Analytical dispersion relation for an infinite axially elongated geometry:

[TABLE]

The roton population grows exponentially when the spectrum becomes imaginary, leading to the appearance of the corresponding momentum peaks at large momenta. From the local instantaneous roton spectrum we may estimate the momentum peak as

[TABLE]

where we assume equal seeding for all unstable modes. We then evaluate the total population of the peak $N_{*}(t)=\int dk_{y}\tilde{n}_{1}(k_{y},t)$, as well as the roton momentum $\langle k\rangle=\int dk_{y}k_{y}\frac{\tilde{n}_{1}(k_{y},t)}{N_{*}(t)}$.

The previously discussed SSM does not describe properly the evolution when the shrinking of the condensate, and hence the increase of its peak density, is too large. As the gas gets dense, quantum fluctuations, which are not considered in the SSM, introduce in our experiments an effective repulsion that crucially prevents large densities or condensate collapse (see Supplementary Information). This limits the use of the SSM to the vicinity of $a_{\rm s}^{*}$. Moreover, for our tightest trap the SSM results are also unreliable, since the theory predicts condensate collapse before the momentum peak develops. That said, as shown in the main text, the SSM provides not only a clear intuitive understanding of our measurements, but to a very large extent also a very good quantitative agreement with both our experimental results and our NS based on a generalized GPE, which is detailed in Supplementary Information.

Data availability: The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.

Supplementary Information

.8 Production of ${}^{166}{\rm Er}$ BECs.

We prepare a ${}^{166}{\rm Er}$ BEC similarly to Refs. Chomaz:2016 ; Aikawa:2012 . From a narrow-line magneto-optical trap with $3\times 10^{7}$ 166Er  atoms, automatically spin-polarized in their absolute lowest Zeeman sub-level  Frisch:2012 , at about $10\mu{\rm K}$, we directly load the atomic gas in a crossed optical dipole trap (ODT) with an efficiency of more than $30\%$. A uniform magnetic field, $B$, is permanently applied along the vertical $z$ axis, fixing the dipole orientation, while its value is varied during the experimental sequence, to tune $a_{\rm s}$ (see Method Quench of the scattering length $a_{\rm s}$). We achieve condensation by means of evaporative cooling in the crossed ODT at $B=1.9\,$G ($a_{\rm s}=80(2)\,a_{0}$). During the evaporation procedure, we first change the power and then the ellipticity of one of the ODT beams (see Section .9). The final atom number, typically $10^{5}$, and condensed fraction, typically $70\%$, are assessed by fitting the time-of-flight (ToF) absorption images of the gas to a bimodal distribution, sum of a TF profile and a broad Gaussian background.

.9 Details on the trapping geometries.

The harmonic potential $V(\mathbf{r})=2m\pi^{2}(\nu_{x}^{2}x^{2}+\nu_{y}^{2}y^{2}+\nu_{z}^{2}z^{2})$ trapping the 166Er atoms results from the crossing of two red-detuned laser beams of $1064\,$nm wavelength at their respective focii; the hODTb propagating along the $y$-axis and the vODTb, propagating nearly collinear to the $z$-axis. The vODTb has a maximum power of 7 W and an elliptical profile with waists of $110$ and $55\,\mu{\rm m}$ along $x$ and $y$ respectively. The hODTb has a maximum power of 24 W, a vertical waist $w_{z}=18\,\mu{\rm m}$, and a controllable horizontal waist, $w_{x}=\lambda^{\prime}w_{z}$. The ellipticity $\lambda^{\prime}$ can be tuned from $1.57$ to $15$ by time averaging the frequency of the first-order deflection of an Acousto-Optic Modulator Ahmadi:2005 . By adjusting independently $\lambda^{\prime}$ and the powers of the vODTb and of the hODTb, we can control the geometry of the trap. Precisely, $\nu_{y}$ is essentially set by the vODTb power, $\nu_{z}$ by that of the hODTb, and $\nu_{x}$ is controlled by both the power and ellipticity of the hODTb, with $\nu_{z}/\nu_{x}\approx\lambda^{\prime}$.

We use the tuning of the hODTb power and ellipticity to perform evaporative cooling to quantum degeneracy (see Section .8). After reaching condensation, we again modify the beam parameters to shape the trap into an axially elongated configuration ($\nu_{y}\ll\nu_{x},\nu_{z}$). The trapping geometries probed in the experiments, $(\lambda,\nu_{x},\nu_{y},\nu_{z})$, are reported in Supplementary Table 1. They are achieved by changing the hODTb power with $\lambda^{\prime}=1.57$ and the vODTb power set to its maximum so that $\nu_{y}$ and $\nu_{z}/\nu_{x}$ are kept roughly constant. Only the green triangle in Fig. 2d is obtained in a distinct configuration, with a vODTb power of 2W, leading to $(\nu_{x},\nu_{y},\nu_{z})=(156,17,198)\,$Hz and $\lambda=11.6$. The $(\lambda,\nu_{x},\nu_{y},\nu_{z})$ are experimentally calibrated via exciting and probing the center-of-mass oscillation of thermal samples. We note that the final atom number $N$, BEC fraction $f$, and temperature $T$ after the shaping procedure depend on the final configurations, as detailed in Supplementary Table 1. $T$ is extracted from the evolution with the ToF duration $t_{\rm f}$ of the size of the background Gaussian in the TF-plus-Gaussian bimodal fit to the corresponding ToF images of the gas. The values of $N_{\rm c}=fN$ and $T$ are used for the initial states $\psi_{i}$ of our real-time simulations (see Section .14).

.10 $\boldsymbol{a_{\rm s}^{*}}$ for various trap geometries.

We estimate $a_{\rm s}^{*}$ for each trap geometry of Fig. 2d by studying the roton population as a function of $a_{\rm s}$. When decreasing $a_{\rm s}$, the roton population and in particular its contrast $\mathcal{C}=A_{*}/A_{0}$ (see Method Fit procedure for the ToF images) exhibits a sharp increase from an essentially zero and flat contrast. We extract $a_{\rm s}^{*,c}$ as the value of the scattering length corresponding to the onset of the increase of the contrast for $t_{\rm h}=3\,$ms; see Supplementary Fig. 1 and Supplementary Table 2. This is a simplified approach with respect to the one in Fig. 4 (see Method Time-rescaling analysis and roton gap estimate.), which we estimate to lead to a maximum underestimate for $a_{\rm s}^{*}$ of about $1.5\,a_{0}$, which lies within our experimental uncertainty on $a_{\rm s}$.

.11 Generalized Non-local Gross-Pitaevskii equation.

Our theory is based on an extended version of the non-linear Gross-Pitaevskii Equation (NLGPE)

[TABLE]

governing the evolution of a macroscopically occupied wavefunction $\psi(\mathbf{r},t)$, with corresponding atomic density $n(\mathbf{r},t)=|\psi(\mathbf{r},t)|^{2}$ at position $\mathbf{r}$ and time $t$. The standard dipolar NLGPE includes the kinetic energy, external trap potential and the mean-field effect of the interactions Pitaevskii:2016 ; Baranov:2008 . These correspond to the three first terms of Eq. (.11), where the mean-field interaction potential takes the form of a convolution of $n$ with the binary interaction potential

[TABLE]

for two particles separated by $\mathbf{r}$ Baranov:2008 . The first term corresponds to contact interactions between the particles with strength $g=\frac{4\pi\hbar^{2}a_{\rm s}}{m}$. The DDI gives rise to the second term, which depends on both the distance and orientation (angle $\theta$) of the vector $\mathbf{r}$ compared to the polarisation axis ($z$ axis) of the dipoles. Most properties of dBECs are well captured by this standard NLGPE (mean-field) Pitaevskii:2016 ; Baranov:2008 .

Recent experimental and theoretical results, however, have established the importance of accounting for quantum fluctuations in dBECs Igor:2016 ; Chomaz:2016 ; Waechtler:2016 ; Waechtler:2016b ; Bisset:2016 . Their effect can be included in the NLGPE in a mean field treatment through a Lee-Huang-Yang correction to the chemical potential, $\Delta\mu[n]=32g(na_{s})^{3/2}(1+3\epsilon_{\rm dd}^{2}/2)/3\sqrt{\pi}$, which is obtained under a local density approximation Pelster:2011 ; Pelster:2012 . The accuracy of this mean field treatment has been established, e.g., in Refs. Waechtler:2016 ; Waechtler:2016b ; Bisset:2016 , and has proven succesful in explaining recent experimental results Igor:2016 ; Chomaz:2016 . The final nonlinear term in the extended NLGPE accounts for three-body losses Kagan:1998 , with an experimentally determined loss parameter $L_{3}$, which is dependent on $a_{\rm s}$ and typically of the order $L_{3}\simeq 10^{-41}\mathrm{m}^{6}\mathrm{s}^{-1}$, as reported in Ref. Chomaz:2016 .

.12 Bogoliubov-de Gennes spectrum.

Collective excitations of the dBEC are obtained by linearising the NLGPE (see Section .11) around a stationary state $\psi_{0}$, which can be obtained by imaginary time propagation (see Section .14). We write $\psi=\mathrm{e}^{-\mathrm{i}\mu t/\hbar}(\psi_{0}+\eta[u\mathrm{e}^{-\mathrm{i}\epsilon t/\hbar}-v^{*}\mathrm{e}^{+\mathrm{i}\epsilon t/\hbar}])$, where $\mu$ is the chemical potential associated with state $\psi_{0}$, and $u,v$ are spatial modes oscillating in time with characteristic frequency $\epsilon/\hbar$ and $\eta\ll 1$ Ronen:2006 . Inserting this ansatz in the NLGPE, and retaining only terms up to linear order in $\eta$ we obtain the Bogoliubov-de Gennes (BdG) equations

[TABLE]

where the operator $A$, acting on a function $f$ and evaluated at point $\mathbf{r}$, is defined as

[TABLE]

The above equations constitute an eigenvalue problem, which we solve numerically using the Arnoldi method to obtain eigenmodes $(u,v)$ and corresponding excitation energies $\epsilon$. The equations presented here are a generalization of the BdG equations for dipolar systems as derived in Ref. Ronen:2006 , to include the LHY correction accounting for quantum fluctuations.

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 $k_{y}^{(\rm eff)}=\langle k_{y}^{2}\rangle^{1/2}$ Bisset:2013 . The spectrum is discrete with phonon-like collective modes at low $k_{y}^{(\rm eff)}$. For higher $k_{y}^{(\rm eff)}$, the spectrum flattens, but eventually bends upwards again due to the dominant kinetic energy. Instead of developing a smooth minimum, roton excitations appear as isolated low-lying modes at intermediate momenta that depart from the overall spectrum Bisset:2013 . These so-called roton fingers are related to confinement of the roton modes in the inhomogeneous BEC of profile $n_{0}(y)$  Lasinio:2013 , see also discussion in the main text.

The confinement is evident from the BdG calculations, in which the lowest roton mode forms a short-wavelength density modulation localized at the trap center (Supplementary Fig. 6e). This contrasts with phonon modes for which the modulation is delocalised over the entire condensate (Supplementary Fig. 6c). The excited states shown in Supplementary Fig. 6b-e correspond to the density $|\psi_{0}+\eta(u-v^{*})|^{2}$, for particular pairs of $(u,v)$ corresponding to phonon and roton modes (Supplementary Fig. 6c and e), and exemplary modes at higher energies (Supplementary Fig. 6b and d). Even while the amplitude $\eta=0.2$ of the excited modes is taken to be equal in Figs. 6b - 6e, the roton excitation (Supplementary Fig. 6e) leads to markedly larger local density modulations than the phonon excitation (Supplementary Fig. 6c). The ToF signatures in Supplementary Fig. 6 are computed by letting the wave function of the excitation, $\eta(u-v^{*})$, expand ballistically for $30\mathrm{ms}$, i.e. neglecting interactions during the expansion. The resulting density $|\eta(u-v^{*})|^{2}$ is then plotted (Supplementary Fig. 6b1 - e1).

Carrying out the numerical BdG spectrum calculation for various trap parameters, scattering lengths and atom numbers confirms the geometrical scaling $k_{\rm rot}\sim 1/\ell_{z}$, as well as its weak dependence on $a_{\rm s}$ close to the instability as expected from the analytical model. We note in particular a good quantitative agreement of $k_{\rm rot}\ell_{z}$ in Supplementary Fig. 6 with the (stationary) analytic factor $\kappa=1.3$ for this same trap geometry (see Method Analytical dispersion relation for an infinite axially elongated geometry). A qualitative agreement of the BdG stationary spectrum with the experimental observations is however not fully reached because of the role of the dynamics. A treatment of this effect is provided by the SSM and NS as discussed in the main text (see also the corresponding Method and Section .14).

Finally, an alternative way of visualizing the excitation spectrum is achieved by computing the dynamical structure factor at $T=0$ for each mode $(u,v)$ Zambelli:2000 ; Blakie:2012 ,

[TABLE]

where $\omega=\epsilon/\hbar$. The structure factor of the axial modes of our finite trapped system (Supplementary Fig. 6a) presents features resembling a roton-maxon spectrum, its continuum limit being a single-valued spectrum as depicted in Fig. 1b. To enhance the visibility, the delta function in Eq. (24) is replaced by a Gaussian with a small finite width in $\omega$. Since the dynamical structure factor determines the response of the system when probed at specific energies and momenta, such as e.g. in Bragg spectroscopy experiments Brunello:2001 ; Blakie:2012 , it is interesting to note the difference in amplitude between the roton modes and other parts of the spectrum. In particular, in our quench experiments where the system is effectively driven at a range of energies and momenta, one would expect the strongest response from the roton modes.

.13 Self-similar dynamics of the BEC.

In our self-similar model we account for the dynamics of the BEC profile and its effect on the roton spectrum (see details in the corresponding Method). To do so, we assume that the condensate preserves its TF shape during the evolution:

[TABLE]

where $X(t)=b_{x}(t)X_{0}$, $Y(t)=b_{y}(t)Y_{0}$, and $Z(t)=b_{z}(t)Z_{0}$ are the re-scaled TF radii, with $b_{x,y,z}(t)$ the scaling coefficients ($b_{x,y,z}(0)=1$). The corresponding hydrodynamic equations reduce to:

[TABLE]

where $\mu({\mathbf{r}},t)=V({\mathbf{r}})+gn({\mathbf{r}},t)+\int d^{3}r^{\prime}V_{dd}({\mathbf{r}}-{\mathbf{r}}^{\prime})n({\mathbf{r}}^{\prime},t)$ is the local chemical potential. The latter acquires the form:

[TABLE]

where $D(t)=1+\epsilon_{dd}F_{1}$, $A(t)=1+\epsilon_{dd}(F_{1}-F_{2}+F_{3})$, $B(t)=1+\epsilon_{dd}(F_{1}-F_{2}-F_{3})$, and $C(t)=1+\epsilon_{dd}(F_{1}+2F_{2})$. Here we have introduced the functions:

[TABLE]

with $\chi_{x}(t)=Z(t)/X(t)$, $\chi_{y}(t)=Z(t)/Y(t)$, $\alpha(\chi_{x},\chi_{y},u)=(\chi_{x}^{2}+\chi_{y}^{2})(1-u^{2})/2+u^{2}$, and $\beta(\chi_{x},\chi_{y},u)=(\chi_{x}^{2}-\chi_{y}^{2})(1-u^{2})/2$. Substituting Eq. (27) into Eq. (26) we obtain a closed set of equations for the scaling parameters:

[TABLE]

Prior to the quench of $a_{\rm s}$, the stationary TF solution is obtained from the self-consistent equations $\chi_{x}(0)=\frac{\nu_{x}}{\nu_{z}}\sqrt{C(0)/A(0)}$ and $\chi_{y}(0)=\frac{\nu_{y}}{\nu_{z}}\sqrt{C(0)/B(0)}$. Solving these equations and using normalisation provides $X(0)$, $Y(0)$, $Z(0)$, and $n(0)$ for known numbers of atoms $N$ and trap frequencies $\nu_{x,y,z}$. We use this stationary solution as the initial condition at the start of the quench ($t=0$), and solve the system of differential equations (31), (32) and (33) to obtain the scaling coefficients. An example for the relevant parameters of Fig. 3 is shown in Supplementary Fig. 3.

.14 Numerical simulations of the evolution.

Our simulations of the NLGPE are performed using a split operator technique. The evolution operator over a time $\Delta t$ ($\Delta t\to i\Delta t$ for imaginary time evolution) may be approximately split as $e^{-i\hat{H}\Delta t/\hbar}=e^{-i\hat{T}\Delta t/\hbar}e^{-i\hat{V}\Delta t/\hbar}+O(\Delta t^{2})$. In this expression, $\hat{T}$ is the kinetic energy term, and $\hat{V}$ the potential energy. The effective potential energy for the evolution is given by the sum of external potential, interparticle interactions, local LHY correction, and three-body losses (see Section .11).

We first evaluate using imaginary time evolution the initial BEC wavefunction, $\psi_{0}(\boldsymbol{r})$, prior to the quench of $a_{\rm s}$. The initial wavefunction for the subsequent real-time evolution is then constructed via $\psi_{i}=\psi_{0}+\Delta\psi$, where $\Delta\psi$ accounts for thermal fluctuations, which we simulate by populating the excited states of the system as described in Ref. Blakie:2008 . The excited states used in this procedure are obtained from the full BdG calculation detailed in Section .12. Starting with this initial wavefunction $\psi_{i}$, we mimic as close as possible the conditions of our experiments, including ramping, holding, and ToF times. In particular, we include the experimentally calibrated $a_{\rm s}(t)$ (Method Quench of the scattering length $a_{\rm s}$). Moreover, for the value of the three-body loss coefficient $L_{3}$ we use a linear fit of the experimentally determined values Chomaz:2016 . From the simulated evolution, we obtain the 3D wavefunction of the gas as a function of $t_{\rm h}$, $\psi(\bf r,t_{\rm h})$, from which we can extract the spatial and momentum distributions.

The simulation of the ToF expansion is performed in two steps. First we use a multi-grid analysis in order to rescale the size of the numerical box as the cloud expands during the ToF expansion. After some expansion time the density drops significantly, and the subsequent evolution can be readily calculated via $e^{-i\hat{T}t/\hbar}$. Our NS show clearly that the effect of nonlinearity is small during the first stages of the evolution, and hence that the ToF expansion indeed may be employed to image the momentum distribution of the condensate at the time in which the trap is opened.

We evaluate the integrated momentum distribution $\tilde{n}(k_{y},t_{\rm h})=\int dk_{x}dk_{z}|\tilde{\psi}({\bf k},t_{\rm h})|^{2}$ with $\tilde{\psi}({\bf k},t_{\rm h})$ the Fourier transform of $\psi(\bf r,t_{\rm h})$. After $t_{\rm h}$ of a few ms, $\tilde{n}(k_{y},t_{\rm h})$ shows clear roton peaks. The exact $t_{\rm h}$ value for the peak emergence depends on the gas characteristics (in particular $T$) and on $a_{\rm s}$. We evaluate the roton momentum as the mean value of the momentum in the roton peak.

From the imaginary-time evolution simulations, we are also able to predict the $a_{\rm s}=a_{\rm s}^{*,\rm st}$ threshold for the mean-field instability of the BEC, which corresponds to the absence of a mean-field stable solution Waechtler:2016b ; Bisset:2016 . To find this instability boundary, we proceed by steps. We start by calculating the ground-state solution for a given $a_{\rm s}$ that we know to be well within the stable regime. We then reduce $a_{\rm s}$ in small steps and successively calculate the corresponding ground-state solution using the solution of the previous step as starting condition. We do so until no mean-field stable solution can be found. The predicted $a_{\rm s}^{*,\rm st}$ are reported in Supplementary Table 2. For the theory predictions shown in Fig. 2d (NS and SSM; see Method Analytical dispersion relation for an infinite axially elongated geometry), we use quenched $a_{\rm s}$ values such that the instability boundary is just slightly crossed, $a_{\rm s}=a_{\rm s}^{*,\rm st}-1a_{0}$.