Observation of vortices and vortex stripes in a dipolar condensate

Quantized vortices are a prototypical feature of superfluidity that have been observed in multiple quantum gas experiments. But the occurrence of vortices in dipolar quantum gases—a class of ultracold gases characterized by long-range anisotropic interactions—has not been reported yet. Here we exploit the anisotropic nature of the dipole–dipole interaction of a dysprosium Bose–Einstein condensate to induce angular symmetry breaking in an otherwise cylindrically symmetric pancake-shaped trap. Tilting the magnetic field towards the radial plane deforms the cloud into an ellipsoid, which is then set into rotation. At stirring frequencies approaching the radial trap frequency, we observe the generation of dynamically unstable surface excitations, which cause angular momentum to be pumped into the system through vortices. Under continuous rotation, the vortices arrange into a stripe configuration along the field, in close agreement with numerical simulations.

Quantized vortices are a prototypical feature of superfluidity that have been observed in multiple quantum gas experiments. But the occurrence of vortices in dipolar quantum gases-a class of ultracold gases characterized by long-range anisotropic interactions-has not been reported yet. Here we exploit the anisotropic nature of the dipole-dipole interaction of a dysprosium Bose-Einstein condensate to induce angular symmetry breaking in an otherwise cylindrically symmetric pancake-shaped trap. Tilting the magnetic field towards the radial plane deforms the cloud into an ellipsoid, which is then set into rotation. At stirring frequencies approaching the radial trap frequency, we observe the generation of dynamically unstable surface excitations, which cause angular momentum to be pumped into the system through vortices. Under continuous rotation, the vortices arrange into a stripe configuration along the field, in close agreement with numerical simulations.
Since the first experiments on gaseous Bose-Einstein condensates (BECs), the observation of quantized vortices has been considered the most fundamental and defining signature of the superfluid nature of such systems. Their very existence sets a unifying concept encompassing a variety of quantum fluids from liquid helium 1 to the core of neutron stars 2 and from superconductors 3 to quantum fluids of light 4 . Their classical counterparts have as well fascinated scientists from different epochs and fields as vortices are found in many scales of physical systems, from tornadoes in the atmosphere to ferrohydrodynamics.
In the quantum realm, a quantized vortex may emerge as a unique response of a superfluid to rotation. It can be understood as a type of topologically protected singularity with a 2π phase winding that preserves the single-valuedness of the superfluid wave function and the irrotational nature of its velocity field. In contact-interacting BECs, vortical singularities have been observed experimentally in the form of single vortices 5,6 , vortex-antivortex pairs 7 , solitonic vortices 8,9 , vortex rings 10 and vortex lattices 6,11 using a number of different techniques. Moreover, vortices play a fundamental role in the description of the Berezinskii-Kosterlitz-Thouless transition in two-dimensional (2D) systems 12 , as well as in the evolution of quantum turbulence 13,14 , and have been observed in interacting Fermi gases along the Bose-Einstein condensate to Bardeen-Cooper-Schrieffer crossover 8,15 .
Recently, a new class of ultracold quantum gases are being created in various laboratories around the world, using strongly magnetic lanthanide atoms 16,17 . Such a system, providing a quantum analogue of classical ferrofluids, enables access to the physics of dipolar BECs, in which atoms feature a strong long-range anisotropic dipole-dipole interaction (DDI) 18,19 on top of the traditional contact-type isotropic one. This intriguing platform provided the key to observe, for example, extended Bose-Hubbard dynamics 20 , roton excitations [21][22][23] , the quantum version of the Rosensweig instability 24 and supersolid states of matter [25][26][27][28] , and is foreseen to host novel phenomena for quantum simulation and metrology 18,19 .
The dipolar interaction is predicted to also intimately change the properties of vortices in quantum gases 29 . For instance, theoretical works predict single vortices to exhibit an elliptic-shaped core for a quasi-2D setting with in-plane dipole orientation [30][31][32][33] or the presence of density oscillations around the vortex core induced by the roton Article https://doi.org/10.1038/s41567-022-01793-8 smoking gun of superfluidity in supersolid states [41][42][43] . However, despite these intriguing predictions, vortices in dipolar quantum gases have not been observed until now.
This Article presents the experimental realization of quantized vortices in a dipolar BEC of highly magnetic dysprosium (Dy) atoms. Following a method proposed in ref. 40 , extended to arbitrary magnetic-field angles in ref. 44 , we show that the many-body phenomenon of magnetostriction 45 , genuinely arising from the anisotropic DDI among atoms, provides a natural route to rotate the systems and nucleate vortices in a dipolar BEC. We carry out studies on the dynamics of the vortex formation, which agree very well with our theoretical predictions. Finally, we observe one of the earliest predictions for vortices in dipolar BECs: the formation of vortex stripes in the system.
In non-dipolar gases, quantized vortices have been produced using several conceptually different techniques, for instance, by rotating non-symmetric optical 6,11 or magnetic 46 potentials, by rapidly shaking the gas 14 , by traversing it with obstacles with large enough velocity 7,47 , by rapidly cooling the gas across the BEC phase transition 48,49 or by directly imprinting the vortex phase pattern 50 . Dipolar quantum gases, while able to form vortices with these same standard procedures 29 , also offer unique opportunities that have no counterpart in contact-interacting gases. Crucially, the DDI gives rise to the phenomenon of magnetostriction in position space 45 . When dipolar BECs are polarized by an external magnetic field B-defining the dipole orientation-the DDI causes an elongation of the cloud along the polarization direction. This is a direct consequence of the system tendency to favour head-to-tail dipole configurations, which effectively reduces the interaction energy 19 .
Such a magnetostrictive effect provides a simple method to induce an elliptic effective potential and drive rotation with a single control parameter. This modification of the effective potential is shown in Fig. 1a for a BEC in an oblate trap with cylindrical symmetry about the z axis. While a non-dipolar BEC takes the same shape as the confining trap ( Fig. 1a(i)), introducing dipolar interactions with polarization axis along z stretches the cloud along this axis yet maintains cylindrical symmetry ( Fig. 1a(ii)). Tilting the magnetic field leads to a breaking of the cylindrical symmetry, resulting in an ellipsoidal deformation of the cloud shape, as seen from the density projection onto the x-y plane ( Fig. 1a(iii)). Finally, under continuous rotation of the magnetic field, which we coin 'magnetostirring', the condensate is predicted to rotate ( Fig. 1a(iv)). This unique approach to stir a dipolar condensate can eventually lead to the nucleation of vortices 40,44 , realizing genuinely interaction-driven vorticity through many-body phenomena.
We explore this protocol using a dipolar BEC of 162 Dy atoms. We create the BEC similar to our previous work 51 with the distinction that here the magnetic-field unit vector, B , is kept tilted at an angle of θ = 35° with respect to the z axis both during evaporative cooling and magnetostirring ( Fig. 1a(iii) and Methods). After preparation, the sample contains about 2 × 10 4 condensed atoms confined within a cylindrically symmetric optical dipole trap (ODT) with typical radial and axial trap frequencies (ω ⊥ , ω z ) = 2π × [50.8(2), 140(1)] Hz. Here, before stirring, the magnetostriction is expected from simulations to increase the cloud aspect ratio (AR) in the horizontal plane from 1 up to 1.03, whereas the trap anisotropy is negligible. We use a vertical (z) absorption imaging to probe the radial (x,y) atomic distribution after a short time-of-flight (TOF) expansion of 3 ms. The atom number is instead measured using horizontal absorption imaging with a TOF of 26 ms.
Similarly to a rotation of a bucket containing superfluid helium or of a smoothly deformed ODT for non-dipolar BECs, magnetostirring is predicted to transfer angular momentum into a dipolar BEC 40,44 . In response to such an imposed rotation, the shape of an irrotational cloud is expected to elongate with an amplitude that increases with the rotation frequency Ω. This phenomenon is clearly visible in our experiments, as shown in Fig. 1b. Here we first revolve the tilted B around the z axis with a linearly increasing rotation frequency (Ω = 2 × 50 Hz s −1 ) and observe that the dipolar BEC starts to rotate minimum in the dispersion relation [30][31][32][33][34] . For vortex pairs, the anisotropic DDI is expected to alter the lifetime and dynamics 33,35 and can even suppress vortex-antivortex annihilation 33 . These interaction properties are predicted to give rise to a vortex lattice structure that can follow a triangular pattern 30,34 , as is typical for non-dipolar BECs 11 , or a square lattice for attractive or zero contact interactions [36][37][38] when the DDI is isotropic (dipoles aligned with the rotation axis). A very striking consequence of the dipoles tilted towards the plane is the formation of vortex stripes 30,39,40 . Moreover, vortices could provide an unambiguous Aspect ratio i ii iii iv Article https://doi.org/10.1038/s41567-022-01793-8 at the same angular speed as the field and deforms with increasing elongation (Fig. 1b, top). We then stop the adiabatic ramp at a given value of Ω and probe the system under continuous rotation. We now find that the cloud continues rotating in the radial plane with an almost constant shape (Fig. 1b, bottom). Note that B is held constant at 5.333(5) G, where we estimate a contact scattering length of about a s = 111a 0 , where a 0 is the Bohr radius (Methods). We further explore the response of our dipolar BEC to magnetostirring by repeating the measurements in Fig. 1b (top), but stopping the ramp at different final values of Ω. The maximum value used for Ω approaches ω ⊥ , corresponding to a ramp duration of 1 s. We quantify the cloud elongation in terms of the aspect ratio AR = σ max /σ min , where the cloud widths σ max and σ min are extracted by fitting a rotated 2D Gaussian function to the density profiles. Figure 1c summarizes our results. We observe that initially the AR slightly deviates from 1 due to magnetostriction. It then slowly grows with increasing Ω, until a rapid increase at around 0.6ω ⊥ occurs, as this allows the angular momentum to increase, which decreases the energy in the rotating frame 52 . Suddenly, at a critical rotation frequency Ω c ≈ 0.74ω ⊥ , the AR abruptly collapses back to AR ≈ 1, showing how the superfluid irrotational nature competes with the imposed rotation. This critical frequency is close to the value found in non-dipolar gases with a rotating elliptical harmonic trap, associated with a resonance at the quadrupole frequency 53 .
To substantiate our observation, we perform numerical simulations of the zero-temperature extended Gross-Pitaevskii equation (eGPE) 54 (Methods). Quantum and thermal fluctuations are added to the initial states, which are important to seed the dynamic instabilities once they emerge at large enough Ω; see later discussion. The lines in Fig. 1c show our results. The dashed line is obtained through the same procedure as the experiment, whereas for the solid line, we halve the ramp rate, spending more time at each frequency. Both ramp procedures show quantitatively the same behaviour up to Ω = 0.8ω ⊥ and are in excellent agreement with the experimental results. The stability of the 1 s ramp exceeds the experimentally observed critical frequency. We partly attribute this discrepancy to asymmetries of the rotation in the experiment that are not present in the simulations, which may lead to an effective speed-up of the dynamical instabilities. However, in all cases, the AR rapidly decreases to about 1.
The growing AR and subsequent collapse to 1 is a signature of the dynamical instability of surface modes, known for being an important mechanism for seeding vortices and allowing them to penetrate into the high-density regions of rotated BECs 52,53,55 , as also predicted for our dipolar system 40 . To search for quantum vortices in our system, we perform a new investigation where we directly set Ω close to Ω c , aiming to trigger the instability at an earlier time when more atoms are condensed. We then hold the magnetic-field rotation fixed at this constant frequency for a time t Ω . As shown in Fig. 2 (bottom), the cloud rapidly elongates, and the density starts to exhibit a spiral pattern, emanating from the tips of the ellipsoid. As early as t Ω = 314 ms, clear holes are observed in the density profile, forming in the density halo around the centre, the first clear indication of vortices in a dipolar gas. These vortices, initially nucleated at the edge of the sample, persist as we continue to stir and eventually migrate towards the central (high-density) region. Vortices are still visible in the experiment after 1 s of magnetostirring, although our atom number decreases throughout this procedure. Article https://doi.org/10.1038/s41567-022-01793-8 Our observations bear a remarkable resemblance to the simulations; Fig. 2 (top) shows the in situ column densities. Taking a fixed atom number of N = 8,000, but otherwise repeating the experimental sequence, we observe many similar features. In the first 100 ms, the system elongates, consistent with Fig. 1, and a spiral density pattern appears before the instability, forming two arms that are filled with vortices close to the central density. Next, turbulent dynamics ensue as the density surface goes unstable and vortices emerge in the central high-density region. For this scattering length and atom number, the relaxation timescale to a stable vortex lattice is longer than the experimentally available (see Extended Data Fig. 3 for more images from this dataset). Note that at angles θ deeper into the plane, more atoms align head-to-tail in the loose radial confinement direction. Thus, when performing the rotation procedure, we find that the BEC is resilient to instability on the timescales of the experiment. The observed evolution of the system under constant rotation shows some concurrence between the appearance of vortices in the absorption images and the formation of a round density pattern in the radial plane with AR ≈ 1 (Fig. 2). Note that the drop in AR observed in Fig. 1 is concurrent with the creation of vortices, but they reside in the low-density regions at this time, and we do not see them. To study this dynamical evolution in more detail, we adopt an analysis protocol for both the experiment and theory that allows us to quantitatively track the evolution of the average number of vortices, v (Methods). The result is shown in Fig. 3a. In brief, for each single image (Fig. 3a, left), we create a blurred reference image by applying a 2D Gaussian filter 56,57 . We then calculate the difference between each single image (Fig. 3a, left) and the corresponding reference (Fig. 3a, middle) to obtain the residual image (Fig. 3a, right), from which we count v by finding local minima below a certain threshold.
For the experimental density profiles, which are affected by both the limited resolution of the imaging system and the weak contrast in the low-density zones (halo) where the vortices initially nest, we expect v to be underestimated relative to the true value and the number expected by theory. However, to carry out a quantitative comparison with the simulations, we apply a blurring filter and add noise to the latter that mimics the actual resolution in the experiment (Methods). Figure 3b shows both the evolution of v and cloud AR as a function of rotation time, t Ω . Solid lines are the results from the eGPE simulations without any adjustable parameters. For t Ω < 200 ms, v is below 1, where vortices, if present, are at the edge of the cloud. For longer times, v increases and saturates to an average value of about three and a maximum of six vortices (see Fig. 3a for an example of five vortices). The observed saturation might be due to the decreased visibility and to the atom-loss-induced shrinking of the BEC size, which is not accounted for in the theory. We also compare the course of the average vortex number with the AR of the cloud. After initial large oscillations, due to the sudden jump in rotation frequency, the AR declines towards ~1 (ref. 53 ). This happens as the vortex number simultaneously increases.
One fascinating prediction with vortices in a strongly dipolar gas under the influence of a rotating magnetic field relates to the structure of the resulting vortex lattice. Due to magnetostriction and the anisotropic vortex cores, the resulting vortex configuration is also anisotropic, producing a stripe phase in the strongly dipolar regime 29,30 , instead of the usual triangular lattice in non-dipolar BECs 6 . The ground state stripe lattice solution for our parameters is shown in Fig. 4a, with a cloud AR = 1.08. In the vortex stripe phase, vertical planes of high-density regions, parallel to the magnetic field, alternate with low-density ones, which host rows of vertical vortex filaments. Such a configuration promotes head-to-tail dipolar attraction within the high-density ridges, and this acts to lower the energy. It should be noted that these states are distinct from the oscillating vortex sheets states, which appear after squeezing a triangular vortex lattice 58 . To explore this prediction, we perform two new surveys. First, we slightly reduce the magnetic-field value, reducing the scattering length to a s ≈ 109a 0 and hence making the system relatively more dipolar. We magnetostir the BEC at a constant rotation frequency Ω = 0.75ω ⊥ for 500 ms, but during TOF, we stop the magnetic-field rotation and keep it in place at θ = 35°. The stripe structure is revealed in Fig. 4b (left) for a single experimental run and is clearly visible in the residual image (Fig. 4b, middle) where the vortices align along three stripes. The spatial structure of the residual image can be assessed through the absolute value of 2D Fourier transform (FT). After taking the FT of each residual image, we then average the result (Fig. 4b, right), finding a clear peak at the wave number k of the inter-stripe spacing. This shows that the stripe spatial structure survives the averaging, implying that the majority of images show stripes with the same spacing, and they also have the same orientation as set by the magnetic field, as evidenced by the example images shown in the right of Fig. 4. Note that these observations do not rely on our ability to resolve individual vortices, as the stripes are an ensemble effect of many aligned vortices. In fact, by comparing with the numerical simulations of the dynamical procedure (Fig. 4c), we expect there are more vortices than detected here that fill in the stripes, forging out this structure. In general, our simulations show that the stripes appear faster when the scattering length is lower and when the atom number is larger. In the long time limit of the scenario presented in Fig. 2, we expect the stationary solution to also be the stripe state, but this is not observable on our timescales.
Remarkably, the stripe structure washes out when we subsequently tilt the magnetic-field orientation to θ = 0° (parallel to the trap symmetry axis), as shown in Fig. 4d (left). Here, after 600 ms of magnetostirring, we add another step in which we spiral up the magnetic field to θ = 0° (with Ω fixed) over 100 ms, before imaging. Under these conditions, all vortex properties are again isotropic within the plane. The non-equilibrium positioning of the vortices is arbitrary, and if we average the FT of the residuals directly, we observe a homogeneous ring in the average FT (Fig. 4d, right). Also, this behaviour is confirmed by the simulations, as shown in Fig. 4e. The vortices survive long after the magnetostirring has stopped (not shown), due to their topological protection.
By exploiting magnetostirring-a novel, robust method of generating angular momentum-we have observed quantized vortices in a dipolar quantum gas and the appearance of the vortex stripe configuration. Future works will focus on investigations of the individual vortex shape and behaviour, such as the anisotropic nature of the vortex cores for in-plane magnetic fields [30][31][32][33] , the interplay between the vortex and roton excitations [30][31][32][33][34] and exotic vortex patterns such as square lattices 29 , and investigations into anisotropic turbulence 59 . This work also opens the door to studying more complex matter under rotation, such as dipolar droplets 60-62 and supersolid states [41][42][43]51 . Such proposals will be challenging due to the intricate density patterns 63 ; however, such observations would provide conclusive evidence of superfluidity in supersolids. Rotating the magnetic field at frequencies far larger than the radial trap frequencies, but smaller than the Larmor frequency, has been observed to tune the sign and magnitude of the dipole-dipole interaction 64 -a method also employed in nuclear magnetic resonance spectroscopybut there remain open questions on the stability of this procedure 65,66 , which if rectifiable would unlock new research directions 64 . Other vortex generation methods, such as thermally activated pairs in quasi-two dimensions to assess the Berezinskii-Kosterlitz-Thouless transition and stochastically generated vortex tangles through temperature quenches to assess the Kibble-Zurek mechanism, remain unexplored in dipolar gases 29 . The technique introduced here is also applicable to a wide range of systems governed by long-range interactions through the manipulation of magnetic or electric fields.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41567-022-01793-8.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Experimental procedure
We prepare an ultracold gas of 162 Dy atoms in an ODT. Three 1,064 nm laser beams, overlapping at their foci, form the ODT. The experimental procedure to BEC is similar to the one followed in our previous work 51 , but the magnetic-field unit vector, B , is tilted by an angle of θ = 35° with respect to the z-axis during the whole sequence. After preparation, the sample contains about 2 × 10 4 condensed atoms. The corresponding trap frequencies are typically (ω ⊥ , ω z ) = 2π × [50.8(2), 140(1)] Hz. For all our measurements, the deviation of the trap AR in the x-y plane AR trap = ω y /ω x from 1 is always smaller than 0.6%. We evaporate the atoms at B = 5.423 (5) G and jump to the final magnetic-field value during the last evaporation ramp. After the preparation of the BEC, the magnetic field is rotated as described in the next section. We use standard absorption imaging to record the atomic distribution. We probe the vortices using the vertical imaging taken along the axis of rotation (z), for which the dark spots within the condensate correspond to the cores of individual vortices. The vertical images are taken with a short TOF of 3 ms and a pulse duration of 3-4 μs. For the data in Figs. 1-3, we let the magnetic-field spinning during TOF, whereas for Fig. 4, we use a static field orientation.

Control of the magnetic field
Calibration. Three pairs of coils-each oriented along a primary axis in the laboratory frame-enable the creation of a homogeneous field with arbitrary orientations. The absolute magnetic-field value B of each pair of coils is independently calibrated using radio frequency (RF) spectroscopy. The RF drives transitions to excited Zeeman states, leading to a resonant dip in the atom number. The long-term stability-measured via the peak position of the RF resonance over the course of several days-is on the order of ΔB = ±1 mG, while shot-to-shot fluctuations, measured via the width of the RF resonance for a single calibration set, is ΔB = ±5 mG.

Rotation.
We drive the rotation of the magnetic field by sinusoidally modulating the magnetic-field value components B x and B y with a phase difference of 90° between them. As we want to keep the absolute magnetic-field value B constant during rotation, we measure it for various values of the azimuthal angle ϕ and fixed θ = 35° by performing Feshbach loss spectroscopy around 5.1 G. We find an average shift of B of about 10 mG from the θ = 0° case, which we take into account.
We also find small deviations as a function of ϕ of ΔB < 20 mG, which might appear due to slightly non-orthogonal alignment of the magnetic fields. We did not correct these deviations for the sake of simplicity.

Scattering length
The scattering length in 162 Dy is currently not known with large accuracy [67][68][69][70] . To estimate the scattering length in the small magnetic-field range around B = 5.3 G, relevant to this work, we use the well-known relation a s = a bg ∏ i [1 − ΔB i /(B − B 0,i )] (ref. 71 ), where B 0,i and ΔB i are the centre position and the width of the i-th feature of the Feshbach loss measurement reported in ref. 70 , respectively. The value of the background scattering length, a bg , is empirically fixed by measuring the magnetic-field value at which the supersolid transition occurs and comparing it with the corresponding critical a s predicted from simulations. Such an approach leads to a s = 111(9)a 0 at B = 5.333 G. Extended Data Fig. 1 shows the resulting scattering lengths for the relevant magnetic fields. Although such an approach gives very good agreement between theory and experiments, future works on a precise determination of a s , similar to the one achieved with erbium 72 , would be desirable.

Magnetostirring
Tilting the magnetic-field vector B away from the symmetry axis of our cylindrical trap leads to an ellipsoidal deformation of the cloud 45 and therefore to a breaking of the cylindrical symmetry. This allows for the transfer of angular momentum to the sample by rotating the magnetic field (magnetostirring). In all our measurements, we use a B tilted with respect to the z axis by 35° and a constant value B. That value is B = 5.333(5) G for the surveys in Figs. 1-3 and B = 5.323(5 )G for Fig. 4. For these parameters, the deformed magnetostricted AR of the cloud is AR − 1 = 0.03. For all our measurements, the measured trap AR trap − 1 < 0.006 is much smaller than the deformation due to magnetostriction. Additionally, we have confirmed with simulations that even with trap asymmetries of up to 10%, for example, (ω x , ω y ) = (55, 50) Hz, this procedure can still generate vortices in a lattice configuration.
At the scattering lengths considered in this work, 35° is an optimal choice to see the vortices within ~500 ms of rotation and anisotropic enough to observe the stripe phase. From the simulations, we find that tilt angles smaller than 35° increase the timescale to vortex nucleation. Similarly, tilting the angle further into the plane increases the number of atoms that are aligned head to tail, making the dipolar interaction dominantly attractive. This attractive force holds the condensate together during the rotation, also increasing the time to vortex nucleation. From the experimental side, increasing the tilt angle reduces the contrast of the absorption imaging, since the magnetic field is not parallel to the imaging axis. As the TOF is only 3 ms, we do not rotate up the magnetic field before imaging to avoid undesired effects, such as losing the anisotropy given by θ ≠ 0°.
Adiabatic frequency ramp. We employ different magnetic-field rotation sequences for the different datasets. For the dataset of Fig. 1c, the rotation frequency of the magnetic field is linearly increased to different final values at a speed of Ω = 2 × 50 Hz s −1 and for a duration of tΩ = 0−1 s. The ramp time is much longer than the period of the rotation Ω −1 for higher rotation frequencies Ω ≳ Ω c , and, therefore, the ramp is adiabatic for the regimes considered, until the onset of dynamical instabilities. After the ramp, the magnetic-field direction is rotated at the target rotation frequency Ω for one final period (as shown in Fig. 1b). We sample ten different final magnetic-field angles during this last rotation, measuring the corresponding AR and averaging the result to remove any potential biases due to latent trap anisotropies. Each data point is then obtained with eight to ten experimental runs.

Constant rotation frequency.
For the dataset of Figs. 2, 3 and 4b, we directly start to rotate at the final rotation frequency Ω without any acceleration phase. The magnetic field is then rotated for a variable time t Ω , after which the atoms are released from the trap and a vertical image is taken. Spiral up magnetic field. For the dataset of Fig. 4d, we employ a similar sequence as described above. However, after constantly rotating the magnetic field at Ω = 0.75ω ⊥ , the magnetic field is spiralled up in 100 ms to θ = 0° by linearly reducing θ while continuing rotating. Afterwards, the atoms are released from the trap and a vertical image is taken.

Theoretical model
We employ an eGP formalism to model our experimental set-up. In this scheme, the inter-particle interactions are described by the two-body pseudo-potential with δ(r) being the Kronecker delta function and r = (x, y, z). The first term describes the short-range interactions governed by the s-wave scattering length a s , with Planck's constant ℏ and particle mass m. The second term represents the anisotropic and long-range dipoledipole interactions, characterized by dipole length a dd = μ 0 μ 2 m m/12 ℏ 2 , with magnetic moment μ m and vacuum permeability μ 0 . We always consider 162 Dy, such that a dd = 129.2a 0 , where a 0 is the Bohr radius.
Beyond-mean-field effects are treated through the inclusion of a Lee-Huang-Yang correction term 73 being the auxiliary function, and the relative dipole strength is given by ϵ dd = a dd /a s . Finally, the full eGPE then reads 54,[74][75][76] iℏ where ω x,y,z are the harmonic trap frequencies. The wave function Ψ is normalized to the total atom number N = ∫d 3 x|Ψ| 2 . The stationary solution for Fig. 4a is found through the imaginary time procedure in the rotating frame, introducing the usual angular momentum operator ΩL z into equation (4). The initial state Ψ(x, 0) of the real-time simulations is always obtained by adding non-interacting noise to the ground state Ψ 0 (x). Given the single-particle eigenstates ϕ n and the complex Gaussian random variables α n sampled with ⟨|α n | 2 ⟩ = (e n /k B T − 1) −1 + 1 2 for a temperature T = 20 nK and Boltzmann's constant k B , the initial state can be described as Ψ (x, 0) = Ψ 0 (x) + ∑ n ′ α n ϕ n (x) , where the sum is restricted only to the modes with ϵ n ≤ 2k B T (ref. 77 ). Throughout, the density images are presented in situ, with a scaling factor to account for the 3 ms TOF for the experimental images.
To obtain the average residual FT images for Fig. 4c,e, we first Fourier transform 115 frames from the simulation between 700 ms and 1.1 s in the rotating frame before averaging the result.

Atom number
Extended Data Fig. 2 shows the condensed atom number N c for the measurement with an adiabatic ramp of the magnetic-field rotational velocity (Ω = 2 × 50 Hz s −1 ), corresponding to the data of Fig. 1c. Three-body losses are negligible in the low-density BEC phase, with losses probably coming from imperfections in the rotation procedure and heating. To extract the atom number, we use the horizontal imaging with 26 ms of TOF. About 3 ms before flashing the imaging resonant light to the atoms, we rotate the magnetic field in the imaging plane and perform standard absorption imaging. From the absorption images, we extract N c from a bimodal fit up to 700 ms. At later times, the system undergoes a dynamic instability (see discussion in the main text), and the density profile deviates from a simple bimodal distribution. During the observation time, we see a slight decrease of N c , and for our theory simulations, we use a constant atom number of N c = 15,000. Note that in all following datasets, in which we abruptly accelerate the magnetic-field rotation to the desired final velocity, we observe a faster decay, and our simulations are performed with either N c = 8,000 or N c = 10,000.

Vortex detection
Vortex detection algorithm. Since vortices appear as dark holes in the density profile of a BEC, which would otherwise have a smooth profile, our approach to extract the number of vortices is to look at deviations between the image and an unmodulated reference image. To extract the vortex number from the raw images, we proceed as follows.
First, we prepare the image n img , the reference image n ref and the residual image n res . The image is normalized such that the maximum density max(n img ) = 1. We create the reference image by blurring the image via applying a 2D Gaussian filter with σ = 5 pixel, corresponding to about 2.1 μm. This blurring smoothens any structure on the lengthscale of the filter width; therefore, any holes in the density profile wash out. We then normalize the atom number of the reference to be the same as from the image N ref = ∫∫n ref ≐ N img = ∫∫n img . The residual image is calculated as the difference between the image and the reference n res = n img − n ref . We additionally mask the region where the density of the reference is below a certain threshold (n res = 0, where n ref ≤ 0.1).
Second, we identify local minima in the residual image and determine whether they are connected to vortices. For this, we create a list of local minima (x min , y min ), defined by the condition that the pixel density n res (x min , y min ) is lower than of all surrounding pixels. Then we remove minima with density values above zero n res (x min , y min ) ≥ 0 or which are within one pixel distance of the mask border. Now we determine a local contrast for each minimum by calculating the difference between its central density value and the mean of the density values ± 2 pixel values away from it n con (x min , y min ) = n res (x min , y min ) − mean(n res (x min ± 2 px, y min ± 2 px)) , and remove minima above a certain threshold n con > −0.11. As a last step, we check the distance d between all remaining minima to avoid double counting of minima too close to each other. In case d is below the threshold d < 5 pixel, the minimum with the higher residual density value n res is discarded.
Preparation of theory density profiles. For the direct comparison of the vortex number shown in Fig. 3b, we apply additional steps to the density profiles obtained from theory. First, we reduce the resolution by a 2 × 2 binning to make the pixel size of the theory density profiles n theo img essentially the same as for the experimental images (sizes are within 5%). After normalizing to max(n theo img ) = 1, we apply Gaussian white noise with zero mean and a variance of 0.01 to each pixel, recreating the noise pattern from empty regions of experimental images. Then we blur the image using a 2D Gaussian filter with σ = 1 pixel (~0.42 μm); this recreates the same resolution condition as our experimental set-up. The resulting density profile is taken as the input image for the vortex detection algorithm described above.
Benchmarking the vortex detection algorithm. As the vortex positions for the simulation images are known a priori due to the available phase map, we can derive the fidelity of the vortex detection algorithm for simulation data. For the theory data shown in Fig. 3b in the time frame between 600 and 700 ms, the average detected vortex number in the simulated density profiles (applying the preparation scheme described above) is about 9, while the real number of vortices present in the same area of the image is about 33 in average. This mismatch is explained by the conservative choice of the thresholds for vortex detection together with the added noise, which results in only counting clear density dips as vortices, throwing out many vortices in the low-density region. This conservative choice of thresholds on the other hand leads to a very high fidelity of >97%, where we define the fidelity as the percentage of detected vortices that correspond to actual present vortices in the data. For raw simulation data (without resolution reduction, added noise and blurring), the vortex detection algorithm would detect up to 80% of the vortices present with a fidelity of >95%.
Note that for the visualization of the vortex positions for Fig. 4, we slightly increased the local threshold n con > −0.08 and decreased the minimum distance between vortices d < 3 pixel, which increases the overall number of vortices detected. For the density distributions obtained from theory, we additionally omit the resolution reduction, addition of noise and blurring steps.