Roton Excitations in an Oblate Dipolar Quantum Gas
Jan-Niklas Schmidt, Jens Hertkorn, Mingyang Guo, Fabian Böttcher, Matthias Schmidt, Kevin S. H. Ng, Sean D. Graham, Tim Langen, Martin Zwierlein, Tilman Pfau
RRoton Excitations in an Oblate Dipolar Quantum Gas
J.-N. Schmidt, ∗ J. Hertkorn, ∗ M. Guo, F. B¨ottcher, M. Schmidt, K.S.H. Ng, S.D. Graham, T. Langen, M. Zwierlein, and T. Pfau †
5. Physikalisches Institut and Center for Integrated Quantum Science and Technology,Universit¨at Stuttgart, Pfaffenwaldring 57, 70569 Stuttgart, Germany MIT-Harvard Center for Ultracold Atoms, Research Laboratory of Electronics,and Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA (Dated: February 3, 2021)We observe signatures of radial and angular roton excitations around a droplet crystallizationtransition in dipolar Bose-Einstein condensates. In situ measurements are used to characterize thedensity fluctuations near this transition. The static structure factor is extracted and used to identifythe radial and angular roton excitations by their characteristic symmetries. These fluctuations peakas a function of interaction strength indicating the crystallization transition of the system. Wecompare our observations to a theoretically calculated excitation spectrum allowing us to connectthe crystallization mechanism with the softening of the angular roton modes.
The roton dispersion relation is essential to understandthe thermodynamics and density fluctuations in super-fluid helium [1–4]. Initially interpreted by Feynman asthe “ghost of a vanishing vortex ring” [5], nowadays theroton is seen as a precursor to the crystallization of asystem [6]. Quantum gases with dipolar interactions fea-ture a similar dispersion relation due to the anisotropicand long-range nature of their interaction [7, 8]. In con-trast to helium, the high tunability of atomic quantumgases allows for systematic studies of roton excitations.Tuning inter-atomic interactions can soften the rotons,which trigger an instability and the formation of a crys-tal of quantum droplets [9–11]. For elongated systemsconfined in cigar-shaped traps, this softening leads tothe emergence of one-dimensional supersolid states thatsimultaneously exhibit superfluid flow and crystalline or-der [12–17].In cylindrically symmetric oblate traps, two types ofroton excitations have been predicted to play a crucialrole in the instability [18–20]. These two modes are theradial and angular roton modes corresponding to the twospatial degrees of freedom in the system. The spectrumof these two-dimensional (2D) modes is more complexthan the spectrum of a previously studied elongated one-dimensional (1D) supersolid [21], making it a challenge todistinguish their individual contributions to the crystal-lization. While the roton modes were directly observedin 1D dipolar systems [22–24], the various angular rotonmodes in 2D and their connection to the crystallizationhave remained elusive.In this Letter, we observe signatures of radial and an-gular roton excitations of an oblate dipolar BEC aroundthe phase transition to a 2D droplet crystal. These ex-citations leave their imprint in the density fluctuationsof the gas, which we measure in situ, giving us directaccess to the static structure factor S ( k ) [21, 25]. Thequantity S ( k ) features a peak at finite momentum and adistinct sixfold angular symmetry upon approaching thecrystallization transition. We use mean-field simulations of the excitation spectrum to interpret the experimentalresults. The observed emergence of angular structure isthereby directly linked to the softening of angular rotonmodes.The excitation spectrum of dipolar BECs in cylin-drically symmetric traps can be theoretically studiedin the Bogoliubov-de Gennes (BdG) framework. Wesolve the BdG-equations including the first beyond mean-field correction [21, 26] term known as the Lee-Huang-Yang (LHY) correction using our experimental param-eters [27]. We calculate the excitation spectrum of thelow-lying modes in the BEC up to the point of the in-stability of 15 × Dy atoms in a trap as a func-tion of scattering length. The trapping frequencies are ω/ π = [35 , . , z . The trap geometry is deliberately made asym-metric to numerically lift the degeneracy in the x-y-direction [28]. The radial and angular roton modes foundin these spectra can soften, by varying the scatteringlength, to energies much smaller than the trap frequen-cies [18–20, 26].The mode patterns of the density fluctuation δn ( r )corresponding to the radial and angular roton modes areshown in Fig. 1(a). We additionally present the spatialpower spectra of the mode patterns given by the squaredmodulus of the Fourier transform |F [ δn ] | ( k ), illustratingthe individual contributions to the static structure factor[24, 30–34]. Radial roton modes are circularly symmet-ric and represent ring-like density modulations at non-zero radial wavevector. Angular rotons have an angularoscillatory structure in addition to the ring-like radialdensity modulation. We describe the angular oscillationwith sin( mφ ), where φ is the azimuthal angle and theinteger m > m nodal lines has a 2 m -fold symmetry resulting ina four- and sixfold symmetric power spectrum for the m = 2 and m = 3 angular roton modes, respectively. Incontrast, the quadrupole mode and higher-lying phonon a r X i v : . [ c ond - m a t . qu a n t - g a s ] F e b m = m = m = m = x y k x k y modepattern δ n ( r ) spectrum | ℱ [ δ n ] ( k ) - m = m = m = m =
77 78 79 80 81 82010203040 scattering length a s ( a ) e x c i t a t i o n e n e r g y ω / π ( H z ) ( a )( b ) FIG. 1. (a) Normalized mode patterns of the density fluctua-tions close to the transition at a s (cid:39) . a and their spatialpower spectra in the x - y -plane. Shown are the lowest radial( m = 0) and angular rotons ( m = 1 , , m angularnodal lines. (b) Corresponding Bogoliubov excitation energiesas a function of scattering length a s , with the insets indicat-ing the emergence of a blood cell shape in the ground-statedensity close to the instability at a s (cid:39) . a . For decreasingscattering lengths, first the m = 2 followed by the m = 1 and m = 3 angular roton modes drops below the radial dipolemodes at the trap frequency and softens further towards theinstability. For information on higher modes (e.g. cyan line)and details see [27, 29]. modes might feature similar azimuthal symmetries butat smaller radial wavevector.In Fig. 1(b) we show the low-lying excitation ener-gies and ground state shapes from the BEC side as afunction of scattering length towards the phase transi-tion point [35]. In the BEC regime at high scatteringlengths a s (cid:39) a , furthest away from the transitionpoint, the dipole mode has the lowest excitation energyand lies at exactly the trap frequency [36]. For decreas-ing scattering lengths, several higher-lying modes rapidlydecrease in energy. These are the radial and twice degen-erate angular roton modes. At around a s (cid:39) . a , the m = 2 mode drops below the trap frequency followed bythe m = 1 and m = 3 mode. Near the phase transitionpoint at a s (cid:39) . a the m = 2 mode is only separatedby a few Hz from the m = 3 mode. The ordering of theangular roton modes close to the phase transition pointdepends on a nontrivial interplay between the trap aspectratio and the interaction strength [37]. Near a s (cid:39) . a the parabola-shaped BEC groundstate transforms into a biconcave blood cell-like shape,which forms as it is energetically favorable in cylindri-cal geometries to push part of the density to the outerrim. This further enhance the softening of angular ro-tons as they become excitations of the ring-shaped re-gion of maximal density [18]. Previous studies on blood-cell-shaped ground states [18–20, 26, 32, 37–41] did notinclude the LHY correction that has since been shownto stabilize the system against collapse. We also findit to enlarge the parameter regime for blood-cell-shapedground states by a range in scattering length of approx-imately 1 a .We experimentally study the emergence of the angu-lar rotons with a dipolar BEC with typically 15 × Dy atoms at a temperature T (cid:39)
20 nK. We ad-just the crossed optical dipole trap after evaporationto an almost cylindrical trap with trapping frequencies ω/ π = [35(1) , , z [27]. In the magnetic field range around 30 Gthat we employ, lower three-body losses lead to dropletcrystal lifetimes on the order of 200 ms after crossing thephase transition, which is an increase by a factor of tencompared to previous experiments at similar densities[13, 27]. We image the cloud in situ with a resolutionof 1 µ m and repeat the experiment around 200 times fora statistical analysis of the atomic densities. We fur-ther quote all scattering lengths relative to a referencescattering length a ref = 91(10) a corresponding to thetransition point because of an overall systematic shift inour scattering length calibration [27].We then extract the static structure factor, which con-nects the spectrum of elementary excitations to the majorcontributing modes in the density fluctuations [24, 30–34]. The intermediate steps of this analysis are shownin Fig. 2 for four distinct scattering lengths: in the BECregime, closer to the transition, in the transition regionand for a droplet crystal.First, we investigate the in situ densities n j ( r )(Fig. 2(a)). To remove residual contributions of thedipole mode the center of mass of the atomic densitydistribution is taken as the origin. In addition, we post-select in an interval of ±
15 % with respect to the meanatom number at each scattering length [27]. We observethat the crystal structure is randomly oriented in indi-vidual images and consequently gets washed out whenaveraging over many images (Fig. 2(b)). This highlightsthe continuous rotational symmetry breaking upon cross-ing the crystalline phase. In order to account for therandomly oriented crystal in the further analysis, we de-termine the individual rotation angles θ in Fourier spaceand align the single-shot images [27]. The rotated images n θj ( r ) are used to create new mean images (cid:104) n θ ( r ) (cid:105) , whichreveal the emergence of the crystal structure (Fig. 2(c)).As expected, the rotation does not affect the mean imagein the BEC regime. (a) − . a +0 . a +2 . a +5 . a ( × . (b) ( × . (c) µ m . . . . . c o l u m nd e n s i t y ( m − ) (d) . µ m − ˆ k x ˆ k y S ( k x , k y ) h n θ ( r ) ih n ( r ) i n θ j ( r ) FIG. 2. (a) Single-shot images for four different relative scat-tering lengths, in the BEC regime (+5 . a ), closer to thetransition point (+2 . a ), in the transition region (+0 . a )and for a droplet crystal ( − . a ). (b) Mean images of theunrotated images showing no clear crystalline structure. (c)Aligned images (see text) indicate the presence of dropletsin the mean image. (d) 2D structure factor showing an in-creasing height of the peaks at finite momentum | k | indicat-ing the approaching transition point. The central area below k min / π (cid:39) . µ m − (see text) was masked out. Second, we calculate the individual fluctuation pat-terns around this mean image δn θj ( r ) = n θj ( r ) − (cid:104) n θ ( r ) (cid:105) given by the deviations of the individual post-selectedand rotated images from the mean (cid:104) n θ ( r ) (cid:105) . We de-termine the mean power spectrum (cid:104)| δn θ ( k ) | (cid:105) from aFourier transform δn θj ( k ) = (cid:82) d r δn θj ( r ) e i k · r of thesefluctuation patterns. This mean power spectrumis closely connected to the static structure factor S ( k ) = (cid:104)| δn θ ( k ) | (cid:105) /N for homogeneous systems [25, 34,42] and also provides valuable insights into the nature ofthe excitations in non-homogeneous systems [43–47].The resulting 2D static structure factor S ( k ) is shownin Fig. 2(d). We restrict the further analysis to momentabetween k min / π (cid:39) . µ m − and k max / π (cid:39) µ m − to account for the finite size of the cloud and the finiteimaging resolution in the experiment [27]. The quantity S ( k ) features several peaks that lie approximately on aring with radius | k | around the origin. The height of theindividual peaks increases when approaching the transi-tion point. Additionally, the angular spreading of thesepeaks across the ring changes with scattering lengths andremains invisible in the case of unrotated images. The en-hancement of the peaks along the alignment axis ( ˆ y -axis)is a result of the rotation algorithm, which always alignsthe images according to their individual most dominantpeak in Fourier space. Exploiting the cylindrical symme- S m a x ( a r b . u . ) (b)BECdroplet crystal − . − . . . . . a s − a ref ) ( a )0 . . k r o t / ( π ) ( µ m − ) (c)0 . . . . . . . k/ (2 π ) ( µ m − )0100200300400 S ( k )( a r b . u . ) +5 . a +3 . a +2 . a +0 . a − . a − . a − . a (a)29.41 29.51 29.61 29.71 29.81Magnetic field (G) FIG. 3. (a) Radial distribution S ( k ) of the two-dimensionalstructure factor after integration over the angular coordinatefor different relative scattering lengths. A clear peak at finitewavevector rises towards the phase transition. For claritythe lines for smaller scattering lengths were shifted vertically.(b) The amplitude of this peak is obtained from a Gaussianfit reaching its maximum at the phase transition. (c) Rotonmomentum shifts towards larger values in the droplet regimewhere it stays constant. The dash-dotted line on the leftindicates the smallest momentum k min / π (cid:39) . µ m − . Thegreen diamond indicates the maximum S ( k ) and the dashedline the roton momentum at that maximum. The gray areain (b) and (c) indicates the transition region. try of the trap we transform S ( k ) to polar coordinates S ( k x , k y ) → S ( k, φ ) and analyze its radial and angularbehaviour separately.The radial behaviour S ( k ) is analyzed by integratingover the angular direction to identify modes at finite mo-mentum | k | independent of their angular symmetry. Theresult is shown in Fig. 3(a) for scattering lengths aroundthe transition. Starting in the BEC regime we find S ( k )to be relatively flat with only a small peak at around k/ π (cid:39) . µ m − . Closer to the transition this peakrises and shifts towards larger momenta. We determinethe center momentum k rot and amplitude S max of thepeak using a Gaussian fit and show it in Fig. 3(b)-(c).The peak of the structure factor S max first increaseswith decreasing scattering lengths and then features amaximum before it decreases again. This is consistentwith the expectation of enhanced fluctuations at thephase transition due to the softening and thermal popu-lation of several roton modes close to the phase transition[24]. Compared to 1D, the increased number of low-lying − . − . . . . . a s − a ref ) ( a )0 . . . . F o u r i e r c o m p o n e n t s ( a r b . u . ) (b) BECdroplet crystal α α − π − π/ π/ π Angle φ (rad)100200300400 S ( φ )( a r b . u . ) +5 . a +3 . a +2 . a +0 . a − . a − . a − . a (a)29.41 29.51 29.61 29.71 29.81Magnetic field (G) FIG. 4. (a) Angular distribution S ( φ ) of the two-dimensional static structure factor integrated over the inter-val k/ π ∈ [0 . , . µ m − around the roton momentum k rot for different relative scattering lengths. (b) Decomposition of S ( φ ) into Fourier components matching the symmetry of thelowest two angular roton modes m = 2 and m = 3 in theBEC. For clarity, the lines in (a) were shifted vertically forsmaller scattering lengths. The gray area in (b) indicates thetransition region. excitations with different symmetries results in highershot-to-shot fluctuations in 2D. From individual imagesit is therefore more challenging to distinguish betweenan angular roton mode and a formed droplet crystal fea-turing the same symmetry. The formation of a dropletcrystal is a direct consequence of the excitation of an an-gular roton mode with the same symmetry assuming adynamical preparation scheme. We account for these un-certainties by marking a transition region (+0 . ± . a rather than a single point [27]. In the droplet regime thepeak amplitude decreases and its position stays approx-imately constant at k rot / π (cid:39) . µ m − which roughlyindicates the inverse droplet distance. The peak alsobroadens further indicating a competition between dif-ferent droplet numbers and spacings.The angular behaviour S ( φ ) is obtained by integrating S ( k, φ ) over the annulus with k ∈ [0 . , . µ m − [27]. S ( φ ) is shown in Fig. 4(a) and allows us to attribute theenhancement of the fluctuations to individual modes.Starting in the BEC regime, only two comparativelysmall peaks at φ = 0 and φ = π are visible which increasefor lower scattering lengths and can mainly be attributedto our rotation algorithm. Closer to the transition, four additional peaks at φ = ± π/ φ = ± π/ m = 3 angular roton mode. In the droplet regime theseintermediate peaks start to wash out, presumably due tothe competition of the three- and four-droplet configura-tions.We further quantify the changing mode populationwith the lowest coefficients of a Fourier expansion of theperiodic function S ( φ ) shown in Fig. 4(b) [27]. On theBEC side, the Fourier coefficients α n give an indicationof the underlying symmetry that then can be connectedto the modes from the simulation of Fig. 1. We focus onthe coefficients α and α which both have a low contri-bution in the BEC regime.An increasing weights indicate a softening of severalangular roton modes. The α and α mode increasetowards the transition, in agreement with the simula-tion shown in Fig. 1(b). The α weight saturates nearthe transition and stays constant in the droplet regime.However, the α weight reaches a maximum after thetransition and becomes smaller in the droplet regime.A stronger α weight than α is not supported by thepresented theory if one assumes the population of thosemodes is in thermal equilibrium, as the m = 2 modehas a slightly lower energy than the m = 3 mode. Thiseffect either hints towards non-equilibrium dynamics ortowards the limits of the LHY approximation in the the-oretical description [48–50].In the crystalline domain the two weights approacheach other again indicating that neither of the two an-gular roton modes are dominant. In this regime S ( k )cannot be viewed as a measurement of excitations ontop of the crystalline ground state because we find stateswith competing droplet numbers broadening the peak inthe radial distribution of the 2D structure factor (seeFig. 3(a) and [27]). Density patterns at single scatter-ing lengths could have two to five droplets present ineach shot (see [27]), indicating that the mean atomic den-sity does not have same symmetry as the ground state.Then, the weights α i reflect the symmetry of the observeddroplet crystal rather than an excitation on top of thiscrystalline ground state. The similarity of the α and α weights therefore indicates similar probabilities to find adroplet crystal with fourfold or sixfold symmetry.In conclusion, we have reported on signatures of radialand angular roton modes by investigating the 2D staticstructure factor S ( k ) of a dipolar BEC. The characteris-tic sixfold symmetry of S ( k ) in the BEC regime can beidentified with the population of the angular roton mode.These observations are supported by simulations of theexcitation spectrum. Our study lays the foundation fora better understanding of the dominant excitations ofdipolar BEC in oblate traps close to the transition to adroplet crystal. It connects the low-lying angular rotonmodes to the crystallization mechanism and the forma-tion of droplets.This work is supported by the German Research Foun-dation (DFG) within FOR2247 under Pf381/16-1 andBu2247/1, Pf381/20-1, FUGG INST41/1056-1 and theQUANT:ERA collaborative project MAQS. M.G. andM.Z. acknowledge funding from the Alexander von Hum-boldt Foundation. ∗ These authors contributed equally to this work. † [email protected][1] L. Landau, Phys. Rev. , 356 (1941).[2] R. P. Feynman, Phys. Rev. , 262 (1954).[3] R. P. Feynman, Phys. Rev. , 1291 (1953).[4] R. P. Feynman and M. Cohen, Prog. Theor. Phys. ,261 (1955).[5] R. P. Feynman, Chapter II Application of Quantum Me-chanics to Liquid Helium , edited by C. Gorter, Progressin Low Temperature Physics, Vol. 1 (Elsevier, 1955) pp.17 – 53.[6] P. Nozi`eres, J. Low Temp. Phys. , 45 (2004).[7] L. Santos, G. V. Shlyapnikov, and M. Lewenstein, Phys.Rev. Lett. , 250403 (2003).[8] M. Wenzel, F. B¨ottcher, J.-N. Schmidt, M. Eisenmann,T. Langen, T. Pfau, and I. Ferrier-Barbut, Phys. Rev.Lett. , 030401 (2018).[9] H. Kadau, M. Schmitt, M. Wenzel, C. Wink, T. Maier,I. Ferrier-Barbut, and T. Pfau, Nature , 194 (2015).[10] I. Ferrier-Barbut, H. Kadau, M. Schmitt, M. Wenzel,and T. Pfau, Phys. Rev. Lett. , 215301 (2016).[11] F. B¨ottcher, J.-N. Schmidt, J. Hertkorn, K. S. H. Ng,S. D. Graham, M. Guo, T. Langen, and T. Pfau, Rep.Prog.Phy. , 012403 (2021).[12] L. Tanzi, E. Lucioni, F. Fam`a, J. Catani, A. Fioretti,C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno,Phys. Rev. Lett. , 130405 (2019).[13] F. B¨ottcher, J.-N. Schmidt, M. Wenzel, J. Hertkorn,M. Guo, T. Langen, and T. Pfau, Phys. Rev. X , 011051(2019).[14] L. Chomaz, D. Petter, P. Ilzh¨ofer, G. Natale, A. Traut-mann, C. Politi, G. Durastante, R. M. W. van Bijnen,A. Patscheider, M. Sohmen, M. J. Mark, and F. Fer-laino, Phys. Rev. X , 021012 (2019).[15] M. Guo, F. B¨ottcher, J. Hertkorn, J.-N. Schmidt,M. Wenzel, H. P. B¨uchler, T. Langen, and T. Pfau,Nature , 386 (2019).[16] G. Natale, R. M. W. van Bijnen, A. Patscheider, D. Pet-ter, M. J. Mark, L. Chomaz, and F. Ferlaino, Phys. Rev.Lett. , 050402 (2019).[17] L. Tanzi, S. Roccuzzo, E. Lucioni, F. Fam`a, A. Fioretti,C. Gabbanini, G. Modugno, A. Recati, and S. Stringari,Nature , 382 (2019).[18] S. Ronen, D. C. Bortolotti, and J. L. Bohn, Phys. Rev.Lett. , 030406 (2007).[19] R. M. Wilson, S. Ronen, J. L. Bohn, and H. Pu, Phys.Rev. Lett. , 245302 (2008).[20] R. N. Bisset, D. Baillie, and P. B. Blakie, Phys. Rev. A , 043606 (2013).[21] J. Hertkorn, F. B¨ottcher, M. Guo, J.-N. Schmidt, T. Lan-gen, H. P. B¨uchler, and T. Pfau, Phys. Rev. Lett. ,193002 (2019).[22] L. Chomaz, R. M. Van Bijnen, D. Petter, G. Faraoni, S. Baier, J. H. Becher, M. J. Mark, F. W¨achtler,L. Santos, and F. Ferlaino, Nat. Phys. , 442 (2018),1705.06914.[23] D. Petter, G. Natale, R. M. W. van Bijnen, A. Patschei-der, M. J. Mark, L. Chomaz, and F. Ferlaino, Phys. Rev.Lett. , 183401 (2019).[24] J. Hertkorn, J.-N. Schmidt, F. B¨ottcher, M. Guo,M. Schmidt, K. Ng, S. Graham, H. P. B¨uchler, T. Lan-gen, M. Zwierlein, and T. Pfau, arXiv:2009.08910(2020).[25] C.-L. Hung, X. Zhang, L.-C. Ha, S.-K. Tung, N. Gemelke,and C. Chin, New J. Phys. , 075019 (2011).[26] S. Ronen, D. C. Bortolotti, and J. L. Bohn, Phys. Rev.A , 013623 (2006).[27] See Supplemental Material at [URL] for further detailson the experiment, analysis and simulation.[28] A. Gallem´ı, S. M. Roccuzzo, S. Stringari, and A. Recati,Phys. Rev. A , 023322 (2020).[29] The m = 1 roton modes weakly lose their degeneracynear the transition due to the slightly asymmetric trapused in the simulations. We show the average of the twoenergies in Fig. 1 and refer to [27] for further details.[30] A. S. Arkhipov, G. E. Astrakharchik, A. V. Belikov, andY. E. Lozovik, J. Exp. Theor. Phys. Lett. , 39 (2005).[31] G. E. Astrakharchik, J. Boronat, I. L. Kurbakov, andY. E. Lozovik, Phys. Rev. Lett. , 060405 (2007).[32] P. B. Blakie, D. Baillie, and R. N. Bisset, Phys. Rev. A , 021604 (2012).[33] P. B. Blakie, D. Baillie, and R. N. Bisset, Phys. Rev. A , 013638 (2013).[34] L. Pitaevskii and S. Stringari, Bose-Einstein Conden-sation and Superfluidity , International series of mono-graphs on physics (Oxford University Press, 2016).[35] The near degeneracy of the roton modes on the BEC sidetranslates into a complex situation on the crystal side. Wefind a number of competing ground states reflecting thedifferent symmetries of the rotons. This precludes a fur-ther meaningful BdG-analysis for the given parameters.[36] W. Kohn, Phys. Rev. , 1242 (1961).[37] R. M. Wilson, S. Ronen, and J. L. Bohn, Phys. Rev. A , 013621 (2009).[38] O. Dutta and P. Meystre, Phys. Rev. A , 053604(2007).[39] R. M. Wilson, S. Ronen, and J. L. Bohn, Phys. Rev. A , 023614 (2009).[40] H. Y. Lu, H. Lu, J. N. Zhang, R. Z. Qiu, H. Pu, andS. Yi, Phys. Rev. A , 023622 (2010).[41] A. D. Martin and P. B. Blakie, Phys. Rev. A , 053623(2012).[42] C.-L. Hung, V. Gurarie, and C. Chin, Science , 1213(2013).[43] J. Esteve, J.-B. Trebbia, T. Schumm, A. Aspect, C. I.Westbrook, and I. Bouchoule, Phys. Rev. Lett. ,130403 (2006).[44] A. Imambekov, I. E. Mazets, D. S. Petrov, V. Gritsev,S. Manz, S. Hofferberth, T. Schumm, E. Demler, andJ. Schmiedmayer, Phys. Rev. A , 033604 (2009).[45] J. Armijo, T. Jacqmin, K. V. Kheruntsyan, and I. Bou-choule, Phys. Rev. Lett. , 230402 (2010).[46] T. Jacqmin, J. Armijo, T. Berrada, K. V. Kheruntsyan,and I. Bouchoule, Phys. Rev. Lett. , 230405 (2011).[47] M. Schemmer, A. Johnson, and I. Bouchoule, Phys. Rev.A , 043604 (2018).[48] R. Sch¨utzhold, M. Uhlmann, Y. Xu, and U. R. Fischer, Int. J. Mod. Phys. B , 3555 (2006).[49] A. R. P. Lima and A. Pelster, Phys. Rev. A , 041604(2011).[50] A. R. P. Lima and A. Pelster, Phys. Rev. A , 063609(2012).[51] K. Baumann, N. Q. Burdick, M. Lu, and B. L. Lev,Phys. Rev. A , 020701 (2014).[52] T. Maier, I. Ferrier-Barbut, H. Kadau, M. Schmitt,M. Wenzel, C. Wink, T. Pfau, K. Jachymski, and P. S.Julienne, Phys. Rev. A , 060702 (2015).[53] T. Maier, H. Kadau, M. Schmitt, M. Wenzel, I. Ferrier-Barbut, T. Pfau, A. Frisch, S. Baier, K. Aikawa,L. Chomaz, M. J. Mark, F. Ferlaino, C. Makrides,E. Tiesinga, A. Petrov, and S. Kotochigova, Phys. Rev.X , 041029 (2015).[54] A. Frisch, M. Mark, K. Aikawa, F. Ferlaino, J. L. Bohn,C. Makrides, A. Petrov, and S. Kotochigova, Nature , 475 (2014).[55] E. Lucioni, L. Tanzi, A. Fregosi, J. Catani, S. Gozzini,M. Inguscio, A. Fioretti, C. Gabbanini, and G. Mod-ugno, Phys. Rev. A , 060701 (2018).[56] Y. Tang, A. Sykes, N. Q. Burdick, J. L. Bohn, and B. L.Lev, Phys. Rev. A , 022703 (2015).[57] Y. Tang, A. G. Sykes, N. Q. Burdick, J. M. DiSciacca,D. S. Petrov, and B. L. Lev, Phys. Rev. Lett. ,155301 (2016).[58] Y. Tang, W. Kao, K.-Y. Li, S. Seo, K. Mallayya,M. Rigol, S. Gopalakrishnan, and B. L. Lev, Phys. Rev.X , 021030 (2018).[59] T. Lindeberg, International Journal of Computer Vision , 283 (1993).[60] T. Lindeberg, International Journal of Computer Vision , 79 (1998).[61] M. Wenzel, F. B¨ottcher, T. Langen, I. Ferrier-Barbut,and T. Pfau, Phys. Rev. A , 053630 (2017).[62] S. M. Roccuzzo and F. Ancilotto, Phys. Rev. A ,041601 (2019).[63] R. N. Bisset, R. M. Wilson, D. Baillie, and P. B. Blakie,Phys. Rev. A , 033619 (2016).[64] D. Baillie, R. M. Wilson, and P. B. Blakie, Phys. Rev.Lett. , 255302 (2017).[65] J. von Neumann and E. Wigner, Phys. Z. , 467 (1929).[66] L. Landau and E. Lifshitz, Quantum Mechanics: Non-Relativistic Theory (Elsevier Science, 1981).
SUPPLEMENTAL MATERIALExperimental protocol and Feshbach resonances
The complete experimental procedure has been de-scribed in detail in Refs [9, 13, 15]. We pre-pare a quasipure BEC of
Dy with T (cid:39)
20 nK ina crossed optical dipole trap at 1064 nm. We thenadiabatically reshape the trap to trapping frequencies ω/ π = [35(1) , , . . z -direction. The magnetic field differencecorresponds to a change in scattering length of about14 . a . This long ramp time and a subsequent waitof additional 80 ms ensures that the system has enoughtime to equilibrate as much as possible before the atomiccloud is imaged in situ with a microscope objective witha resolution of about 1 µ m.We use two broader Feshbach resonances withinthe rich spectrum of resonances in Dy [51–54].We determine the position of these resonances to beat B = 22 . B = 26 . = 2 . = 100(10) mG, respectively. Thecharacterization of the broad resonance is challenging asit overlaps with a dense group of narrow resonances [55].These two resonances dominate the scattering length cal-ibration and allow us to work several Gauss away tosignificantly reduce three-body losses. The lifetime ofdroplets several a away from the phase transition isabout 200 ms in the used range of magnetic fields, whichis longer than in our previous measurements [13]. Thescattering length calibration is perturbed further by twosmall resonances near 30 G, the magnetic field we workat. These resonances are located at B = 29 . B = 29 . = 3(1) mG and∆ = 1(1) mG, respectively. The relatively large error inthe width of the broad resonances comes from these nar-row resonances being near the zero-crossing of the broadresonance. By changing the width of the broad resonancewithin the stated uncertainty the scattering length cali-bration can easily be shifted by several a . We work farfrom the resonances so that this uncertainty only leadsto an overall shift in the scattering length. Thereforewe can determine differences in scattering length with ahigher precision than the overall shift. We quote all scat-tering lengths throughout this work relative to a refer-ence scattering length a ref = 91(10) a corresponding tothe center of the transition region. The error includesthe mentioned uncertainty in the width of the broadFeshbach resonances. The background scattering length a bg = 140(20) a also suffers from a large uncertainty [56–58].For our configuration, lower magnetic fields correspond − . − . . . . . a S − a ref ) ( a )11121314 a t o m nu m b e r N ( ) BECdroplet crystal29.41 29.51 29.61 29.71 29.81Magnetic Field (G)
Supplementary Figure S1. Average atom number for the scat-tering length range in the experiment. Crossing the phasetransition from BEC to droplet crystal the increasing densityleads to larger three-body losses causing lower atom numbersat smaller scattering lengths. Error bars indicate the standarderror of the mean. to lower scattering length. The densities and three-bodylosses increase at lower scattering lengths, leading to alower mean atom number in the droplet. Figure S1 showsthe mean atom number for the complete scattering lengthrange. The mean atom number is about 20 % lower inthe droplet regime than the BEC regime.
Analysis method and angular-averaged staticstructure factor
The procedure to extract the static structure factorfrom the individual in situ image is similar to the onedescribed in [24]. After centering the images to theircenter-of-mass to remove the effect of the dipole mode,we perform a post-selection on atom number for everyscattering length. We only consider images with atomnumbers of ±
15 % around the mean atom number at ev-ery scattering length.The rotational symmetry of the trap generates rota-tionally symmetric averaged density distributions thatwashes out the droplet structure found in individual im-ages. We therefore remove the rotational symmetry andalign the droplet crystals from each experimental realiza-tion to the same angle. After calculating the 2D Fouriertransform we transform the image to polar coordinatesand integrate over radial momenta larger than the small- k threshold k > k min . We then find the angle θ thatcorresponds to the maximum and rotate the original im-age by − θ . Minor imperfections in the imaging systemsuch as astigmatism and distortion can lead to a non-uniform distribution of the rotation angles. This effectbecomes stronger in the BEC where there is no strongfeature at finite momentum. These arbitrary rotationsshould also not affect the evaluation because there is noangular structure.The calculation of the static structure factor con-tains an ensemble average of the power spectrum ofthe fluctuations which is sensitive to single-shot atomnumber fluctuations. We reduce this effect by normal-izing each rotated image to its specific atom number˜ n θj = n θj /N j and calculate the normalized fluctuations δ ˜ n θj ( r ) = ˜ n θj ( r ) − (cid:104) ˜ n θ ( r ) (cid:105) as the deviation of the nor-malized image ˜ n θj ( r ) from the mean normalized image (cid:104) ˜ n θ ( r ) (cid:105) . The structure factor is given by S ( k ) = ¯ N (cid:104)| δ ˜ n θ ( k ) | (cid:105) , (1)where (cid:104)| δ ˜ n θ ( k ) | (cid:105) is the mean power spectrum of thesefluctuations and ¯ N is the mean atom number. The2D Fourier transform of the normalized fluctuations δ ˜ n θj ( k ) = F [ δ ˜ n θj ]( k ) = (cid:82) d r δ ˜ n θj ( r ) e i k · r is obtained fromthe line integrated images. The presented measurementof the static structure factor is limited at small- k and atlarge- k . The finite size of the atom cloud L (cid:39) µ m setsthe small- k limit at k min / π (cid:39) . µ m − . The large- k limit is set by the finite imaging resolution of 1 µ m or k max / π (cid:39) µ m − .We transform the two-dimensional static structurefactor to polar coordinates S ( k x , k y ) → S ( k, φ ) and thenanalyse the radial and angular dependence separately.The angular distribution of the static structure factor S ( φ ) is a periodic function by definition, and it is possi-ble to write it as a series expansion in harmonic function f ( x ) = (cid:80) n α n cos( nx + ϕ n ) with x ∈ [0 , π ). The coeffi-cients α n = (cid:112) a n + b n and the phase ϕ n = arctan( a n /b n )are given by a n = (1 /π ) (cid:82) π f ( x ) cos( nx )d x and b n = (1 /π ) (cid:82) π f ( x ) sin( nx )d x respectively. The weights α i can be related to the weight of the correspondingmode in the Bogoliubov-de Gennes calculation becausethese harmonic functions have exactly the same sym-metry as the density fluctuation patterns of the variousangular roton modes. The angular structure factor wasnormalized to its maximum value before calculating theFourier components, minimizing the influence of theabsolute scale at each scattering length. Determination of the phase transition region
The transition from a BEC to a droplet crystal oc-curs approximately where the peak of the static structurefactor at finite momentum S max reaches a maximum as afunction of scattering length (see Fig. 3 of the main text).We extract this amplitude S max by fitting a S ( k ) using aGaussian fit after subtracting the background. In Fig. S2we compare S max with two other measures to define thephase transition: the average peak density ¯ n max and theaverage spectral weight SW . The spectral weight is de-fined as the integral of the spatial power spectrum | n ( k ) | − . − . . . . . a s − a ref ) ( a )0 . . . . . n o r m a li ze d m e a s u r e BECdroplet crystal
SWS max /S max , ¯ n max / ¯ n max , Supplementary Figure S2. Three different measures to char-acterize the phase transition are shown. The peak amplitudeof the structure factor S max (blue), the spectral weight SW (green), and the average peak density n max (black) are allnormalized to their maxima. While the SW and S max in-crease around +1 . a , the peak density reaches values closeto its maximum at higher relative scattering lengths. Theshaded gray area indicates the region where the phase tran-sition occurs and was determined by these measures (seetext). The maximum values are ¯ n max , = 7 . × m − and S max , = 119. over the annulus with k ∈ [0 . , . µ m − . We normalizeall quantities to their maximum for comparison.Starting in the BEC both the spectral weight SW and the peak amplitude of the structure factor S max show a strong increase towards lower scattering lengths. S max seems to have a extended peak between − . a and +0 . a and decreases again for even lower scatter-ing length. The spectral weight SW saturates around+0 . a , indicating the transition to a droplet crystal.By contrast, the peak density increases at higher scat-tering lengths in the BEC and reaches its maximum at+1 . a .The maximum peak amplitude in the static structurefactor hints towards enhanced fluctuations close to thetransition point due to a softening and thermal popu-lation of multiple roton modes. The peak density andthe spectral weight are useful indicators of this transi-tion but they cannot distinguish between an angular ro-ton mode of a certain symmetry and a formed dropletcrystal with the same symmetry. The increased numberof low-lying modes with different symmetries in oblatetrap geometries might lead to a higher shot-to-shot fluc-tuations. These enhanced fluctuations might result insmall differences between the three measures. We there-fore identify a transition region centered at +0 . a witha width of ± . a (shaded gray area). − − − a s − a ref ) ( a )246 N u m b e r o f d r o p l e t s . . . . . P r o b a b ili t y Supplementary Figure S3. Droplet number as a function ofthe rel. scattering length in the droplet regime. Startingwith a maximum probability for three droplets right after thephase transition, the peak of the distribution shifts towardsfour droplets for lower scattering length.
Droplet number probability distribution
The angular static structure factor S ( φ ) for scatter-ing lengths in the droplet regime suggests the competi-tion between a fourfold and sixfold symmetric dropletcrystal. Both symmetries can be seen in the Fouriercomponent analysis in Fig. 4(b) of the main text. Toidentify the two symmetries, we determine the numberof droplets per image and compare the droplet numberprobability distribution for different scattering lengths.We compute the number of droplets per image by us-ing a droplet number detection algorithm based on theLaplacian-of-Gaussian technique [59, 60]. The resultinghistograms for several scattering lengths in the dropletregime are shown in Fig. S3 as a color plot.We are most likely to observe three droplets directlyafter entering the droplet regime, but four droplets aremore probable deeper in the droplet regime. This changein the average droplet number indicates a change of thedroplet crystal lattice structure. The lattice structureis a triangle for three droplets and a square for fourdroplets. The varying probability of droplet number indi-cates the competition between the sixfold (three droplets)and fourfold (four droplets) symmetry in the static struc-ture factor. Simulation details and higher-lying excitation modes
We describe our system of dipolar atoms in BEC closeto the crystallization transition by the extended Gross-Pitaevskii equation (eGPE) [26, 61, 62] and refer toRef. [21] for details on the simulations. i ¯ h∂ t ψ = H GP ψ, (2) dipolequadrupolehexapole m = m = m = m = m = m =
78 80 82 84102030405060 scattering length a s ( a ) e x c i t a t i o n e n e r g y ω / π ( H z ) Supplementary Figure S4. The excitation spectrum shown inFig. 1 of the main text for higher excitation energies and largerscattering lengths. The quasi-periodic boundary conditionalong the angular direction results in multiple angular rotonmodes that all soften towards the transition point. The firstfive are shown here including their coupling to the dipole andquadrupole mode of the BEC. As an inset we show the avoidedcrossing of the m = 2 angular roton mode with the quadrupolemode in greater detail. where we define H GP = H + g | ψ | + Φ dip [ ψ ] + g qf | ψ | and ψ is normalized to the atom number N = (cid:82) d r | ψ ( r ) | . The term H = − ¯ h ∇ / M + V ext contains the kinetic energy and trap confinement V ext ( r ) = M ( ω x x + ω y y + ω z z ) /
2, where M is themass. The contact interaction strength g = 4 π ¯ h a s /M isgiven by the scattering length a s . The dipolar mean fieldpotential is Φ dip = (cid:82) d r (cid:48) V dd ( r − r (cid:48) ) | ψ ( r (cid:48) , t ) | where V dd ( r ) = g dd π − ϑ | r | is the dipolar interaction foraligned dipoles with an angle ϑ between r and the mag-netic field axis. The strength of the dipolar interaction isgiven by the parameter g dd = 4 π ¯ h a dd /M characterizedby the dipolar length a dd = µ µ m M/ (12 π ¯ h ) andthe magnetic moment µ m . Furthermore the quantumfluctuations within the local density approximationare given by the quantity g qf = 32 / (3 √ π ) ga / Q ( (cid:15) dd ),where (cid:15) dd = g dd /g = a dd /a s is the relative dipolarstrength. In our simulations, we use the approximation Q ( (cid:15) dd ) (cid:39) (cid:15) [10, 50, 61, 63, 64].We use the Bogoliubov-de Gennes (BdG) formalism asdescribed in Ref. [21] and linearly expand the wavefunc-tion ψ ( r , t ) = ψ ( r ) + λ [ u ( r ) e − iωt + v ∗ ( r ) e iωt ] e − iµt/ ¯ h around the ground state ψ with the chemical potential µ . This ansatz together with equation (2) leads to a sys-tem of linear equations that can be expressed in matrixform. See Refs. [21, 22, 64] for the complete form of theBdG matrix representation. We numerically solve theseequations to obtain the modes u and v corresponding tothe lowest excitation energies ¯ hω . From the solutions u and v we obtain the density fluctuation patterns δn ∝ ( u + v ) ψ [21, 26] discussed in the main text (see0Fig. 1). We obtain a discrete spectrum of elementaryexcitations because of the finite-sized system.In Fig. S4 we show a wider view of the excitation spec-trum discussed in the main text (Fig. 1) with higher exci-tation energies and larger scattering lengths. The shownexcitation spectrums of a 2D and a 1D BEC are simi-lar [21], but the inclusion of more modes complicates the2D BEC spectrum. Due to the quasi-periodic boundaryconditions along the angular direction, multiple (angu-lar) roton modes exist, which all soften near the transi-tion point. Which mode softens fastest is determined byan interplay between the interaction strengths and theaspect ratio [37].Generally, modes in the same m -subspace can cou-ple to each other [20, 26, 37, 65, 66]. While the angu-lar roton modes are softening, they couple to the BECphonon modes in the same m -subspace. This can lead toavoided level crossings and hybridization between angu-lar roton modes and BEC multipole (dipole, quadrupole,hexapole, ...) modes [20]. As an example we show aninset of the avoided crossing between the m = 2 angu- lar roton mode and the quadrupole mode of the BEC.The hybridization between the angular roton mode andthe quadrupole mode leads to an additional radial nodein the mode pattern of the quadrupole mode after theanticrossing ( a s < ∼ . m = 0) shows such an avoided crossing with thebreathing mode at around a s (cid:39) . a and the m = 3angular roton mode with the hexapole mode at around a s (cid:39) a . Here, the deformation of the branches islarger compared to the avoided crossing of the m = 2 an-gular roton mode, suggesting a larger coupling betweenthe two participating modes.We note that labelling the angular roton modes withtheir number of nodal lines m is only a good choice fora circularly symmetric ground state. With an increasingasymmetry in the radial trap aspect ratio this degener-acy is lifted and m ceases to uniquely classify the modepatterns. The splitting of the m = 1 mode close to the in-stability is a first signature of this effect, even though thedeviation from a perfectly circular trap is about 0 ..