Vortices and vortex lattices in quantum ferrofluids

At low energies, and far from any two-body bound states, the interatomic interactions can be described by an effective pseudo-potential which is the sum of a contact term originating from the van der Waals interactions and a bare dipolar term [32, 79],

$$\begin{eqnarray}\begin{array}{*{35}{l}} U\left(\mathbf{r}\mathbf{-}{{\mathbf{r}}^{\prime}}\,\right) & = & {{U}_{\text{vdW}}}\left(\mathbf{r}-{{\mathbf{r}}^{\prime}}\,\right)+{{U}_{\text{dd}}}\left(\mathbf{r}\mathbf{-}{{\mathbf{r}}^{\prime}}\,\right) \\ {} & = & g\delta \left(\mathbf{r}\mathbf{-}{{\mathbf{r}}^{\prime}}\right)+\frac{{{C}_{\text{dd}}}}{4\pi}\frac{1-3{{\cos}^{2}}\theta}{|\mathbf{r}\mathbf{-}{{\mathbf{r}}^{\prime}}{{|}^{3}}}, \end{array}\end{eqnarray} \tag{ 1 }$$

where we have assumed the dipoles are polarised by an external field with θ being the angle between the polarisation direction and the inter-atom vector $\mathbf{r}-{{\mathbf{r}}^{\prime}}$. The short-range interactions are characterized by the coupling constant $g=4\pi {{\hbar}^{2}}{{a}_{s}}/m$, where a_s is the s-wave scattering length and m is the mass of the atoms. The strength of the dipolar interactions is set by C_dd. For magnetic dipoles ${{C}_{\text{dd}}}={{\mu}_{0}}{{d}^{2}}$ [81], where ${{\mu}_{0}}$ is the permeability of free space and d is the magnetic dipole moment of the atoms. Equation (1) also holds for electric dipoles induced by a static electric field $\mathbf{E}=\mathbf{\hat{k}}E$, for which the coupling constant is ${{C}_{\text{dd}}}={{E}^{2}}\alpha _{\text{p}}^{2}/{{\epsilon}_{0}}$ [82, 83], where ${{\alpha}_{\text{p}}}$ is the static polarizability and ${{\epsilon}_{0}}$ is the permittivity of free space. The formation and cooling to degeneracy of polar molecules, with significant electric dipole moments, is proving to be challenging. However, progress is ongoing with ⁴⁰K⁸⁷Rb [84] and ¹³³Cs⁸⁷Rb [85], which are expected to have a significant coupling constant for electric fields on the order of several hundred V cm⁻¹. The dipolar interaction U_dd, illustrated in figure 2, is negative for $\theta =0$, representing the attraction of head-to-tail dipoles, and positive for $\theta =\pi /2$, representing the repulsion of side-by-side dipoles. At the 'magic angle', ${{\theta}_{\text{m}}}=\arccos \left(1/\sqrt{3}\right)\approx {{54.7}^{\circ}}$, the dipolar interaction is zero.

Zoom In Zoom Out Reset image size

It is often convenient to work in momentum space. The Fourier transform ${{\tilde{U}}_{\text{dd}}}\left(\mathbf{k}\right)=\mathcal{F}\left[{{U}_{\text{dd}}}\right]={\int}^{}{{U}_{\text{dd}}}\left(\mathbf{r}\right){{\text{e}}^{-\text{i}\mathbf{k}\cdot \mathbf{r}}}\text{d}\mathbf{r}$ of the dipolar interaction is [81],

$$\begin{eqnarray}&&{{\tilde{U}}_{\text{dd}}}\left(\mathbf{k}\right)={{C}_{\text{dd}}}\left({{\cos}^{2}}\alpha -\frac{1}{3}\right),\end{eqnarray} \tag{ 2 }$$

where α is the angle between $\mathbf{k}$ and the polarization direction.

The strength of the dipolar interactions is conveniently parameterized by the ratio [32],

$$\begin{eqnarray}&&{{\varepsilon}_{\text{dd}}}=\frac{{{C}_{\text{dd}}}}{3g},\end{eqnarray} \tag{ 3 }$$

where g can be tuned between $-\infty$ and $+\infty$ via a Feshbach resonance [86, 87]. In effect ${{\varepsilon}_{\text{dd}}}$ gives the relative importance of the anisotropic, long-range dipole–dipole interactions to the isotropic, short-range van der Waals interactions. It is defined with a factor of 3 in the denominator so that the homogeneous dipolar condensate is unstable when ${{\varepsilon}_{\text{dd}}}>1$—see section 4.1.1. For ⁵²Cr, ¹⁶⁸Er and ¹⁶⁴Dy, the natural value of ${{\varepsilon}_{\text{dd}}}$ is 0.16 [88], 0.4 [11] and 1.45 [9, 89], respectively.

While C_dd is conventionally positive and set to the natural value of the given atom, it is predicted to be possible to reduce C_dd below its natural value, including to negative values, by tilting the polarization direction off-axis and rotating it rapidly [90]. Hence it is feasible to consider $-\infty <{{\varepsilon}_{\text{dd}}}<\infty$, with both negative and positive C_dd. Note that for ${{C}_{\text{dd}}}<0$ the dipole–dipole interaction becomes repulsive for head-to-tail dipoles and attractive for side-by-side dipoles.

In the mean-field limit, at zero temperature, a single wavefunction, $\Psi \left(\mathbf{r},t\right)$, can be used to describe the condensate. The condensate density $n\left(\mathbf{r},t\right)=\,| \Psi \left(\mathbf{r},t\right){{|}^{2}}$ is normalized such that

$$\begin{eqnarray}&&N={\int}^{}| \Psi {{|}^{2}}\,\text{d}\mathbf{r},\end{eqnarray} \tag{ 4 }$$

where N is the number of atoms in the condensate. The wavefunction obeys the dipolar GPE [78, 81, 82],

$$\begin{eqnarray}&&\text{i}\hbar \frac{\partial \Psi}{\partial t}=\left[-\frac{{{\hbar}^{2}}}{2m}{{\nabla}^{2}}+V+g| \Psi {{|}^{2}}\,+\, \Phi \right] \Psi,\end{eqnarray} \tag{ 5 }$$

where $V=V\left(\mathbf{r}\right)$ is the external potential acting on the condensate (which in principle may also be time-dependent, but here we consider it to be static). The local term, $g| \Psi {{|}^{2}}$, arises from the van der Waals interactions and the non-local term, $\Phi$, arises from the dipolar interactions [79],

$$\begin{eqnarray} &&\Phi \left(\mathbf{r},t\right)={\int}^{}{{U}_{\text{dd}}}\left(\mathbf{r}-{{\mathbf{r}}^{\prime}}\right)n\left({{\mathbf{r}}^{\prime}},t\right)\,\text{d}{{\mathbf{r}}^{\prime}}.\end{eqnarray} \tag{ 6 }$$

If we take the dipoles to be polarized along the z-direction, then using identities from potential theory the dipolar potential can be expressed as [91, 92],

$$\begin{eqnarray} &&\Phi \left(\mathbf{r},t\right)=-3g{{\varepsilon}_{\text{dd}}}\left(\frac{{{\partial}^{2}}}{\partial {{z}^{2}}}\phi \left(\mathbf{r},t\right)+\frac{1}{3}n\left(\mathbf{r},t\right)\right),\end{eqnarray} \tag{ 7 }$$

where ϕ is a fictitious 'electrostatic' potential defined as

$$\begin{eqnarray}&&\phi \left(\mathbf{r},t\right)=\frac{1}{4\pi}{\int}^{}\frac{n\left(\mathbf{r}\mathbf{,}{{\mathbf{t}}^{\prime}}\right)}{|\mathbf{r}-{{\mathbf{r}}^{\prime}}|}\,\text{d}{{\mathbf{r}}^{\prime}}.\end{eqnarray} \tag{ 8 }$$

This effectively reduces the problem of calculating the dipolar potential $\Phi$ to one of calculating an electrostatic potential of the form (8) which is easier to compute because the Green's function $1/|\mathbf{r}-{{\mathbf{r}}^{\prime}}|$ has no angular dependence. Furthermore, hundreds of years of literature exists providing analytic methods for solving electrostatic and gravitational problems with this form of interaction [93–96]. Alternatively, $\Phi$ can be evaluated in momentum space by exploiting the convolution theorem,

$$\begin{eqnarray} &&\Phi (\mathbf{r},t)={{\mathcal{F}}^{-1}}[{{\tilde{U}}_{\text{dd}}}(\mathbf{k})\tilde{n}(\mathbf{k},t)],\end{eqnarray} \tag{ 9 }$$

where $\tilde{{n}}(\mathbf{k},t)=\mathcal{F}\,[n(\mathbf{r},t)]$.

In condensate experiments the external potential V is typically harmonic with the general form,

$$\begin{eqnarray}&&V(\mathbf{r})=\frac{1}{2}m(\omega _{x}^{2}{{x}^{2}}+\omega _{y}^{2}{{y}^{2}}+\omega _{z}^{2}{{z}^{2}}),\end{eqnarray} \tag{ 10 }$$

where ${{\omega}_{j}}$ (j  =  x, y, z) are the trap's angular frequencies. In general, this gives rise to three harmonic oscillator lengths, ${{\ell}_{j}}=\sqrt{\hbar /m{{\omega}_{j}}}$, which are the characteristic length scales imposed by the trap on the wavefunction in the three different directions. However, cylindrically-symmetric traps are common, defined as,

$$\begin{eqnarray}&&V(\mathbf{r})=\frac{1}{2}m(\omega _{\bot}^{2}{{\rho}^{2}}+\omega _{z}^{2}{{z}^{2}})=\frac{1}{2}m\omega _{\bot}^{2}({{\rho}^{2}}+{{\gamma}^{2}}{{z}^{2}}),\end{eqnarray} \tag{ 11 }$$

where $\gamma ={{\omega}_{z}}/{{\omega}_{\bot}}$ is the so-called trap ratio. When $\gamma \gg 1$ the BEC shape will typically be oblate (flattened) while for $\gamma \ll 1$ it will typically be prolate (elongated).

Time-independent solutions of the GPE satisfy,

$$\begin{eqnarray} &&\Psi (\mathbf{r},t)=\psi (\mathbf{r}){{\text{e}}^{-\text{i}\mu t/\hbar}},\end{eqnarray} \tag{ 12 }$$

where μ is the chemical potential⁵. Inserting this into equation (5), the time-independent dipolar GPE for the time-independent wavefunction $\psi (\mathbf{r})$ is

$$\begin{eqnarray}&&\mu \psi =-\frac{{{\hbar}^{2}}}{2m}{{\nabla}^{2}}\psi +V\psi +g|\psi {{|}^{2}}\psi + \Phi \psi.\end{eqnarray} \tag{ 13 }$$

Solutions of the time-independent GPE are stationary solutions of the system, and the lowest energy of these is the ground state.

The energy of the condensate is given by,

$$\begin{eqnarray}\begin{array}{*{35}{l}} E & = & {\int}^{}\left[\frac{{{\hbar}^{2}}}{2m}|\nabla \Psi {{|}^{2}}\,+\,V| \Psi {{|}^{2}}\,+\,\frac{g}{2}| \Psi {{|}^{4}}\,+\,\frac{\Phi}{2}| \Psi {{|}^{2}}\right]\,\text{d}\mathbf{r} \\ {} & = & {{E}_{\text{kin}}}+{{E}_{\text{pot}}}+{{E}_{\text{vdW}}}+{{E}_{\text{dd}}}. \end{array}\end{eqnarray} \tag{ 14 }$$

The terms represent (from left to right) kinetic energy E_kin, potential energy E_pot, the van der Waals interaction energy E_vdW and the dipolar interaction energy E_dd. Provided that the potential V is independent of time, then the total energy E is conserved during the time evolution of the GPE.

Comparing the relative size of the kinetic term and the net interaction term in the dipolar GPE in an untrapped (V  =  0) system defines a length scale termed the healing length,

$$\begin{eqnarray}&&\xi =\frac{\hbar}{\sqrt{m\mu}},\end{eqnarray} \tag{ 15 }$$

which may be interpreted as the minimum length-scale over which the wavefunction changes appreciably.

Efficient numerical methods for solving the dipolar GPE are available [97–101] and progress has been made on extending this treatment to include finite temperature effects and quantum fluctuations [102–105].

There is a deep link between the GPE and fluid dynamics. Indeed, the condensate can be thought of as a fluid, characterised by its density and velocity distributions. This is revealed by writing the condensate wavefunction in the Madelung form $\Psi \left(\mathbf{r},t\right)=\sqrt{n\left(\mathbf{r},t\right)}{{\text{e}}^{\text{i}S\left(\mathbf{r},t\right)}}$, where the local phase, $S\left(\mathbf{r},t\right)$, defines the fluid velocity field $\mathbf{v}\left(\mathbf{r},t\right)$,

$$\begin{eqnarray}&&\mathbf{v}\left(\mathbf{r},t\right)=\frac{\hbar}{m}\boldsymbol{\nabla} S\left(\mathbf{r},t\right).\end{eqnarray} \tag{ 16 }$$

Inserting the Madelung form into the GPE, and separating real and imaginary terms, yields two equations which together are exactly equivalent to the GPE. The first is the continuity equation,

$$\begin{eqnarray}&&\frac{\partial n}{\partial t}+\boldsymbol{\nabla} \cdot \left(n\mathbf{v}\right)=0.\end{eqnarray} \tag{ 17 }$$

This embodies the conservation of the number of atoms. The second equation is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} m\frac{\partial \mathbf{v}}{\partial t}=-\boldsymbol{\nabla} \left(\frac{1}{2}m{{v}^{2}}+V+gn+ \Phi -\frac{{{\hbar}^{2}}}{2m}\frac{{{\nabla}^{2}}\sqrt{n}}{\sqrt{n}}\right). \end{array}\end{eqnarray} \tag{ 18 }$$

The ${{\nabla}^{2}}\sqrt{n}/\sqrt{n}$-term is the quantum pressure, arising from the zero-point kinetic energy of the atoms. It can be dropped when the interactions and external potential dominate the zero-point motion, leading to the Thomas–Fermi approximation. In this regime, and in the absence of dipolar interactions, equations (17) and (18) resemble the continuity and Euler hydrodynamical equations for inviscid fluids. As such they are often referred to as the superfluid hydrodynamic equations [106–110]. Equations (17) and (18) have been extended to include dipolar interactions and are referred to as the dipolar superfluid hydrodynamic equations.

4.1.1. Three-dimensional case.

For $V\left(\mathbf{r}\right)=0$ (uniform condensate of infinite extent), the stationary solution is,

$$\begin{eqnarray}&&\psi =\sqrt{{{n}_{0}}},\quad \quad \mu ={{n}_{0}}g\left(1-{{\varepsilon}_{\text{dd}}}\right),\end{eqnarray} \tag{ 19 }$$

i.e. a state of uniform density n₀. The two contributions to the chemical potential μ are the uniform mean-field potentials generated by the van der Waals and the dipolar interactions, respectively. In the absence of dipolar interactions (${{\varepsilon}_{\text{dd}}}=0$), the corresponding solution has chemical potential $\mu ={{n}_{0}}g$. By comparison, the homogeneous dipolar system is akin to a non-dipolar system but with an effective coupling

$$\begin{eqnarray}&&{{g}_{\text{eff}}}=g\left(1-{{\varepsilon}_{\text{dd}}}\right).\end{eqnarray} \tag{ 20 }$$

For a three-dimensional homogeneous dipolar condensate, the Bogoliubov dispersion relation between the energy E_B and momentum $\mathbf{p}$ of a perturbation is given by,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{E}_{\text{B}}}\left(\mathbf{p}\right) & = & \sqrt{{{c}^{2}}\left(\theta \right)\,{{p}^{2}}+{{\left(\frac{{{p}^{2}}}{2m}\right)}^{2}}}, \end{array}\end{eqnarray} \tag{ 21 }$$

where $c\left(\theta \right)$ is the speed of sound,

$$\begin{eqnarray}&&c\left(\theta \right)=\sqrt{\frac{g{{n}_{0}}}{m}\left[1+{{\varepsilon}_{\text{dd}}}\left(3{{\cos}^{2}}\theta -1\right)\right]}.\end{eqnarray} \tag{ 22 }$$

The angle θ is that between the excitation momentum $\mathbf{p}$ and the polarization direction. For low momenta the spectrum is linear ${{E}_{\text{B}}}\approx c\left(\theta \right)\,p$ which is characteristic of phonons with a phase velocity $c\left(\theta \right)$ that depends on direction. For higher momenta the relation becomes quadratic in p which is characteristic of free-particle excitations.

The amplitude of a mode specified by momentum $\mathbf{p}$ evolves in time as $\exp \left(-i{{E}_{\text{B}}}\left(\mathbf{p}\right)t/\hbar \right)$. If ${{E}_{\text{B}}}\left(\mathbf{p}\right)$ should become imaginary the relevant amplitude grows exponentially, signifying a dynamical instability. In the case of the three-dimensional homogeneous dipolar BEC considered in this section, with ${{E}_{\text{B}}}\left(\mathbf{p}\right)$ provided by equation (21), such an instability arises for small p, i.e. long wavelengths, and this is known as the phonon instability, familiar from non-dipolar attractive (g  <  0) condensates [108]. Examining the parameter space over which equation (21) is real-valued indicates that the three-dimensional homogeneous system is stable against the phonon instability in the range $-0.5\leqslant {{\varepsilon}_{\text{dd}}}\leqslant 1$ for g  >  0, and ${{\varepsilon}_{\text{dd}}}\leqslant -0.5,{{\varepsilon}_{\text{dd}}}>1$ for g  <  0.

A full theoretical treatment of a trapped BEC involves solving the dipolar GPE, given in equation (5) [108, 109]. The non-local nature of the mean-field potential describing dipolar interactions means that this task is more challenging than for purely s-wave BECs. Moreover, the stability of the condensate becomes non-trivial, becoming dependent on the geometry of the trap and the number of atoms (in addition to the dipole strength). Additionally, a dipolar condensate can suffer from a density-wave instability associated with a novel type of excitation called a roton in analogy with a similar type of excitation in superfluid helium [111–113]. To characterise the stability of a dipolar condensate we first derive and examine the Thomas–Fermi ground state solutions for a dipolar BEC.

4.2.1. Thomas–Fermi solutions.

The problem of finding the ground state solution (as well as low-energy dynamics) is greatly simplified by making use of the Thomas–Fermi approximation, whereby density gradients in the GPE (or, equivalently, the hydrodynamic equations) are ignored, allowing analytic solutions [106]. For a non-dipolar condensate, with repulsive van der Waals interactions, this is valid for $N{{a}_{s}}/\bar{\ell}\gg 1$, where $\bar{\ell}={{\left({{\ell}_{x}}{{\ell}_{y}}{{\ell}_{z}}\right)}^{1/3}}$ is the geometric mean of the harmonic oscillator lengths [108, 109]. This regime is relevant to many experiments.

In the dipolar case, the Thomas–Fermi approximation is valid when the net interactions are repulsive and the number of atoms is large; rigorous criteria have been established for certain geometries in [114]. Although the governing equations for a dipolar BEC contain the non-local potential $\Phi \left(\mathbf{r}\right)$, exact solutions known from the pure s-wave case hold, in modified form, in the dipolar case too [91, 92], and we make extensive use of them throughout this review.

Consider a dipolar condensate polarized in the z-direction, with repulsive van der Waals interactions (a_s  >  0), and confined by a cylindrically-symmetric trap of the form of equation (11). We limit the analysis to the regime of $-0.5\leqslant {{\varepsilon}_{\text{dd}}}\leqslant 1$, where the Thomas–Fermi approach predicts that stationary solutions are stable [91]. Outside of this regime the condensate becomes prone to collapse [115, 116]. Under the Thomas–Fermi approximation the time-independent GPE (13) reduces to,

$$\begin{eqnarray}&&\frac{1}{2}m\omega _{\bot}^{2}\left({{\rho}^{2}}+{{\gamma}^{2}}{{z}^{2}}\right)+ \Phi \left(\mathbf{r}\right)+gn\left(\mathbf{r}\right)=\mu.\end{eqnarray} \tag{ 23 }$$

Making use of the electrostatic formulation given in equations (7) and (8), exact solutions of equation (23) can be obtained for any general parabolic trap, as proven in appendix A of [92]. In particular, the solutions for the density profile take the form,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} n\left(\mathbf{r}\right)={{n}_{\text{cd}}}\left(1-\frac{{{\rho}^{2}}}{R_{\bot}^{2}}-\frac{{{z}^{2}}}{R_{z}^{2}}\right)\,\,\,\,\text{for}\,\,\,n\left(\mathbf{r}\right)\geqslant 0, \end{array}\end{eqnarray} \tag{ 24 }$$

where ${{n}_{\text{cd}}}=15N/\left(8\pi R_{\bot}^{2}{{R}_{z}}\right)$⁶ is the central density, and R_z and ${{R}_{\bot}}$ are the Thomas–Fermi radii of the condensate in the axial and transverse directions.

Remarkably, the inverted parabolic density profile, equation (24), is of the same form as that found in non-dipolar BECs [108, 109]. However, whereas non-dipolar BECs have the same aspect ratio, $\kappa ={{R}_{\bot}}/{{R}_{z}}$, as the trap, for dipolar BECs $\kappa \ne \gamma$ (in general) and must be evaluated using the following transcendental equation [91, 92],

$$\begin{eqnarray}\begin{array}{*{35}{l}} 0 & = & 3{{\kappa}^{2}}{{\varepsilon}_{\text{dd}}}\left[\left(\frac{{{\gamma}^{2}}}{2}+1\right)\frac{f\left(\kappa \right)}{1-{{\kappa}^{2}}}-1\right] \\ {} & {} & +\,\,\left({{\varepsilon}_{\text{dd}}}-1\right)\left({{\kappa}^{2}}-{{\gamma}^{2}}\right), \end{array}\end{eqnarray} \tag{ 25 }$$

where,

$$\begin{eqnarray}&&f\left(\kappa \right)=\frac{1+2{{\kappa}^{2}}}{1-{{\kappa}^{2}}}-\frac{3{{\kappa}^{2}}\text{arctanh}\sqrt{1-{{\kappa}^{2}}}}{{{\left(1-{{\kappa}^{2}}\right)}^{3/2}}},\end{eqnarray} \tag{ 26 }$$

which takes the value f  =  1 at $\kappa =0$, and monotonically decreases towards f  =  −2 as $\kappa \to \infty$, passing through zero at $\kappa =1$. This is a robust feature: the same transcendental equation is recovered using a variational approach based on a gaussian ansatz for the condensate wave function [79, 80]. For a non-dipolar (${{\varepsilon}_{\text{dd}}}=0$) condensate one finds the expected result that $\kappa =\gamma$. However, the presence of dipolar interactions leads to magnetostriction of the condensate, such that $\kappa <\gamma$ for ${{\varepsilon}_{\text{dd}}}>0$ and $\kappa >\gamma$ for ${{\varepsilon}_{\text{dd}}}<0$. This behaviour is shown in figure 3 (top) [117]. Note that, within the range $-0.5\leqslant {{\varepsilon}_{\text{dd}}}\leqslant 1$ these are global solutions; elsewhere the solutions are either metastable (light grey shading) or unstable (dark grey shading). For conventional dipoles (${{C}_{\text{dd}}}>0$, ${{\varepsilon}_{\text{dd}}}>0$), the condensate is least stable in prolate ($\gamma <1$) traps; here the dipoles lie predominantly in the attractive head-to-tail configuration and undergo collapse when ${{\varepsilon}_{\text{dd}}}$ becomes too large. By contrast, in oblate ($\gamma >1$) traps stability is enhanced since the dipoles lie predominantly in the repulsive side-by-side configuration. Meanwhile the opposite is true for anti-dipoles (${{C}_{\text{dd}}}<0$, ${{\varepsilon}_{\text{dd}}}<0$). Away from the instabilities, these solutions agree well with numerical solutions of the full dipolar GPE in the Thomas–Fermi regime [115]. Close to the instabilities, zero-point kinetic energy (neglected within the Thomas–Fermi approach) can enhance the stability of the solutions.

Zoom In Zoom Out Reset image size

Once the BEC aspect ratio κ is found from the transcendental equation, the Thomas–Fermi radii are determined by the expressions,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{R}_{\bot}} & = & {{\left[\frac{15Ng\kappa}{4\pi m\omega _{\bot}^{2}}\left\{1+{{\varepsilon}_{\text{dd}}}\left(\frac{3}{2}\frac{{{\kappa}^{2}}f\left(\kappa \right)}{1-{{\kappa}^{2}}}-1\right)\right\}\right]}^{1/5}}, \end{array}\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{R}_{z}} & = & \frac{{{R}_{\bot}}}{\kappa}, \end{array}\end{eqnarray} \tag{ 28 }$$

and the total energy is given by,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{E}_{\text{TF}}} & = & \frac{N}{14}m\omega _{x}^{2}R_{x}^{2}\left(2+\frac{{{\gamma}^{2}}}{{{\kappa}^{2}}}\right) \\ {} & {} & +\,\,\frac{15}{28\pi}\frac{{{N}^{2}}g}{R_{x}^{2}{{R}_{z}}}\left[1-{{\varepsilon}_{\text{dd}}}\,f\left(\kappa \right)\right]. \end{array}\end{eqnarray} \tag{ 29 }$$

The first term corresponds to the trapping energy and the second to the s-wave and dipolar interaction energies. Finally, the dipolar potential inside the condensate can be explicitly obtained as [92],

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\Phi}_{\text{TF}}}\left(\rho,z\right) & = & \frac{{{n}_{\text{cd}}}{{C}_{\text{dd}}}}{3}\left[\frac{{{\rho}^{2}}}{R_{x}^{2}}-\frac{2{{z}^{2}}}{R_{z}^{2}}-f\left(\kappa \right)\left(1-\frac{3}{2}\frac{{{\rho}^{2}}-2{{z}^{2}}}{R_{x}^{2}-R_{z}^{2}}\right)\right]. \end{array}\end{eqnarray} \tag{ 30 }$$

This is generally a saddle-shaped function that reflects the anisotropic nature of the dipolar interactions and drives the elongation of the BEC along the polarization direction. A more general version of ${{\Phi}_{\text{TF}}}\left(\mathbf{r}\right)$ for the case of a dipolar BEC without cylindrical symmetry is given later in equation (71).

4.2.2. Outside the Thomas–Fermi regime: rotons and density oscillations.

For a condensate strongly confined in one dimension it is possible to reduce the effective dimensionality of the system to form a quasi-two-dimensional condensate. This offers a simplified platform to study vortices and vortex lattices in dipolar condensates, while still retaining the key physics.

Consider the dipoles to be polarized at an angle α to the z-axis, lying in the x  −  z plane, and strong harmonic confinement $V(z)=\frac{1}{2}m\omega _{z}^{2}{{z}^{2}}$ in the z-direction which satisfies $\hbar {{\omega}_{z}}\gg \mu$, i.e. the trapping energy dominates over the condensate energy scale. This set-up is illustrated in figure 4 [129]. In this regime, one can approximate the wavefunction by the ansatz,

$$\begin{eqnarray} &&\Psi \left(\boldsymbol{\rho},z,t\right)={{\Psi}_{\bot}}\left(\boldsymbol{\rho},t\right){{\psi}_{z}}(z).\end{eqnarray} \tag{ 31 }$$

Zoom In Zoom Out Reset image size

Axially, the condensate is taken to be frozen into the axial ground harmonic oscillator state ${{\psi}_{z}}(z)={{\left(\pi \ell _{z}^{2}\right)}^{-1/4}}{{\text{e}}^{-{{z}^{2}}/2\ell _{z}^{2}}}$. The dynamics then become planar, parametrised by the two-dimensional time-dependent wavefunction, ${{\Psi}_{\bot}}\left(\boldsymbol{\rho},t\right)$. Note that ${{\Psi}_{\bot}}$ is normalized to the number of atoms, i.e. $N={\int}^{}|{{\Psi}_{\bot}}{{|}^{2}}\text{d}\boldsymbol{\rho}$. Inserting this ansatz into the dipolar GPE, equation (5), and integrating out the axial direction then leads to the effective two-dimensional dipolar GPE [130],

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{i}\hbar \frac{\partial {{\Psi}_{\bot}}}{\partial t}=\left[-\frac{{{\hbar}^{2}}}{2m}\nabla _{\bot}^{2}+V\left(\boldsymbol{\rho}\right)+\frac{g}{\sqrt{2\pi}{{l}_{z}}}|{{\Psi}_{\bot}}{{|}^{2}}\,+\,{{\Phi}_{\bot}}\right]{{\Psi}_{\bot}}. \end{array}\end{eqnarray} \tag{ 32 }$$

The $g/\sqrt{2\pi}{{l}_{z}}$ coefficient characterises the effective van der Waals interactions in the plane, and ${{\Phi}_{\bot}}$ is the effective planar dipolar potential,

$$\begin{eqnarray}&&{{\Phi}_{\bot}}\left(\boldsymbol{\rho},t\right)={\int}^{}U_{\text{dd}}^{\bot}\left(\boldsymbol{\rho}-\boldsymbol{\rho}^{{\prime}}\right)~{{n}_{\bot}}\left(\boldsymbol{\rho}^{{\prime}},t\right)\,\text{d}\boldsymbol{\rho}^{{\prime}},\end{eqnarray} \tag{ 33 }$$

where ${{n}_{\bot}}=|{{\Psi}_{\bot}}{{|}^{2}}$ is the two-dimensional density. The real-space form of the effective two-dimensional dipolar interaction potential $U_{\text{dd}}^{\bot}$ is given elsewhere [131], while in this review its Fourier transform is used [130, 132],

$$\begin{eqnarray}&&\tilde{U}_{\text{dd}}^{\bot}\left({{\mathbf{\tilde{q}}}}\right)=\frac{4\pi {{C}_{\text{dd}}}}{9\sqrt{2\pi}{{l}_{z}}}\left[{{F}_{\parallel}}\left({{\mathbf{\tilde{q}}}}\right){{\sin}^{2}}\alpha +{{F}_{\bot}}\left({{\mathbf{\tilde{q}}}}\right){{\cos}^{2}}\alpha \right],\end{eqnarray} \tag{ 34 }$$

where ${{F}_{\parallel}}\left({{\mathbf{\tilde{q}}}}\right)=-1+3\sqrt{\pi}\frac{\tilde{{q}}_{x}^{2}}{\tilde{{q}}}{{\text{e}}^{{{\tilde{{q}}}^{2}}}}\text{erfc}\left(\tilde{{q}}\right)$, ${{F}_{\bot}}\left({{\mathbf{\tilde{q}}}}\right)=2-3\sqrt{\pi}\tilde{{q}}{{\text{e}}^{{{\tilde{{q}}}^{2}}}}\text{erfc}\left(\tilde{{q}}\right)$ and ${{\mathbf{\tilde{q}}}}=\mathbf{q}{{l}_{z}}/\sqrt{2}$ with $\mathbf{q}$ being the projection of $\mathbf{k}$ onto the x  −  y plane, i.e. the reciprocal space analogue of $\boldsymbol{\rho}$⁸. From this momentum space representation ${{\Phi}_{\bot}}$ can then be evaluated using the convolution theorem as in equation (9). An important parameter is the ratio $\sigma ={{l}_{z}}/\xi$, where ξ is the healing length. The two-dimensional approximation requires $\sigma <1$.

Under a cylindrically-symmetric harmonic trap and for dipoles polarized along z, the above 'two-dimensional mean-field regime' is formally entered when $N{{a}_{s}}\left(1+2{{\varepsilon}_{\text{dd}}}\right)l_{z}^{3}/l_{\bot}^{4}\ll 1$. In the opposing regime, when $N{{a}_{s}}\left(1+2{{\varepsilon}_{\text{dd}}}\right)l_{z}^{3}/l_{\bot}^{4}\gg 1$, the system enters the three-dimensional Thomas–Fermi regime [114]. A more general analysis of flattened condensates in [133] has established the validity of the two-dimensional mean-field regime for arbitrary polarization direction.

In the absence of any planar trapping potential $\left[V\left(\boldsymbol{\rho}\right)=0\right]$, the stationary solution of the quasi-two-dimensional dipolar condensate is the homogeneous state [130],

$$\begin{eqnarray}&&{{\psi}_{\bot}}=\sqrt{{{n}_{0}}},\quad \mu =\frac{{{n}_{0}}}{\sqrt{2\pi}{{l}_{z}}}\left[g+\frac{{{C}_{\text{dd}}}}{3}\left(3{{\cos}^{2}}\alpha -1\right)\right],\end{eqnarray} \tag{ 35 }$$

where n₀ is the uniform two-dimensional density. This system undergoes the phonon instability when the net local interactions become attractive in the plane, i.e. when $g+{{C}_{\text{dd}}}\left[3{{\cos}^{2}}\alpha -1\right]/3<0$. The phonon unstable regions in the ${{\varepsilon}_{\text{dd}}}-\alpha$ plane are shown in figure 5. These can be understood by considering the trade-off between the van der Waals and the dipolar interactions [129, 134]. Note the divergent behaviour at the magic angle ${{\alpha}_{\text{m}}}$, across which the planar dipolar interactions switch between repulsive and attractive.

Zoom In Zoom Out Reset image size

The roton instability also arises in this quasi-two-dimensional BEC [112], as indicated in figures 5(a) and (b) (blue shaded regions). For g  >  0 the roton instability is induced by the attractive part of the dipolar interaction and is only possible for $\alpha \ne 0$; for $\alpha =0$ the condensate cannot probe the attractive part of $U_{\text{dd}}^{\bot}$ [132] (this is true only in the strict quasi-two-dimensional limit). For g  <  0 the roton exists for all α [135] (excluding the magic angle). For small α it is induced by the attractive van der Waals interactions, while for larger α it is driven by the attractive axial component of the dipolar interactions.

The extent of the dipolar BEC in the z-direction (σ) effects the stability of the roton. Specifically, as the condensate becomes narrower (σ decreases) the out-of-plane component of $U_{\text{dd}}^{\bot}$ decreases [129] and the regimes of roton instability (blue regions in figure 5) shrink.

Even single vortices are giant excitations involving a considerable fraction of the entire BEC. The energy associated with the formation of a vortex ${{E}_{\text{v}}}\equiv E-{{E}_{0}}$, where E₀ is the energy of the non-rotating (vortex-free) state, is generically much larger than the energy of elementary excitations described by the Bogoliubov spectrum given in equation (21). In a frame rotating at angular frequency $\Omega$, the total energy of the system is shifted to $E-\mathbf{\Omega}\cdot \mathbf{L}$, where $\mathbf{L}$ is the angular momentum in the laboratory frame, and hence it only becomes energetically favourable to form a vortex if the angular momentum is such that $|\mathbf{\Omega}\cdot \mathbf{L}|>{{E}_{\text{v}}}$ leading to a critical rotation frequency,

$$\begin{eqnarray}&&{{\Omega}_{\text{v}}}=\frac{{{E}_{\text{v}}}}{N\hbar {{q}_{\text{v}}}}.\end{eqnarray} \tag{ 37 }$$

E_v can be computed analytically in certain situations as we discuss below. Before we do so, it is important to point out that equation (37) considers the energetics but not the kinetics of vortex formation. Both theory and experiment reveal that the true value of ${{\Omega}_{\text{v}}}$ is often considerably higher than predicted by equation (37). This is because the vortex-free state can remain a local energy minimum separated by an energy barrier from the global minimum corresponding to the vortex state. The kinetics of vortex formation are examined in section 7.

In the simplest case of an infinite system with a vortex the condensate wavefunction can be written,

$$\begin{eqnarray}&&\psi \left(\rho,\theta,z\right)=f\left(\rho \right){{\text{e}}^{\text{i}{{q}_{\text{v}}}\theta}}.\end{eqnarray} \tag{ 38 }$$

Substituting into the dipolar GPE (5) an equation for the amplitude f about the vortex is obtained,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mu f=-\frac{{{\hbar}^{2}}}{2m}\frac{1}{\rho}\frac{\partial}{\partial \rho}\left(\rho \frac{\partial f}{\partial \rho}\right)+\frac{{{\hbar}^{2}}q_{\text{v}}^{2}}{2m{{\rho}^{2}}}\,f+g\,{{f}^{3}}+ \Phi f, \end{array}\end{eqnarray} \tag{ 39 }$$

where the Laplacian has been expressed in its cylindrically symmetric form,

$$\begin{eqnarray}&&{{\nabla}^{2}}=\frac{1}{\rho}\frac{\partial}{\partial \rho}\left(\rho \frac{\partial}{\partial \rho}\right)+\frac{1}{{{\rho}^{2}}}\frac{{{\partial}^{2}}}{\partial {{\theta}^{2}}}+\frac{{{\partial}^{2}}}{\partial {{z}^{2}}}.\end{eqnarray} \tag{ 40 }$$

The second term on the right-hand side of equation (39) is the only difference from the non-vortex case, and is associated with the kinetic energy of the circulating flow, giving rise to a centrifugal barrier. Note that $|{{q}_{\text{v}}}|>1$ vortices are energetically unstable compared to multiple singly-charged vortices, and rarely arise unless engineered; for this reason we will consider only $|{{q}_{\text{v}}}|\,=1$.

In any system uniform along the polarization direction (z) the dipolar potential reduces to a local potential $\Phi \left(\mathbf{r}\right)\to -g{{\varepsilon}_{\text{dd}}}n\left(\mathbf{r}\right)$. This can be seen from equation (7), which gives the relationship between the fictitious electrostatic potential $\phi \left(\mathbf{r}\right)$ and the dipolar potential $\Phi \left(\mathbf{r}\right)$, and noting that ${{\partial}^{2}}\phi /\partial {{z}^{2}}$ must equal zero⁹. Thus, the last two terms in equation (39) can be combined into a single contact term $g\left(1-{{\varepsilon}_{\text{dd}}}\right){{f}^{3}}$. Results that hold for the usual s-wave case in this context therefore also hold for the dipolar case provided one replaces g by $g\left(1-{{\varepsilon}_{\text{dd}}}\right)$. For example, analysis of equation (39) for the s-wave case reveals that the centrifugal barrier term dictates that the density relaxes as $1/{{\rho}^{2}}$ to the asymptotic background value $\sqrt{{{n}_{0}}}$ as $\rho \to \infty$, and for $\rho \to 0$ the density tends to zero as ${{\rho}^{2|{{q}_{\text{v}}}|}}$ [108]. Furthermore, although equation (39) can not be solved in terms of known functions, it can be solved numerically, and with appropriate scaling the result for $f\left(\rho \right)$ is universal. This solution for $f\left(\rho \right)$ can then be used in the energy functional given in equation (14) and the extra energy per unit length due to the introduction of a vortex evaluated. In the pure s-wave case one obtains [108],

$$\begin{eqnarray}&&{{\epsilon}_{\text{v}}}=\pi {{n}_{0}}\frac{{{\hbar}^{2}}}{m}\log \left(1.464\frac{b}{\xi}\right).\end{eqnarray} \tag{ 41 }$$

This result was originally obtained for superfluid ⁴He by Ginzburg and Pitaevskii in 1958 [136]. In this expression the healing length ξ, which gives the size of the vortex core, forms a lower cutoff and the length b, which could be the system radius, is the upper cutoff. Although $b\gg \xi$, their finite values avoid a logarithmic singularity that originates from the centrifugal barrier term. The only place that interactions enter this expression for ${{\epsilon}_{\text{v}}}$ is through the healing length $\xi =\hbar /\sqrt{m\mu}$ which can immediately be adapted to the dipolar case using $\mu ={{n}_{0}}g\left(1-{{\varepsilon}_{\text{dd}}}\right)$. The critical rotation frequency for this case can now be evaluated by replacing E_v by ${{\epsilon}_{v}}$ and N by N/L (number of atoms per unit length) in equation (37).

In order to obtain analytic results for the untrapped system it can sometimes be useful to use the following approximate solution for $f\left(\rho \right)$ which incorporates the correct behaviour for $\rho \to 0$ and $\rho \to \infty$,

$$\begin{eqnarray}&&{{f}^{2}}\left(\rho \right)=n\left(\rho,z\right)={{n}_{0}}\frac{{{\rho}^{\prime 2}}}{1+{{\rho}^{\prime 2}}},\end{eqnarray} \tag{ 42 }$$

where ${{\rho}^{\prime}}=\rho /\xi$ and n₀ is the density at infinity. It will be convenient later to write this in the form of a homogeneous background density n₀ and a negative vortex density n_v,

$$\begin{eqnarray}&&n\left(\rho,z\right)={{n}_{0}}+{{n}_{\text{v}}},\end{eqnarray} \tag{ 43 }$$

where ${{n}_{\text{v}}}=-{{n}_{0}}/\left(1+{{\rho}^{\prime 2}}\right)$.

A schematic of a straight singly-charged vortex line through a non-dipolar harmonically-trapped condensate is shown in figure 6. The vortex has a well defined core, with zero density and a phase singularity at its centre, relaxing to the background condensate density over a distance given by the conventional healing length $\xi =\hbar /\sqrt{m\mu}$. This is typically of the order of $0.1-1\mu$m but can be tuned by means of a Feshbach resonance.

Zoom In Zoom Out Reset image size

In three dimensions the vortices may bend, e.g. into tangles and rings, carry linear or helical Kelvin wave excitations and undergo reconnections. However, under strong axial confinement of the condensate, the dynamics become effectively two-dimensional. Being topological defects, vortices can only disappear via annihilation with an oppositely-charged vortex (the two-dimensional analog of a reconnection) or by exiting the condensate at a boundary. In trapped condensates, an off-centre vortex precesses about the trap centre; this can be interpreted in terms of a Magnus force [36]. Thermal dissipation causes a precessing vortex to spiral out of a trapped condensate [137–140]. Acceleration of a vortex (or an element of a three-dimensional vortex line) leads to emission of phonons, analogous to the Larmor radiation from an accelerating charge, although under suitable confinement these phonons can be reabsorbed to prevent net decay of the vortex [141].

Optical absorption imaging of the vortices is typically preceded by expansion of the cloud to enlarge the cores [35, 142]. This method has been extended to provide real-time imaging of vortex dynamics [36]. While this imaging approach detects density only, the vortex circulation is detectable via gyroscopic techniques [143].

Yi and Pu [144] performed the first study of vortices in a dipolar BEC, obtaining numerical solutions for a quasi-two-dimensional trapped dipolar condensate featuring a vortex. For dipoles polarized perpendicular to the plane they found the striking result that density ripples form about the vortex core for trap ratios $\gamma \sim 100$ and attractive van der Waals interactions. These ripples are not contained in the simple ansatz given above in equation (42) and seem to be a rather special feature associated with non-local interactions. Indeed, they had previously been seen in numerical simulations of vortices in superfluid ⁴He where non-local potentials are employed [145–148]. For purely dipolar (g  =  0), oblate condensates, Wilson et al [120, 122] numerically found vortex ripples for moderate trap ratios $\gamma \sim 17$, see figure 7 (top) [120], and established the link to roton mixing into the vortex solution, similar to the biconcave structure that they found was induced in vortex-free dipolar condensates (energetic favourability of dipoles aligning head-to-tail). Vortex ripples have since been studied in other works [123, 129, 134], and similar ripples arise in the presence of other localized density depletions, such as due to localized repulsive potentials [122, 130] and dark solitons [149, 150]. The presence of the vortex slightly reduces the stable parameter space for the condensate relative to the vortex-free condensate [120, 144]. For dipoles tilted perpendicular to the axis of the vortex, the vortex core becomes elliptical, see figure 7 (bottom) [144], due to the anisotropic dipolar potential in the plane.

Zoom In Zoom Out Reset image size

The properties of an off-axis straight vortex line in a trapped dipolar condensate have been considered in [151, 152] in the Thomas–Fermi regime, showing that the dipolar interaction lowers (raises) the precession speed in an oblate (prolate) trap. In the presence of thermal dissipation, making the dipolar interactions partially attractive by changing the polarization direction leads to a reduction in the condensate size and a faster decay rate of the precessing vortex [153].

For a general vortex line in three dimensions, the vortex elements interact with each other at long-range via the dipolar interactions [154], as well as the usual hydrodynamic interaction [155, 156]. This modifies the Kelvin-wave (transverse) modes of the vortex line, and can support a roton minimum in their dispersion relation. For large dipolar interactions, the Kelvin waves can undergo a roton instability, leading to novel helical or snake-like configurations [157].

With ${{C}_{\text{dd}}}<0$ and tight axial trapping, stable two-dimensional bright solitons have been predicted [158, 159], i.e. wavepackets which are self-trapped by interactions in two dimensions. This idea was extended by Tikhonenkov et al [160] to predict stable two-dimensional vortex solitons, which may be considered as a two-dimensional bright soliton carrying a central vortex.

In the hydrodynamic (Thomas–Fermi) regime appropriate for large condensates, the problem of a dipolar BEC with single vortex in a three-dimensional trap can be tackled semi-analytically [161]. In this case the energy associated with the curvature of the density due both to the trapping and the vortex core can be ignored in comparison to each of the rotational, interaction, and trapping energies. These remaining energies can be evaluated analytically by assuming a density profile much like the one given in equation (42), i.e. an unperturbed background piece plus a negative vortex 'density'. The difference is that we now take the unperturbed background density to be the inverted parabola ${{n}_{\text{TF}}}\left(\rho,z\right)$ given in equation (24) which is an exact solution of the vortex-free Thomas–Fermi problem. Thus, the density profile reads,

$$\begin{eqnarray}&&n\left(\rho,z\right)={{n}_{\text{TF}}}\left(\rho,z\right)+{{n}_{\text{v}}}\left(\rho,z\right)\end{eqnarray} \tag{ 44 }$$

where,

$$\begin{eqnarray}&&{{n}_{\text{v}}}\left(\rho,z\right)=-{{n}_{\text{TF}}}\left(\rho,z\right)\frac{{{\beta}^{2}}}{{{\beta}^{2}}+{{\rho}^{2}}}.\end{eqnarray} \tag{ 45 }$$

The length scale β parameterizes the size of the vortex core and is one of three variational parameters $\left\{\beta,{{R}_{\bot}},\kappa \right\}$ with respect to which the total energy functional must be minimized in order to find their stationary values. Notice that this ansatz does not include ripples which are beyond the Thomas–Fermi approximation.

The energy functionals for the rotational, s-wave interaction, and trapping energies can all be evaluated analytically, albeit laboriously, using the above density profile. The results are given in [161] and will not be repeated here. Obtaining an analytic result for the dipolar interaction energy ${{E}_{\text{dd}}}=(1/2){\int}^{} \Phi \left(\mathbf{r}\right)n\left(\mathbf{r}\right)\text{d}\mathbf{r}$ is more difficult. However, in the hydrodynamic regime we have $\beta \ll {{R}_{\bot}}$, implying that the contribution to $\Phi \left(\mathbf{r}\right)$ from n_v is negligible in comparison to that from ${{n}_{\text{TF}}}$. Thus, to a very good approximation we can write,

$$\begin{eqnarray}&&{{E}_{\text{dd}}}\approx \frac{1}{2}{\int}^{}{{\Phi}_{\text{TF}}}\left(\mathbf{r}\right)~\left[{{n}_{\text{TF}}}\left(\mathbf{r}\right)+{{n}_{\text{v}}}\left(\mathbf{r}\right)\right]\text{d}\mathbf{r}\end{eqnarray} \tag{ 46 }$$

i.e. replace the true $\Phi \left(\mathbf{r}\right)$ by that purely due to the unperturbed background ${{\Phi}_{\text{TF}}}\left(\mathbf{r}\right)$ which is known analytically and is given in equation (30). Since ${{\Phi}_{\text{TF}}}\left(\mathbf{r}\right)$ is a quadratic function of the coordinates this integral can be done exactly [161]. Finally, to find the energy E_v associated with exciting a vortex it is necessary to subtract from E the energy E₀ of the vortex-free state, but this latter energy is also known analytically and is given as ${{E}_{\text{TF}}}$ in equation (29).

In this way E_v, and hence ${{\Omega}_{\text{v}}}$, can be computed for the trapped dipolar BEC. It is found that dipolar interactions increase ${{\Omega}_{\text{v}}}$ in prolate traps and lower it in oblate traps when compared to the pure s-wave case [161, 162]. Intuitively, this makes sense because in the prolate case dipolar interactions tend to reduce ${{R}_{\bot}}$ but in the oblate case they increase it. The rotational energy density $(1/2)n\left(\rho,z\right){{v}^{2}}\left(\rho \right)$ is lower at larger radii because $v\propto 1/\rho$ and hence E_v is lowered if atoms are moved to larger radii, like in the oblate case, and vice-versa in the prolate case. This interpretation is backed up by the following expression for the critical rotation frequency derived in the pure s-wave case in the Thomas–Fermi limit [163],

$$\begin{eqnarray}&&{{\Omega}_{\text{v}}}=\frac{5}{2}\frac{\hbar}{mR_{\bot}^{2}}\ln \frac{0.67{{R}_{\bot}}}{\xi}.\end{eqnarray} \tag{ 47 }$$

The numerical factors 5/2 and 0.67 arise from the inverted parabola of the Thomas–Fermi density profile: since the parabolic profile is maintained in the dipolar case we expect a similar expression to hold there. If ${{R}_{\bot}}$ in equation (47) is replaced by its dipolar version as given in equation (27), the resulting prediction for ${{\Omega}_{\text{v}}}$ is in very close agreement with the variational calculation described above [161]. In principle, the healing length should also be changed to its dipolar version but this only leads to a logarithmic correction.

We now review the vortex solutions in the simpler context of the homogeneous quasi-two-dimensional dipolar condensate, following the work of [129, 134]. Figure 8 [129] plots these solutions, found by numerical solution of the quasi-two-dimensional dipolar GPE, as a function of ${{\varepsilon}_{\text{dd}}}$.

Zoom In Zoom Out Reset image size

The vortex profile has a non-trivial dependence on the polarization angle and ${{\varepsilon}_{\text{dd}}}$. For $\alpha =0$, the vortex density is axisymmetric. For ${{\varepsilon}_{\text{dd}}}=0$ the condensate is non-dipolar and the vortex takes the standard form (see left inset of figure 8(a)) of a circularly-symmetric core of vanishing density at the centre that monotonically returns to its background value over a healing length $\xi =\hbar /\sqrt{m\mu}$ [108]. Since the length scaling applied in figure 8 is the dipolar healing length, the vortex core structure shown in figure 8(a) for ${{\varepsilon}_{\text{dd}}}\ne 0$ remains, for the most part, similar to that for ${{\varepsilon}_{\text{dd}}}=0$, i.e. the main effect of dipolar interactions is to rescale the size of the vortex core, but not change its structure. The exceptions to this are close to the phonon and roton instabilities. As the phonon instability boundary at ${{\varepsilon}_{\text{dd}}}=-0.5$ is approached from the stable side, the vortex core takes on a funnel-like profile (middle inset of figure 8(a)). This is associated with the cancellation of explicit s-wave van der Waals interactions in the system, i.e. the van der Waals interactions cancel the contact contribution from the dipolar interactions. The right inset of figure 8(a) also shows that as the roton instability is approached density ripples emerge around the vortex core. Moving away from the vortex core these ripples decay. The maximum amplitude of the density ripples is  ∼$20 \%$ of n₀, and their wavelength is of order the roton wavelength $\approx 4\xi$.

For $\alpha \ne 0$, see figures 8(b) and (c), the vortex profile becomes anisotropic. In particular, as the roton instability is approached, density ripples again form, but now aligned in the direction of the attractive dipolar interactions (along the polarization direction for ${{C}_{\text{dd}}}>0$ and perpendicular for ${{C}_{\text{dd}}}<0$). These anisotropic ripples are related to the anisotropic mixing of the roton into the ground state [130].

5.4.1. Dipolar mean-field potential due to a vortex: giant anti-dipoles.

Considered as a density defect in a homogenous background, the vortex gives rise to its own dipolar mean-field potential. Furthermore, because the creation of a vortex core displaces a large number of atomic dipoles, the vortex can be treated as a single giant anti-dipole [154, 161]. Take the case of a vortex within an otherwise uniform background of density n₀, with the density expressed using the decomposition as in equation (43). Then, by noting that the dipolar potential $\Phi$ is a linear functional of density, the total dipolar potential can be written as,

$$\begin{eqnarray}\begin{array}{*{35}{l}} \Phi [n]\left(\mathbf{r}\right) & = & \Phi \left[{{n}_{0}}+{{n}_{\text{v}}}\right]\left(\mathbf{r}\right) \\ {} & = & \Phi \left[{{n}_{0}}\right]\left(\mathbf{r}\right)+ \Phi \left[{{n}_{\text{v}}}\right]\left(\mathbf{r}\right)={{\Phi}_{0}}+ \Phi \left[{{n}_{\text{v}}}\right]\left(\mathbf{r}\right), \end{array}\end{eqnarray} \tag{ 48 }$$

i.e. a contribution ${{\Phi}_{0}}$, which is just a constant, from the uniform background n₀ and a spatially dependent contribution $\Phi \left[{{n}_{\text{v}}}\right]$ from the hole created by the vortex core, i.e. the giant anti-dipole. This decomposition, illustrated in figure 9, assists in understanding the interaction between vortices in dipolar systems.

Zoom In Zoom Out Reset image size

Consider a vortex in the quasi-two-dimensional dipolar condensate. For $\alpha =0$ and away from the roton and phonon instabilities, the vortex profile is well-approximated by the non-dipolar ansatz given by equations (42) and (43), where the healing length is taken to be the dipolar healing length $\xi =\hbar /\sqrt{m\mu}$. The dipolar potential generated by the density defect can be determined via the convolution result ${{\Phi}_{\bot}}(\boldsymbol{\rho},t)={{\mathcal{F}}^{-1}}[\tilde{U}_{\text{dd}}^{\bot}(\mathbf{q}){{\tilde{n}}_{\bot}}(\mathbf{q},t)]$. Due to the cylindrical symmetry of the $\alpha =0$ case, one can perform the Fourier transforms through Hankel transforms. The Hankel transform of the two-dimensional equivalent of equation (43) is,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{{\tilde{n}}}_{\bot}}\left(\mathbf{q}\right)={{n}_{0}}\left(\frac{\delta (q)}{q}+\frac{{{K}_{0}}\left(q\xi /\sqrt{2}\right)}{2}\right), \end{array}\end{eqnarray} \tag{ 49 }$$

where ${{K}_{0}}\left(\cdot \right)$ is a modified Bessel function of the second kind. Expanding the dipolar interaction potential $\tilde{U}_{\text{dd}}^{\bot}\left(\mathbf{q}\right)$ given in equation (34) with $\alpha =0$ as a series in the condensate width parameter, σ, gives

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \tilde{U}_{\text{dd}}^{\bot}\left(\mathbf{q}\right)=\frac{4\pi {{C}_{\text{dd}}}}{9\sqrt{2\pi}{{l}_{z}}}\left(2-\sqrt{\frac{9\pi}{2}}q{{l}_{z}}\sigma \right)+\mathcal{O}\left({{\sigma}^{2}}\right). \end{array}\end{eqnarray} \tag{ 50 }$$

Then, to first order in $\sigma ={{l}_{z}}/\xi$ and third order in $1/{{\rho}^{\prime}}=\xi /\rho$ ¹⁰, the dipolar potential due to a vortex is [129, 134],

$$\begin{eqnarray}&&{{\Phi}_{\bot}}\left(\boldsymbol{\rho}\right)={{\Phi}_{0}}\left[1-\frac{1}{{{\rho}^{\prime 2}}}+\left(\frac{A\ln {{\rho}^{\prime}}+B}{{{\rho}^{\prime 3}}}\right)\sigma \right],\end{eqnarray} \tag{ 51 }$$

with constants $A=-\sqrt{9\pi /8}\approx -1.88$, $B=\left(\ln 2-1\right)A\approx 0.577$, and ${{\Phi}_{0}}$ the dipolar potential at infinity.

This result indicates that the vortex causes a reduction in the mean-field dipolar potential which is consistent with the reduced density of dipoles in the vicinity of the vortex. One also sees from equation (51) that the dipolar potential generated by the vortex relaxes predominantly as $1/{{\rho}^{\prime 2}}$ to the background value ${{\Phi}_{0}}$, and this is confirmed by numerical solutions, as shown in figure 10(c) [129]. This dependence arises because the vortex density itself relaxes as $1/{{\rho}^{\prime 2}}$, and the leading contribution to ${{\Phi}_{\text{v}}}$ is from the local density. Indeed, the mean-field potential due to van der Waals interactions from the vortex also scales in proportion to the local density with a $1/{{\rho}^{\prime 2}}$ dependence [164]. The long-range contribution to ${{\Phi}_{\text{v}}}$ can be interpreted as arising from effective anti-dipoles in the vortex core, see figure 9. This non-local contribution is represented in equation (51) by terms linear in σ. In the limit $\sigma \to 0$ the volume of the anti-dipoles in the vortex core vanishes and hence this long-range contribution also vanishes. Unlike the topological potential associated with quantised superfluid flow around a vortex core the dipolar potential due to the vortex core is not topological, i.e. it depends on the volume of the vortex core. The dominant contribution to the non-local vortex potential scales as $\ln {{\rho}^{\prime}}/{{\rho}^{\prime 3}}$, i.e. the absence of dipoles in the vortex core can not be considered to be point-like. This is unsurprising since in the vicinity of the vortex core the density scales as a power law [108].

Zoom In Zoom Out Reset image size

For dipoles tilted away from the vertical ($\alpha \ne 0$), the vortex core and its dipolar potential become anisotropic and an analytic treatment is challenging. Figures 10(a) and (b) shows an example numerical solution of the vortex density and dipolar potential for $\alpha \ne 0$. The dipolar potential is indeed anisotropic about the vortex. Remarkably, the modification to the dipolar potential induced by the vortex mimics the dipole–dipole interaction itself, with an angular dependence which resembles $1-3{{\cos}^{2}}\theta$, being reduced (attractive) along y and increased (repulsive) along x. Thus, at least in its angular dependence, the vortex shares qualitatively the characteristics of a mesoscopic dipole. Like the above $\alpha =0$ case, the dipolar potential is found to decay at long-range as $1/{{\rho}^{\prime 2}}$ to the background value, as seen in figure 10(c). At short range, $\rho \lesssim 10\xi$, the dipolar potential is dominated by the core structure.

Two vortices (or indeed two elements of the same three-dimensional vortex line [154–156]) have a well-known hydrodynamic interaction due to the kinetic energy associated with the mutual cancellation/reinforcement of their velocity fields. Consider a cylindrical condensate of radius R and height L, featuring two straight vortices at planar positions $\boldsymbol{\rho}_{{1}}$ and $\boldsymbol{\rho}_{{2}}$, a distance d apart. The vortices have charge q₁ and q₂, and individual velocity fields ${{\mathbf{v}}_{1}}$ and ${{\mathbf{v}}_{2}}$, respectively. The net velocity field of the two vortices is ${{\mathbf{v}}_{1}}+{{\mathbf{v}}_{2}}$. The energy of the vortices can be estimated by integrating the total kinetic energy across the system,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{E}_{\text{kin}}}=L{\int}^{}\frac{1}{2}mn\left(\boldsymbol{\rho}\right)|{{\mathbf{v}}_{\mathbf{1}}}\left(\boldsymbol{\rho}\right)+{{\mathbf{v}}_{\mathbf{2}}}\left(\boldsymbol{\rho}\right){{|}^{2}}\,\text{d}\boldsymbol{\rho}. \end{array}\end{eqnarray} \tag{ 52 }$$

For simplicity, one can ignore the vortex core density and set $n\left(\boldsymbol{\rho}\right)=n$. Assuming $\xi \ll d\ll R$ then the (kinetic) energy of the pair is,

$$\begin{eqnarray}&&{{E}_{\text{kin}}}=\frac{L\pi n{{\hbar}^{2}}}{m}\left[q_{1}^{2}\ln \frac{R}{\xi}+q_{2}^{2}\ln \frac{R}{\xi}+2{{q}_{1}}{{q}_{2}}\ln \frac{R}{d}\right].\end{eqnarray} \tag{ 53 }$$

Note that, to avoid singularities in the velocity field the integration region excludes a disc of radius of one healing length about each vortex centre. The first two terms are the energies of the individual vortices if they were isolated, and the third term is the pair interaction energy. For a vortex–antivortex pair (${{q}_{1}}=-{{q}_{2}}$) the interaction energy is negative, whilst for a co-rotating pair (${{q}_{1}}={{q}_{2}}$) it is positive. This is explained physically by the fact that for a vortex–antivortex pair the flow fields tend cancel out in the bulk, reducing the net kinetic energy in the bulk, whilst for a co-rotating pair the flow fields tend to reinforce, increasing the total kinetic energy.

In the presence of dipolar interactions the vortices feature an additional long-range interaction. This interaction can be pictured as the interaction between two lumps of anti-dipoles in empty space, as illustrated in figure 11. Before we review how two vortices interact in the presence of dipolar interactions, we first make more precise the definition of the vortex energy introduced in section 5, and hence allow the identification of vortex–vortex interaction energy. In non-dipolar condensates, the vortex energy is conventionally defined as the energy difference between a system with and without a vortex, where both systems have the same number of particles [108]. Imagine first a system (quasi-two-dimensional) with a vortex and N particles covering an area A. If the asymptotic density is n₀, then the number of particles in this system can be expressed as,

$$\begin{eqnarray}&&N=A{{n}_{0}}-{{{\int}^{}}_{A}}\left({{n}_{0}}\,-\,|{{\psi}_{\bot}}{{|}^{2}}\right)\,\text{d}\boldsymbol{\rho}.\end{eqnarray} \tag{ 54 }$$

Zoom In Zoom Out Reset image size

Now consider the system without a vortex, but with the same number of particles. It has constant density n₀  =  N/A and its energy is

$$\begin{eqnarray}&&{{E}_{0}}=\frac{{{N}^{2}}}{2\sqrt{2\pi}{{l}_{z}}A}{{g}_{\text{eff}}},\end{eqnarray} \tag{ 55 }$$

where ${{g}_{\text{eff}}}=g+\left({{C}_{\text{dd}}}/3\right)\left[3{{\cos}^{2}}\alpha -1\right]$. Inserting equation (54) into equation (55) gives,

$$\begin{eqnarray}&&{{E}_{0}}\approx \frac{{{g}_{\text{eff}}}}{2\sqrt{2\pi}{{l}_{z}}}\left(An_{0}^{2}-2{{n}_{0}}\underset{A}{\int}\left({{n}_{0}}\,-\,|{{\psi}_{\bot}}{{|}^{2}}\right)\,\text{d}\boldsymbol{\rho}\right),\end{eqnarray} \tag{ 56 }$$

where a term, negligible in the limit $A\gg {{\xi}^{2}}$, has been omitted. Defining the energy of a single vortex as ${{E}_{1}}=E-{{E}_{0}}$, the interaction energy between two vortices must be [108],

$$\begin{eqnarray}&&{{E}_{12}}\left(\boldsymbol{\rho}_{{1}}-\boldsymbol{\rho}_{{2}}\right)={{E}_{2}}\left(\boldsymbol{\rho}_{{1}},\boldsymbol{\rho}_{{2}}\right)-{{E}_{1}}\left(\boldsymbol{\rho}_{{1}}\right)-{{E}_{1}}\left(\boldsymbol{\rho}_{{2}}\right),\end{eqnarray} \tag{ 57 }$$

where ${{E}_{2}}\left(\boldsymbol{\rho}_{{1}},\boldsymbol{\rho}_{{2}}\right)$ is the energy of the 2-vortex system. This is plotted in figure 12 for $\alpha =0$ and for (a) vortex–antivortex (VA) pairs and (b) vortex–vortex (VV) pairs, as a function of their separation d, based on numerical dipolar GPE two-vortex solutions [134].

Zoom In Zoom Out Reset image size

In the absence of dipoles, E₁₂ increases with d for the VA pair and decreases for the VV pair. For $d\gg \xi$, E₁₂ closely follows the logarithmic scaling of the hydrodynamic prediction, while for $d\lesssim \xi$ the overlap of the cores causes a breakdown of the logarithmic behaviour. With dipoles, E₁₂ is significantly modified at short and intermediate length scales up to $d\approx 5\xi$, but at larger scales the effects of the dipoles are small in comparison to the hydrodynamic effects. The modification due to the dipoles arises from a non-trivial combination of the modified density profile and non-local interactions.

When the dipoles are tilted in the plane, E₁₂ becomes dependent on the in-plane angle of the pair relative to the polarization direction, η. This is illustrated for VA and VV pairs by the examples in figure 13 [134]. For small separations, the angular dependence is dominated by local effects arising from the density profile of the pairs, particularly by any density ripples. However, for $d\gg \xi$, $\Delta E_{12}^{\prime}={{E}_{12}}\left(\eta \right)-{{E}_{12}}\left(\eta =0\right)$ approaches a sinusoidal dependence on η, analogous to the dipolar interaction itself.

Zoom In Zoom Out Reset image size

In a two-dimensional system, a vortex co-moves with the local fluid velocity. Thus, for a vortex pair, each vortex is carried along by the flow field of the other vortex. For a VA pair this means that the vortices move in the same direction, perpendicular to the inter-vortex axis. This solitary wave has speed $v=\hbar /md$ for well-separated vortices.

When the VA separation is small, $d\sim \xi$, the vortex and anti-vortex are susceptible to annihilation, an event which results in a burst of density waves. Numerical simulations [134] show that dipoles modify this separation threshold (figure 14(a)). Moreover, since for $\alpha \ne 0$ the speed of sound varies with angle [130], one can expect the pair speed to be anisotropic in space.

Zoom In Zoom Out Reset image size

For a VV pair, the flow which carries each vortex now acts in opposite directions (again, perpendicular to the line separating the vortices and with the above speed), resulting in the co-rotation of the vortices about their mid-point. Viewed another way, the vortices follow a path of constant energy; since the interaction energy in the absence of dipoles depends only on d, this path is circular. The same is true for axisymmetric dipoles, $\alpha =0$. However, for $\alpha \ne 0$, the vortices co-rotate on an anisotropic path, as shown in figure 14(b) [134], due to the anisotropic interaction energy of the pair (refer to figure 13). Moreover, for extreme cases where the vortex is highly elliptical with significant ripples, co-rotation can be prevented altogether, with the vortices being localized (figure 14(b), red solid lines). In this limit, the vortices act as extended, highly-elongated objects, with effective geometrical restrictions on their motion past each other, reminiscent of the smectic phase of liquid crystals.

Gautam [153] numerically considered the dynamics of a corotating VV pair in a dipolar BEC in the presence of dissipation. For symmetric configurations in the trap, the vortices decay with equal decay times, while for asymmetric initial configurations the decay is modified with one vortex decaying slower at the expense of the other.

In conventional condensates, vortices have been generated through several mechanisms. Below we list the main ones, as well as relevant considerations in the presence of dipoles.

• The most common and intuitive approach to generate vortices is via mechanical rotation of the system [35, 40, 41], analogous to the 'rotating bucket' experiments in Helium II [118]. Both the thermodynamic threshold for vortices to be favoured, as well as the process by which vortices nucleate [165–167] into the condensate, are sensitive to dipolar interactions; this will be analysed in detail below.
• Motion of a localized obstacle or potential (as generated by a tightly-focussed blue-detuned laser beam) through a condensate (or, equivalently, motion of the condensate relative to a static obstacle) leads to the nucleation of vortices above a critical relative speed [37, 38, 142], forming a quantum wake downstream of the obstacle. The critical speed is related to the Landau criterion which predicts the formation of elementary excitations in the fluid for relative speeds exceeding ${{v}_{\text{c}}}=\text{min}\left[\omega \left(\mathbf{k}\right)/k\right]$ [108, 109]. Ticknor et al [130] examined this process in a quasi-two-dimensional dipolar condensate. For dipoles tilted into the plane ($\alpha >0$), the critical speed becomes anisotropic, a consequence of Landau's criterion and the anisotropic dispersion relation in the plane. The critical velocity for vortex nucleation can also be derived by considering the energetics of a vortex–antivortex pair [107], implying that the aniostropic critical velocity is directly related to the anisotropic vortex interaction energy.
• The phase of the condensate can be directly engineered via optical imprinting to produce vortex phase singularities, as employed to generate both singly- and multiply-charged vortices [168]. This mechanism is independent of the dipoles themselves.
• Following a rapid quench through the transition temperature for the onset of Bose–Einstein condensation, the growth of local phase-coherent domains leads to the entrapment of phase singularities and hence vortices (i.e. the Kibble–Zurek mechanism [169, 170]) [36, 46, 171]. A relevant consideration is the effect of the dipoles on the critical temperature. This shift is sensitive to the shape of the trap, relative to the polarization direction, but is only up to a few percent for ⁵²Cr [172]. However, this shift may be more significant in ¹⁶⁸Er and ¹⁶⁴Dy condensates.
• Dark solitons are dimensionally unstable to decay into vortex pairs or vortex rings [39, 173, 174] via the so-called snake instability, previously established in nonlinear optics [175]. As shown by Nath et al [176] the nonlocal character of dipolar interactions can stabilize the dark soliton against this instability.
• Instead of rotating the system it is possible to introduce a time-reversal symmetry breaking synthetic magnetic field [177]. Since the atoms in an atomic BEC are charge-neutral it may at first seem counter-intuitive to consider the effects of a synthetic magnetic field in these systems. However, time-reversal symmetry-breaking in charge-neutral systems can be overcome through mechanical rotation and exploiting the equivalence of the Coriolis force and the Lorentz force to create a synthetic vector potential which gives rise to a synthetic magnetic field. This case of mechanical rotation is analysed in detail below. It is also possible to realise an optically synthesised vector potential field for BECs using spatially dependent optical coupling between internal states of the atoms in the condensate. This spatially dependent coupling can yield Berry phases [178] sufficient to create large synthetic magnetic fields. As in the case of mechanical rotation, vortex nucleation is dependent on the properties of the stationary solutions and can be analyzed in the Thomas–Fermi approximation [179]. Numerical investigations of the dipolar GPE, carried out by Zhao and Gu [180], find that the nucleation of vortices depends on the dipole strength, the strength of the synthetic magnetic field, the potential geometry, and the orientation of the dipoles, with anisotropic interactions significantly altering vortex nucleation.

The most common approach for generating vortices and vortex lattices in trapped condensates is via rotation. Since a cylindrically-symmetric trap set into rotation applies no torque to the condensate, the trap is made anisotropic in the plane of rotation. In the simplest case, this leads to a trap which is weakly elliptical in the plane of rotation [35, 40], with a potential of the form,

$$\begin{eqnarray}&&V\left(\mathbf{r}\right)=\frac{1}{2}m\omega _{\bot}^{2}\left[\left(1-\epsilon \right){{x}^{2}}+\left(1+\epsilon \right){{y}^{2}}+{{\gamma}^{2}}{{z}^{2}}\right],\end{eqnarray} \tag{ 58 }$$

where rotation is performed about the z-axis.

For typical parameters in the absence of dipolar interactions, vortices become energetically favourable in harmonically-trapped condensates for rotation frequencies $\Omega \sim 0.3{{\omega}_{\bot}}$. Surprisingly, in non-dipolar BEC experiments the observed nucleation of vortices occurs at considerably larger rotation frequencies $\Omega \sim 0.7{{\omega}_{\bot}}$. Theoretical analysis based on the hydrodynamic equations reveals the important role of collective modes. Specifically, for $\Omega \lesssim 0.7{{\omega}_{\bot}}$ low-lying collective modes are excited via elliptical deformation. The seeding of vortices, at higher rotation frequencies, arises when one or more of these modes becomes unstable [181, 182]. Evidence for this comes from comparison between experiments [40, 183] and full numerical simulations of the GPE [184–187].

The hydrodynamic description of condensates in rotating elliptical traps can be extended to include dipolar interactions. For rotation about the z-axis, described by the rotation vector $\mathbf{\Omega}$ where $\Omega =|\mathbf{\Omega}|$ is the rotation frequency, the Hamiltonian in the rotating frame is given by,

$$\begin{eqnarray}&&{{H}_{\text{eff}}}={{H}_{0}}-\mathbf{\Omega}\cdot \hat{L},\end{eqnarray} \tag{ 59 }$$

where the Hamiltonian in the absence of rotation is H₀ and the quantum mechanical angular momentum operator is $\hat{L}=-\text{i}\hbar \left(\mathbf{r}\times \boldsymbol{\nabla} \right)$. Applying this result for $\mathbf{\Omega}=\left(0,0, \Omega \right)$, the dipolar GPE in the rotating frame is [188, 189],

$$\begin{eqnarray}\begin{array}{*{35}{l}} \text{i}\hbar \frac{\partial \Psi}{\partial t} & = & \left[-\frac{{{\hbar}^{2}}}{2m}{{\nabla}^{2}}+V+ \Phi +g| \Psi {{|}^{2}} \right. \\ {} && \left. \, - \Omega \frac{\hbar}{i}\left(x\frac{\partial}{\partial y}-y\frac{\partial}{\partial x}\right)\right] \Psi, \end{array}\end{eqnarray} \tag{ 60 }$$

where in the rotating frame the trapping potential V, given by equation (58), is stationary. The spatial coordinates $\mathbf{r}$ are those of the rotating frame, with the momentum coordinates expressed in the laboratory frame [188–190].

Using the Madelung transform, as per section 3.3, leads to the following dipolar hydrodynamical equations in the rotating-frame

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{\partial n}{\partial t}=-\boldsymbol{\nabla} \cdot \left[n\left(\mathbf{v}-\mathbf{\Omega}\times \mathbf{r}\right)\right], \end{array}\end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} m\frac{\partial \mathbf{v}}{\partial t}=-\boldsymbol{\nabla} \left(\frac{1}{2}m{{v}^{2}}+V+gn+ \Phi -m\mathbf{v}\cdot \left[\mathbf{\Omega}\times \mathbf{r}\right]\right), \end{array}\end{eqnarray} \tag{ 62 }$$

where the quantum pressure is assumed to be small and is neglected, i.e. the Thomas–Fermi limit.

Stationary solutions of equations (61) and (62) satisfy the equilibrium conditions,

$$\begin{eqnarray}&&\frac{\partial n}{\partial t}=0~~\text{and}~~\frac{\partial \mathbf{v}}{\partial t}=0.\end{eqnarray} \tag{ 63 }$$

Assuming the irrotational ($\nabla \times \mathbf{v}=0$) velocity field ansatz [181]

$$\begin{eqnarray}&&\mathbf{v}={{\alpha}_{\text{v}}}\nabla (xy)\end{eqnarray} \tag{ 64 }$$

permits us to examine the rotating solutions in terms of the velocity field amplitude ${{\alpha}_{\text{v}}}$. A physical interpretation of velocity field amplitude, ${{\alpha}_{\text{v}}}$, can be gained from the continuity equation (61). Specifically, it can be written as [109]

$$\begin{eqnarray}&&{{\alpha}_{\text{v}}}=-\mathcal{D} \Omega,\end{eqnarray} \tag{ 65 }$$

where $\mathcal{D}$ is the deformation of the BEC in the x  −  y plane,

$$\begin{eqnarray}&&\mathcal{D}=\frac{\langle \,{{y}^{2}}-{{x}^{2}}\rangle}{\langle \,{{y}^{2}}+{{x}^{2}}\rangle}=\frac{\kappa _{y}^{2}-\kappa _{x}^{2}}{\kappa _{y}^{2}+\kappa _{x}^{2}},\end{eqnarray} \tag{ 66 }$$

where $\langle \ldots \rangle$ denotes the expectation value in the stationary state and ${{\kappa}_{x}}={{R}_{x}}/{{R}_{z}}$ and ${{\kappa}_{y}}={{R}_{y}}/{{R}_{z}}$ represent the aspect ratios of the BEC along x and y with respect to z.

Combining equations (62) and (64) gives,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mu =\frac{m}{2}(\tilde{\omega}_{x}^{2}{{x}^{2}}+\tilde{\omega}_{y}^{2}{{y}^{2}}+\omega _{z}^{2}{{z}^{2}})+gn(\mathbf{r})+ \Phi (\mathbf{r}), \end{array}\end{eqnarray} \tag{ 67 }$$

where the effective trap frequencies ${{\tilde{\omega}}_{x}}$ and ${{\tilde{\omega}}_{y}}$ are given by,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \tilde{\omega}_{x}^{2}=\omega _{\bot}^{2}(1-\epsilon )+\alpha _{\text{v}}^{2}-2{{\alpha}_{\text{v}}} \Omega, \end{array}\end{eqnarray} \tag{ 68 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \tilde{\omega}_{y}^{2}=\omega _{\bot}^{2}(1+\epsilon )+\alpha _{\text{v}}^{2}+2{{\alpha}_{\text{v}}} \Omega. \end{array}\end{eqnarray} \tag{ 69 }$$

The breaking of cylindrical symmetry means that the BEC has an ellipsoidal shape and in the Thomas–Fermi approximation it adopts an inverted parabolic density profile of the form,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} n\left(\mathbf{r}\right)={{n}_{\text{cd}}}\left(1-\frac{{{x}^{2}}}{R_{x}^{2}}-\frac{{{y}^{2}}}{R_{y}^{2}}-\frac{{{z}^{2}}}{R_{z}^{2}}\right)\,\,\,\,\text{for}\,\,\,n\left(\mathbf{r}\right)\geqslant 0, \end{array}\end{eqnarray} \tag{ 70 }$$

where ${{n}_{\text{cd}}}=15N/\left(8\pi {{R}_{x}}{{R}_{y}}{{R}_{z}}\right)$. This is an exact solution of the stationary dipolar hydrodynamic equations given in equations (61) and (62) with velocity field (64), and it only remains to find the radii $\left\{{{R}_{x}},{{R}_{y}},{{R}_{z}}\right\}$ and velocity amplitude ${{\alpha}_{\text{v}}}$. The exact dipolar potential due to a parabolic density distribution of general ellipsoidal symmetry, with dipoles aligned in the z-direction, i.e. the generalized version of that given in equation (30), is derived in the Appendices of [92] and [117], and is given by,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & \Phi =3g{{\varepsilon}_{\text{dd}}} \\ {} & \times \,\left(\frac{{{n}_{\text{cd}}}{{\kappa}_{x}}{{\kappa}_{y}}}{2}\left[{{\beta}_{001}}-\frac{{{x}^{2}}{{\beta}_{101}}+{{y}^{2}}{{\beta}_{011}}+3{{z}^{2}}{{\beta}_{002}}}{R_{z}^{2}}\right]-\frac{n}{3}\right), \end{array}\end{eqnarray} \tag{ 71 }$$

where the coefficients ${{\beta}_{ijk}}$ are,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\beta}_{ijk}}={\int}_{0}^{\infty}\frac{1}{{{\left(\kappa _{x}^{2}+s\right)}^{i+\frac{1}{2}}}{{\left(\kappa _{y}^{2}+s\right)}^{j+\frac{1}{2}}}{{(1+s)}^{k+\frac{1}{2}}}}\,\text{d}s, \end{array}\end{eqnarray} \tag{ 72 }$$

for integer-valued i, j and k. Thus, equation (67) can be rearranged to obtain an expression for the density profile [165–167],

$$\begin{eqnarray}\begin{array}{*{35}{l}} n & = & \frac{\mu -\frac{m}{2}(\tilde{\omega}_{x}^{2}{{x}^{2}}+\tilde{\omega}_{y}^{2}{{y}^{2}}+\omega _{z}^{2}{{z}^{2}})}{g(1-{{\varepsilon}_{\text{dd}}})} \\ {} & + & \frac{3g{{\varepsilon}_{\text{dd}}}\frac{{{n}_{\text{cd}}}{{\kappa}_{x}}{{\kappa}_{y}}}{2R_{z}^{2}}[{{x}^{2}}{{\beta}_{101}}+{{y}^{2}}{{\beta}_{011}}+3{{z}^{2}}{{\beta}_{002}}-R_{z}^{2}{{\beta}_{001}}]}{g(1-{{\varepsilon}_{\text{dd}}})}. \end{array}\end{eqnarray} \tag{ 73 }$$

By equating the x², y² and z² terms in equations (70) and (73) three self-consistency conditions are found. These conditions define the size and shape of the condensate,

$$\begin{eqnarray}\begin{array}{*{35}{l}} \kappa _{x}^{2} & = & {{\left(\frac{{{\omega}_{z}}}{{{{\tilde{\omega}}}_{x}}}\right)}^{2}}\frac{1+{{\varepsilon}_{\text{dd}}}\left(\frac{3}{2}\kappa _{x}^{3}{{\kappa}_{y}}{{\beta}_{101}}-1\right)}{\zeta}, \end{array}\end{eqnarray} \tag{ 74 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \kappa _{y}^{2} & = & {{\left(\frac{{{\omega}_{z}}}{{{{\tilde{\omega}}}_{y}}}\right)}^{2}}\frac{1+{{\varepsilon}_{\text{dd}}}\left(\frac{3}{2}\kappa _{y}^{3}{{\kappa}_{x}}{{\beta}_{011}}-1\right)}{\zeta}, \end{array}\end{eqnarray} \tag{ 75 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} R_{z}^{2} & = & \frac{2g{{n}_{\text{cd}}}}{m\omega _{z}^{2}}\zeta, \end{array}\end{eqnarray} \tag{ 76 }$$

where $\zeta =1-{{\varepsilon}_{\text{dd}}}\left[1-\frac{9{{\kappa}_{x}}{{\kappa}_{y}}}{2}{{\beta}_{002}}\right]$. Furthermore, by inserting equation (73) into equation (61), the stationary solutions are seen to satisfy the condition [165–167],

$$\begin{eqnarray}\begin{array}{*{35}{l}} 0 & = & \left({{\alpha}_{\text{v}}}+ \Omega \right)\left(\tilde{{\omega}}_{x}^{2}-\frac{3}{2}{{\varepsilon}_{\text{dd}}}\frac{\omega _{\bot}^{2}{{\kappa}_{x}}{{\kappa}_{y}}{{\gamma}^{2}}}{\zeta}{{\beta}_{101}}\right) \\ {} & {} & +\,\,\left({{\alpha}_{\text{v}}}- \Omega \right)\left(\tilde{{\omega}}_{y}^{2}-\frac{3}{2}{{\varepsilon}_{\text{dd}}}\frac{\omega _{\bot}^{2}{{\kappa}_{x}}{{\kappa}_{y}}{{\gamma}^{2}}}{\zeta}{{\beta}_{011}}\right). \end{array}\end{eqnarray} \tag{ 77 }$$

Equation (77) gives the velocity field amplitude ${{\alpha}_{\text{v}}}$ for a given ${{\varepsilon}_{\text{dd}}}$, $\Omega$ and trap geometry. For ${{\varepsilon}_{\text{dd}}}=0$, ${{\alpha}_{\text{v}}}$ is independent of the s-wave interaction strength g and the trap ratio γ. Dipolar interactions qualitatively alter this scenario with ${{\alpha}_{\text{v}}}$ becoming dependent on ${{\varepsilon}_{\text{dd}}}$ and γ. The solutions to equation (77), as a function of $\Omega$, have significantly different properties depending on whether the traps are circular ($\epsilon =0$) or elliptical ($\epsilon >0$) in the x  −  y plane. We restrict our analysis to $\Omega <{{\omega}_{\bot}}$, since the static solutions are known to disappear for $\Omega \sim {{\omega}_{\bot}}$ due to centre of mass instabilities [181]. Below, the circular ($\epsilon =0$) and elliptical ($\epsilon >0$) cases are considered.

7.2.1. Circular trapping in the x  −  y plane: $\epsilon =0$.

In figure 15(a) [167] the solutions of equation (77) are plotted as a function of rotation frequency $\Omega$ for a spherically-symmetric trap, $\gamma =1$ and $\epsilon =0$, for various values of ${{\varepsilon}_{\text{dd}}}$. For a given value of ${{\varepsilon}_{\text{dd}}}$ the solutions have the same qualitative structure. Specifically, only one solution exists (${{\alpha}_{\text{v}}}=0$) up to some critical frequency. Two additional solutions (${{\alpha}_{\text{v}}}>0$ and ${{\alpha}_{\text{v}}}<0$) bifurcate from this single solution at the critical rotation frequency, denoted as the bifurcation frequency ${{\Omega}_{\text{b}}}$.

Zoom In Zoom Out Reset image size

As expected, when ${{\varepsilon}_{\text{dd}}}=0$ the results of [181, 182] are reproduced with ${{\Omega}_{\text{b}}}={{\omega}_{\bot}}/\sqrt{2}$ and ${{\alpha}_{\text{v}}}=\pm \sqrt{2{{\Omega}^{2}}-\omega _{\bot}^{2}}$ for the two additional solutions when $\Omega >{{\Omega}_{\text{b}}}$. The critical frequency ${{\Omega}_{\text{b}}}$ is associated with the spontaneous excitation of quadrupole modes. Specifically, in the Thomas–Fermi regime, the surface excitation dispersion is given by (see p 183 of [109])

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \omega _{l}^{2}=\left({{q}_{l}}/m\right){{\nabla}_{R}}V=\left({{q}_{l}}\omega _{\bot}^{2}/2\right){{\nabla}_{R}}\left({{R}^{2}}\right), \end{array}\end{eqnarray} \tag{ 78 }$$

where R is the radius of the BEC. For a surface excitation with angular momentum $\hbar l=\hbar {{q}_{l}}R$, in the absence of rotation, equation (78) reduces to ${{\omega}_{l}}=\sqrt{l}{{\omega}_{\bot}}$. Rotation shifts the mode frequency by $-l \Omega$, see equation (59). Hence for a rotating BEC the quadrupole surface excitation (l  =  2) frequency is ${{\omega}_{2}}\left(\Omega \right)=\sqrt{2}{{\omega}_{\bot}}-2 \Omega$ [109]. The bifurcation frequency, ${{\Omega}_{\text{b}}}={{\omega}_{\bot}}/\sqrt{2}$, occurs at the same rotation frequency at which the energy of the quadrupole mode is zero. Within the context of the Thomas–Fermi approximation, this critical frequency is independent of the strength of the contact interactions (g) and for $\Omega \geqslant {{\omega}_{\bot}}/\sqrt{2}$ the two additional solutions arise from the excitation of the quadrupole mode.

Referring back to equation (78), the inclusion of non-local dipolar interactions implies that the force $-\nabla V$ no longer has a simple dependance on R [91]. Hence, there is no reason to suspect that the condition to excite the quadrupole mode will be independent of ${{\varepsilon}_{\text{dd}}}$. Indeed, the dependence on ${{\varepsilon}_{\text{dd}}}$ can be seen in figure 15(a) where the introduction of dipolar interactions leads to a shift in ${{\Omega}_{\text{b}}}$. Specifically, for ${{\varepsilon}_{\text{dd}}}>0$ (${{\varepsilon}_{\text{dd}}}<0$) the bifurcation frequency is increased (decreased). It is possible to evaluate ${{\Omega}_{\text{b}}}$ analytically by realising that ${{\kappa}_{x}}={{\kappa}_{y}}=\kappa$, for ${{\alpha}_{\text{v}}}=0$. In this limit the aspect ratio κ is determined by the transcendental equation (25) [91, 92]. As ${{\alpha}_{\text{v}}}\to {{0}_{+}}$, the first order corrections to ${{\kappa}_{x}}$ and ${{\kappa}_{y}}$ with respect to κ from equations (74) and (75) can be calculated and inserted into equation (77). Solving for $\Omega$ (noting that $\Omega \to {{\Omega}_{\text{b}}}$ as ${{\alpha}_{\text{v}}}\to {{0}_{+}}$) gives,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{{{\Omega}_{\text{b}}}}{{{\omega}_{\bot}}}=\sqrt{\frac{1}{2}+\frac{3}{4}{{\kappa}^{2}}{{\varepsilon}_{\text{dd}}}{{\gamma}^{2}}\frac{{{\kappa}^{2}}{{\beta}_{201}}-{{\beta}_{101}}}{1-{{\varepsilon}_{\text{dd}}}\left(1-\frac{9}{2}{{\kappa}^{2}}{{\beta}_{002}}\right)}}, \end{array}\end{eqnarray} \tag{ 79 }$$

which is plotted in figure 15(b) as a function of γ for various values of ${{\varepsilon}_{\text{dd}}}$. When ${{\varepsilon}_{\text{dd}}}=0$ then ${{\Omega}_{\text{b}}}={{\omega}_{x}}/\sqrt{2}$, which is independent of $\gamma ={{\omega}_{z}}/{{\omega}_{x}}$ [181, 182]. However, for ${{\varepsilon}_{\text{dd}}}\ne 0$ the value of γ for which ${{\Omega}_{\text{b}}}$ reaches a minimum changes. Specifically, as ${{\varepsilon}_{\text{dd}}}$ is increased from  −0.5 the minimum value of ${{\Omega}_{\text{b}}}$ changes from occurring at values of γ where the trap shape is oblate ($\gamma >1$) to ones where it is prolate ($\gamma <1$). Fixing γ and increasing ${{\varepsilon}_{\text{dd}}}$ leads to a monotonic decrease in ${{\Omega}_{\text{b}}}$. Increasing ${{\varepsilon}_{\text{dd}}}$ can lead to a significant reduction in the bifurcation frequency, i.e. for ${{\varepsilon}_{\text{dd}}}=0.99$, ${{\Omega}_{\text{b}}}$ is reduced by $20 \%$ compared to its non-dipolar value to $\approx 0.55{{\omega}_{\bot}}$. It is tempting to consider the case where ${{\varepsilon}_{\text{dd}}}>1$ to induce shifts in ${{\Omega}_{\text{b}}}$ to even lower values, however, in these regimes the premise of the calculation, i.e. the Thomas–Fermi approximation, may not be valid.

7.2.2. Elliptical trapping in the x  −  y plane: $\epsilon >0$.

Rotating elliptical traps have been created experimentally with lasers and magnetic fields [35, 40]. Following the experiment by Madison et al [35], below we consider a trap with a weak ellipticity of $\epsilon =0.025$. Figure 16(a) [167] shows the solutions to equation (77) for various values of ${{\varepsilon}_{\text{dd}}}$ for $\gamma =1$. As in the non-dipolar case [181, 182], the solutions become heavily modified for $\epsilon >0$. There is an upper branch solution (${{\alpha}_{\text{v}}}>0$) which extends over the range $0\leqslant \Omega \leqslant {{\omega}_{\bot}}$ and a lower solution (${{\alpha}_{\text{v}}}<0$) which is doubled valued and exists above some critical rotation frequency. The critical rotation frequency for the lower brach solution is denoted as the back-bending frequency ${{\Omega}_{\text{b}}}$, which in the limit $\epsilon =0$ can be regarded as the limiting case of the bifurcation frequency and hence we use the same notation for both. Unlike the circular trap, considered in section 7.2.1, there are no ${{\alpha}_{\text{v}}}=0$ solutions for $\Omega >0$. In the absence of dipolar interactions increasing the trap ellipticity results in an increase of the back-bending frequency ${{\Omega}_{\text{b}}}$. As shown in figure 16(b) dipolar interactions, as in the case of $\epsilon =0$, reduce (increase) ${{\Omega}_{\text{b}}}$ for ${{\varepsilon}_{\text{dd}}}>0$ (${{\varepsilon}_{\text{dd}}}<0$).

Zoom In Zoom Out Reset image size

Due to the anisotropy of the dipolar interactions increasing ${{\varepsilon}_{\text{dd}}}$ decreases both ${{\kappa}_{x}}$ and ${{\kappa}_{y}}$, i.e. the BEC becomes more prolate. Under rotation dipolar interactions also increase the deformation of the BEC in the x  −  y plane, as can be deduced from figure 16(a). Specifically, as ${{\varepsilon}_{\text{dd}}}$ is increased, for $\Omega >0$, ${{\alpha}_{\text{v}}}$ increases and hence, see equation (65), the deformation ($\mathcal{D}$) of the dipolar BEC increases.

The static solutions derived above are stationary but not necessarily stable. In this section the dynamical stability of the stationary solutions is analyzed. This is done by considering small perturbations in the BEC density and phase of the form $n={{n}_{\text{eq}}}+\delta n$ and $S={{S}_{\text{eq}}}+\delta S$. By linearizing the dipolar hydrodynamic equations (61) and (62), the dynamics of such perturbations can be described as [165–167],

$$\begin{eqnarray}&&\frac{\partial}{\partial t}\left[\begin{array}{c} \delta S \\ \delta n \end{array}\right]\,=-\left[\begin{array}{c} {{\mathbf{v}}_{\mathbf{c}}}\cdot \boldsymbol{\nabla} & \frac{g}{m}\left(1+{{\varepsilon}_{\text{dd}}}K\right) \\ \boldsymbol{\nabla} \cdot {{n}_{\text{eq}}}\boldsymbol{\nabla} & \boldsymbol{\nabla} \cdot \mathbf{v}+{{\mathbf{v}}_{\mathbf{c}}}\cdot \boldsymbol{\nabla} \end{array}\right]\left[\begin{array}{c} \delta S \\ \delta n \end{array}\right]\end{eqnarray} \tag{ 80 }$$

where ${{\mathbf{v}}_{\mathbf{c}}}=\mathbf{v}-\mathbf{\Omega}\times \mathbf{r}$ and the integral operator K is defined as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \left(K\delta n\right)\left(\mathbf{r}\right)=-3\frac{{{\partial}^{2}}}{\partial {{z}^{2}}}{\int}^{}\frac{\delta n\left({{\mathbf{r}}^{\prime}}\right)}{4\pi |\mathbf{r}-{{\mathbf{r}}^{\prime}}|}\,\text{d}{{\mathbf{r}}^{\prime}}-\delta n\left(\mathbf{r}\right). \end{array}\end{eqnarray} \tag{ 81 }$$

To investigate the stability of the BEC the eigenfunctions and eigenvalues of the operator in equation (80) can be examined. Dynamical instability arises when one or more eigenvalues λ possess a positive real part. While the imaginary parts of the eigenvalues characterise the frequencies of collective modes of the BEC [191], the magnitudes of any real parts control the rate of growth of unstable modes. If the density and phase fluctuations are expressed as polynomials of degree N, the operators in equation (80), including K [93, 94, 96], result in polynomials which are of degree N or less. This enables equation (80) to be recast as a scalar matrix operator which acts on vectors of polynomial coefficients. Finding the eigenvalues (determining stability) and eigenvectors (characterising modes) is then a simple computational task [165–167, 182].

Below, the stability of the rotating solutions are considered for $\epsilon >0$ (figure 16(a)). For ${{\alpha}_{\text{v}}}<0$ there are two static solutions for $\Omega >{{\Omega}_{\text{b}}}$. The solution nearest the ${{\alpha}_{\text{v}}}=0$ axis is always dynamically stable except when $\Omega \to {{\omega}_{\bot}}$, where the BEC is susceptible to a centre-of-mass instability [192]. The other solution, in the ${{\alpha}_{\text{v}}}<0$ half-plane, is always dynamically unstable and hence experimentally irrelevant. The dynamical stability of solutions in the upper half-plane (${{\alpha}_{\text{v}}}>0$) is more interesting. In figure 17 [167] the maximum positive real eigenvalues of the upper-branch solution is plotted as a function of $\Omega$ for a fixed trap geometry and various values of ${{\varepsilon}_{\text{dd}}}$. To obtain these results a maximum polynomial perturbation of $\delta n={{x}^{p}}{{y}^{q}}{{z}^{r}}$ with $p+q+r\leqslant 3=N$ was considered.

Zoom In Zoom Out Reset image size

In figure 17(inset) an example of the unstable region (for the upper-branch solution) is shown, in the $\epsilon - \Omega$ plane, for ${{\varepsilon}_{\text{dd}}}=0$ and $\gamma =1$. The shaded crescents [182] denote the regimes of dynamical instability ($\text{Re}\left(\lambda \right)>0$). Each crescent corresponds to a single value of N, the total degree of the polynomial perturbation. Each higher value of N adds another crescent from above. The crescents merge for large rotation frequencies, with the eigenvalues being comparatively large. At lower rotation frequencies the crescents become vanishingly thin with comparatively small eignevalues [187] (at least an order of magnitude smaller than the main instability region). The relative smallness of the eigenvalues in this region indicates that instabilities grow over a much longer time-scale as compared to the main region of instability. For non-dipolar BECs this was numerically investigated by solving the GPE [187]. These numerical results showed that the narrow instability regions have negligible effect when ramping $\Omega$ at rates greater than $\text{d} \Omega /\text{d}t=2\times {{10}^{-4}}\omega _{\bot}^{2}$. Hence these narrow regions of instability have minimal consequence over the time-scales typically considered in an experiment and can be ignored, with the main region of instability denoted by the dashed line in figure 17(inset). Experimental trap ellipticities are usually $\epsilon \leqslant 0.1$ and the unstable regime can be quantified with N  =  3 for the perturbations. Denoting the lower bound for the dynamical instability as ${{\Omega}_{\text{i}}}$, figure 17 indicates that as ${{\varepsilon}_{\text{dd}}}$ is increased the lower bound for ${{\Omega}_{\text{i}}}$ decreases.

As the size of the scalar matrix operator (80) is increased to $N=4,5,\ldots$, the higher lying modes develop real eigenvalues as $\Omega$ is increased. These higher lying modes fall within the region of instability already shown in figure 17 for N  =  3 and so do not alter the range of parameters where instability occurs. As we shall now explain, this result implies that the Thomas–Fermi spectrum does not contain a roton minimum as a function of angular momentum L. Ultimately, this is because a roton minimum means, somewhat counter-intuitively, that higher lying modes can have lower energy than lower lying modes. This has important consequences in the presence of fluid flow. For example, Pitaevskii [193] considered the case of superfluid ⁴He flowing through a pipe with rough walls. If the excitation spectrum of the superfluid at rest is E(p), when it flows at speed v the spectrum is Galilean-shifted $E\left(\,p\right)\to E\left(\,p\right)-pv$. There is, therefore, a critical velocity ${{v}_{c}}=\text{min}~E\left(\,p\right)/p$, where min means the minimum with respect to p, at which excitations have zero energy and can be freely produced, depleting the superfluid. However, because superfluid ⁴He has a roton minimum in its spectrum, the excitations created at v_c have a specific momentum $p\approx {{p}_{r}}$, where p_r is the roton minimum, triggering an instability to the formation of a density wave with wavelength $\propto p_{r}^{-1}$. For a rotating system similar arguments can be made with the Galilean-shifted energy taking the form $E\to E-L \Omega$ with angular roton modes [119] becoming unstable at some critical rotation frequency. Crucially, in the dynamical instability analysis presented above it was found that the modes become unstable in order as $\Omega$ is increased. This implies that for the regimes considered, an angular-roton minimum at finite angular momentum does not exist. Within the context of the analysis presented ($-0.5\leqslant {{\varepsilon}_{\text{dd}}}\leqslant 1$) this is not surprising since roton minima in dipolar BECs do not occur for $-0.5\leqslant {{\varepsilon}_{\text{dd}}}\leqslant 1$. To access regimes where a roton minimum occurs requires a treatment beyond the Thomas–Fermi approximation since the zero-point energy plays a crucial role in determining the properties of a BEC in the presence of strong dipolar interactions.

For a non-dipolar BEC in the Thomas–Fermi limit, only the rotation frequency and trap ellipticity determine the stability of the rotating frame solutions. As such adiabatic changes in and $\Omega$ can be utilised to trigger dynamical instability. For non-dipolar BECs this has been realised both experimentally [40, 183] and numerically [184–186], with excellent agreement with the hydrodynamic predictions.

For dipolar BECs the static solutions and their instability depend on the rotation frequency $\Omega$, the trap ellipticity , the trap ratio γ and the interaction parameter ${{\varepsilon}_{\text{dd}}}$, see sections 7.2 and 7.3. In principle all of these parameters can be adiabatically tuned in time. This provides many routes through parameter space to induce instability in the system.

Figure 18 [167] shows examples of these routes. Specifically, figure 18 shows the static solutions ${{\alpha}_{\text{v}}}$ of equation (77) as a function of $\Omega$ (figure 18(a)) (with $\epsilon =0.025$, $\gamma =1$ and ${{\varepsilon}_{\text{dd}}}=0.99$), (figure 18(b)) (with $\gamma =1$, $\Omega =0.7{{\omega}_{\bot}}$ and ${{\epsilon}_{\text{dd}}}=0.99$), ${{\epsilon}_{\text{dd}}}$ (figure 18(c)) (with $\epsilon =0.025$, $\gamma =1$ and $\Omega =0.7{{\omega}_{\bot}}$) and γ (figure 18(d)) (with $\epsilon =0.025$, $\Omega =0.7{{\omega}_{\bot}}$ and ${{\epsilon}_{\text{dd}}}=0.99$). In each figure the routes towards instability, due to an adiabatic change in the free parameters $\Omega$, , ${{\epsilon}_{\text{dd}}}$ or γ, are denoted by grey arrows, with the onset of instability highlighted by asterisks. Two types of instability can arise: (i) a dynamical instability (dashed grey arrows), or (ii) back-bending of the stationary solution such that it ceases to exist (solid grey arrows). For ${{\alpha}_{\text{v}}}>0$, the instability always arises from dynamical instabilities in the rotating frame stationary solution. For ${{\alpha}_{\text{v}}}<0$ the instability always arises from the rotating frame stationary solution ceasing to exist.

Zoom In Zoom Out Reset image size

7.4.1. Does the final state of the system contain vortices?

In section 7.4 we considered the routes to instability in a rotating dipolar BEC. This section addresses the question of whether such instabilities lead to the seeding of a vortex or vortex lattice in the BEC. We need to consider both the hydrodynamical surface instability of the rotating BEC and the energetic favourability for a vortex to be supported in the BEC.

For a non-dipolar BEC the admittance of a vortex into the BEC becomes energetically favourable when the rotation frequency exceeds the critical frequency ${{\Omega}_{\text{v}}}$ defined in equation (37), where in the Thomas–Fermi limit (assuming cylindrical symmetry) ${{\Omega}_{\text{v}}}$ is given by equation (47). For typical condensate parameters ${{\Omega}_{\text{v}}}\sim 0.3{{\omega}_{\bot}}$. This result is inconsistent with experimental observations of a much higher threshold for vortex lattice formation, $\Omega \sim 0.7{{\omega}_{\bot}}$. This discrepancy can be understood by considering the mechanism for vortices to be seeded into the BEC. For example, consider the case where $\Omega$ is increased adiabatically from zero (figure 18(a)). Then the pertinent stationary solution is in the ${{\alpha}_{\text{v}}}>0$ half-plane, which becomes dynamically unstable when $\Omega \approx 0.7{{\omega}_{\bot}}$. This instability provides a mechanism for vortices to be seeded into the BEC and ultimately relax into a vortex lattice. Alternatively, for $\Omega >{{\Omega}_{\text{v}}}$ the global energy minimum of the BEC is one which contains a vortex or vortex lattice, with the vortex-free solution being a local energy minimum. However, as $\Omega$ is adiabatically increased from zero, for $\Omega <{{\Omega}_{\text{v}}}$ the global energy minimum is defined by the rotating frame stationary solutions. As such as $\Omega$ passes through ${{\Omega}_{\text{v}}}$ there is no adiabatic path from the rotating frame stationary solution to the vortex solution, i.e. to move from the local energy minimum of the rotating frame stationary solution to the global energy minimum of the vortex state requires overcoming a significant energy barrier, due to the change in topology of the BEC between the two states. The dynamical instabilities discussed in section 7.4 provide a mechanism to overcome this topological barrier and hence seed the formation of vortices and vortex lattices at rotation frequencies $\Omega \approx 0.7{{\omega}_{\bot}}$.

Of course, the details of how a vortex or a vortex lattice forms once a hydrodynamical instability occurs is non-trivial, however the instability is the first step in the process [184–186]. For example, consider a BEC undergoing an adiabatic introduction of $\Omega$ (from zero), see dashed grey arrow in figure 18(a). At a critical rotation frequency, ${{\Omega}_{\text{i}}}$, the BEC becomes dynamically unstable. This leads to the exponential growth of surface ripples in the BEC [35, 183, 186]. Alternatively, consider a BEC undergoing an adiabatic introduction of (from zero), see solid grey arrow in figure 18(b). Above some critical trap ellipticity the stationary solutions no longer exist and the BEC undergoes large shape oscillations. In each case, and also for the adiabatic introduction of ${{\varepsilon}_{\text{dd}}}$ (figure 18(c)) or γ (figure 18(d)), these instabilities provide a mechanism for vortices to nucleate into the condensate. From this a turbulent vortex state emerges which then relaxes to a vortex lattice [186, 194].

The rotation frequency at which it becomes energetically favourable for a vortex to be admitted into a dipolar condensate depends crucially on the trap geometry γ and the strength of the dipolar interactions ${{\varepsilon}_{\text{dd}}}$ as discussed in section 5.3 based upon the results in [161]. There it was assumed that the system was cylindrically-symmetric, however, if we assume a very weak ellipticity $\epsilon =0.025$, it is expected that the correction to the critical frequency will be correspondingly small.

Below we consider two regimes for admittance of a vortex or vortex lattice into a dipolar BEC, as a function of rotation frequency and ${{\varepsilon}_{\text{dd}}}$. Initially an oblate trap is considered with $\gamma =10$, see figure 19(a) which plots the instability frequencies ${{\Omega}_{\text{i}}}$ (short dashed curve) and ${{\Omega}_{\text{b}}}$ (long dashed red curve) as a function of the dipolar interactions ${{\varepsilon}_{\text{dd}}}$. For adiabatic changes in $\Omega$ (vertical path) or ${{\varepsilon}_{\text{dd}}}$ (horizontal path), the system becomes unstable when it reaches one of the instability lines. The key feature to note is that the instability frequencies decrease weakly as the dipolar interactions are increased and have the approximate value ${{\Omega}_{\text{i}}}\approx {{\Omega}_{\text{b}}}\approx 0.75{{\omega}_{\bot}}$. Also shown in figure 19(a) is the rotation frequency ${{\Omega}_{\text{v}}}$ (solid curve) at which it becomes energetically favourable for a single vortex to reside in the dipolar BEC. Dipolar interactions, in this oblate system, lead to a weak decrease in ${{\Omega}_{\text{v}}}$, with ${{\Omega}_{\text{v}}}\approx 0.1{{\omega}_{\bot}}$ for the parameters considered [161]. These results show us that when the condensate becomes dynamically unstable a vortex state is already energetically favoured, i.e. the rotation frequency at which it is energetically preferable to have a vortex in the BEC is much lower than the rotation frequency required to induce an instability. As such, it is expected that in an oblate dipolar BEC a vortex lattice will ultimately form when these instabilities are reached.

Zoom In Zoom Out Reset image size

For a prolate trap with $\gamma =0.1$, figure 19(b) plots the instability frequencies ${{\Omega}_{\text{i}}}$ (short dashed curve) and ${{\Omega}_{\text{b}}}$ (long dashed red curve) and ${{\Omega}_{\text{v}}}$ (solid curve). The instability frequencies follow a similar trend to the oblate case. However, ${{\Omega}_{\text{v}}}$ is drastically different, increasing significantly with ${{\varepsilon}_{\text{dd}}}$ [161]. Depending on the strength of the dipolar interactions, two regimes are predicted. For ${{\varepsilon}_{\text{dd}}}\leqslant 0.8$ both ${{\Omega}_{\text{i}}}$ and ${{\Omega}_{\text{b}}}$ are larger than ${{\Omega}_{\text{v}}}$, and hence it is expected that an instability will lead to the admission of vortices into the condensate (and subsequently a vortex lattice). However, for ${{\varepsilon}_{\text{dd}}}\geqslant 0.8$, both ${{\Omega}_{\text{i}}}$ and ${{\Omega}_{\text{b}}}$ are lower than ${{\Omega}_{\text{v}}}$, implying that although an instability in the vortex free solution occurs, a vortex state is not energetically favourable. In this scenario the final state is not clear due to the net attractive dipolar interactions in this prolate configuration, and similar behaviour to non-dipolar BECs with attractive contact interactions (g  <  0) may arise. For non-dipolar BECs with attractive interactions the formation of vortices is also unfavourable with the possible final states including centre-of-mass motion, quadrupole oscillations and higher angular momentum-carrying shape excitations [195–197]. However, the true nature of the final state warrants further investigation.

Numerical results of the dipolar GPE in the quasi-two-dimensional regime [198, 199] show that, for an adiabatic introduction of the rotation frequency, the strength of the dipolar interaction influences the rotation frequency at which vortices are admitted into the condensate. In agreement with analysis presented above, it is also found that as the dipolar interaction is increased the rotation frequency required to nucleate vortices is reduced.

In principle the analysis presented in sections 7.2–7.4 can be modified to consider vortex generation through the application of a synthetic magnetic field and extended to investigate how off-axis alignment of the dipoles affects the stationary solutions and their dynamical stability.

7.5.1. Synthetic magnetic fields.

In a remarkable experiment in non-dipolar BECs, synthetic magnetic fields have been used to nucleate vortices [177]. The nucleation of vortices in such systems can be analysed via the general methods presented sections 7.2 and 7.3 [179], with favourable agreement being found with experimental results. Although the extension of this calculation to the dipolar case has not yet been pursued in the literature, the analysis presented in sections 7.2–7.4 for dipolar BECs can be modified to include synthetic magnetic fields. Let us outline the calculation: The starting point is,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{\partial n}{\partial t}=-\boldsymbol{\nabla} \cdot \left[n\boldsymbol{\nu}\right], \end{array}\end{eqnarray} \tag{ 82 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} m\frac{\partial \mathbf{v}}{\partial t}=-\boldsymbol{\nabla} \left(\frac{1}{2}m{{\nu}^{2}}+V+gn+ \Phi \right), \end{array}\end{eqnarray} \tag{ 83 }$$

where the generalized velocity is $\boldsymbol{\nu}=\mathbf{v}-{{q}^{\ast}}{{\mathbf{A}}^{\mathbf{\ast}}}/m$, with ${{\mathbf{A}}^{\mathbf{\ast}}}$ being the synthetic vector potential associated with the synthetic magnetic field. From this it is possible to derive an equivalent set of stationary solutions and determine their stability via,

$$\begin{eqnarray}&&\frac{\partial}{\partial t}\left[\begin{array}{c} \delta S \\ \delta n \end{array}\right]\,=-\left[\begin{array}{c} \boldsymbol{\nu}\cdot \boldsymbol{\nabla} & \frac{g}{m}\left(1+{{\varepsilon}_{\text{dd}}}K\right) \\ \boldsymbol{\nabla} \cdot {{n}_{(0)}}\boldsymbol{\nabla} & \boldsymbol{\nabla} \cdot \boldsymbol{\nu}+\boldsymbol{\nu}\cdot \boldsymbol{\nabla} \end{array}\right]\left[\begin{array}{c} \delta S \\ \delta n \end{array}\right].\end{eqnarray} \tag{ 84 }$$

Such an analysis would provide the basis to quantify under what regimes it is favourable for a vortex to be seeded into a dipolar BEC in the presence of a synthetic magnetic field.

7.5.2. Off-axis dipole orientation.

Assuming that a condensate is confined in a cylindrically-symmetric harmonic trap, the energy functional in the rotating frame is given by,

$$\begin{eqnarray}\begin{array}{*{35}{l}} E\left[\psi \right] & = & {\int}^{}\psi \left(\mathbf{r}\right){{H}^{\prime}}\psi \left(\mathbf{r}\right)\,\text{d}\mathbf{r} \\ {} & {} & +\,\,\frac{1}{2}{\int}^{}n\left({{\mathbf{r}}_{\mathbf{1}}}\right)U\left({{\mathbf{r}}_{\mathbf{1}}}-{{\mathbf{r}}_{\mathbf{2}}}\right)n\left({{\mathbf{r}}_{\mathbf{2}}}\right)\,\text{d}{{\mathbf{r}}_{\mathbf{1}}}\text{d}{{\mathbf{r}}_{\mathbf{2}}}, \end{array}\end{eqnarray} \tag{ 85 }$$

where,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{H}^{\prime}} & = & \frac{1}{2m}{{\left(\text{i}\hbar \nabla + \Omega m{{\mathbf{\hat{z}}}}\times \boldsymbol{\rho}\right)}^{2}} \\ {} & {} & +\,\,\frac{m}{2}\left(\omega _{\bot}^{2}-{{\Omega}^{2}}\right){{\rho}^{2}}+\frac{m}{2}\omega _{z}^{2}{{z}^{2}}. \end{array}\end{eqnarray} \tag{ 86 }$$

There are two distinct contributions to the total energy: a single-particle energy contribution, E₀ (the first term in equation (85)) and an interaction energy contribution, which in the absence of dipolar interactions is E_vdW (the second term in equation (85)).

To study the properties of vortex lattices it is appropriate to consider the quasi-two-dimensional regime. Physically, this corresponds to a situation where the condensate is rapidly rotating so that centrifugal spreading is significant, and where the longitudinal trapping is strong ($\hbar {{\omega}_{z}}\gg gn(0)$). In such circumstances it is appropriate to assume that the longitudinal motion is described by the ground state of the z-confinement, so that the condensate wavefunction can be factorized according to equation (31). The resulting form for the energy functional, in the absence of dipolar interactions is,

$$\begin{eqnarray}\begin{array}{*{35}{l}} E & = & \frac{N}{2}\hbar {{\omega}_{z}}+{\int}^{}\psi _{\bot}^{\star}\left(\boldsymbol{\rho}\right)H_{\bot}^{\prime}\psi _{\bot}^{\star}\left(\boldsymbol{\rho}\right)\,\text{d}\boldsymbol{\rho} \\ {} & {} & +\,\,\frac{g}{2\sqrt{2\pi}{{l}_{z}}}{\int}^{}n_{\bot}^{2}\left(\boldsymbol{\rho}\right)\,\text{d}\boldsymbol{\rho}, \end{array}\end{eqnarray} \tag{ 87 }$$

where,

$$\begin{eqnarray}\begin{array}{*{35}{l}} H_{\bot}^{\prime} & = & -\frac{1}{2m}{{\left(\text{i}\hbar {{\nabla}_{\bot}}+ \Omega m{{\mathbf{\hat{z}}}}\times \boldsymbol{\rho}\right)}^{2}} \\ {} & {} & +\,\,\frac{m}{2}\left(\omega _{\bot}^{2}-{{\Omega}^{2}}\right){{\rho}^{2}}. \end{array}\end{eqnarray} \tag{ 88 }$$

Consider the fast-rotating limit where $\Omega \to {{\omega}_{\bot}}$. This is the point at which the centrifugal spreading due to rotation almost overwhelms the confinement due to the radial trap. In this limit, the single-particle Hamiltonian in quasi-two dimensions tends to $-{{\left(\text{i}\hbar {{\nabla}_{\bot}}+ \Omega m{{\mathbf{\hat{z}}}}\times \boldsymbol{\rho}\right)}^{2}}/2m$, up to a constant. The eigenfunctions of this Hamiltonian are well known: they are the Landau level orbitals ${{u}_{{{m}^{\prime}},{{n}^{\prime}}}}(x,y)$ with corresponding eigenenergies ${{\epsilon}_{{{n}^{\prime}}}}=\hbar \Omega \left({{n}^{\prime}}+1/2\right)$. The n quantum number labels the Landau level, and may take on any non-negative integer value. Each Landau level is infinitely degenerate since ${{\epsilon}_{{{n}^{\prime}}}}$ does not depend on ${{m}^{\prime}}$ (which also takes on non-negative integer values).

Although the Landau level orbitals do not necessarily form a complete basis for the full quasi-two-dimensional Hamiltonian, in the limit of weak interactions they are a good approximate basis choice [195]. In searching for the ground state of the condensate it is assumed that it is adequately described by a superposition of ${{n}^{\prime}}=0$ Landau level orbitals, i.e. the lowest Landau level (LLL) approximation. Using an unrestricted minimization this assumption implies that [52],

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \left(1-\frac{\Omega}{{{\omega}_{\bot}}}\right)\lesssim \frac{\sqrt{2\pi}{{l}_{z}}}{2N{{a}_{s}}}. \end{array}\end{eqnarray} \tag{ 89 }$$

In this limit the two-dimensional ground state wavefunction of the condensate can be written as,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\psi}_{\bot}}\left(\boldsymbol{\rho}\right) & = & \underset{{{m}^{\prime}}=0}{\overset{\infty}{\sum}}{{c}_{{{m}^{\prime}}}}{{u}_{{{m}^{\prime}},0}}\left(\boldsymbol{\rho}\right) \\ {} & = & \underset{{{m}^{\prime}}=0}{\overset{\infty}{\sum}}\frac{{{c}_{{{m}^{\prime}}}}}{\sqrt{2\pi {{m}^{\prime}}!}}{{\left(\frac{x+\text{i}y}{{{l}_{\bot}}}\right)}^{{{m}^{\prime}}}}\exp \left(-\frac{{{x}^{2}}+{{y}^{2}}}{2l_{\bot}^{2}}\right), \end{array}\end{eqnarray} \tag{ 90 }$$

where the right-hand side follows from the explicit form of ${{u}_{{{m}^{\prime}},0}}(x,y)$. The length scale ${{l}_{\bot}}=\sqrt{\hbar /m \Omega}$ characterises the radial extent of the condensate and is effectively equal to the transverse trap length since $\Omega \to {{\omega}_{\bot}}$.

At this point it turns out to be convenient to convert to a complex number representation, rather than working with the components of a two-dimensional vector. To this end $\boldsymbol{\rho}=x{{\mathbf{\hat{x}}}}+y{{\mathbf{\hat{y}}}}$ is mapped onto the complex number $w=x+\text{i}y$. Then the ground state wavefunction of equation (90) may be written as,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\psi}_{\bot}}(w)=h(w)\exp \left(-\frac{|w{{|}^{2}}}{2l_{\bot}^{2}}\right), \end{array}\end{eqnarray} \tag{ 91 }$$

where h is an analytic function of w. With this definition, the coefficients of the superposition ${{c}_{{{m}^{\prime}}}}$ have been absorbed into h. It is possible to fully specify h in terms of its roots since it is an analytic function, where each root specifies the location of a vortex core at the corresponding x and y coordinates in the condensate [208].

A vortex lattice ground state corresponds to a situation where the roots of h lie on a lattice. The next step is to construct an analytic function with roots that satisfy this property. Fortunately, there is a well-studied function in complex analysis which will be useful here: the Jacobi theta function ${{\theta}_{1}}\left(z,\zeta \right)$. The roots of ${{\theta}_{1}}\left(z,\zeta \right)$ are able to describe any regular lattice up to a rotation by making an appropriate choice for the parameter $\zeta$. This means that h may be expressed in terms of ${{\theta}_{1}}$ to obtain a vortex lattice ground state wavefunction. However, in order to make this connection we must introduce a way to describe a regular lattice mathematically.

A two-dimensional lattice may be fully specified by a pair of basis vectors: ${{\mathbf{b}}_{\mathbf{1}}}$ and ${{\mathbf{b}}_{\mathbf{2}}}$. The points of the lattice are obtained from these basis vectors by constructing all possible linear combinations of the form ${{m}_{1}}{{\mathbf{b}}_{\mathbf{1}}}+{{m}_{2}}{{\mathbf{b}}_{\mathbf{2}}}$ where m₁ and m₂ are integers, i.e. a two-dimensional Bravais lattice. In general, four real parameters are required to fully specify a lattice: two real components for each of the two basis vectors. However, one of these parameters may be fixed since there is no need to distinguish between lattices which are equivalent up to a rotation. This is justified because the Hamiltonian which describes the condensate is cylindrically symmetric. In order to fix one of the parameters the first basis vector ${{\mathbf{b}}_{\mathbf{1}}}$ is chosen to be orientated along the x-axis.

The three remaining parameters which describe the lattice are depicted in figure 20, along with the lattice basis vectors. The first basis vector ${{\mathbf{b}}_{\mathbf{1}}}$ is specified solely in terms of its magnitude, which is denoted by the parameter b₁, since its direction is fixed along ${{\mathbf{\hat{x}}}}$. The second basis vector ${{\mathbf{b}}_{\mathbf{2}}}$ is defined with reference to the first basis vector through a rotation and rescaling. This requires two additional parameters: a rotation angle η and a scaling factor τ. Writing out the basis vectors explicitly in terms of the parameters b₁, τ and η, results in ${{\mathbf{b}}_{\mathbf{1}}}={{b}_{1}}{{\mathbf{\hat{x}}}}$ and ${{\mathbf{b}}_{\mathbf{2}}}=\tau {{b}_{1}}\left(\cos \eta {{\mathbf{\hat{x}}}}+\sin \eta {{\mathbf{\hat{y}}}}\right)$.

Zoom In Zoom Out Reset image size

The other parameter that appears in figure 20 is the area of the unit cell, given by ${{v}_{c}}=b_{1}^{2}\tau \sin \eta$. This parameter is redundant in specifying the lattice because it depends on the other parameters: b₁, τ and η. However, it is useful to mention it because it appears in a number of places in the subsequent analysis.

Another useful concept regarding our mathematical description of the vortex lattice is the reciprocal lattice. It is needed to represent a function which is defined on a lattice as a Fourier series. The reciprocal lattice is obtained from the original lattice by a transformation of the lattice basis vectors. Denoting the basis vectors of the reciprocal lattice as ${{\mathbf{q}}_{\mathbf{1}}}$ and ${{\mathbf{q}}_{\mathbf{2}}}$, we have ${{\mathbf{q}}_{\mathbf{1}}}=2\pi {{\mathbf{b}}_{\mathbf{2}}}\times {{\mathbf{\hat{z}}}}/{{v}_{c}}=2\pi {{b}_{1}}\tau \left(\sin \eta {{\mathbf{\hat{x}}}}-\cos \eta {{\mathbf{\hat{y}}}}\right)/{{v}_{c}}$ and ${{\mathbf{q}}_{\mathbf{2}}}=2\pi {{\mathbf{\hat{z}}}}\times {{\mathbf{b}}_{\mathbf{1}}}/{{v}_{c}}=2\pi {{b}_{1}}{{\mathbf{\hat{y}}}}/{{v}_{c}}$. As for the case of the original lattice, the reciprocal lattice is constructed from its basis vectors by considering all linear combinations of the form ${{m}_{1}}{{\mathbf{q}}_{\mathbf{1}}}+{{m}_{2}}{{\mathbf{q}}_{\mathbf{2}}}$ where m₁ and m₂ are integers. Each one of these linear combinations is a reciprocal lattice vector, which is denoted in general by $\mathbf{q}$.

From this, it is possible to describe the two-dimensional density as,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{n}_{\bot}}\left(\boldsymbol{\rho}\right)=\frac{N}{\pi {{\chi}^{2}}}{{\text{e}}^{-\frac{{{\rho}^{2}}}{{{\chi}^{2}}}}}g\left(\boldsymbol{\rho}\right), \end{array}\end{eqnarray} \tag{ 92 }$$

where $g\left(\boldsymbol{\rho}\right)=\sum\nolimits_{\mathbf{q}}{{\tilde{g}}_{\mathbf{q}}}\exp \left(i\mathbf{q}\cdot \boldsymbol{\rho}\right)$ and,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{g}_{\mathbf{q}}}=\frac{{{(-1)}^{{{m}_{1}}+{{m}_{2}}}}{{\text{e}}^{-{{v}_{c}}{{q}^{2}}/8\pi}}}{\sqrt{\tau \sin \eta}}, \end{array}\end{eqnarray} \tag{ 93 }$$

with $\mathbf{q}={{m}_{1}}{{\mathbf{q}}_{\mathbf{1}}}+{{m}_{2}}{{\mathbf{q}}_{\mathbf{2}}}$. In other words, the two-dimensional condensate density is the product of a Gaussian envelope and a function $g\left(\boldsymbol{\rho}\right)=\sum\nolimits_{\mathbf{q}}{{\tilde{g}}_{\mathbf{q}}}\exp \left(i\mathbf{q}\cdot \boldsymbol{\rho}\right)$ which is periodic on the vortex lattice. The Fourier coefficients of $g\left(\boldsymbol{\rho}\right)$ are rescaled compared to those in equation (93) so that ${{\tilde{g}}_{\mathbf{q}}}={{g}_{\mathbf{q}}}/\sum\nolimits_{\mathbf{v}}{{\tilde{g}}_{\mathbf{v}}}\exp \left(-i{{v}^{2}}{{\chi}^{2}}/4\right)$. The radial extent of the condensate cloud is quantified by the length-scale $\chi ={{\left(l_{\bot}^{-2}-\pi v_{c}^{-1}\right)}^{-1/2}}$ which is related to the number of vortices in the system.

From this ansatz for the vortex lattice ground state, it is possible to calculate the energy of the condensate as a function of the lattice parameters. Using the fact that ${{\psi}_{\bot}}\left(\boldsymbol{\rho}\right)$ is in the LLL, the single-particle energy contribution, in the limit of large vortex number, is independent of the vortex lattice parameters. The interaction energy contribution may be rewritten in the following form,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{E}_{\text{vdW}}}=\frac{{{N}^{2}}g}{{{\left(2\pi \right)}^{3/2}}{{l}_{z}}{{\chi}^{2}}}\mathcal{I}[n]. \end{array}\end{eqnarray} \tag{ 94 }$$

Here $\mathcal{I}[n]={{\pi}^{2}}{{\chi}^{2}}/{{N}^{2}}{\int}^{}n_{\bot}^{2}\left(\boldsymbol{\rho}\right)\text{d}\boldsymbol{\rho}$ is a dimensionless analogue of the interaction energy contribution. Substituting the ansatz for the vortex lattice ground state from equation (92) gives,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathcal{I}\left(\tau,\eta \right)=\frac{1}{2}\underset{\mathbf{q},\mathbf{v}}{\sum}{{{\tilde{g}}}_{\mathbf{q}}}{{{\tilde{g}}}_{\mathbf{v}}}{{\text{e}}^{-\frac{{{\chi}^{2}}|\mathbf{q}+\mathbf{v}{{|}^{2}}}{8}}}. \end{array}\end{eqnarray} \tag{ 95 }$$

This quantity depends on τ and η implicitly through its dependence on the reciprocal lattice vectors. Considering the limit of large vortex number where ${{\chi}^{2}}{{\mathbf{q}}^{2}}\gg 1$, it is possible to simplify this expression significantly. In fact, one may assume that the exponential ${{\text{e}}^{-\frac{{{\chi}^{2}}|\mathbf{q}+\mathbf{v}{{|}^{2}}}{8}}}$ is so sharply peaked at $\mathbf{q}=-\mathbf{v}$ that it may be approximated by a Kronecker delta, ${{\delta}_{\mathbf{q},-\mathbf{v}}}$. This collapses the double sum to a single sum over $\mathbf{q}$, which is much simpler to evaluate,

$$\begin{eqnarray}\begin{array}{*{35}{l}} \mathcal{I}\left(\tau,\eta \right) & \approx & \frac{1}{2}\underset{\mathbf{q}}{\sum}{{\left({{{\tilde{g}}}_{\mathbf{q}}}\right)}^{2}} \\ {} & = & \frac{1}{2}\underset{{{m}_{1}},{{m}_{2}}=-\infty}{\overset{\infty}{\sum}}{{\text{e}}^{\pi \left(2{{m}_{1}}{{m}_{2}}\cot \eta -m_{1}^{2}\tau \csc \eta -\frac{m_{2}^{2}\csc \eta}{\tau}\right)}}. \end{array}\end{eqnarray} \tag{ 96 }$$

In order to minimise $\mathcal{I}\left(\tau,\eta \right)$ as defined in equation (96) a numerical treatment is required, which means cutting off the sums over m₁ and m₂ at some upper and lower bounds. Defining the upper and lower bounds to be at M and  −M respectively, with M set to 15 ensures that $\mathcal{I}\left(\tau,\eta \right)$ is accurate to double precision for values of τ and η greater than about 0.05. To explore regions in which τ or η is less than 0.05, M needs to be increased above 15 to include higher frequency terms in the Fourier series. Fortunately, it is reasonable to exclude these regions, in the absence of dipolar interactions, because they correspond to an unphysical situation where the vortex lattice begins to collapse onto a line.

It is illustrative to perform the minimisation of $\mathcal{I}\left(\tau,\eta \right)$ graphically, by generating a contour plot of $\mathcal{I}\left(\tau,\eta \right)$ as a function of τ and η. The result is shown in figure 21. In this plot, the dark purple shading corresponds to regions where $\mathcal{I}\left(\tau,\eta \right)$ approaches its minimum value. The white areas correspond to regions where $\mathcal{I}\left(\tau,\eta \right)$ is greater than the cut-off value, which is in this case set to 0.6. In figure 21 each local minimum has the same value of $\mathcal{I}$ equal to 0.5797. However, a careful consideration shows that each of these solutions in fact corresponds to the same type of lattice; just at a different scale and orientation. The solution labelled '1' in figure 21 is the standard parametrisation of the triangular lattice. It is illustrated in figure 22, along with the standard parametrisations of the three other types of lattice: square, rectangular and parallelogrammic. In general, the standard parametrisation is defined to be the one for which τ is closest to 1. The alternative parametrisations of the same lattice have smaller values of τ compared to the standard one. For example, the second solution, denoted as '2' in figure 21 is an equally valid parametrisation of the triangular lattice.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The above calculation verifies that a triangular vortex lattice geometry is always favoured in non-dipolar BECs. Of course, this was to be expected based on the results of numerous experiments and previous theoretical studies. It is interesting to note that the triangular lattice geometry is not significantly favoured over other possible lattice geometries. In particular, the energy corresponding to the square lattice geometry at $\tau =1$ and $\eta =\pi /2$ (the saddle point in figure 21) is only 1.8% larger than that of the triangular lattice geometry. It is therefore conceivable that the energy minimum may shift to a non-triangular lattice geometry if the functional form of the interaction energy contribution is altered. With this motivation the above calculation is generalized below to include dipolar interactions.

For the case of dipolar BECs the dipolar interaction potential, ${{U}_{\text{dd}}}\left(\mathbf{r}\right)$, needs to be included. This is the only modification to the theory that is required to account for the effect of the dipoles. Assuming that the dipoles are aligned perpendicular to the plane of rotation, i.e. along the z-axis, implies that the cylindrical symmetry of the Hamiltonian is maintained. In this limit the criterion for being in the LLL becomes

$$\begin{eqnarray}\begin{array}{*{35}{l}} \left(1-\frac{\Omega}{{{\omega}_{\bot}}}\right) & \lesssim & \frac{\sqrt{2\pi}{{l}_{z}}}{2N{{a}_{s}}\left(1+\frac{16\pi {{\varepsilon}_{\text{dd}}}}{3}\left[1-\overline{\sigma}\sqrt{\frac{9\pi}{8}}\right]\right)} \\ {} & + & \frac{{{\pi}^{2}}{{l}_{z}}{{\varepsilon}_{\text{dd}}}\overline{\sigma}}{{{a}_{s}}N{{\left(1+\frac{16\pi {{\varepsilon}_{\text{dd}}}}{3}\left[1-\overline{\sigma}\sqrt{\frac{9\pi}{8}}\right]\right)}^{2}}}, \end{array}\end{eqnarray} \tag{ 97 }$$

where $\overline{\sigma}={{l}_{z}}/{{R}_{\bot}}$. The above has been obtained by considering the expansion of the $\tilde{U}_{\text{dd}}^{\bot}\left(\mathbf{q}\right)$ in terms of the width of the BEC to first order, i.e. equation (50), and then calculating the dipolar potential to second order in $\rho /{{R}_{\bot}}$.

The additional contribution to the energy functional arises due to the dipolar interactions,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{E}_{\text{dd}}}=\frac{1}{2}{\int}^{}n\left({{\mathbf{r}}_{\mathbf{1}}}\right){{U}_{\text{dd}}}\left({{\mathbf{r}}_{\mathbf{1}}}-{{\mathbf{r}}_{\mathbf{2}}}\right)n\left({{\mathbf{r}}_{\mathbf{2}}}\right)\,\text{d}{{\mathbf{r}}_{\mathbf{1}}}\text{d}{{\mathbf{r}}_{\mathbf{2}}}. \end{array}\end{eqnarray} \tag{ 98 }$$

Adding this to the single-particle and contact interaction energy contributions gives the total energy in the rotating frame for a dipolar BEC,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} E={{E}_{0}}+{{E}_{\text{vdW}}}+{{E}_{\text{dd}}}. \end{array}\end{eqnarray} \tag{ 99 }$$

The results obtained in section 8.1 for the single-particle (E₀) and contact interaction (E_vdW) energy contributions are unchanged for the dipolar case. All that remains then, is to perform the calculation for ${{E}_{\text{dd}}}$.

Rewriting the dipolar interaction energy contribution in reciprocal space by applying the Fourier convolution theorem leads to,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{E}_{\text{dd}}}=\frac{1}{2}\frac{1}{{{\left(2\pi \right)}^{3}}}{\int}^{}\tilde{n}\left(\mathbf{k}\right)\tilde{n}\left(-\mathbf{k}\right){{{\tilde{U}}}_{\text{dd}}}\left(\mathbf{k}\right)\,\text{d}\mathbf{k}. \end{array}\end{eqnarray} \tag{ 100 }$$

To evaluate this in quasi-two dimensions, we substitute the Fourier transform of the quasi-two-dimensional condensate density, $\tilde{n}\left(\mathbf{k}\right)={{\text{e}}^{-k_{z}^{2}l_{z}^{2}/4}}{{\tilde{n}}_{\bot}}\left(\mathbf{q}\right)$, such that

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{E}_{\text{dd}}}=\frac{1}{2}\frac{1}{{{\left(2\pi \right)}^{2}}}{\int}^{}{{{\tilde{n}}}_{\bot}}\left(\mathbf{q}\right){{{\tilde{n}}}_{\bot}}\left(-\mathbf{q}\right)\tilde{U}_{\text{dd}}^{\bot}\left(\mathbf{q}\right)\,\text{d}\mathbf{q}, \end{array}\end{eqnarray} \tag{ 101 }$$

where $\tilde{U}_{\text{dd}}^{\bot}\left(\mathbf{q}\right)$ is given by equation (34), with $\alpha =0$.

Now that an expression for the dipolar interaction energy in quasi-two dimensions has been derived, it is possible to evaluate the energy assuming that the condensate is in the vortex lattice ground state. Performing the Fourier transform on the vortex lattice condensate density specified in equation (92) results in,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \tilde{n}\left(\mathbf{q}\right)=N\underset{{{\mathbf{q}}^{\prime}}}{\sum}{{{\tilde{g}}}_{{{\mathbf{q}}^{\prime}}}}{{\text{e}}^{-{{\chi}^{2}}|{{\mathbf{q}}^{\prime}}-\mathbf{q}{{|}^{2}}/4}}. \end{array}\end{eqnarray} \tag{ 102 }$$

This enters into the expression for the dipolar interaction energy, leading to,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{E}_{\text{dd}}}\left(\tau,\eta \right)=\frac{{{N}^{2}}{{C}_{\text{dd}}}}{3{{\left(2\pi \right)}^{\frac{3}{2}}}l_{z}^{3}}\underset{\mathbf{q},\mathbf{v}}{\sum}{{{\tilde{g}}}_{\mathbf{q}}}{{{\tilde{g}}}_{\mathbf{v}}}{{\text{e}}^{-{{\chi}^{2}}\left({{q}^{2}}+{{v}^{2}}\right)/4}}A\left(\mathbf{v}-\mathbf{q}\right), \end{array}\end{eqnarray} \tag{ 103 }$$

where,

$$\begin{eqnarray}\begin{array}{*{35}{l}} A\left(\mathbf{q}\right) & = & \frac{1}{2\pi}{\int}^{}{{\text{e}}^{-\frac{{{\chi}^{2}}}{l_{z}^{2}}\left({{u}^{2}}+\frac{{{l}_{z}}\mathbf{q}\cdot \mathbf{u}}{\sqrt{2}}\right)}}\left[2-3\sqrt{\pi}u{{e}^{u}}\text{erfc}(u)\right]\,\text{d}\mathbf{u} \\ {} & = & 2{\int}_{0}^{\infty}{{\text{e}}^{-\frac{{{\chi}^{2}}{{u}^{2}}}{l_{z}^{2}}}}u{{I}_{0}}\left(\frac{{{\chi}^{2}}qu}{\sqrt{2}l_{z}^{2}}\right)\,\text{d}u \\ {} & - & 3\sqrt{\pi}{\int}_{0}^{\infty}{{\text{e}}^{-\frac{{{\chi}^{2}}{{u}^{2}}}{l_{z}^{2}}}}{{u}^{2}}{{e}^{u}}\text{erfc}(u){{I}_{0}}\left(\frac{{{\chi}^{2}}qu}{\sqrt{2}l_{z}^{2}}\right)\,\text{d}u, \end{array}\end{eqnarray} \tag{ 104 }$$

and ${{I}_{0}}\left(\cdot \right)$ is the modified Bessel function of the first kind. At this point, the integral may be separated into two terms ${{A}_{1}}\left(\mathbf{q}\right)+{{A}_{2}}\left(\mathbf{q}\right)$. The first term has a simple solution,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{A}_{1}}\left(\mathbf{q}\right)=2{\int}_{0}^{\infty}u{{\text{e}}^{-\frac{{{\chi}^{2}}{{u}^{2}}}{l_{z}^{2}}}}{{I}_{0}}\left(\frac{{{\chi}^{2}}qu}{\sqrt{2}l_{z}^{2}}\right)\,\text{d}u={{\left(\frac{{{l}_{z}}}{\chi}\right)}^{2}}{{\text{e}}^{\frac{{{\chi}^{2}}{{q}^{2}}}{8}}}. \end{array}\end{eqnarray} \tag{ 105 }$$

Returning to the expression for the dipolar interaction energy given in equation (103), and substituting the simplified result for ${{A}_{1}}\left(\mathbf{q}\right)$ gives

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{E}_{\text{dd}}}\left(\tau,\eta \right)=\frac{2{{N}^{2}}{{C}_{\text{dd}}}}{3{{\left(2\pi \right)}^{\frac{3}{2}}}{{l}_{z}}{{\chi}^{2}}}\mathcal{I}\left(\tau,\eta \right)+\frac{{{N}^{2}}{{C}_{\text{dd}}}}{3{{\left(2\pi \right)}^{\frac{3}{2}}}l_{z}^{3}}\mathcal{W}\left(\tau,\eta \right), \end{array}\end{eqnarray} \tag{ 106 }$$

where $\mathcal{I}\left(\tau,\eta \right)$ is as defined in equation (95) and,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathcal{W}\left(\tau,\eta \right)=\underset{\mathbf{q},\mathbf{v}}{\sum}{{{\tilde{g}}}_{\mathbf{q}}}{{{\tilde{g}}}_{\mathbf{v}}}{{\text{e}}^{-{{\chi}^{2}}{{\left({{\mathbf{q}}^{2}}+\mathbf{v}\right)}^{2}}/4}}{{A}_{2}}\left(\mathbf{v}-\mathbf{q}\right). \end{array}\end{eqnarray} \tag{ 107 }$$

This expression shows that the dipolar energy has been separated into two distinct contributions: a local contribution which is proportional to $\mathcal{I}\left(\tau,\eta \right)$ and a non-local contribution which is proportional to $\mathcal{W}\left(\tau,\eta \right)$. In principle, it is now possible to calculate the dipolar interaction energy as a function of the lattice parameters, however the function for the non-local contribution, $\mathcal{W}\left(\tau,\eta \right)$ is not analytically tractable and it must be evaluated numerically. In order to resolve this the difficult integral is expanded in a series of simpler integrals, each of which can be performed analytically. This can be done by expressing the complementary error function, which appears in the integrand, as a power series,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{erfc}(u)=1-\frac{1}{\sqrt{\pi}}\underset{n=0}{\overset{\infty}{\sum}}\frac{{{(-1)}^{n}}{{u}^{2n+1}}}{n!(2n+1)} \end{array}\end{eqnarray} \tag{ 108 }$$

and bringing the sum outside of the integral. Additionally, as in the case of the contact interactions it is possible to reduce the double sum, in equation (107) to a single sum in the limit of large vortex number, resulting in,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathcal{W}\left(\tau,\eta \right)\approx \underset{\mathbf{q}}{\sum}{{\left({{{\tilde{g}}}_{\mathbf{q}}}\right)}^{2}}{{\text{e}}^{-\frac{{{\chi}^{2}}{{q}^{2}}}{2}}}\left(A_{2}^{\text{a}}\left(2\mathbf{q}\right)+A_{2}^{\text{b}}\left(2\mathbf{q}\right)\right), \end{array}\end{eqnarray} \tag{ 109 }$$

where,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\text{e}}^{-\frac{{{\chi}^{2}}{{q}^{2}}}{2}}}A_{2}^{\text{a}}(2\mathbf{q}) & = & \frac{3\pi}{8{{({{\beta}^{2}}-1)}^{\frac{5}{2}}}}{{\text{e}}^{-\frac{({{\beta}^{2}}-2){{q}^{2}}{{\chi}^{2}}}{4({{\beta}^{2}}-1)}}} \\ {} & {} & \times \,\,\left[{{\beta}^{4}}l_{z}^{2}{{q}^{2}}{{I}_{1}}\left(\frac{{{\beta}^{4}}l_{z}^{2}{{q}^{2}}}{4-4{{\beta}^{2}}}\right) \right. \\ {} & {} & \left. -\,\,\left(2{{\beta}^{2}}+{{\beta}^{4}}l_{z}^{2}{{q}^{2}}-2\right){{I}_{0}}\left(\frac{{{\beta}^{4}}l_{z}^{2}{{q}^{2}}}{4-4{{\beta}^{2}}}\right)\right] \end{array}\end{eqnarray} \tag{ 110 }$$

and,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & {{\text{e}}^{-\frac{{{\chi}^{2}}{{q}^{2}}}{2}}}A_{2}^{\text{b}}\left(2\mathbf{q}\right)=3{{\text{e}}^{-\frac{{{q}^{2}}{{\chi}^{2}}}{2({{\beta}^{2}}-1)}}} \\ {} & ~ & \times \underset{n=0}{\overset{\infty}{\sum}}\frac{{{(-1)}^{n}}(n+1)}{{{({{\beta}^{2}}-1)}^{2+n}}(2n+1)}{{L}_{n+1}}\left(\frac{{{\beta}^{4}}l_{z}^{2}{{q}^{2}}}{2-2{{\beta}^{2}}}\right), \end{array}\end{eqnarray} \tag{ 111 }$$

where $\beta =\chi /{{l}_{z}}$ and ${{L}_{n}}\left(\cdot \right)$ is the nth Laguerre polynomial.

With the results of section 8.1 it is now feasible, from a computational perspective, to numerically minimise the condensate energy, $E\left(\tau,\eta \right)={{E}_{0}}+{{E}_{\text{int}}}\left(\tau,\eta \right)$, with respect to τ and η to determine the optimal vortex lattice geometry. In minimising the condensate energy, as in section 8.1, only the interaction energy contribution needs to be considered,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{E}_{\text{int}}}\left(\tau,\eta \right)=\frac{{{N}^{2}}{{C}_{\text{dd}}}}{3{{\left(2\pi \right)}^{\frac{3}{2}}}l_{z}^{3}} \\ {} & \times \,\left[\left(2+\frac{1}{{{\varepsilon}_{\text{dd}}}}\right){{\left(\frac{{{l}_{z}}}{\chi}\right)}^{2}}\mathcal{I}\left(\tau,\eta \right)+\mathcal{W}\left(\tau,\eta \right)\right], \end{array}\end{eqnarray} \tag{ 112 }$$

since the single-particle contribution does not depend on τ and η. Assuming ${{C}_{\text{dd}}}>0$ so that the factor outside the square brackets is positive the vortex lattice configuration is determined by the the minimization of

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \left(2+\frac{1}{{{\varepsilon}_{\text{dd}}}}\right){{\left(\frac{{{l}_{z}}}{\chi}\right)}^{2}}\mathcal{I}\left(\tau,\eta \right)+\mathcal{W}\left(\tau,\eta \right), \end{array}\end{eqnarray} \tag{ 113 }$$

with respect to τ and η. The results of this minimization are shown in figure 23 for particular choices of length scale parameters which are given in the figure caption. In both figures, the optimal values of τ and η are plotted against $\varepsilon _{\text{dd}}^{-1}$. The optimal value of τ is represented by a solid line with reference to the scale on the left vertical axis, while the optimal value of η is represented by a dotted line with reference to the scale on the right vertical axis. Taken together, the optimal values of τ and η describe the optimal vortex lattice geometry.

Zoom In Zoom Out Reset image size

By comparing the optimal values of τ and η to the standard lattice parameterisations given in figure 22, it is possible to classify the geometry of the vortex lattice. For example, at the point $\varepsilon _{\text{dd}}^{-1}=-1.82$ in figure 23(a), a square lattice is favourable since the optimal values of τ and η at that point are 1 and $\pi /2$ respectively. Continuing in this way, three types of lattice emerge depending on the value of $\varepsilon _{\text{dd}}^{-1}$: triangular, square or rectangular. Each of these lattice geometries occurs in a distinct region of the phase space, indicated by distinct shading in figure 23. The fact that the transition to different lattice geometries occurs when $\varepsilon _{\text{dd}}^{-1}\approx -2$ indicates that the transition arises in the regime where the long-range contribution to the interactions $\left[\mathcal{W}\left(\tau,\eta \right)\right]$ dominates, see equation (113). Comparing figures 23(a), with ${{l}_{\bot}}/{{l}_{z}}=40$, and (b), with ${{l}_{\bot}}/{{l}_{z}}=80$, the relative size of each region is comparable for both condensate aspect ratios. There is however an overall translation and scaling difference between the two cases: the regions in figure 23(b), where the condensate is more oblate, are contracted and shifted to the left compared to those in figure 23(a), where the condensate is less oblate. A more accurate description of the phase regions is given in table 1 in terms of inequalities.

Table 1. Definition of the four regions in phase space shown in figure 23.

Phase	Figure 23(a)	Figure 23(b)
Collapse	$\varepsilon _{\text{dd}}^{-1}<-1.895$	$\varepsilon _{\text{dd}}^{-1}<-1.946$
Rectangular	$-1.895<\varepsilon _{\text{dd}}^{-1}<-1.825$	$-1.946<\varepsilon _{\text{dd}}^{-1}<-1.911$
Square	$-1.825<\varepsilon _{\text{dd}}^{-1}<-1.815$	$-1.911<\varepsilon _{\text{dd}}^{-1}<-1.906$
Triangular	$\varepsilon _{\text{dd}}^{-1}>-1.815$	$\varepsilon _{\text{dd}}^{-1}>-1.906$

In addition to the three pattern-shaded regions, there is also a solid grey-shaded region for which the condensate is in the so-called collapse phase. In the collapse phase, the optimal value of τ tends to zero and the vortex lattice analysis begins to break down. Physically, this phase corresponds to a situation where the vortices are arranged in densely-packed lines, with the spacing between the lines being larger than the extent of the condensate in the x  −  y plane. Since the unit cell of the lattice extends beyond the boundaries of the condensate in this situation, the vortex lattice can be considered to have collapsed. When the system enters the collapse phase, there are also signs that the analysis becomes invalid. Since τ approaches zero in this phase, the expression for interaction energy becomes inaccurate because only enough terms in the Fourier decomposition to consider values of τ greater than about 0.05 are included (as in section 8.1 M  =  15).

Apart from the stability of the vortex lattice, we may also assess the stability of the condensate itself by looking at the sign of the interaction energy. If the interaction energy is negative, then it can approximately be regarded that the condensate as being prone to collapse. The region of phase space for which interaction energy is negative is the area to the left of the red line in figures 23(a) and (b). Interestingly, the interaction energy becomes negative at roughly the value of $\varepsilon _{\text{dd}}^{-1}$ where the optimal value of τ approaches zero. This suggests that there may be a link between the collapse of the condensate and the collapse of the vortex lattice.

In this section the vortex lattice geometry in dipolar BECs for the special case of on-axis polarisation has been analysed. The results show that three lattice geometries are possible, depending on the value of $\varepsilon _{\text{dd}}^{-1}$ and the values of the length scales ${{l}_{\bot}}$, l_z and χ. In general, a triangular lattice geometry is favoured in regions where the local interaction contribution dominates, as was seen in the non-dipolar case. However, in regions where the non-local interaction contribution becomes significant, the favoured lattice geometry changes from triangular to square or rectangular. Below a certain value of $\varepsilon _{\text{dd}}^{-1}$ (corresponding to reasonably strong, attractive contact interactions) the vortex lattice and the condensate appear to collapse concurrently.

The results obtained above qualitatively agree with previous results of Zhang and Zhai [202] and Cooper et al [201, 203]. Zhang and Zhai also find that the lattice geometry undergoes a transition from triangular  →  square  →  rectangular  →  collapse as the value of $\varepsilon _{\text{dd}}^{-1}$ decreases. Cooper et al also find the same transitions between lattice geometries. However, they do not find that the lattice collapses after passing through the rectangular lattice phase. Instead, they observe a bubble phase—a different kind of periodic vortex structure in which the vortices are arranged around bubbles of high particle density. Such states do not occur in the above analysis since they fall outside the scope of the analytic treatment used.

Yi and Pu [144] also conducted a similar study of vortex lattice geometry based on numerical simulations of the GPE, however their results are not in agreement with those obtained above, nor with those of Zhang and Zhai [202] and Cooper et al [201, 203]. They only observe triangular lattice geometries in their simulations, and conclude that the square and rectangular lattice geometries do not exist. Possible explanations for this discrepancy include that the particular parameter values chosen for the simulations do not fall in the square and rectangular lattice regions and/or the simulations were not in the LLL regime.

Although it is no longer appropriate to assume that the wavefunction envelope has cylindrical symmetry for the case of off-axis polarisation, there is still another useful symmetry that can be exploited: reflection symmetry. This reflection symmetry occurs about the x  −  z plane—the plane which contains both the polarisation vector and the axis of rotation. In order to derive a new ansatz for the vortex lattice ground state which assumes reflection symmetry, only minor modifications need to be made to the derivation given in section 8.2. Specifically, it is necessary to introduce two new variational parameters: λ and ζ.

The parameter λ is required to describe the deviation of the condensate cloud from cylindrical symmetry. It is the ratio of the width of the condensate cloud along the y-direction divided by the width along the x-direction. If the density profile of the cloud is expressed in the form $\exp \left[-{{x}^{2}}/l_{x}^{2}-{{y}^{2}}/l_{y}^{2}\right]$, then the aspect ratio is written as $\lambda ={{l}_{y}}/{{l}_{x}}$. For a cylindrically symmetric BEC, the width of the cloud along the x- and y-directions must be the same, which implies that $\lambda =1$. For $\alpha >0$ (see figure 4), since dipolar BECs elongate along the direction of polarisation, it is expected that ${{l}_{x}}>{{l}_{y}}$ and hence $0<\lambda <1$.

The other new parameter, ζ, is required to allow the vortex lattice to adopt any orientation with respect to the polarisation direction. It is defined to be the angle between the first lattice basis vector, ${{\mathbf{b}}_{\mathbf{1}}}$, and the projection of the dipole polarisation onto the plane of rotation. For a cylindrically symmetric BEC, the energy will be independent of ζ. However, for non-cylindrically symmetric dipolar BECs ($\alpha >0$), it is conceivable that the energy may depend on ζ.

By modifying the derivation of the ansatz for the vortex lattice ground state to incorporate the new parameters, it is possible to show that the condensate density must be of the following form,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{n}_{\bot}}\left(\boldsymbol{\rho}\right) & = & \frac{2N\lambda}{\pi {{\chi}^{2}}\left(1+{{\lambda}^{2}}\right)}{{\text{e}}^{-\frac{2}{1+{{\lambda}^{2}}}\left(\frac{{{\lambda}^{2}}{{x}^{2}}+{{y}^{2}}}{{{\chi}^{2}}}\right)}} \\ {} & {} & \times \,\,\underset{\mathbf{q}}{\sum}{{{\overline{g}}}_{\mathbf{q}}}{{\text{e}}^{\text{i}\mathbf{q}\cdot {{{\hat{R}}}_{\zeta}}\boldsymbol{\rho}}}, \end{array}\end{eqnarray} \tag{ 114 }$$

where $\hat{R}$ represents the standard two-dimensional rotation operator, $\chi ={{[{{({{l}_{x}}{{l}_{y}})}^{-1}}-\pi v_{c}^{-1}]}^{-1/2}}$ and,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{{\overline{g}}}_{\mathbf{q}}}=\frac{{{g}_{\mathbf{q}}}}{\sum\nolimits_{\mathbf{v}}{{g}_{\mathbf{v}}}\exp \left(-\frac{{{\chi}^{2}}\left(1+{{\lambda}^{2}}\right)}{8}[\mathbf{v}\cdot {{{\hat{R}}}_{\zeta}}{{\mathbf{\hat{x}}}}/{{\lambda}^{2}}+\mathbf{v}\cdot {{{\hat{R}}}_{\zeta}}{{\mathbf{\hat{y}}}}]\right)}. \end{array}\end{eqnarray} \tag{ 115 }$$

From this starting point is possible to generalize the approach presented in section 8.2 to the case where $\alpha >0$. The total energy can be broken into two components,

$$\begin{eqnarray}\begin{array}{*{35}{l}} E\left(\tau,\eta,\lambda,\zeta \right) & = & {{E}_{0}}\left(\lambda \right)+{{E}_{\text{int}}}\left(\tau,\eta,\lambda,\zeta \right). \end{array}\end{eqnarray} \tag{ 116 }$$

The non-interacting component of the energy is given by,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{E}_{0}}\left(\lambda \right) & = & N\left(\hbar \Omega +\frac{1}{2}\hbar {{\omega}_{z}}\right) \\ {} & {} & +\,\,\frac{N\hbar {{\chi}^{2}}\left(\omega _{\bot}^{2}-{{\Omega}^{2}}\right)}{2 \Omega \left(l_{x}^{2}+l_{y}^{2}\right)}{{\left(\frac{1+{{\lambda}^{2}}}{2\lambda}\right)}^{2}}, \end{array}\end{eqnarray} \tag{ 117 }$$

which in the limit $\Omega \to {{\omega}_{\bot}}$ tends to a constant. Following the same procedure as in section 8.2, using equation (92), the interaction energy is given by,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{E}_{\text{int}}}\left(\tau,\eta,\lambda,\zeta \right) & = & \frac{{{N}^{2}}{{C}_{\text{dd}}}}{6{{\left(2\pi \right)}^{3/2}}{{l}_{z}}{{\chi}^{2}}}\left[2{{\left(\frac{\chi}{{{l}_{z}}}\right)}^{2}}\mathcal{W}\left(\tau,\eta,\lambda,\zeta \right) \right. \\ {} & {} & \left. +\,\,\left(1+3\cos \left(2\alpha \right)+\frac{2}{{{\varepsilon}_{\text{dd}}}}\right)\mathcal{I}\left(\tau,\eta,\lambda \right)\right], \end{array}\end{eqnarray} \tag{ 118 }$$

where,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathcal{I}\left(\tau,\eta,\lambda \right)\approx \frac{{{\lambda}^{2}}}{1+{{\lambda}^{2}}}\underset{\mathbf{q}}{\sum}{{\left({{{\tilde{g}}}_{\mathbf{q}}}\right)}^{2}} \end{array}\end{eqnarray} \tag{ 119 }$$

and,

$$\begin{eqnarray}\begin{array}{*{35}{l}} \mathcal{W}\left(\tau,\eta,\lambda,\zeta \right) & \approx & \underset{\mathbf{q}}{\sum}{{\left({{\tilde{{g}}}_{\mathbf{q}}}\right)}^{2}}{{\text{e}}^{-{{\chi}^{2}}\frac{\left(1+{{\lambda}^{2}}\right)}{2}\left[{{\left(\frac{q_{x}^{\left(-\zeta \right)}}{\lambda}\right)}^{2}}+q_{y}^{\left(-2\zeta \right)}\right]}} \\ {} & {} & \times \,\,{{A}_{2}}\left(\frac{2q_{x}^{\left(-\zeta \right)}}{\lambda},2q_{y}^{\left(-\zeta \right)}\right), \end{array}\end{eqnarray} \tag{ 120 }$$

where $q_{x}^{\left(-\zeta \right)}$ ($q_{y}^{\left(-\zeta \right)}$) represents the x-component (y-component) of ${{\hat{R}}_{-\zeta}}\mathbf{q}$ and,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{A}_{2}}\left({{q}_{x}},{{q}_{y}}\right) & = & \frac{3}{2\sqrt{\pi}}{\int}_{0}^{\infty}\text{d}u\,{\int}_{0}^{2\pi}\text{d}{{u}_{\phi}}\,{{u}^{2}}{{\text{e}}^{{{u}^{2}}}}\text{erfc}(u) \\ {} & {} & \times \,\left({{\sin}^{2}}\alpha {{\cos}^{2}}{{u}_{\phi}}-{{\cos}^{2}}\alpha \right) \\ {} & {} & \times \,\exp \left[-\frac{{{\chi}^{2}}\left(1+{{\lambda}^{2}}\right)}{2l_{z}^{2}}\left(\frac{{{u}^{2}}{{\cos}^{2}}{{u}_{\phi}}}{2}\right)\right] \\ {} & {} & \times \,\exp \left[-\frac{{{\chi}^{2}}\left(1+{{\lambda}^{2}}\right)}{2l_{z}^{2}}{{u}^{2}}{{\sin}^{2}}{{u}_{\phi}}\right] \\ {} & {} & \times \,\,\exp \left[-\frac{{{\chi}^{2}}\left(1+{{\lambda}^{2}}\right)}{2l_{z}^{2}}\left(\frac{{{l}_{z}}{{q}_{x}}u\cos {{u}_{\phi}}}{\sqrt{2}\lambda}\right)\right] \\ {} & {} & \times \,\,\exp \left[-\frac{{{\chi}^{2}}\left(1+{{\lambda}^{2}}\right)}{2l_{z}^{2}}\left(\frac{{{l}_{z}}{{q}_{y}}u\sin {{u}_{\phi}}}{\sqrt{2}}\right)\right]. \end{array}\end{eqnarray} \tag{ 121 }$$

Although there are severe computational limitations surrounding the minimisation of the above expression for the condensate energy, it is still possible to make some meaningful calculations if only triangular and square lattices are considered. For example, it is possible to address the question of whether there is a transition between triangular and square lattices. By considering the minimization at two points in parameter space: $\left(\varepsilon _{\text{dd}}^{-1}=0.9,\alpha =\pi /2\right)$ and $\left(\varepsilon _{\text{dd}}^{-1}=0.95,\alpha =\pi /2\right)$ we find evidence for a transition. In the limit $\Omega \to {{\omega}_{\bot}}$ (${{E}_{0}}\left(\lambda \right)\to N\left[\hbar \Omega +\hbar {{\omega}_{z}}/2\right]$) table 2 shows ${{E}_{\text{int}}}\left(\tau,\eta,\lambda,\zeta \right)$ for the two vortex lattices at the two points in parameter space considered. Looking at the results, there is a phase transition in the vortex lattice geometry from triangular to square as a function of ${{\varepsilon}_{\text{dd}}}$. It is also found that the variational parameter ζ is essentially irrelevant, since the minimum value of ${{E}_{\text{int}}}\left(\tau,\eta,\lambda,\zeta \right)$ is found to be the same for any choice of ζ. This suggests that the orientation of the vortex lattice is unaffected by the broken cylindrical symmetry due to the off-axis polarisation. However, the optimal value of λ does deviate from the cylindrically symmetric value of 1. At both points considered, the optimal value is less than 1, which indicates, as expected, that the dipolar BEC is elongated along the direction of polarisation. Interestingly, the optimal value of λ does not depend on the lattice geometry.

Table 2. Results of the minimisation of ${{\tilde{{E}}}_{\text{int}}}=(3{{(2\pi )}^{\frac{3}{2}}}l_{z}^{3})/({{N}^{2}}{{C}_{\text{dd}}}){{E}_{\text{int}}}$ with respect to λ and ζ, for triangular and square lattices, at two values of ${{\varepsilon}_{\text{dd}}}$ with $\alpha =\pi /2$. The specific parameters used are ${{v}_{c}}/(\pi {{l}_{x}}{{l}_{y}})=1.0191$, $\sqrt{{{l}_{x}}{{l}_{y}}}/{{l}_{z}}=40$ and $\chi /{{l}_{z}}=292.33$. The optimal value of ζ is not included because ${{\tilde{E}}_{\text{int}}}$ was found to be independent of ζ. The results show that the optimal value of λ is the same for both triangular and square lattice geometries.

$\varepsilon _{\text{dd}}^{-1}$	Phase	Optimal λ	${{\tilde{E}}_{\text{int}}}$
0.9	Square	0.93	−4979.49
0.9	Triangular	0.93	−4788.31
0.95	Square	0.78	−1377.22
0.95	Triangular	0.78	−1426.81

Numerical studies [144], based on the dipolar GPE, suggest that for $\alpha >0$ the vortex lattice can undergo a phase transition from triangular structure to a non-triangular structure as ${{\varepsilon}_{\text{dd}}}$ is increased. This is consistent with analysis presented above, however, to our knowledge there has not, to date, been a thorough study of vortex lattice structures for the regime of $\alpha >0$.

The theoretical study of non-dipolar two-component BECs [208, 209] has shown how interspecies interactions g₁₂ can influence the vortex lattice structure of the two components. The analysis presented in section 8 can be adapted to analyse the vortex lattice structure of two component condensates. Doing this Mueller and Ho [208] were able to quantify various regimes through the parameter $\beta ={{g}_{12}}/\sqrt{{{g}_{1}}{{g}_{2}}}$, where g₁₂ is the interspecies interaction parameter and g₁₍₂₎ is the intraspecies interaction parameter for component 1 (2) of the two-component BEC. This treatment and subsequent numerical analysis [209] shows that for $\beta <0$ the two components overlap and a single triangular vortex lattice arises. For $\beta >0$ the two components separate to form interlaced triangular vortex lattices. As $\beta \sim 1$ the triangular lattice distorts to form square or rectangular arrays. This theoretical analysis is consistent with experiment [210].

It is also possible to consider two-component dipolar BEC systems, with both interspecies and intraspecies contact and dipolar interactions. Work by Shirley et al [211] showed that such systems (where one of the components has zero dipolar interactions), under rotation, exhibit a rich phase diagram, which includes triangular vortex lattices, square vortex lattices, vortex sheets (where half quantum vortices of one component align in a winding sheet, which is interwoven with a sheet in the other component [212]), half quantum vortex chains (where vortices, alternating between each component, line up along a chain) and half quantum vortex molecules (where a vortex in a given component pairs up with a vortex in the other component). This analysis is consistent with further studies [213–215] and has been extended to consider how dipole alignment in the plane of rotation [215] and component-dependent optical lattices [216] influences the phase diagram.

The aim of this review was to take the reader on a journey, starting with the fundamental concepts and methodologies used to understand the properties of dipolar BECs and then show how these have been applied to understand the properties of vortices and vortex lattices in these systems. Throughout, dipolar interactions are seen to enrich the physical properties of the system and vortices therein. The journey started in earnest in section 3 where we met the dipolar interactions; these introduce a long-range and anisotropic component to the interactions, making a significant departure from conventional s-wave interactions which appear in the theory as isotropic contact interactions. In section 4 we saw that the dipolar interaction significantly modifies the fundamental stationary solutions of a dipolar condensate, including the introduction of collapse instabilities dependent on the shape of the boundary relative to the polarization direction. In section 5 we found that dipolar interactions can alter the energy and structure of a single vortex. Specifically, in quasi-two-dimensional systems when the dipole alignment is in the plane of the condensate, the vortex core is no longer circularly symmetric. Additionally, density ripples appear in the vicinity of the vortex core, as the roton instability is approached, due to the roton mixing with the ground state of the system. In section 6 we found that the interaction between vortices can be altered by the absence of dipoles in the vortex core, introducing an additional long-range and anisotropic contribution to the vortex–vortex interaction. This can, for instance, lead to the suppression of the annihilation of vortex–antivortex pairs and induce the co-rotational dynamics of vortex–vortex pairs to become anisotropic. In section 7 we summarised the methods for generating vortices in condensates, and discussed the role of dipolar interactions in these processes. Concentrating on the properties and instabilities of rotating condensates, dipolar interactions were shown to significantly alter the regimes of stability and the critical rotation frequencies for vortices to be nucleated. This also allowed us to identifiy routes to vortex formation under rotation. Finally, in section 8 we found that dipolar interactions lead to new and exotic vortex lattice phases; whereas lattices in non-dipolar condensates are well-known to follow a triangular pattern, dipolar interactions can support rectangular, square and bubble phases.

Of course any journey is a compromise between taking an efficient route and a scenic path, i.e. in this case a compromise between completing the review, in a timely manner, and detailing every contribution to the field. Unfortunately our path has been fairly efficient and as such we have omitted several other aspects relating to vortices and vortex lattices in quantum ferrofluids. Below we, all too briefly, provide a snapshot of some of the scenery we have missed along the way and avenues for further exploration.

9.2.1. Dipolar Bose–Einstein condensates in toroidal traps.

9.2.2. Fractional quantum Hall physics in dipolar Bose–Einstein condensates.

9.2.3. Dipolar fermions.

9.2.4. Berezinskii–Kosterlitz–Thouless transition.

In a strictly homogeneous two-dimensional system at finite temperature, long range phase coherence cannot be established and condensation will not occur. Superfluidity can still occur at very low temperatures where quasi-long-range order exists [266, 267], but even this is destroyed when the temperature is raised through the Berezinskii–Kosterlitz–Thouless transition [268, 269] at which point spatial correlations change from power law to exponential decay. Physically, this transition can be thought of as occurring when virtual vortex–antivortex pairs unbind and there is a proliferation of free vortices. The Berezinskii–Kosterlitz–Thouless and BEC transition have been studied in trapped ultracold gases via observations of phase defects [270], vortices [271] and changes in the density profile due to the onset of superfluidity [272]. The Berezinskii–Kosterlitz–Thouless transition in the presence of dipolar interactions has been studied using, for a homogeneous system, Monte Carlo methods [273], the mean-field Hartree–Fock–Bogoliubov–Popov model [274], the dipolar XY-model [275] and most recently renormalized Hartree–Fock theory [77]. A simple picture underpinning our understanding of the Berezinskii–Kosterlitz–Thouless transition comes from asking the question: what is the energy required to separate a vortex–antivortex pair? For vortex–antivortex pairs in a uniform two-dimensional non-dipolar BEC (of two-dimensional background density n₀), the energy of a pair, separated by a distance d and calculated from hydrodynamical arguments, is approximately given by $V(d)=2\pi {{\hbar}^{2}}{{n}_{0}}\ln \left(d/\xi \right)/m$. The critical temperature associated with the transition is given by the relation $2\pi {{\hbar}^{2}}{{n}_{0}}/\left(m{{k}_{\text{B}}}T\right)=4$. This is calculated by determining the average distance between the pairs,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle {{d}^{2}}\rangle =\frac{{\int}_{\xi}^{\infty}{{d}^{3}}{{\text{e}}^{-\beta V(d)}}\,\text{d}d}{{\int}_{\xi}^{\infty}d{{\text{e}}^{-\beta V(d)}}\,\text{d}d}={{\xi}^{2}}\frac{\Gamma -2}{\Gamma -4}, \end{array}\end{eqnarray} \tag{ 122 }$$

where we have used the short-hand notation $\Gamma =2\pi {{\hbar}^{2}}{{n}_{0}}/\left(m{{k}_{\text{B}}}T\right)$. This result diverges at $\Gamma =4$, signifying the Berezinskii–Kosterlitz–Thouless transition.

In this simple model the inclusion of dipolar interactions adds three complications. The first complication is that as dipolar interactions are introduced ξ will change. The second complication is, as seen in section 6, the interaction between a vortex and an antivortex is modified in the presence of dipolar interactions due to the absence of dipoles in the vortex cores. As such the energy scaling of the pair with separation, i.e. V(d), is changed. The third complication arises if the dipole alignment has some component in the two-dimensional plane of the gas. In this case the interaction between the vortex–antivortex pair is no longer just a function of the distance between them, it also depends on the in-plane angle of the pair relative to the polarization direction, i.e. $V(d)\to V\left(d,\eta \right)$ (see figure 13(b)). The above generalisations can be incorporated into equation (122) by considering the following analysis [276]. To quantify if dipolar interactions have a significant effect on the Berezinskii–Kosterlitz–Thouless transition temperature we consider a simple model for the interaction between a vortex–antivortex pair. Specifically, the vortices carry a dipole moment $\propto \left(0,\sin \alpha,\cos \alpha \right)$. If the separation between the vortices is $\mathbf{d}=d\left(\cos \eta,\sin \eta,0\right)$ then the interaction between the vortex and the anti-vortex can be written as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{V\left(d,\eta \right)}{{{k}_{\text{B}}}T}= \Gamma \ln \left(\frac{d}{\xi}\right)-\lambda {{\left(\frac{\xi}{d}\right)}^{3}}f\left(\alpha,\eta \right). \end{array}\end{eqnarray} \tag{ 123 }$$

For simplicity, we have assumed that the interaction between vortices arising from the absence of dipoles in their cores is $\propto 1/{{d}^{3}}$. Given the results in section 5.4.1, which finds that due to the absence of dipoles in a vortex core ${{\Phi}_{\text{dd}}}$ has a term $\propto \ln (d)/{{d}^{3}}$, this represents something of a simplification. From this

$$\begin{eqnarray}\begin{array}{*{35}{l}} \langle {{d}^{2}}\rangle & = & \frac{{\int}_{\xi}^{\infty}{{d}^{3- \Gamma}}\,\text{d}d{\int}_{0}^{2\pi}\exp \left[\lambda f\left(\alpha,\eta \right){{\left(\frac{\xi}{d}\right)}^{3}}\right]\text{d}\eta}{{\int}_{\xi}^{\infty}{{d}^{1- \Gamma}}\,\text{d}d{\int}_{0}^{2\pi}\exp \left[\lambda f\left(\alpha,\eta \right){{\left(\frac{\xi}{d}\right)}^{3}}\right]\text{d}\eta} \end{array}\end{eqnarray} \tag{ 124 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & = & {{\xi}^{2}}\left[\frac{\sum\nolimits_{n=0}^{\infty}{{a}_{n}}{{\lambda}^{n}}{{\left(\Gamma +3n-4\right)}^{-1}}}{\sum\nolimits_{n=0}^{\infty}{{a}_{n}}{{\lambda}^{n}}{{\left(\Gamma +3n-2\right)}^{-1}}}\right], \end{array}\end{eqnarray} \tag{ 125 }$$

where

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{a}_{n}}=\frac{1}{n!}{\int}_{0}^{2\pi}{{f}^{n}}\left(\alpha,\eta \right)\text{d}\eta. \end{array}\end{eqnarray} \tag{ 126 }$$

When $\lambda =0$ the result obtained from equation (122) is regained, i.e. $\langle {{d}^{2}}\rangle ={{\xi}^{2}}\left(\Gamma -2\right)/\left(\Gamma -4\right)$, with the Berezinskii–Kosterlitz–Thouless transition arising from the divergence at $\Gamma =4$. However, the inclusion of the additional power law interaction between the vortex and antivortex, arising from the absence of dipoles in the region of the cores, does not fundamentally change this result. Specifically, the mean-square separation, given by equation (125), diverges first at $\Gamma =4$, i.e. this simple analysis implies that although the dipolar interaction can change the interaction between vortex–antivortex pairs the long-range hydrodynamic interaction always wins, suggesting that the Berezinskii–Kosterlitz–Thouless transition is unaffected by dipolar interactions. This is consistent with more detailed analysis presented in [273, 275] which found, at best, a weak dependence of the Berezinskii–Kosterlitz–Thouless transition on the strength of the dipolar interactions. For example, working within the XY-model Vasiliev et al [275] concluded that polarization of the system (i.e. the generation of additional vortex–antivortex pairs between any two points) screens out the long-range interactions. However, in our opinion there is still room for more work in this area such as considering finite-size effects and the true form of the dipolar interaction between a vortex–antivortex pair in order to firmly establish the role of dipolar interactions on this phase transition in physical systems.

9.2.5. Vortex lattices in the supersolid phase.

9.2.6. The Onsager vortex phase transition.