Rotational tuning of the dipole-dipole interaction in a Bose gas of magnetic atoms

TL;DR

This paper explores how rotating magnetic fields can tune dipole-dipole interactions in a Bose-Einstein condensate of magnetic atoms, analyzing the emergence of instabilities and their dependence on rotation frequency.

Contribution

It provides simulations of realistic tuning sequences and characterizes the development of dynamic instabilities, offering insights for experimental feasibility.

Findings

Instabilities develop more slowly at higher rotation frequencies

Tuned dipole-dipole interactions are experimentally feasible

Dynamic instabilities can be characterized through simulations

Abstract

We investigate the dynamics of a Bose-Einstein condensate of magnetic atoms in which the dipoles are rotated by an external magnetic field. The time-averaged dipole-dipole interaction between the atoms is effectively tuned by this rotation, however recent experimental and theoretical developments show that dynamic instabilities emerge that may cause heating. We present simulations of a realistic tuning sequence in this system, and characterize the system behavior and the emergence of instabilities. Our results indicate that the instabilities develop more slowly as the rotation frequency increases, and indicate that experiments with tuned dipole-dipole interactions should be feasible.

Tables

Table 1. Table 1: Details of Fig.  5 cases. Here N 𝑁 N is the expected average total number of atoms in the field ψ 𝜓 \psi and E / N 𝐸 𝑁 E/N is the initial average energy per particle including noise up to the cutoff energy. The chemical potential μ 𝜇 \mu is for the initial state ψ 0 subscript 𝜓 0 \psi_{0} . The instability time is the average over five trajectories with Ω / 2 ​ π = 3 ​ kHz Ω 2 𝜋 3 kHz \Omega/2\pi=3\>\mathrm{kHz} .

$f$	${ϵ}_{dd}$	$N_c$	$N$	$E/Nh$	$\mu /h$	$t_{\mathrm{inst}}$
(Hz)				(Hz)	(Hz)	(ms)
50	0.50	1000	1057	170	224	$>300$
50	0.75	1000	1034	149	192	$>300$
50	0.90	1000	1028	140	178	$>300$
70	0.50	1000	1060	251	333	139
70	0.75	1000	1042	219	284	175
70	0.90	1000	1028	205	263	108
50	0.50	2500	2642	231	313	103
50	0.75	2500	2588	198	264	115
50	0.90	2500	2568	184	243	90
70	0.50	2500	2677	344	467	52
70	0.75	2500	2607	294	393	42
70	0.90	2500	2582	272	361	38
50	0.50	5000	5326	297	407	54
50	0.75	5000	5194	252	342	42
50	0.90	5000	5142	232	312	36
70	0.50	5000	5394	443	608	27
70	0.75	5000	5233	374	510	12
70	0.90	5000	5177	344	466	12
50	0.50	10000	10730	386	532	28
50	0.75	10000	10429	324	445	13
50	0.90	10000	10322	297	406	11
70	0.50	10000	10884	576	796	11
70	0.75	10000	10523	483	665	9
70	0.90	10000	10390	442	606	8

Equations

$$\mathbf{e}(t)=(\hat{x}\cos \Omega t+\hat{y}\sin \Omega t)\sin \phi (t)+\hat{z}\cos \phi (t).$$

$${\mathcal{L}}_{\mathrm{GP}}=H_{sp}+g_s|\psi {|}^2+{\Phi }_D(\mathbf{x},t),$$

$$H_{sp}=-\frac{{\hslash }^2}{2m}{\nabla }^2+\frac{m}{2}({\omega }_x^2{x}^2+{\omega }_y^2{y}^2+{\omega }_z^2{z}^2),$$

$$U_{dd}(\mathbf{r},t)$$

$$E[\psi ]=\int d\mathbf{x}{\psi }^{\ast }\left[H_{sp}+\frac{1}{2}{g}_s|\psi {|}^2+\frac{1}{2}{\Phi }_D(\mathbf{x},t)\right]\psi .$$

$${\bar{U}}_{dd}(\mathbf{r})$$

$${\bar{g}}_{dd}(\phi )$$

$$\bar{\mu }\bar{\psi }={\bar{\mathcal{L}}}_{\mathrm{GP}}\bar{\psi },$$

$$E_K=-\frac{{\hslash }^2}{2m}\int d\mathbf{x}{\psi }^{\ast }{\nabla }^2\psi .$$

$$\psi (\mathbf{x},0)={\psi }_0(\mathbf{x})+\sum _j{}^{{}^{\prime }}{\alpha }_j{\varphi }_j(\mathbf{x}),$$

$${\tilde{U}}_{dd}^R(\mathbf{k},t)={\tilde{U}}_{dd}(\mathbf{k},t)\left[1+3\frac{\cos (kR)}{(kR{)}^2}-3\frac{\sin (kR)}{(kR{)}^3}\right],$$

$${\tilde{U}}_{dd}(\mathbf{k},t)=g_{dd}\left(3[\mathbf{e}(t)\cdot \hat{k}{]}^2-1\right).$$

Full text

I Introduction

The ability to manipulate the short-ranged interactions in ultra-cold quantum gases using Feshbach resonances Chin et al. (2010) has enabled a wide range of quantum many-body physics to be investigated Bloch et al. (2008). With the production of quantum degenerate gases of highly magnetic atoms (e.g. Cr Griesmaier et al. (2005), Dy Lu et al. (2011) and Er Aikawa et al. (2012)) the field of quantum gases now has access to systems with long ranged dipole-dipole interactions (DDIs). It is desirable to tune the magnitude and sign of the DDI independently of the short-ranged interaction. Being able to do so would allow experiments to explore new phenomena, for example, stabilizing two-dimensional bright solitons Pedri and Santos (2005), realizing rotonic Kelvin-wave excitations on vortex lines Klawunn et al. (2008), or allowing the interactions between vortices to be manipulated Mulkerin et al. (2013) (also see Lahaye et al. (2009)). In 2002, Giovanazzi et al. Giovanazzi et al. (2002) proposed a scheme for tuning the DDIs by rotating the dipoles fast enough that the interaction can be time-averaged. This scheme was recently implemented in an experiment with a Bose-Einstein condensate (BEC) of 162Dy atoms by Tang et al. Tang et al. (2018). In these experiments, the anisotropic expansion dynamics of the condensate (which depends on the magnitude and sign of the DDI) was used to reveal the effects of the tuned interactions.

In practice, DDI tuning is implemented by rotating the magnetic field (used to polarize the atoms) about an axis at frequency of $\Omega$. If $\Omega$ is smaller than the Larmor frequency $\omega_{\mathrm{L}}$ (typically $\omega_{\mathrm{L}}\gtrsim 10^{6}\,$Hz), then the dipoles will follow the magnetic field. It has been assumed that if this rotation is fast compared to the trap frequencies $\omega_{\mathrm{tr}}$ (typically $\omega_{\mathrm{tr}}\sim 10^{2}\,$Hz), and hence the typical timescale of atomic motion, then the rotating dipoles can be time-averaged to describe their effect on the condensate. The angle $\varphi$ of the magnetic field with respect to its axis of rotation, determines the strength of the time-averaged DDI: the time-averaged DDI takes the form of a static DDI for dipoles polarized along the axis of rotation, but with the bare DDI coupling constant $g_{dd}$ scaled by a factor of $\frac{1}{2}(3\cos^{2}\varphi-1)$. Thus, by varying $\varphi$, the time-averaged DDI can be reduced in strength, made zero, or even negative Giovanazzi et al. (2002); Lahaye et al. (2009).

The experiments of Tang et al. Tang et al. (2018) observed that when tuning the DDI, the lifetime of their condensate was reduced by more than an order of magnitude. They speculated that the configuration of coils used in the experiment, which generated residual magnetic field gradients, could have contributed to this lifetime reduction. However, with a condensate lifetime of $\sim 160\,$ms, they concluded it would be sufficiently long for many types of experiments with tuned DDIs.

A recent theoretical treatment of the meanfield dynamics of a condensate with a rotating DDI was presented by Prasad et al. Prasad et al. (2019). Their analysis, performed in the Thomas Fermi limit (where the zero-point kinetic energy terms are neglected), predicted that with rotating dipoles the condensate was dynamically unstable in a broad parameter regime. Significantly, they concluded that this instability prevents the formation of a stable long-lived rotationally-tuned BEC. Their work was performed in the co-rotating frame of the magnetic dipoles and did not directly model the kind of dynamical scenario used in experiments to tune the DDIs. They present a dynamical simulation demonstrating the instability, but for a rather slow rotation frequency case ($\Omega=3\omega_{\mathrm{tr}}$), notably with $\Omega<\mu/\hbar$, where $\mu$ is the condensate chemical potential that usefully quantifies the interaction energy scale. So, important questions remain about the instability time-scale for experimentally relevant tuning processes, and how the instability time depends on the rotation frequency relative to characteristic frequencies of the condensate.

Here we address these questions by simulating the instability dynamics of a dipolar BEC for a realistic tuning scenario. We consider the case of magnetic dipoles taken into a tuned DDI configuration from an initially stationary (untuned) state using an appropriate ramping procedure. Our simulations are based on the truncated Wigner formalism using a nonlocal Gross-Pitaevskii equation (GPE) with noise added to simulate quantum fluctuation effects. We use measures of width and the DDI energy to quantify the tuning of the dipole interactions in our simulations. In general we find that this ramp excites collective dynamics of the condensate, but this can be minimized for ramp times exceeding the characteristic trap period. We identify the kinetic energy as a good observable to reveal the onset of dynamic instabilities arising from the DDI tuning. We find that increasing rotation frequency of the dipoles delays the onset of instability. Our results indicate that DDI tuning with minimal heating over long time scales ($\gtrsim 100\,$ms) should be feasible in experiments with sufficiently fast rotation frequencies.

II Theory

II.1 Formalism

We consider a uniform external magnetic field aligned at polar angle $\varphi(t)$ to the $z$ axis and rotating at rate $\Omega$, i.e. $\mathbf{B}(t)=B\mathbf{e}(t)$ where

[TABLE]

We assume that $\Omega\ll\omega_{\mathrm{L}}=\mu_{m}B/\hbar$, where $\mu_{m}$ is the magnetic moment of the atoms, so that magnetic dipoles follow the external field, i.e. the dipoles are aligned along $\mathbf{e}(t)$. The time-dependent GPE for the wavefunction $\psi$ is $i\hbar\dot{\psi}=\mathcal{L}_{\mathrm{GP}}\psi$, with

[TABLE]

where

[TABLE]

is the single particle Hamiltonian. The atoms interact by a contact interaction with coupling constant $g_{s}=4\pi\hbar^{2}a_{s}/m$, where $a_{s}$ is the $s$-wave scattering length (which we take to be positive). They also interact by a DDI described by the potential Lahaye et al. (2009)

[TABLE]

where $g_{dd}=\mu_{0}\mu_{m}^{2}/3\geq 0$ is the DDI coupling constant. In Eq. (2) the DDI appears via the effective dipolar interaction potential $\Phi_{D}(\mathbf{x},t)=\int d\mathbf{x}^{\prime}\,U_{dd}(\mathbf{x}-\mathbf{x}^{\prime},t)|\psi(\mathbf{x}^{\prime})|^{2}$ (we omit the $t$ dependence of $\psi$). The energy of the system is

[TABLE]

The strength of the DDIs is also conveniently characterized in terms of the dipole length $a_{dd}\equiv m\mu_{0}\mu_{m}^{2}/12\pi\hbar^{2}$, such that $g_{dd}=4\pi a_{dd}\hbar^{2}/m$. The dimensionless number $\epsilon_{dd}\equiv{g_{dd}}/{g_{s}}={a_{dd}}/{a_{s}}$ describes the ratio of the (untuned) DDIs to the contact interactions.

II.2 Time-averaged DDI

For $\varphi(t)=\varphi$ fixed and $\Omega$ sufficiently high we can time-average the dipole motion to obtain the averaged DDI

[TABLE]

where $\theta$ is the angle between $\mathbf{r}$ and the $z$ axis and the effective time-averaged coupling constant is Giovanazzi et al. (2002)

[TABLE]

By choice of the angle $\varphi$, we can vary the time-averaged coupling constant from $\bar{g}_{dd}(0)=g_{dd}$ to $\bar{g}_{dd}(\frac{\pi}{2})=-\frac{1}{2}g_{dd}$.

Stationary states $\bar{\psi}$ of the time-averaged Hamiltonian satisfy

[TABLE]

with $\bar{\mathcal{L}}_{\mathrm{GP}}=H_{sp}+g_{s}|\bar{\psi}|^{2}+\bar{\Phi}_{D}(\mathbf{x},t)$, and $\bar{\Phi}_{D}=\int d\mathbf{x}^{\prime}\,\bar{U}_{dd}(\mathbf{x}-\mathbf{x}^{\prime})|\bar{\psi}(\mathbf{x}^{\prime})|^{2}$.

III Results

We consider the dynamics of a dipolar BEC taken from an initial condition of a static dipole oriented along $z$, into a configuration with the dipoles rotating in the $xy$-plane $\varphi=\tfrac{\pi}{2}$, for which $\bar{g}_{dd}=-\tfrac{1}{2}g_{dd}$. We do this by maintaining a constant angular frequency of rotation $\Omega$ about the $z$ axis, and tilt the dipole by linearly increasing $\varphi$ from [math] to $\pi/2$ over the time period $t_{\mathrm{on}}$ (see Fig. 1). We work in the regime that $\Omega\gg 1/t_{\mathrm{on}}$, so that it is reasonable for analysis to time average over the rotation to obtain an effective interaction $\bar{g}_{dd}$ that varies with $\varphi$ on the slower timescale of $t_{\mathrm{on}}$.

For our simulations the initial BEC is taken as the ground state $\psi_{0}$ of the static configuration, i.e. $\mu\psi_{0}=\mathcal{L}_{\mathrm{GP}}\psi_{0}$ with $\varphi=0$, and $N_{c}=\int d\mathbf{x}\,|\psi_{0}|^{2}$ condensate atoms. We add noise to the ground state to account for the effects of quantum fluctuations, which can seed instabilities in the system dynamics. This noise is added as half an atom per mode (on average) in the single particle harmonic oscillator basis up to a single particle energy cutoff of $\epsilon_{\mathrm{cut}}$. Formally this method of adding noise corresponds to the truncated Wigner prescription (see Blakie et al. (2008)). This state is then evolved according to the GPE $i\hbar\dot{\psi}=\mathcal{L}_{\mathrm{GP}}\psi$ using the time-dependent dipole polarization outlined above. More details about the initial state preparation and evolution are given in Appendix A.

We can monitor the effect of the dipole dynamics on the density distribution through the evolution of its widths, which we characterize by the second moment of the position coordinates, i.e. $\sigma_{\nu}=(\int d\mathbf{x}\,\nu^{2}|\psi|^{2}/N)^{1/2}$, for $\nu=\{x,y,z\}$, where $N$ is the norm of $\psi$. The results in Fig. 2(b) show the width dynamics for a ramp with $t_{\mathrm{on}}=20\,$ms and $\Omega/2\pi=2\,$kHz. The widths change significantly during the dipole ramp in a manner that is consistent with the system experiencing an effective dipole interaction of $\bar{g}_{dd}(\varphi)$. To see this, we compare the widths from the dynamic simulation to the those of the ground state $\bar{\psi}$ of the time-averaged interaction (dashed lines), and find good quantitative agreement. Because the ramp time is comparable to the trap period we find that collective modes are excited during the ramp, as revealed by the oscillating widths that persist after the ramp has concluded.

The width behavior can be understood by considering the DDI energy $E_{D}=\frac{1}{2}\int d\mathbf{x}\,\Phi_{D}|\psi|^{2}$. Because the DDI is anisotropic, the density distribution, and importantly its aspect ratio, can change to reduce the DDI energy – a phenomenon generally referred to as magnetostriction Tang et al. (2018). The dynamics of $E_{D}$ are shown in Fig. 2(c). We note that since $E_{D}$ is evaluated using the instantaneous DDI potential it has small and rapid oscillation at a frequency of $2\Omega$ (see inset) about a mean value that varies on the time-scale of the ramp. For comparison we also show the dipole energy for the ground state of the time-averaged interaction [i.e. $\bar{E}_{D}=\frac{1}{2}\int d\mathbf{x}\,\bar{\Phi}_{D}|\bar{\psi}|^{2}$, dashed line in Fig. 2(c)] and find that it is similar to the dynamical result. This confirms that rotating the dipoles can effectively tune the DDIs in the system.

At the start of the ramp $\varphi=0$, so the bare and time-averaged interactions are identical, $\bar{g}_{dd}=g_{dd}$, and the density distribution is prolate ($\sigma_{z}>\sigma_{x,y}$), with a negative value of $E_{D}$. At $t\!\approx\!12.2\,$ms the polar angle is at the magic value $\varphi_{M}=\cos^{-1}\!\frac{1}{\sqrt{3}}\approx 54.7^{\circ}$, where the time-averaged DDI is zero and the density distribution is close to isotropic, with $E_{D}\approx 0$.111In our results the time-averaged ground state becomes isotropic, and $E_{D}$ is maximized, prior to the magic angle because the harmonic trap is oblate. At the conclusion of the ramp $\bar{g}_{dd}=-\frac{1}{2}g_{dd}$ and the energy is minimized by the density distribution being oblate ($\sigma_{z}<\sigma_{x,y}$).

In Fig. 2(d)-(f) we consider a faster ramp of $t_{\mathrm{on}}=10\,$ms that is shorter than the trap period. In this case, the collective modes are more strongly excited. The out of phase oscillation of $\sigma_{z}$ and $\sigma_{x,y}$ [see Fig. 2(e)] reveals that the dominant mode excited has a quadrupolar character. This also manifests as an oscillation in the dipole energy. For this ramp time we consider the condensate to be strongly excited as the amplitude of the energy oscillation in $E_{D}$ is comparable to the expected mean value $|\bar{E}_{D}|$ for the time-averaged ground state. From hereon we focus on the slower $t_{\mathrm{on}}=20\,$ms ramp to ensure that the state prepared at the conclusion of the ramp is close to the time averaged ground state.

We now consider the emergence of dynamical instabilities after the dipole ramp. We find that these instabilities are revealed by monitoring the kinetic energy

[TABLE]

In Fig. 3(a) we show the evolution of $E_{K}$ for various rotation frequencies $\Omega$. Initially $E_{K}$ displays oscillatory dynamics at a frequency comparable to the trap frequency, arising from the collective modes excited. However, at later times $E_{K}$ suddenly starts rapid growth, corresponding to the instability identified in Prasad et al. (2019) and marking where the condensate begins to heat. For the $\Omega/2\pi=2\,$kHz case, corresponding to the same simulation shown in Fig. 2(a)-(c), $E_{K}$ is seen to start growing at $t\approx 50\,$ms and appears to strongly diverge by $60\,$ms, even though little evidence of the instability is apparent at these times in either the widths or dipole energy [see Fig. 2(b)-(c)]. For simulations with faster rotation frequencies the initial phase of oscillatory $E_{K}$ dynamics is almost identical (independent of $\Omega$), however the transition to unstable behavior is seen to be delayed to later times. For $\Omega/2\pi=3\,$kHz, the instability does not emerge until $t\approx 110\,$ms.

We repeat the above analysis for a condensate with twice as many atoms ($N_{c}=5\times 10^{3}$) in Fig. 3(b). Here we see that for the same value of $\Omega$ the dynamic instabilities tend to manifest at earlier times than for the lower atom number case of Fig. 3(a). Nevertheless, for sufficiently high rotation frequencies (here $\Omega/2\pi>8\,$kHz), the instability can be pushed out to $t\gtrsim 100\,$ms.

To quantify the instability time, we define $t_{\mathrm{inst}}$ as the time when $E_{K}$ increases to twice its initial value. To tune the DDIs successfully requires that $t_{\mathrm{inst}}$ is appreciably greater than $t_{\rm{on}}$. The dynamic instabilities are dependent on the sampling of noise in the initial state, and different trajectories will have different instability times. We show some examples of variance of $t_{\mathrm{inst}}$ for different trajectories in Fig. 4(a).

We also examine the dependence of our results on the single particle energy cutoff $\epsilon_{\mathrm{cut}}$ used to define the modes in which the noise is added. Notably as $\epsilon_{\mathrm{cut}}$ increases the number of modes with noise added grows as $N_{\mathrm{cut}}\propto\epsilon_{\mathrm{cut}}^{2}$, and hence the total noise added to the simulations increases. Our results in Fig. 4(b) show that the instability time decreases as we add noise to higher energy modes (i.e. as $N_{\mathrm{cut}}$ increases). However, once modes with a single particle energy of up to $\epsilon_{\mathrm{cut}}\approx 2\mu$ are populated with noise, the instability time levels off and has a weak dependence on the noise. A case with twice as many atoms ($N_{c}=5\times 10^{3}$) and a higher chemical potential (see Appendix B) is shown in the inset to Fig. 4(b), also confirming that the results converge for $\epsilon_{\mathrm{cut}}\approx 2\mu$. This dependence on noise suggests that the most important modes for the instability occur in the energy range $\epsilon\in[0,2\mu]$, since populating higher energy modes with noise has negligible effect on the instabilities.

Finally, we explore the instability time over a wide parameter regime, varying the condensate atom number $N_{c}$, the ratio of the DDI to the contact interaction $\epsilon_{dd}$, the trap frequency $f$ (as defined in the caption of Fig. 5), and the rotation frequency $\Omega$. For each system case $(N_{c},\epsilon_{dd},f)$ we have simulated the tuning dynamics using five trajectories for each value of $\Omega$.

We have checked that the results we show are independent of the choice of numerical simulation grid, as detailed in Appendix A. For sufficiently low rotation frequency (all of the results shown), the results are independent of the grid, i.e. for different grids the changes in $t_{\mathrm{inst}}$ are smaller than 1 ms. For higher rotation frequencies our simulations results become grid dependent and are not shown. In these high frequency simulations we necessarily must use small numerical time-steps (see Appendix A) and individual simulations in these cases take up to four days to run using Titan-V graphics processor unit (GPU) hardware. For faster rotation rates grid-converged results may be possible with larger grids, but the limited GPU memory means we are unable to thoroughly explore such cases at this time.

The results for the instability times of our simulations, shown in Fig. 5(a), are seen to vary strongly with parameters. The shortest instability times are $t_{\mathrm{inst}}\approx 10\,$ms (i.e. $<t_{\mathrm{on}}$) indicating the system becomes unstable before the conclusion of the tuning ramp. This case tends to occur for large atom numbers at low rotation frequencies. The instability time of the $N_{c}=5\times 10^{3}$ and $f=70\,$Hz cases are seen to increase with $\Omega$, and we extrapolate that they would have a $100\,$ms instability time for $\Omega/2\pi\sim 15\times 10^{3}\,$Hz. The $N_{c}=10^{4}$ cases follow a similar trend, but increase more slowly with increasing $\Omega$.

We find that our data can be scaled to better reveal the general trends. We utilize our earlier observation (related to the sensitivity of the instability time to the noise) that the relevant instabilities appear to scale with an energy $\epsilon\sim\mu$, to identify $\mu$ as a key energy scaling parameter. Scaling the rotation frequency by $\mu N_{c}$, and the instability time by the trap frequency, the rescaled data is shown Fig. 5(b). While the data does not completely collapse to a universal line, the general trend of instability time increasing with $\Omega$ is clearly revealed. Interestingly the $\mu N_{c}$ scaling of $\Omega$ suggests that the required rotation frequencies scale extensively with $N_{c}$, indicating that rotational tuning for large condensates will require large rotation frequencies, that may be challenging to obtain in experiments.

IV Discussion and Conclusions

In this paper we have simulated the process of rotational DDI tuning implemented using a realistic tuning ramp. We have examined the properties of the resulting state in dynamical simulations to verify the effects of the tuned interactions, and the ramp time scale needed to avoid excessively exciting collective modes. We observe the dynamic instabilities predicted by Prasad et al. Prasad et al. (2019), and find that the kinetic energy provides a useful observable of the instability and allows us to identify an instability time. Our results indicate that the time scale of the instabilities is sensitive to the rotation frequency of the magnetic field, and we find that in many cases the instability can be delayed for sufficiently high frequency. We emphasize that our results here are for the maximally tuned case, i.e. with $\varphi=\pi/2$ where $\bar{g}_{dd}=-\frac{1}{2}g_{dd}$. In general we expect that the instabilities should be weaker for smaller values of $\varphi$.

The first experiments examining rotational tuning have observed a lifetime reduction, evidenced through the loss of atoms in the system. Our work suggests that a systematic study in terms of the kinetic energy evolution (or equivalently momentum width) will allow the instabilities to be systematically quantified. This may also help to understand to what extent the lifetime in experiment was due to magnetic field gradients versus the rotationally induced instability.

A deeper understanding of the microscopic origin of the instabilities is clearly needed. Prasad et al. Prasad et al. (2019) presented the results of a polynomial basis approach that they used to quantify the eigenvalues of the linear excitations of the tuned condensate in the Thomas-Fermi approximation. A non-zero imaginary part of an eigenvalue indicates that the associated mode is dynamically unstable, with the magnitude of the imaginary part relating to the growth rate. Our results motivate the need for a better understanding of the instabilities beyond the Thomas-Fermi approximation, and of the spatial character of the unstable modes. This may help identify factors that affect the instability timescale, and suggest methods for extending the lifetime of condensates with rotationally tuned DDIs.

Another important direction to consider is the effect of tuned DDIs in the regime of quantum droplets. Such droplets can occur when $\epsilon_{dd}>1$ (e.g., see Kadau et al. (2016); Chomaz et al. (2016); Schmitt et al. (2016); Wächtler and Santos (2016a, b); Baillie et al. (2016); Bisset et al. (2016)) and the condensate becomes unstable to mechanical collapse, but is stabilized by beyond meanfield terms, i.e. leading order repulsive effect of quantum fluctuations. We have not included these quantum fluctuation effects in our results presented here, but they should have a minor effect in the $\epsilon_{dd}<1$ regime we have considered. In the case of a self-bound droplet, the chemical potential $\mu$ is typically negative (e.g. see Baillie et al. (2017)), and the scaling we predict in Fig. 5(b) clearly cannot hold.

The results shown in this paper for dysprosium apply to other species with the same $\epsilon_{dd}$ and atom number, but with $t_{\mathrm{on}}$, $t_{\mathrm{inst}}$, $\Omega$, and $f$ scaled by $ma_{dd}^{2}$. For 164Dy to 168Er, times must be divided by (and frequencies multiplied by) approximately $3.8$. The results of Fig. 5(b) are unchanged, but $f=50$ and $f=70\>\mathrm{Hz}$ for 164Dy correspond to $f\approx 190$ and $f\approx 267\>\mathrm{Hz}$ respectively for 168Er.

Acknowledgements.

We thank S. Prasad, B. Mulkerin, A. M. Martin, Y. Tang, L. Chomaz, and B. Lev for helpful discussions. We acknowledge the contribution of NZ eScience Infrastructure (NeSI) high-performance computing facilities, and support from the Marsden Fund of the Royal Society of New Zealand.

Appendix A Numerical Treatment

The initial state used in simulations is

[TABLE]

where $H_{sp}\phi_{j}=\epsilon_{j}\phi_{j}$, and $\alpha_{j}$ is a complex Gaussian noise with $\langle\alpha_{j}\rangle=0$ and $\langle|\alpha_{j}|^{2}\rangle=\tfrac{1}{2}$. We restrict the sum to modes with energies $\epsilon_{j}\leq\epsilon_{\mathrm{cut}}$, exclude the ground state, and denote the total number of modes within the sum as $N_{\mathrm{cut}}$.

We start with the dipoles aligned along $\hat{z}$ and tilt into the $xy$-plane linearly while the dipoles rotate at frequency $\Omega$ using the ramp outlined in Section III [see Fig. 2(a)]. For calculating the dipolar interaction, we use a cutoff which is given by (4) for $r<R$ and zero otherwise. The resulting Fourier transform is Ronen et al. (2006)

[TABLE]

where

[TABLE]

For calculations shown in Section III, we used numerical grids with $192$ gridpoints in each direction and a spatial grid resolution of $0.15\sqrt{\hbar/m\omega_{x}}$. For our grid checks, for each set of parameters and each rotation frequency, we also ran one of the trajectories on four additional grids, with (i) 50% more gridpoints in the radial directions, (ii) 50% more gridpoints in the axial direction, (iii) 50% greater grid range in the radial directions, (iv) 50% greater grid range in the axial direction. For one case, we also checked that if the results are independent of the simulation grid for one trajectory, then they are independent for all five trajectories. Our time step was $\Omega\Delta t=0.01$.

Appendix B Simulation properties

Here we provide additional details about the states and simulations used in Fig. 5. Table 1 gives the trap parameter $f$ [recall $(\omega_{x},\omega_{y},\omega_{z})=2\pi f\times(1,1,1.1)$], the relative dipole strength $\epsilon_{dd}$ (note we fix $a_{dd}=130.8\,a_{0}$), and condensate number $N_{c}$. Also, the chemical potential $\mu$ for the initial condensate (without noise), the total number of atoms $N$ when half-an-atom of noise is added to single particle modes with energies $\epsilon_{j}\leq 2\mu$ and $E/N$, where $E$ is the average total energy of the initial state including noise.