Vortex Lattice Formation in Dipolar Bose-Einstein Condensates via Rotation of the Polarization

TL;DR

This paper investigates how rotating the polarization of a dipolar Bose-Einstein condensate induces vortex lattice formation, combining analytical solutions, stability analysis, and numerical simulations to reveal a novel vortex creation mechanism.

Contribution

It introduces a new method of vortex formation in dipolar BECs through polarization rotation, supported by semi-analytical and numerical analysis.

Findings

Rotation of polarization leads to vortex lattice formation.

Dynamical instabilities seed vortex creation.

Vortices relax into a triangular lattice configuration.

Abstract

The behaviour of a harmonically trapped dipolar Bose-Einstein condensate with its dipole moments rotating at angular frequencies lower than the transverse harmonic trapping frequency is explored in the co-rotating frame. We obtain semi-analytical solutions for the stationary states in the Thomas-Fermi limit of the corresponding dipolar Gross-Pitaevskii equation and utilise linear stability analysis to elucidate a phase diagram for the dynamical stability of these stationary solutions with respect to collective modes. These results are verified via direct numerical simulations of the dipolar Gross-Pitaevskii equation, which demonstrate that dynamical instabilities of the co-rotating stationary solutions lead to the seeding of vortices that eventually relax into a triangular lattice configuration. Our results illustrate that rotation of the dipole polarization represents a new route to vortex formation in dipolar Bose-Einstein condensates.

Figures

Equations

$$i\hslash \frac{\partial \psi }{\partial t}=\left[-\frac{{\hslash }^2{\nabla }^2}{2m}+V_{\text{T}}+V_{\text{int}}+i\hslash \Omega \left(x\frac{\partial }{\partial y}-y\frac{\partial }{\partial x}\right)\right]\psi .$$

$$V_{\text{T}}\left(\mathbf{r}\right)=\frac{1}{2}m\left[{\omega }_{\perp }^2{\rho }^2+{\omega }_z^2{z}^2\right],$$

$$V_{\text{int}}(\mathbf{r},t)=g|\psi {|}^2+\frac{C_{\text{dd}}}{4\pi }\int \mathrm{d}{\mathbf{r}}^{\prime }\frac{1-3{\cos }^2{\theta }_{\mathbf{r},{\mathbf{r}}^{\prime }}(t)}{|\mathbf{r}-{\mathbf{r}}^{\prime }{|}^3}{\left|\psi ({\mathbf{r}}^{\prime },t)\right|}^2.$$

$$V_{\text{int}}\left(\mathbf{r},t\right)=g(1-{ϵ}_{\text{dd}})n\left(\mathbf{r},t\right)-3g{ϵ}_{\text{dd}}{\left(\hat{\mathbf{e}}\cdot \nabla \right)}^2{\varphi }_{\text{dd}}\left(\mathbf{r},t\right),$$

$${\varphi }_{\text{dd}}\left(\mathbf{r},t\right)=\frac{1}{4\pi }\int {\mathrm{d}}^3{r}^{\prime }\frac{{\left|\psi \left({\mathbf{r}}^{\prime },t\right)\right|}^2}{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|},$$

$$\frac{\partial n}{\partial t}$$

$$m\frac{\partial \mathbf{v}}{\partial t}$$

$$+\nabla \left[\frac{{\hslash }^2}{2m\sqrt{n}}{\nabla }^2(\sqrt{n})\right].$$

$$\mathbf{v}=\frac{\hslash \nabla S}{m}.$$

$$m\frac{\partial \mathbf{v}}{\partial t}=-\nabla \left[\frac{1}{2}m{\mathbf{v}}^2-m\mathbf{v}\cdot \left(\mathbf{\Omega }\times \mathbf{r}\right)+V_{\text{T}}+V_{\text{int}}\right].$$

$$\psi (\mathbf{r},t)=\psi (\mathbf{r},t=0)\exp \left(\frac{-i\mu t}{\hslash }\right).$$

$$0$$

$$\mu$$

$$+g(1-{ϵ}_{\text{dd}})n-3g{ϵ}_{\text{dd}}\frac{{\partial }^2{\varphi }_{\text{dd}}}{\partial x^2}.$$

$$n_{\text{TF}}(\mathbf{r})$$

$$S_{\text{TF}}(\mathbf{r},t)$$

$$\alpha =\left(\frac{{\kappa }_x^2-{\kappa }_y^2}{{\kappa }_x^2+{\kappa }_y^2}\right)\Omega .$$

$${\varphi }_{\text{dd}}(x,y,z)$$

$$+\frac{n_0{\kappa }_x{\kappa }_y}{8R_z^2}(x^4{\beta }_{200}+y^4{\beta }_{020}+z^4{\beta }_{002})$$

$$+\frac{n_0{\kappa }_x{\kappa }_y}{4R_z^2}(x^2{y}^2{\beta }_{110}+y^2{z}^2{\beta }_{011}+x^2{z}^2{\beta }_{101}),$$

$${\beta }_{ijk}\left({\kappa }_x,{\kappa }_y\right)={\int }_0^{\infty }\frac{\mathrm{d}s}{({\kappa }_x^2+s{)}^{i+\frac{1}{2}}({\kappa }_y^2+s{)}^{j+\frac{1}{2}}(1+s{)}^{k+\frac{1}{2}}}.$$

$${\tilde{\omega }}_x^2$$

$${\tilde{\omega }}_y^2$$

$${\kappa }_x^2=\frac{1}{\zeta }{\left(\frac{{\omega }_{\perp }\gamma }{{\tilde{\omega }}_x}\right)}^2\left[1+{ϵ}_{\text{dd}}\left(\frac{9}{2}{\kappa }_x^3{\kappa }_y{\beta }_{200}-1\right)\right],$$

$${\kappa }_y^2=\frac{1}{\zeta }{\left(\frac{{\omega }_{\perp }\gamma }{{\tilde{\omega }}_y}\right)}^2\left[1+{ϵ}_{\text{dd}}\left(\frac{3}{2}{\kappa }_y^3{\kappa }_x{\beta }_{110}-1\right)\right],$$

$$R_z^2=\frac{2g{n}_0}{m{\gamma }^2{\omega }_{\perp }^2}\zeta ,$$

$$\zeta =1+{ϵ}_{\text{dd}}\left(\frac{3}{2}{\kappa }_x{\kappa }_y{\beta }_{101}-1\right),$$

$$0$$

$$+(\alpha -\Omega )\left[{\tilde{\omega }}_y^2-\frac{3}{2}{ϵ}_{\text{dd}}\frac{{\omega }_{\perp }^2{\kappa }_x{\kappa }_y{\gamma }^2}{\zeta }{\beta }_{110}\right].$$

$$V_{\text{int}}(\text{r},t)=g|\psi {|}^2+\frac{C_{\text{dd}}}{3}{\mathcal{F}}^{-1}\left[{\tilde{U}}_{\text{dd}}^{R_c}(\text{k})\mathcal{F}\left[|\psi {|}^2\right]\right],$$

$${\tilde{U}}_{\text{dd}}^{R_c}(\text{k})=\left[1+3\frac{\cos \left(R_{c}k\right)}{R_c^2{k}^2}-3\frac{\sin \left(R_{c}k\right)}{R_c^3{k}^3}\right]\left(3{\cos }^2{\theta }_k-1\right).$$

$$\frac{\mathrm{d}\Omega }{\mathrm{d}t}=5\times 10^{-4}\Theta ({\Omega }_{\text{f}}-\Omega ){\omega }_{\perp }^2,$$

$$n(\mathbf{r},t)$$

$$S(\mathbf{r},t)$$

$$\frac{\partial }{\partial t}\left(\begin{array}{c}\delta S\\ \delta n\end{array}\right)=\mathbb{M}\left(\begin{array}{c}\delta S\\ \delta n\end{array}\right),$$

$$\mathbb{M}=-\left(\begin{array}{cc}{\mathbf{v}}_c\cdot \nabla & \frac{g}{\hslash }\left(1-{ϵ}_{\text{dd}}\hat{K}\right)\\ \frac{\hslash }{m}\nabla \cdot \left(n_{\text{TF}}\nabla \right)& {\mathbf{v}}_c\cdot \nabla \end{array}\right),$$

$${\mathbf{v}}_c=\frac{\hslash }{m}\nabla S_{\text{TF}}-\mathbf{\Omega }\times \mathbf{r},$$

$$\hat{K}\left[\delta n\right]=\delta n+3\frac{{\partial }^2}{\partial x^2}{\int }_{{\Gamma }_{\text{TF}}}{\mathrm{d}}^3{r}^{\prime }\frac{\delta n\left({\mathbf{r}}^{\prime },t\right)}{4\pi \left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}.$$

$$\left(\begin{array}{c}\delta S(\mathbf{r},t)\\ \delta n(\mathbf{r},t)\end{array}\right)=\sum _{\nu }\left(\begin{array}{c}\delta S_{\nu }(\mathbf{r})\\ \delta n_{\nu }(\mathbf{r})\end{array}\right)\exp ({\lambda }_{\nu }t),$$

Full text

I Introduction

Interacting quantum gases have proved to be a fertile testing ground for theories of many-body physics in recent years. In particular, dipolar quantum gases such as dipolar Bose-Einstein condensates (BECs) and degenerate dipolar fermionic gases offer a unique route to novel many-body phenomena where the subtle interplay between long- and short-range interactions plays a significant role Baranov (2008); Lahaye et al. (2009); Baranov et al. (2012), and considerable progress has been made in the study of dipolar BECs in the last two decades. After the first dipolar BEC was produced using chromium in 2005 Griesmaier et al. (2005), the effects of the dipolar interaction on the stability and excitations of the condensate were of much interest, with particular focus on the anisotropy of the condensate’s superfluidity Lahaye et al. (2008); Koch et al. (2008); Wenzel et al. (2018) and the theoretical possibility of roton modes and supersolidity Santos et al. (2003); Ronen et al. (2007). More recently the production of dysprosium and erbium BECs, which feature considerably stronger dipolar interactions than Cr BECs, has led to the discovery of an increasingly diverse range of exotic phenomena Lu et al. (2011); Aikawa et al. (2012). Examples of these include self-bound dipolar droplets Kadau et al. (2016); Schmitt et al. (2016); Ferrier-Barbut et al. (2016); Chomaz et al. (2016) and the experimental confirmation of the presence of roton modes Chomaz et al. (2018); Petter et al. (2019) and the supersolid phase Tanzi et al. (2019); Böttcher et al. (2019); Chomaz et al. (2019). Many of these have been found to be explicitly dependent on the role of beyond-mean-field quantum fluctuations in stabilizing a strongly dipolar BEC Lima and Pelster (2012); Bisset et al. (2016).

One area of considerable research in BECs has centred on the emergence and properties of vortices, whose circulations are necessarily quantised due to the condensate wavefunction being locally single-valued Abo-Shaeer et al. (2001); Coddington et al. (2003); Henn et al. (2009); Kwon et al. (2016); Gauthier et al. (2019); Johnstone et al. (2019). Examples of the methods that experimentalists have used to create vortices include: direct imprinting of phase defects into the condensate Leanhardt et al. (2002), rotation of either a laser beam stirrer or the external trapping potential of the condensate itself Madison et al. (2000); Hodby et al. (2001), dragging a barrier through the condensate Neely et al. (2010); Gauthier et al. (2019); Johnstone et al. (2019), applying a rapidly oscillating perturbation to the trapping potential Henn et al. (2009), Bose-condensing a rotating normal Bose gas Haljan et al. (2001), and utilising the Kibble-Zurek mechanism to trigger the formation of topological defects Weller et al. (2008). While vortices have not yet been experimentally observed in dipolar BECs, there exists an extensive body of theoretical research regarding vortex structure, vortex lattice structure, and vortex-vortex interactions in these systems Cooper et al. (2005); Zhang and Zhai (2005); Yi and Pu (2006); Komineas and Cooper (2007); Klawunn et al. (2008); Abad et al. (2009); Mulkerin et al. (2013); Kumar et al. (2016); Cai et al. (2018); Cidrim et al. (2018); Bland et al. (2018); Martin et al. (2017).

In this work we propose a new mechanism for producing vortices in dipolar BECs that involves the direct rotation of the polarizing field of the dipoles. For a dipolar BEC in a harmonic trap with an initially static dipole polarization we consider the effects of slowly increasing the rotation frequency of the dipole polarization about an orthogonal axis. Using a semi-analytical model based on dipolar Gross-Pitaevskii theory, it is shown that the co-rotating stationary state of the condensate becomes dynamically unstable against symmetry-breaking collective modes as the rotation frequency approaches the transverse trapping frequency. Numerically, it is demonstrated that this dynamical instability results in the formation of vortices in the previously vorticity-free condensate, and that these vortices ultimately self-order into a triangular Abrikosov vortex lattice. Experimentally, the rotation of the polarization of a dipolar BEC has already been achieved in the context of tuning the effective strength of the time-averaged dipole-dipole interaction Giovanazzi et al. (2002); Tang et al. (2018). Crucially, the rotation frequency in such studies is at least two orders of magnitude larger than that investigated in our work, suggesting that the mechanism we propose for the production of vortices in dipolar BECs is readily achievable via the use of existing technologies.

This paper is structured as follows. In Section II we discuss a semi-analytical formalism, based on the dipolar Gross-Pitaevskii equation and its reformulation as the equations of dipolar superfluid hydrodynamics, which we use to solve self-consistently for the Thomas-Fermi stationary solutions of a dipolar BEC with its dipoles rotating perpendicularly to their alignment. These solutions are presented in detail in Section III. In Section IV, we use the dipolar Gross-Pitaevskii equation to numerically simulate a dipolar BEC with the dipole rotation frequency being slowly increased from zero. and demonstrate that while the condensate initially obeys the Thomas-Fermi stationary solution, it develops an instability that results in vortex formation as the rotation frequency approaches the transverse trapping frequency. Finally in Section V, it is shown via linearization of the time-dependent superfluid hydrodynamic equations that the Thomas-Fermi stationary solutions are dynamically unstable against collective modes in the regime in which vortex formation is observed in Section IV.

II Thomas-Fermi Theory for the Rotating Dipolar BEC

Our semi-analytical description of the dipolar BEC with rotating dipole moments begins with a summary of the mean-field theory that we shall employ throughout this paper. We consider a dilute condensate whose bosons have intrinsic magnetic or electric dipoles that are aligned along an externally applied polarizing field. We utilise the notion of the mean-field order parameter, $\psi\equiv\langle\hat{\psi}_{0}\rangle$, which is defined as the vacuum expectation value of the bosonic annihilation operator and is normalized such that $\int\mathrm{d}\mathbf{r}|\psi(\mathbf{r},t)|^{2}$ equals $N$, the number of condensed bosons. At the lowest order of approximation, the dipolar BEC is described by the dipolar Gross-Pitaevskii theory Baranov (2008); Lahaye et al. (2009); Baranov et al. (2012) – a classical field theory for $\psi$. In a reference frame that rotates at an angular frequency $\mathbf{\Omega}=\Omega\hat{z}$ with respect to the laboratory frame, the dipolar Gross-Pitaevskii equation (dGPE) is given by

[TABLE]

where $m$ is the mass of each dipole. Here we consider the external trapping potential $V_{\text{T}}$ to be a harmonic trap of the form

[TABLE]

The axial symmetry of the trapping about $\hat{z}$ is imposed so that the trapping is identical in both the rotating and laboratory frames, and we also define the axial trapping ratio $\gamma=\omega_{z}/\omega_{\perp}$. The relevant interactions between the dipoles are the short-range, isotropic, van der Waals interaction and the anisotropic, long-range, dipole-dipole interaction. These are accounted for in the dGPE by a two-body interaction term, $V_{\text{int}}$ Góral et al. (2000); Santos et al. (2000); Yi and You (2001):

[TABLE]

Here the interaction strength of the local contribution is given by $g=4\pi\hbar^{2}a_{\text{s}}/m$, where $a_{\text{s}}$ is the $s$-wave scattering length of the relevant two-body scattering potential, while the angle $\theta_{\mathbf{r},\mathbf{r}^{\prime}}$ is defined through the relation $\cos^{2}\theta_{\mathbf{r},\mathbf{r}^{\prime}}(t)=\left[\hat{e}(t)\cdot(\hat{\mathbf{r}}-\hat{\mathbf{r}}^{\prime})\right]^{2}$, where $\hat{\mathbf{e}}(t)$ is the polarization direction of the dipoles. We also define the interaction ratio $\epsilon_{\text{dd}}=C_{\text{dd}}/(3g)$ and for the later ease of use, we also state an equivalent definition of $V_{\text{int}}$ in terms of a dipolar pseudopotential $\phi_{\text{dd}}$,

[TABLE]

such that $\nabla^{2}\phi_{\text{dd}}=-4\pi|\psi|^{2}$ O’Dell et al. (2004); Eberlein et al. (2005). In this work, we assume that $\hat{\mathbf{e}}$ rotates in the $x$-$y$ plane about $\hat{z}$ at the angular frequency $\Omega$ and, as such, we fix $\hat{\mathbf{e}}=\hat{x}\,\forall\,t$ in the co-rotating frame without loss of generality.

The form of the interaction strength $C_{\text{dd}}$ is dependent on whether the polarizing field is electric or magnetic. For a species with electric dipole moment $d$, $C_{\text{dd}}=d^{2}/\epsilon_{0}$, where $\epsilon_{0}$ is the permittivity of free space, while the expression for a species with magnetic dipole moment $\mu$ is $C_{\text{dd}}=\mu_{0}\mu^{2}$, where $\mu_{0}$ is the permeability of free space. In terms of commonly used dipolar BEC species, the bare value of $\epsilon_{\text{dd}}=C_{\text{dd}}/(3g)$ ranges from $0.16$ for 52Cr Griesmaier et al. (2006) to as high as $1.42$ for 164Dy Tang et al. (2015), both of which are orders of magnitude higher than those of ‘nondipolar’ species such as 87Rb Lahaye et al. (2009). However, $\epsilon_{\text{dd}}$ represents an additional tuneable parameter alongside $\Omega$ and $\gamma$ in this system due to the ability to tune $g$ experimentally by exploiting Feshbach resonances Lahaye et al. (2008); Koch et al. (2008). Recent theoretical studies of dipolar BECs, motivated by experimental findings highlighting the role of beyond-mean-field effects in the behaviour of strongly dipolar BECs ($\epsilon_{\text{dd}}\gtrsim 1$) have included an additional term proportional to $|\psi|^{5}$ on the right hand side of Eq. (1) Lima and Pelster (2012); Bisset et al. (2016); Baillie and Blakie (2018). However, our analysis is restricted to values of $\epsilon_{\text{dd}}$ considerably less than $1$, and as such, we expect that Eq. (1) accurately describes the behaviour of such a system in the dilute limit.

To solve Eq. (1) we shall re-express the order parameter $\psi$ as $\psi=\sqrt{n}\exp{(iS)}$, where $n=|\psi|^{2}$ and $S$ are interpreted as the condensate’s number density and phase respectively. By applying this transformation to Eq. (1) and separating out the resulting real and imaginary terms, we obtain the dipolar superfluid hydrodynamic equations Recati et al. (2001); Sinha and Castin (2001); O’Dell et al. (2004):

[TABLE]

Here the condensate’s velocity field, $\mathbf{v}$,is defined as

[TABLE]

Equation (6) is analogous to the continuity equation for a classical fluid, while Eq. (7) is analogous to the Euler equation for inviscid fluid flow. In the Thomas-Fermi (TF) regime, which generally exists when $Na_{s}$ is sufficiently high, the effects of zero-point kinetic energy fluctuations of the condensate are negligible Edwards and Burnett (1995); Parker and O’Dell (2008). This is equivalent to ignoring the ‘quantum pressure’ term in Eq. (7) that is proportional to $\nabla[\nabla^{2}(\sqrt{n})/\sqrt{n}]$, with the resulting simplified Euler equation given by

[TABLE]

When $\Omega=0$, the rotating frame is equivalent to the laboratory frame, with Eqs. (6) and (9) reducing to the description of a dipolar BEC in the TF regime with a time-independent polarization axis, a limit that has been documented extensively in previous studies O’Dell et al. (2004); Eberlein et al. (2005).

Initially, we focus on the stationary solutions of Eqs. (1) and (9), whose time dependence is specified by the condensate’s chemical potential, $\mu$, as

[TABLE]

Via the definitions of $n$ and $v$, and Eqs. (4), (5), (6), (9) and (10), the stationary solutions in the TF limit obey

[TABLE]

Let us impose Ansätze for $n$ and $S$ of the form

[TABLE]

and define the ratios $\kappa_{x}=R_{x}/R_{z}$ and $\kappa_{y}=R_{y}/R_{z}$. Here the TF density is defined as nonzero only when the quantity expressed in Eq. (13) is positive, while $n_{0}=15N/(8\pi R_{x}R_{y}R_{z})$ is a factor that normalises $n_{\text{TF}}$ to $N$. Exact solutions within the TF approximation may be determined via a set of self-consistency relations that are satisfied by the parameters $\kappa_{x}$, $\kappa_{y}$, $R_{z}$ and $\alpha$, which we shall proceed to outline briefly. Firstly, via Eq. (8), the choice of $S_{\text{TF}}$ implies that the resulting velocity field, $\mathbf{v}=\alpha\nabla(xy)$, has zero vorticity. Substituting Eqs. (13) and (14) into Eq. (11) results in the condition Recati et al. (2001)

[TABLE]

Equation (15) provides an intuitive meaning for the amplitude of the velocity field, $\alpha$: for positive nonzero $\Omega$ the cross-section of the condensate’s density profile, in the $x$-$y$ plane, is elongated along the co-rotating $x$-axis ($y$-axis) when $\alpha$ is positive (negative). Thus the condensate is axially symmetric about $\hat{z}$ when $\alpha=0$.

The dipolar pseudopotential, $\phi_{\text{dd}}$, that corresponds to the TF density specified in Eq. (13) is exactly given by Eberlein et al. (2005); van Bijnen et al. (2010); Prasad et al. (2019)

[TABLE]

where $\beta_{ijk}\left(\kappa_{x},\kappa_{y}\right)$ denotes the following integral:

[TABLE]

If we define the effective $x$ and $y$ trapping frequencies as

[TABLE]

evaluating Eq. (12) gives us the self-consistency relations we seek Prasad et al. (2019):

[TABLE]

[TABLE]

To obtain the velocity field as well as the shape of the density it is sufficient to solve Eqs. (20), (21) and (24) self-consistently.

III Stationary Solutions of the Thomas-Fermi Problem

In Fig. 1, we plot $\alpha/\omega_{\perp}$ as a function of $\Omega/\omega_{\perp}$ at fixed $\gamma=1$ for $\epsilon_{\text{dd}}\in\{0,0.1,0.5,0.8\}$ (a) and at fixed $\epsilon_{\text{dd}}=0.5$ for $\gamma\in\{0.1,1,10\}$ (b). In Ref. Prasad et al. (2019), it was found that a bifurcation exists in this system at a given rotation frequency $\Omega_{\text{b}}$, a quantity that is dependent on $\epsilon_{\text{dd}}$ and $\gamma$, at which the number of stationary solutions increases from 1 to 3 Prasad et al. (2019). When $\epsilon_{\text{dd}}=0$, the bifurcation diagram is symmetric about the $\Omega$ axis and two additional symmetric branches emerge when $\Omega=\omega_{\perp}/\sqrt{2}$. Let us denote the $\alpha>0$ and $\alpha<0$ branches as Branches I and III respectively; both of these subcritical branches terminate at $\Omega=\omega_{\perp}$, a limit which is characterised by $\alpha\rightarrow\pm\omega_{\perp}$. However the $\alpha=0$ branch, which we shall denote as Branch II, persists for all $\Omega$. Conversely, when $\epsilon_{\text{dd}}>0$, Branch I starts from $\alpha=0$ at $\Omega=0$ and exhibits a monotonic increase of $\alpha$ until the limit $\alpha\rightarrow\omega_{\perp}$ as $\Omega\rightarrow\omega_{\perp}$. Instead of the symmetric bifurcation that is present when $\epsilon_{\text{dd}}=0$, a nonzero $\epsilon_{\text{dd}}$ results in a symmetry-breaking bifurcation where Branches II and III obey $\alpha<0$ and are simply connected at $\Omega=\Omega_{\text{b}}(\epsilon_{\text{dd}},\gamma)$ but are disconnected from Branch I. While Branch III terminates at $\Omega=\omega_{\perp}$ with $\alpha\rightarrow-\omega_{\perp}$ as $\Omega\rightarrow\omega_{\perp}$, the overcritical Branch II persists for $\Omega\rightarrow\infty$. Though we do not consider this limit, it has been found that Branch II monotonically approaches the limit $\alpha\rightarrow 0^{-}$ as $\Omega\rightarrow\infty$ Prasad et al. (2019). In Fig. 1 (b), it is seen that the effect of the trapping aspect ratio $\gamma$ is qualitative and leads to shifts in the bifurcation frequency relative to the case of the spherically symmetric trap. We also note that for Branch I, $\frac{\mathrm{d}\alpha}{\mathrm{d}\Omega}$ is higher for smaller $\gamma$ and higher $\epsilon_{\text{dd}}$.

The qualitative features of the bifurcation diagram shown in Fig. 1 are similar to the diagram corresponding to the TF limit of a BEC, with or without $z$-polarised dipoles, subject to a harmonic trap that is anistropic in the $x$-$y$ plane and rotating about $\hat{z}$ Recati et al. (2001); van Bijnen et al. (2009). However in such systems the symmetry-breaking parameter is not $\epsilon_{\text{dd}}$ but is instead the $x$-$y$ trapping ellipticity. It is also possible to obtain stationary solutions and a corresponding phase diagram in the presence of the rotation of both an anisotropic trap and the dipole polarization. In this regime, the interplay between a positive $\epsilon_{\text{dd}}$ with the dipoles polarized along $\hat{x}$ in the co-rotating frame and a negative trapping ellipticity – that is, an elongation of the trap along the co-rotating $y$-axis – may lead to interesting effects. However, as a minimal model of vortex nucleation in dipolar BECs we assume that the trapping is cylindrically symmetric about $\hat{z}$.

IV Dipolar Gross-Pitaevskii Equation Simulations

In this system the properties of Branch II ($\alpha<0\,\forall\,\Omega\in[\Omega_{\text{b}},\infty)$) have been studied in previous investigations into the effective rotational tuning of $\epsilon_{\text{dd}}$ Prasad et al. (2019); Baillie and Blakie (2019). Instead we are interested in Branch I of the stationary solutions ($\alpha\geq 0$) in the presence of a spherically-symmetric trap, that is, $\gamma=1$. Via Eq. (15) we see that a nonzero $\alpha$ implies that the condensate density profile ‘rotates’ in the laboratory frame despite the velocity field being irrotational and free of vorticity. When the harmonic trapping of a (non)dipolar condensate in the TF limit is subject to a quasi-adiabatic angular acceleration from $\Omega=0$, the condensate deviates from the expected TF solution at a critical rotation frequency $\Omega=\Omega_{\text{v}}$ Sinha and Castin (2001); Madison et al. (2001); Parker et al. (2006); van Bijnen et al. (2007, 2009) where vortices enter the system and the condensate acquires a nonzero angular momentum. It has also been established in such systems that this vortex state, while energetically favourable to the TF profile at nonzero rotation frequencies, will emerge from the TF solution only in the presence of a dynamical instability that forces the system away from the stationary solution. Since the bifurcation diagram shown in Fig. 1(a) is qualitatively similar to that of a (non)dipolar condensate subject to a rotating harmonic trap with a nonzero ellipticity in the $x$-$y$ plane Recati et al. (2001); van Bijnen et al. (2009), we expect that slowly increasing the rotation frequency of the dipole polarization from zero will result in the triggering of a dynamical instability and a subsequent transition to a state with vortices. We explore this idea via two complementary methods, the first being simulations of the dGPE in the TF regime of a dipolar BEC being subject to an acceleration of the rotation of the dipole polarization about $\hat{z}$, and the second being a linearization of Eqs. (6) and (9) about the corresponding TF stationary states along Branch I over a range of values of $\epsilon_{\text{dd}}$.

The numerical integration of Eq. (1) is performed with the ADI-TSSP method Bao and Wang (2006), which is an extension of the split-step Fourier method incorporating rotation. All simulations are undertaken with a 1923 grid with the spatial step $d=0.15\sqrt{\hbar/(m\omega_{\perp})}\equiv 0.15l_{\perp}$ and temporal step $\Delta t=0.004\omega_{\perp}^{-1}$. In order to make direct comparisons with the TF analysis, we require a large number of bosons Parker and O’Dell (2008) and thus we fix $N=10^{5}$. Since fast Fourier transform algorithms are naturally periodic, we employ a spherical cut-off to the dipolar potential, the second term of Eq. (3), and restrict the range of the DDI to a sphere of radius $R_{c}$ in real-space in order to reduce the effects of alias copies. Therefore, using the convolution theorem, the interaction potential can be re-expressed as Ronen et al. (2006)

[TABLE]

where $\mathcal{F}$ ($\mathcal{F}^{-1}$) denotes the (inverse) Fourier transform and the k-space dipolar pseudo-potential is

[TABLE]

Here, $\theta_{k}$ is the angle between k and the direction of the dipoles, and the cut-off radius $R_{c}$ is chosen to be larger than the system size.

In order to verify the results presented in Fig. 1(a) we take the imaginary time stationary solution for $\Omega=0$ and a fixed $\epsilon_{\text{dd}}$ and very slowly ramp up $\Omega$ in real time, traversing through the stationary solutions. We employ the following procedure for the real-time ramp rate of $\Omega$,

[TABLE]

where $\Theta(\cdot)$ is the Heaviside function and $\Omega_{\text{f}}>0$ is a final constant choice of rotation frequency, at which the acceleration is halted. In this procedure we define $t=0$ to be the time when $\Omega=\Omega_{\text{f}}$. To assess the stability of these stationary states we employ the above ramp procedure, terminating the angular acceleration at the specified $\Omega_{\text{f}}$, then allowing the system to evolve for $t=500\omega_{\perp}^{-1}$. Any slow growing instabilities will have had time to seed vortices into the system, allowing us to precisely define a critical $\Omega=\Omega_{\text{v}}$ at the boundary between vorticity and vorticity-free solutions. We also model the the random external symmetry-breaking perturbations that may shift the condensate state away from the stationary state in an experimental scenario by modifying the condensate density at the initial timestep; at each spatial grid point, the condensate density is subjected to a random, local perturbation of up to $5\%$ of the local density. To ensure that the ramp rate is as close to being adiabatic as possible, we have performed several test simulations for half and double the ramp rate, and found a negligible difference for the onset of instability.

Figure 2 depicts the density with fixed $\epsilon_{\text{dd}}=0.1$ (panels (a)-(f)) and $\epsilon_{\text{dd}}=0.5$ (panels (g)-(i)), as cross-sections at $z=0$ ramped up to rotation frequencies $\Omega_{\text{f}}=0.55\omega_{\perp}$, $0.85\omega_{\perp}$, and $0.95\omega_{\perp}$. Odd numbered rows show snapshots taken to $t=0$, i.e. the first moment when $\Omega=\Omega_{\text{f}}$, and even numbered rows show snapshots after $t=500\omega_{\perp}^{-1}$ evolution. Overlaid on these density profiles is the predicted TF solution, found through the following procedure: first, the value of $R_{z}$ is found through a numerical fit of the simulation to a 3D TF profile, then the $(x,y)$ TF radii are generated from Eq. (20) as $R_{x}=\kappa_{x}R_{z}$ and Eq. (21) as $R_{y}=\kappa_{y}R_{z}$ from the solutions for $\alpha$ for Branch I. For low $\Omega$, the condensate density is consistent with the TF stationary solution: the density is smooth and approximates the paraboloid profile of Eq. (13). However, at some critical $\Omega=\Omega_{\text{v}}$ the solution goes through a dynamical instability. This is first visible through a rippling of the density and then the formation of a spiral-like structure that later evolves towards a turbulent vortex state.

With the dGPE simulations in hand, it is also possible to compare the results from these simulations to the predicted density profiles from the TF self-consistency relations in Eqs. (20), (21) and (24). Figure 3 compares the numerically obtained solutions for $\alpha$ (markers) with those from the TF analysis (solid lines), for the same parameters as those shown in Fig. 2. All values of $\alpha$ are measured at the end of the ramp procedure, at $t=0$, however the red crosses show simulations where solutions are unstable at long times (c.f. panels (h) and (k)). Quite remarkably, the larger $\epsilon_{\text{dd}}$ solution is dynamically stable for faster $\Omega$. For $\epsilon_{\text{dd}}=0.1$ and $\Omega>0.8\omega_{\perp}$ vortices enter the condensate and there is no sensible measure of $\alpha$. The boundary between the red crosses and the blue triangles or green circles is indicative of the value of $\Omega_{\text{v}}$ for each $\epsilon_{\text{dd}}$.

We also find that, given sufficient time to evolve, the vortices we observe in the dGPE simulation form a vortex lattice. This is shown in Fig. 4 for $\epsilon_{\text{dd}}=0.5$ and $\Omega_{\text{f}}=0.850\omega_{\perp}$, taking the initial condition from Fig. 2(h), after a total temporal evolution of $t=5000\omega_{\perp}^{-1}$. After $t=100\omega_{\perp}^{-1}$ the condensate density shows small fluctuations on its surface and at $t=350\omega_{\perp}^{-1}$ large quantities of vortices have entered the condensate, but are still confined to the peripheries of the density. Finally, by $t=5000\omega_{\perp}^{-1}$ the vortices have entered the condensate and have relaxed into a triangular Abrikosov vortex lattice. Several theoretical studies have indicated that this may be the ground state for a system of vortices in (non-)dipolar BECs at finite rotation frequencies Cooper et al. (2005); Zhang and Zhai (2005); Komineas and Cooper (2007); Fetter (2009); Cai et al. (2018); Kumar et al. (2016); Martin et al. (2017). In nondipolar BECs subject to a ramping up of the trapping rotation frequency, a vortex lattice is seen to form following a period of evolution of the nondipolar GPE at constant rotation frequency after the vortex instability has been triggered Lobo et al. (2004); Parker et al. (2006). In the context of dipolar BECs, the results of $3$D study with no contact interactions, i.e. $a_{s}=0,\epsilon_{\text{dd}}\rightarrow\infty$, also predict a square condensate density like that of Fig. 2(l) Kumar et al. (2016).

However, a comprehensive numerical dGPE study of the ground state configurations of vortex lattices in quasi-$2$D trapping geometries, i.e. $\gamma\gg 1$, finds that when the dipole orientation is in the $x$-$y$ plane, a triangular lattice configuration is favoured for $\epsilon_{\text{dd}}\lesssim 0.8$ and a striped condensate phase with a square lattice configuration is favourable for $\epsilon_{\text{dd}}\gtrsim 0.8$ in the $\Omega\rightarrow\omega_{\perp}$ limit Cai et al. (2018). Given that a triangular vortex lattice is predicted for the regimes that we consider in our work, albeit in a quasi-$2$D geometry, it is possible that a striped condensate phase with a square vortex lattice geometry may be seen if we allow $\Omega_{\text{f}}$ to approach $\omega_{\perp}$ and if $\epsilon_{\text{dd}}\rightarrow 1$. However, the dimensionality of the system would be expected to play a significant role in the nature of the vortex lattice configuration, and beyond-mean-field effects play a significant role in the behaviour of dipolar BECs with large $\epsilon_{\text{dd}}$ Lima and Pelster (2012); Bisset et al. (2016). As such, a direct comparison of our methodology with these results cannot be made. A thorough consideration of beyond-mean-field effects as well as the eventual vortex lattice configurations across different regimes of $\epsilon_{\text{dd}}$ and $\Omega_{\text{f}}$ is beyond the scope of this work, but still warrants further study.

V Linearization of the Thomas-Fermi Stationary Solutions

The solutions of the dGPE equation clearly depict the emergence of vortices in the vorticity-free TF stationary states after a certain critical rotation frequency, $\Omega=\Omega_{\text{v}}(\epsilon_{\text{dd}})$. As each real-time evolution of the dGPE begins with a random symmetry-breaking perturbation of the density, a natural conclusion would be that this perturbation induces a dynamical instability of the condensate that allows for a dynamical route from the TF stationary states to lower energy states such as vortex states at high rotation frequencies Sinha and Castin (2001); van Bijnen et al. (2009). In order to understand when the stationary states become dynamically unstable, we linearize the fully time-dependent hydrodynamic equations, Eqs. (6) and (9), about the TF stationary states as predicted by Eq. (11) and (12). This yields a formulation of the fully time-dependent state in terms of collective modes about the stationary states, one or more of which may be dynamically unstable against external perturbations. Initially, we express the time-dependent solutions to Eq. (6) and (9) as

[TABLE]

Here, $n_{\text{TF}}(\mathbf{r})$ and $S_{\text{TF}}(\mathbf{r},t)$ are given by Eqs. (13) and (14), thus satisfying Eqs. (11) and (12), while $\delta n$ and $\delta S$ are considered to be ‘small’ perturbations about the stationary states. The linearization of the time-dependent problem is subsequently carried out by substituting Eqs. (28) and (29) into Eqs. (6) and (9), utilising Eqs. (11) and (12) to simplify the result, and neglecting contributions from terms that are of higher order than linear in the fluctuations $\delta n$ and $\delta S$. The resulting system of coupled first-order ordinary differential equations is given by Sinha and Castin (2001); Parker et al. (2006); van Bijnen et al. (2007, 2009):

[TABLE]

Here $\Gamma_{\text{TF}}$ is defined as $\Gamma_{\text{TF}}=\{\mathbf{r}\in\mathbb{R}^{3}:n_{\text{TF}}(\mathbf{r})>0\}$. The linearized fluctuations are given by a sum of collective modes, indexed by $\nu$, such that

[TABLE]

To diagonalize Eq. (35), we use a monomial basis for $\delta n_{\nu}$ and $\delta S_{\nu}$ of the form $\{x^{i}y^{j}z^{k}\}$ Sinha and Castin (2001). By inspection of $\mathbb{M}$, a given mode will feature the same value of $l=\max{[i+j+k]}$ for both $\delta n$ and $\delta S$, a number which we refer to as the order of the mode. While it is relatively simple to calculate how a given monomial is transformed by the nondipolar components of $\mathbb{M}$, the evaluation of the dipolar contribution is quite involved, but by using methods originally developed for the study of Newtonian potentials inside classical self-gravitating ellipsoidal fluids, it is possible to calculate $\widehat{K}\left[\delta n\right]$ Levin and Muratov (1971); van Bijnen et al. (2010); Sapina et al. (2010). For a polynomial $x^{i}y^{j}z^{k}$, we rewrite the exponents as

[TABLE]

where $\{\delta_{\lambda},\delta_{\mu},\delta_{\nu}\}\in\{0,1\}$ and $\{\lambda,\mu,\nu\}$ are integers. We also introduce the definition $\sigma=\lambda+\mu+\nu+1$. The integral in $\mathbb{M}$ is then given by

[TABLE]

where the expression is simplified by the definitions

[TABLE]

and $\beta_{ijk}$ represents the integrals defined in Eq. (17). The dipolar contribution to $\mathbb{M}$ is obtained by taking the second partial derivative of the expression on the RHS of Eq. (37) with respect to $x$. We also note that, formally, Eq. (16) may be derived via use of Eqs. (37) – (39).

After diagonalizing Eq. (35), the TF stationary state is said to be dynamically unstable if any mode $\nu$ is characterised by $\text{Re}[{\lambda_{\nu}}]>0$, since the magnitude of the mode grows exponentially with time. Within this framework, we may be certain that a stationary state is stable only if there are no eigenvalues $\lambda_{\nu}$ possessing a positive real component. As there are an infinite-number of possible collective modes, it is necessary to fix a truncation parameter, $N_{\text{max}}$, such that $l=\max{[i+j+k]}\leq N_{\text{max}}$. For a choice of $N_{\text{max}}$, it is possible that there are modes of a higher order at a given point in parameter space that are unstable and, as such, the absence of unstable modes for a given $N_{\text{max}}$ is not a guarantee of dynamical stability. However, if a sensible choice of $N_{\text{max}}$ yields at least one unstable mode, it is sufficient to claim that the corresponding stationary state is dynamically unstable. Due to computational constraints we fix $N_{\text{max}}=10$, a choice which we consider to be sufficient for a qualitative analysis. We also note that the $l=0$ and $l=1$ modes are always dynamically stable. Firstly, by inspection of $\mathbb{M}$, it is clear that truncating the basis to $N_{\text{max}}=0$ yields two null eigenvalues. As for the $l=1$ dipole modes, which represent centre-of-mass rigid-body oscillations of the condensate, these are exactly decoupled from the two-body interactions of the condensate as per Kohn’s theorem: the eigenvalues of the six dipole modes are purely imaginary and are given by $\text{Im}(\lambda_{l=1})\in\{\omega_{\perp}\pm\Omega,\gamma\omega_{\perp}\}$, along with their complex conjugates Recati et al. (2001).

Let us consider the domain $\Omega/\omega_{\perp}\in[0.5,1]$ and $\epsilon_{\text{dd}}\in[0,0.5]$, which roughly corresponds to the domain in which the vortex instability was observed in the dGPE simulations. In Fig. 5, points in parameter space are shaded if we find at least one eigenvalue of $\mathbb{M}$ with a positive real component greater than $10^{-4}$ (a) and $10^{-2}$ (b). Figure 5 demonstrates that, at higher values of $\Omega<\omega_{\perp}$, the Branch I TF stationary solutions become dynamically unstable. The resulting region of dynamical instability is found to be continuous as $\Omega\rightarrow\omega_{\perp}$ but separating into distinct arcs for lower $\Omega$. Each of these arcs corresponds to distinct polynomial orders $l\geq 3$ and is connected to the bulk of the dynamically unstable region at successively higher values of both $\Omega$ and $\epsilon_{\text{dd}}$. There is also an instability of an $l=2$ mode, corresponding to a narrow arc that is disconnected from the connected region of instability and emerges at $\epsilon_{\text{dd}}=0$ at $\Omega/\omega_{\perp}\approx 0.74$. In Fig. 5 (a), we also observe additional narrow arcs of instability corresponding to slowly growing high-order unstable modes. These are not present in Fig. 5(b) due to the higher eigenvalue cutoff. We note that for larger $\epsilon_{\text{dd}}$, Figs. 5(a) and (b) imply dynamical stability near $\Omega=\omega_{\perp}$, but this is merely due to the choice of $N_{\text{max}}$ being insufficiently large in order to accurately probe the dynamical stability of the solutions in this regime. However, Fig. 5 is qualitatively similar to the corresponding stability diagrams for both nondipolar BECs and $z$-polarized dipolar BECs in asymmetrical rotating harmonic traps, where the transverse trapping ellipticity plays an analogous role to $\epsilon_{\text{dd}}$ in our system Sinha and Castin (2001); van Bijnen et al. (2007, 2009). In these analogous systems, for a given ellipticity, the stationary solutions are always unstable above a critical rotation frequency and so we expect that the solutions are dynamically unstable when $\Omega\rightarrow\omega_{\perp}$, a property that would be clear if a larger value of $N_{\text{max}}$ had been chosen.

We now proceed to compare the predictions of the dGPE simulations and the linearization of the TF hydrodynamic equations regarding the dynamical instability of the condensate. Overlain on Fig. 5 are red dashed lines at constant $\epsilon_{\text{dd}}$, each representing the trajectories of dGPE simulations utilising ramp procedure $(2)$ as described in Section IV. Each dashed line terminates in a red star at $\Omega=\Omega_{\text{v}}(\epsilon_{\text{dd}})$, the rotation frequency which we estimate to be the lowest rotation frequency at which the dynamical instability develops from an average over $10$ simulations. For low $\epsilon_{\text{dd}}$ we find that $\Omega_{\text{v}}$ corresponds roughly to the narrow, high-order modes seen in Fig. 5(a), while the bulk of the unstable region of parameter space that survives in Fig. 5(b) seems to be responsible for the dynamical instability when $\epsilon_{\text{dd}}$ equals $0.4$ or $0.5$. This is a somewhat unusual departure from previous analyses of similar systems, such as the rotating asymmetrically trapped condensate, as the instabilities in these systems seem to be driven exclusively by the bulk of the instability region rather than the narrow instability arcs in the corresponding eigenvalue plots Sinha and Castin (2001); van Bijnen et al. (2007, 2009). This tendency for the condensate to be more unstable than expected warrants further study.

VI Conclusion

We have identified a new mechanism for the dynamical formation of vortex lattices in dipolar BECs, a phenomenon which has previously been considered in the context of time-dependent rotating traps. Using the dGPE and the dipolar superfluid hydrodynamic equations in the Thomas-Fermi limit, we have investigated the effects of a slow increase of the rotation frequency of the dipole moments about the cylindrical axis of the applied harmonic trap. Simulations of the dGPE demonstrate that the stationary solutions in the TF regime suffer a dynamical instability that causes the dipolar condensate to spontaneously reduce its energy and acquire a nonzero angular momentum via the seeding of vortices. Further time evolution of the dGPE reveals that the the vortices self-order into a triangular Abrikosov vortex lattice, a ground state that has been predicted for a dipolar BEC in numerous prior studies. The domain of parameter space that supports such a dynamical instability is further explored via linearization of the dipolar superfluid hydrodynamic equations about the TF stationary states at a given rotation frequency and dipolar coupling strength. By characterizing the linear stability of the collective modes obtained via this procedure, we find that a dynamical route to the vortex lattice ground state is available for any value of the dipolar interaction strength as the rotation frequency approaches the transverse trapping frequency. This result paves the way to studying vortex lattices in dipolar BECs in experiments without any modification to the externally applied trap, but instead via manipulation of the orientation of the dipole moments themselves. In particular, such systems may harbor novel phases such as a square vortex lattice phase Cooper et al. (2005); Zhang and Zhai (2005); Martin et al. (2017). The technology to realise this scheme experimentally already exists at the time of writing, as the polarization of a dipolar BEC has been rotated at considerably higher rotation frequencies in an experiment investigating the effective rotational tuning of the dipolar interaction in the rapid rotation limit Tang et al. (2018). While we predict the formation of a triangular vortex lattice in the parameter regimes we have explored, it may be possible to obtain square lattice configurations for $\epsilon_{\text{dd}}\sim 1$ and this, as well as the effects of modifying the trapping aspect ratio, $\gamma$, warrants further study. Finally, we note that the mechanism detailed here may be generalized by considering the rotation of dipole moments tilted away from the $x$-$y$ axis, which might alter the nature of the vortex instability and the properties of the resulting vortex lattice.

Acknowledgements.

S.B.P. is supported by an Australian Government Research Training Program Scholarship and by the University of Melbourne. A.M.M. would like to thank the Institute of Advanced Study (Durham University, U.K.) for hosting him during the initial stages of developing this collaborative research project and the Australian Research Council (Grant No. LE180100142) for support. T.B. and N.G.P. thank the Engineering and Physical Sciences Research Council of the UK (Grant No. EP/M005127/1) for support.