Vortices and vortex lattices in quantum ferrofluids

The experimental realization of quantum-degenerate Bose gases made of atoms with sizeable magnetic dipole moments has created a new type of fluid, known as a quantum ferrofluid, which combines the extraordinary properties of superfluidity and ferrofluidity. A hallmark of superfluids is that they are constrained to rotate through vortices with quantized circulation. In quantum ferrofluids the long-range dipolar interactions add new ingredients by inducing magnetostriction and instabilities, and also affect the structural properties of vortices and vortex lattices. Here we give a review of the theory of vortices in dipolar Bose-Einstein condensates, exploring the interplay of magnetism with vorticity and contrasting this with the established behaviour in non-dipolar condensates. We cover single vortex solutions, including structure, energy and stability, vortex pairs, including interactions and dynamics, and also vortex lattices. Our discussion is founded on the mean-field theory provided by the dipolar Gross-Pitaevskii equation, ranging from analytic treatments based on the Thomas-Fermi (hydrodynamic) and variational approaches to full numerical simulations. Routes for generating vortices in dipolar condensates are discussed, with particular attention paid to rotating condensates, where surface instabilities drive the nucleation of vortices, and lead to the emergence of rich and varied vortex lattice structures. We also present an outlook, including potential extensions to degenerate Fermi gases, quantum Hall physics, toroidal systems and the Berezinskii-Kosterlitz-Thouless transition.

Ferrohydrodynamics describes the motion of fluids comprised of particles with significant magnetic (or electric) dipole moments [1,2].The dipole-dipole (from henceforth: dipolar) inter-particle interaction causes magnetostriction (or electrostriction) and gives rise to spectacular instabilities such as the normal field instability [3] that can lead to complex pattern formation.Classical ferrofluids have been investigated since the 1960s, the first ferrofluid having been invented at NASA with the intention of making a jet fuel whose flow could be directed in a zero-gravity environment using a magnetic field [4].They have subsequently found a broad range of applications reaching from liquid seals around rotating shafts (such as in hard disks), where the fluid is held in place using magnets, to magnetically targeted drugs in medicine [5].However, the focus of this review is upon quantum ferrofluids as first realized in 2005 with the creation of a Bose-Einstein condensate (BEC) in a vapour of 52 Cr atoms by the Stuttgart group [6].These atoms have a magnetic dipole moment of 6µ B , six times larger than that found in the alkalis which are used in the majority of BEC experiments. 52Cr atoms therefore have dipolar interactions which are 36 times larger than in standard BECs.Other groups have also studied 52 Cr BECs [7,8], as well as BECs made of atoms with even larger magnetic dipoles such as 164 Dy [9,10], and 168 Er [11].Many of the signatures of ferrohydrodynamic behaviour have now been observed in these gases, including magnetostriction [12] ‡, collapse due to dipolar interactions [13,14], and the quantum analog of the Rosensweig instability [10,15].Recently, fully self-bound dipolar droplets have been reported [16,17].Additionally, the production of ultracold fermionic 40 K- 87 Rb [18] polar molecules and the cooling of fermonic 161 Dy [19] and 167 Er [20], all with significant dipole moments, pave the way for a new generation of quantum degenerate Fermi gas experiments, where dipolar interactions dominate.In Fermi gas systems the partially attractive nature of the dipolar interaction opens up the possibility of Bardeen-Cooper-Schrieffer (BCS) pairing at sufficiently low temperatures [21,22,23,24,25,26,27,28,29].Excellent reviews of the field of ultracold dipolar gases can be found in Refs.[30,31,32,33,34].
In geometries approaching the onedimensional limit, so-called solitonic vortices have been formed [45,46] which share properties between vortices and their one-dimensional analogs: dark solitons.‡ This paper coined the term "quantum ferrofluid".
Several reviews exist which summarise the significant experimental and theoretical aspects of vortices and vortex lattices in non-dipolar BECs [47,48,49,50,51,52].Vortices have yet to be observed in quantum ferrofluids, although numerical simulations suggest the formation of vortex rings in the dipolar collapse experiment of Ref. [13], and the formation of vortexantivortex pairs [53] in the droplet experiment of Ref. [10].
Here we establish the properties of vortices and vortex lattices in quantum ferrofluids, reviewing the theoretical progress that has been made over the last decade.Whilst it is possible to also consider the properties of vortices and vortex lattices in dipolar Fermi gases, this review is confined to the bosonic case and only a brief discussion of fermionic systems will be given in the Summary and Outlook (Section 9).The structure of this review is built upon the philosophy of taking the reader on a journey.This journey starts in Section 2 where the properties of classical and strongly correlated quantum ferrofluids are briefly discussed.Sections 3 and 4 provide a brief introduction to the properties of dipolar BECs in the absence of vortices.In Section 3 we examine the mathematical form of the dipolar interaction in quantum ferrofluids, and present the most widely used model for quantum ferrofluids -the dipolar Gross-Pitaevskii equation (GPE) -along with its hydrodynamical interpretation.Section 4 builds on this theory to consider the stability of dipolar BECs.
Specifically, we look at stability in the Thomas-Fermi (hydrodynamic) regime, where interactions dominate, and more general dipolar GPE solutions, in three-dimensional and quasi-two-dimensional systems.
Section 5 focuses on the properties of single vortex lines in threedimensional condensates and single vortices in quasitwo-dimensional systems.In Section 6 we consider vortex-vortex and vortex-antivortex dynamics in quasitwo-dimensional dipolar BECs, primarily focusing on solutions of the dipolar GPE.Section 7 addresses the routes to vortex and vortex lattice formation.This focuses primarily on stationary solutions (in the rotating frame) and their dynamical stability, enabling us to ascertain under what conditions it might be expected that vortices will nucleate into the dipolar BEC.Section 8 analyses how dipolar interactions can induce changes to vortex lattice structures.This revisits previous work and presents some new variational calculations that elucidate the properties of vortex lattice structures in dipolar condensates.In Section 9 we give a brief summary and provide an outlook to several topical aspects for future development which are not covered in the main body of the review.Prospects for quantum turbulence with quantum ferrofluids, which are not covered here, are discussed elsewhere [54].

Classical Ferrofluids and Strongly Correlated Quantum Ferrofluids
A classical ferrofluid can be formed from a suspension of small permanently magnetized particles, with a typical size of 0.01 µm, in a non-magnetic solvent [1].
Their most closely related electrical counterparts are known as electrorheological fluids which are suspensions of electrically polarizable particles, typically 1 − 100 µm in size, in an insulating solvent [55,56,57] (there are also magnetorheological fluids where small micelles of magnetizable fluid are suspended in a nonmagnetizable fluid [58,59].) In the presence of an external field, a classical ferrofluid can form a zoo of different patterns including hexagonal cells [3], columns [60], stripe and bubble phases [61], and disordered stripe phases producing labyrinthine structures [1].Some of these patterns also occur in quantum ferrofluids: stripe phases (density wave modulations) will be discussed in the quantum case in the absence of vortices in Section 4.2.2 and in the presence of vortices beginning in Section 5.2.Vortex lattices in quantum fluids tend to form hexagonal patterns (Abrikosov lattices) even in the absence of dipolar interactions where the longrange logarithmic hydrodynamic interaction between vortices plays an important role.However, as we shall see in Sections 8.2 and 8.3, the presence of dipolar interactions can change the lattice configuration to square and bubble geometries.
The above patterns can all be captured to some degree by mean-field theory where no attempt is made to describe the fluid at the molecular level.However, it has long been predicted that with strong dipolar interactions ferrofluid molecules can form chains and rings [62,63,64].In order to describe such strong correlation effects, which lie beyond the mean-field description, it is necessary to use computationally intensive Monte Carlo and Molecular Dynamics techniques [65,66,67,68].These methods suggest a complex and rich phase diagram.In zero applied field at relatively low temperatures, and as the density is increased from low to high, the molecules initially form rings which unbind into chains which then cluster into networks, which break down into a normal liquid, then form a ferroelectrically ordered liquid, followed by a possible ferroelectric columnar ordering, and finally form a ferroelectric solid.In the presence of an external field chains, columns, sheets, bent walls, lamellar, labyrinthine or wormlike structures, and hexagonal structures all appear [68].Transitions from single chains to double chains (zig-zag) can also occur if initially strong transverse confinement is reduced [69].
In the rest of this review we restrict ourselves to mean-field phenomena.It is important to mention, however, that there is a sizeable body of theoretical work in strongly correlated quantum dipolar systems in two dimensions.Strongly correlated dipolar gases do not yet exist in the laboratory but ideas to realize them include using ultracold molecules with very large dipole moments dressed by microwaves [70] and Rydberg dressed ultracold gases [71].One of the main interests in these systems is the formation of so-called supersolids which are both crystalline and yet also have superfluid properties [70,72,73,74,75,76].Unlike the density wave structures we shall study later in this review, which have many atoms per wavelength, in the strongly correlated case the periodicity can be at the single atom or few atom length scale.Liquid crystal phases have also been identified [77].

Quantum Ferrofluids: Theory and Basic Properties
The successful Bose-Einstein condensation of gases of 52 Cr atoms [6,7], 164 Dy [9, 10] and 168 Er [11] have realized BECs with significant dipolar interactions.A basic property of these interactions is that their net effect depends on the shape of the BEC, as illustrated in Figure 1.For dipoles polarised along the long axis of a prolate (elongated) dipolar gas [Figure 1(a)] the net contribution to the dipolar interaction is attractive.By contrast, for dipoles polarised along the the short axis of an oblate (flattened) dipolar gas [Figure 1(b)] the net contribution to the dipolar interaction is repulsive.Compared to BECs with interactions which are dominated by isotropic s-wave scattering, a dipolar BEC will be elongated along the direction of the polarising field (magnetostriction) [78,79,80].

The dipolar interaction
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], 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 r − r .The short-range interactions are characterized by the coupling constant g = 4π 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 where µ 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 E = kE, for which the coupling constant is [82,83], where α p is the static polarizability and 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 40 K 87 Rb [84] and 133 Cs 87 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 θ = 0, representing the attraction of headto-tail dipoles, and positive for θ = π/2, representing the repulsion of side-by-side dipoles.At the "magic angle", θ m = arccos(1/ √ 3) ≈ 54.7 • , the dipolar interaction is zero.
It is often convenient to work in momentum space.The Fourier transform Ũdd (k) = F[U dd ] = U dd (r)e −ik•r dr of the dipolar interaction is [81], where α is the angle between k and the polarization direction.
The strength of the dipolar interactions is conveniently parameterized by the ratio [32], where g can be tuned between −∞ and +∞ via a Feshbach resonance [86,87].In effect ε 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 ε dd > 1 -see Section 4.1.1.For 52 Cr, 168 Er and 164 Dy, the natural value of ε 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 −∞ < ε dd < ∞, with both negative and positive C dd .Note that for C dd < 0 the dipole-dipole interaction becomes repulsive for head-to-tail dipoles and attractive for side-by-side dipoles.

The dipolar Gross-Pitaevskii equation
In the mean-field limit, at zero temperature, a single wavefunction, Ψ(r, t), can be used to describe the condensate.The condensate density n(r, t) where N is the number of atoms in the condensate.The wavefunction obeys the dipolar GPE [78,81,82], where V = V (r) 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|Ψ| 2 , arises from the van der Waals interactions and the non-local term, Φ, arises from the dipolar interactions [79]: If we take the dipoles to be polarized along the zdirection, then using identities from potential theory the dipolar potential can be expressed as [91,92], where φ is a fictitious 'electrostatic' potential defined as This effectively reduces the problem of calculating the dipolar potential Φ to one of calculating an electrostatic potential of the form (8) which is easier to compute because the Green's function 1/ |r − r | 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,94,95,96].Alternatively, Φ can be evaluated in momentum space by exploiting the convolution theorem, where ñ(k, t) = F[n(r, t)].
In condensate experiments the external potential V is typically harmonic with the general form, where ω j (j = x, y, z) are the trap's angular frequencies.
In general, this gives rise to three harmonic oscillator lengths, j = /mω 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, where γ = ω z /ω ⊥ is the so-called trap ratio.When γ 1 the BEC shape will typically be oblate (flattened) while for γ 1 it will typically be prolate (elongated).Time-independent solutions of the GPE satisfy, where µ is the chemical potential §.Inserting this into Eq.( 5), the time-independent dipolar GPE for the time-independent wavefunction ψ(r) is 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, 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, which may be interpreted as the minimum lengthscale over which the wavefunction changes appreciably.
Efficient numerical methods for solving the dipolar GPE are available [97,98,99,100,101] and progress has been made on extending this treatment to include finite temperature effects and quantum fluctuations [102,103,104,105].

Dipolar hydrodynamic equations
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 Ψ(r, t) = n(r, t)e iS(r,t) , where the local phase, S(r, t), defines the fluid velocity field v(r, t), Inserting the Madelung form into the GPE, and separating real and imaginary terms, yields two § Throughout this review Ψ (ψ) denotes the time-(in)dependent condensate wavefunction.
equations which together are exactly equivalent to the GPE.The first is the continuity equation, This embodies the conservation of the number of atoms.The second equation is (18) The ∇ 2 √ n/ √ 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, Eqs. ( 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,107,108,109,110].Equations ( 17) and (18) have been extended to include dipolar interactions and are referred to as the dipolar superfluid hydrodynamic equations.

Vortex-Free Solutions and Stability
Before discussing vortices, we next describe the solutions and stability of the dipolar condensates themselves, in homogeneous and trapped systems, and introduce some key analytical tools and physical concepts.

Three-dimensional case
For V (r) = 0 (uniform condensate of infinite extent), the stationary solution is, i.e. a state of uniform density n 0 .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 (ε dd = 0), the corresponding solution has chemical potential µ = n 0 g.By comparison, the homogeneous dipolar system is akin to a non-dipolar system but with an effective coupling For a three-dimensional homogeneous dipolar condensate, the Bogoliubov dispersion relation between the energy E B and momentum p of a perturbation is given by, where c(θ) is the speed of sound, The angle θ is that between the excitation momentum p and the polarization direction.For low momenta the spectrum is linear E B ≈ c(θ) p and characteristic of phonons with a phase velocity c(θ) 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 p evolves in time as exp(−iE B (p)t/ ).If E B (p) 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 B (p) provided by Eq. ( 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 Eq. ( 21) is real-valued indicates that the three-dimensional homogeneous system is stable against the phonon instability in the range −0.5 ≤ ε dd ≤ 1 for g > 0, and ε dd ≤ −0.5, ε dd > 1 for g < 0.

Trapped dipolar condensates
A full theoretical treatment of a trapped BEC involves solving the dipolar GPE, given in Eq. ( 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 nontrivial, 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,112,113].To characterise the stability of a dipolar condensate we first derive and examine the Thomas-Fermi ground state solutions for a dipolar BEC.

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 / ¯ 1, where ¯ = ( x y z ) 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 Ref. [114].
Although the governing equations for a dipolar BEC contain the non-local potential Φ(r), 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 zdirection, with repulsive van der Waals interactions (a s > 0), and confined by a cylindrically-symmetric trap of the form of Eq. ( 11).We limit the analysis to the regime of −0.5 ≤ ε dd ≤ 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, Making use of the electrostatic formulation given in Eqs. ( 7) and ( 8), exact solutions of Eq. ( 23) can be obtained for any general parabolic trap, as proven in Appendix A of Ref. [92].In particular, the solutions for the density profile take the form, where n cd = 15N/(8πR 2 ⊥ R z ) is the central density, and R z and R ⊥ are the Thomas-Fermi radii of the condensate in the axial and transverse directions.
Remarkably, the inverted parabolic density profile, Eq. ( 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, κ = R ⊥ /R z , as the trap, for dipolar BECs κ = γ (in general) and must be evaluated using the following transcendental equation [91,92], Usually the central density of the Thomas-Fermi profile is denoted n 0 .However to avoid confusion, later in the review, we have used n cd to define the central density of the Thomas-Fermi profile.
which takes the value f = 1 at κ = 0, and monotonically decreases towards f = −2 as κ → ∞, passing through zero at κ = 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 (ε dd = 0) condensate one finds the expected result that κ = γ.However, the presence of dipolar interactions leads to magnetostriction of the condensate, such that κ < γ for ε dd > 0 and κ > γ for ε dd < 0. This behaviour is shown in Figure 3 (top) [117].Note that, within the range −0.5 ≤ ε dd ≤ 1 these are global solutions; elsewhere the solutions are either metastable (light grey shading) or unstable (dark grey shading).For conventional dipoles (C dd > 0, ε dd > 0), the condensate is least stable in prolate (γ < 1) traps; here the dipoles lie predominantly in the attractive head-to-tail configuration and undergo collapse when ε dd becomes too large.By contrast, in oblate (γ > 1) traps stability is enhanced since the dipoles lie predominantly in the repulsive sideby-side configuration.Meanwhile the opposite is true for anti-dipoles (C dd < 0, ε 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, zeropoint kinetic energy (neglected within the Thomas-Fermi approach) can enhance the stability of the solutions.
Once the BEC aspect ratio κ is found from the transcendental equation, the Thomas-Fermi radii are determined by the expressions, and the total energy is given by, 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], 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 Φ TF (r) for the case of a dipolar BEC without cylindrical symmetry is given later in Eq. (71).

Outside the Thomas-Fermi regime: rotons and density oscillations
According to the Thomas-Fermi approach, a trap which is sufficiently oblate (γ > ∼ 5.2) is stable to collapse even in the limit ε dd → ∞.However, numerical solutions reveal a different fate, whereby the condensate undergoes instability even for γ → ∞ [112].This is associated with the development of a roton minimum in the dispersion relation [111,112,132], reminiscent of rotons in superfluid helium [118].For certain parameters, this minimum can approach zero energy, triggering an instability at finite k known as the roton instability.The Thomas-Fermi approach, which is limited to the class of inverted parabolic solutions, is unable to account for this phenomenon.The roton is a strict consequence of the non-local interactions, and does not arise for conventional condensates.
The effect of this in trapped, purely dipolar condensates was revealed by Ronen et al., with the stability diagram shown in Figure 3 (bottom) [119].When the dipolar interaction parameter D = N mC dd /(4π 2 l ⊥ ) ¶ exceeds a critical value, for any trap ratio, the system is unstable to collapse.The condensate becomes unstable to modes with increasingly large number of radial and angular modes as the trap aspect ratio increases, signifying that collapse proceeds on a local, rather than global, scale [120].Of particular interest is the appearance, close to the instability boundary and under oblate traps, of ground state solutions with a biconcave, red blood cell-like, shape [see Figure 3 (bottom)] [119,121].Subsequent works confirmed these density oscillations as being due to the roton, which, for certain parameters, mixes with the ground state of the system [122].More generally, when van der Waals interactions are included [123,124], both biconcave and dumbbell shapes can arise [125].Under box-like potentials, which have been realized in recent years [126,127], density oscillations associated with the roton can arise at the condensate edge [79].
An intuitive interpretation of the roton in an oblate trap was put forward by Bohn et al. [128].As the dipole strength is increased, it is energetically favourable for the dipoles to locally move out of the plane and align head-to-tail perpendicular to the plane, thereby taking advantage of this attractive configuration.This leads to a periodic density in the plane, with a wavenumber corresponding to that of the roton minimum.
In this geometry it is also interesting to note that quantum depletion of the condensate is predicted to diverge at the roton instability [132].

Quasi-two-dimensional dipolar Bose-Einstein condensate
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 ¶ In the original work by Ronen et al. [119] D was defined as D = (N − 1)mC dd /(4π 2 l ⊥ ).However, the derivation of the dipolar GPE requires N 1.Hence, in this review, we have defined D = N mC dd /(4π 2 l ⊥ ). 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) = 1 2 mω 2 z z 2 in the z-direction which satisfies ω z µ, 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, Axially, the condensate is taken to be frozen into the axial ground harmonic oscillator state ψ z (z) = (π 2 z ) −1/4 e −z 2 /2 2 z .The dynamics then become planar, parametrised by the two-dimensional time-dependent wavefunction, Ψ ⊥ (ρ, t).Note that Ψ ⊥ is normalized to the number of atoms, i.e.N = |Ψ ⊥ | 2 dρ.Inserting this ansatz into the dipolar GPE, Eq. ( 5), and integrating out the axial direction then leads to the effective two-dimensional dipolar GPE [130], The g/ √ 2πl z coefficient characterises the effective van der Waals interactions in the plane, and Φ ⊥ is the effective planar dipolar potential, where n ⊥ = |Ψ ⊥ | 2 is the two-dimensional density.
The real-space form of the effective two-dimensional dipolar interaction potential U ⊥ dd is given elsewhere [131], while in this review its Fourier transform is used [130,132], where F (q) = −1 + 3 √ π q2 x q e q2 erfc(q), F ⊥ (q) = 2 − 3 √ π qe q2 erfc(q) and q = ql z / √ 2 with q being the projection of k onto the x − y plane, i.e. the reciprocal space analogue of ρ + .From this momentum space representation Φ ⊥ can then be evaluated using the convolution theorem as in Eq. ( 9).An important parameter is the ratio σ = l z /ξ, where ξ is the healing length.The two-dimensional approximation requires σ < 1.
Under a cylindrically-symmetric harmonic trap and for dipoles polarized along z, the above "twodimensional mean-field regime" is formally entered when N a s (1 + 2ε dd )l 3 z /l 4 ⊥ 1.In the opposing regime, when N a s (1 + 2ε dd )l 3 z /l 4 ⊥ 1, the system enters the three-dimensional Thomas-Fermi regime [114].A more general analysis of flattened condensates in Ref. [133], has established the validity of the two-dimensional mean-field regime for arbitrary polarization direction.
In the absence of any planar trapping potential [V (ρ) = 0], the stationary solution of the quasi-twodimensional dipolar condensate is the homogeneous state [130], where n 0 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 dd [3 cos 2 α −1]/3 < 0. The phonon unstable regions in the ε dd − α 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 α m , across which the planar dipolar interactions switch between repulsive and attractive.
+ Throughout this review k is the reciprocal lattice vector in three spatial dimensions r and q is the reciprocal lattice vector in two spatial dimensions ρ.The roton instability also arises in this quasi-twodimensional BEC [112], as indicated in Figure 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 α = 0; for α = 0 the condensate cannot probe the attractive part of U ⊥ dd [132] (this is true only in the strict quasi-twodimensional 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 outof-plane component of U ⊥ dd decreases [129] and the regimes of roton instability (blue regions in Figure 5) shrink.

Single Vortices
Quantized vortices are a consequence of the condensate's phase coherence.To preserve the singlevaluedness of the condensate wavefunction, the total change in phase around some closed path C must be 2πq v , with q v = 0, ±1, ±2, ....If q v = 0 then there exists one or more phase singularities within C. These singularities are the quantized vortices, and q v is the vortex charge.
Consider an isolated vortex at the origin in a uniform system, which is straight along z.The condensate phase about the vortex follows the azimuthal angle, S(ρ, θ, z) = θ = arctan(y/x).From Eq. ( 16) this gives rise to a circulating azimuthal flow with speed, Within the fluid the flow is irrotational with zero vorticity, i.e. ∇ × v = 0.This can also be seen directly from the definition of the fluid velocity, which is curlfree ∇ × v = 0 and thus very different to the classical solid-body rotation, for which v ∝ ρ.At the point of the singularity, however, the vorticity takes the finite value q v h/m.Also, at this point the density is zero, preventing the unphysical scenario of infinite mass current.The cylindrically-symmetric flow associated with the velocity field given in Eq. ( 36) carries a total angular momentum L z = N q v .

Energetics of vortex formation
Even single vortices are giant excitations involving a considerable fraction of the entire BEC.The energy associated with the formation of a vortex E v ≡ E −E 0 , where E 0 is the energy of the non-rotating (vortexfree) state, is generically much larger than the energy of elementary excitations described by the Bogoliubov spectrum given in Eq. (21).In a frame rotating at angular frequency Ω, the total energy of the system is shifted to E − Ω • L, where 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 |Ω•L| > E v leading to a critical rotation frequency, E v can be computed analytically in certain situations as we discuss below.Before we do so, it is important to point out that Eq. ( 37) considers the energetics but not the kinetics of vortex formation.Both theory and experiment reveal that the true value of Ω v is often considerably higher than predicted by Eq. ( 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, Substituting into the dipolar GPE (5) an equation for the amplitude f about the vortex is obtained, ) where the Laplacian has been expressed in its cylindrically symmetric form, The second term on the right-hand side of Eq. ( 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 v | > 1 vortices are energetically unstable compared to multiple singly-charged vortices, and rarely arise unless engineered; for this reason we will consider only In any system uniform along the polarization direction (z) the dipolar potential reduces to a local potential Φ(r) → −gε dd n(r).This can be seen from Eq. ( 7), which gives the relationship between the fictitious electrostatic potential φ(r) and the dipolar potential Φ(r), and noting that ∂ 2 φ/∂z 2 must equal zero * .
Thus, the last two terms in Eq. ( 39) can be combined into a single contact term g(1 − ε dd )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(1 − ε dd ).For example, analysis of Eq. ( 39) for the s-wave case reveals that the centrifugal barrier term dictates that the density relaxes as 1/ρ 2 to the asymptotic background value √ n 0 as ρ → ∞, and for ρ → 0 the density tends to zero as ρ 2|qv| [108].Furthermore, although Eq. ( 39) can not be solved in terms of known functions, it can be solved numerically, and with appropriate scaling the result for f (ρ) is universal.This solution for f (ρ) can then be used in the energy functional given in Eq. ( 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], This result was originally obtained for superfluid 4 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 ξ, their finite values avoid a logarithmic singularity that originates from the centrifugal barrier term.The only place that interactions enter this expression for v is through the healing length ξ = / √ mµ which can immediately be adapted to the dipolar case using µ = n 0 g(1 − ε dd ).The critical rotation frequency for this case can now be evaluated by replacing E v by v and N by N/L (number of atoms per unit length) in Eq. (37).

General features of a vortex in a quantum ferrofluid
In order to obtain analytic results for the untrapped system it can sometimes be useful to use the following approximate solution for f (ρ) which incorporates the correct behaviour for ρ → 0 and ρ → ∞, where ρ = ρ/ξ and n 0 is the density at infinity.It will be convenient later to write this in the form of a homogeneous background density n 0 and a negative vortex density n v , where 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 ξ = / √ mµ.This is typically of the order of 0.1 − 1µm but can be tuned by means of a Feshbach resonance.
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,138,139,140].Acceleration of a vortex (or an element of a threedimensional 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 γ ∼ 100 and attractive van der Waals interactions.These ripples are not contained in the simple ansatz given above in Eq. ( 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 4 He where nonlocal potentials are employed [145,146,147,148].For purely dipolar (g = 0), oblate condensates, Wilson et al. [120,122] numerically found vortex ripples for moderate trap ratios γ ∼ 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.
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 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 twodimensional vortex solitons, which may be considered as a two-dimensional bright soliton carrying a central vortex.

Vortex in a trapped dipolar Bose-Einstein condensate in the hydrodynamic regime
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 Eq. ( 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 TF (ρ, z) given in Eq. ( 24) which is an exact solution of the vortex-free Thomas-Fermi problem.Thus, the density profile reads, where, The length scale β parameterizes the size of the vortex core and is one of three variational parameters {β, R ⊥ , κ} 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 Ref. [161] and will not be repeated here.Obtaining an analytic result for the dipolar interaction energy E dd = (1/2) Φ(r)n(r)dr is more difficult.However, in the hydrodynamic regime we have β R ⊥ , implying that the contribution to Φ(r) from n v is negligible in comparison to that from n TF .Thus, to a very good approximation we can write, i.e. replace the true Φ(r) by that purely due to the unperturbed background Φ TF (r) which is known analytically and is given in Eq. ( 30).Since Φ TF (r) 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 0 of the vortex-free state, but this latter energy is also known analytically and is given as E TF in Eq. ( 29).
In this way E v , and hence Ω v , can be computed for the trapped dipolar BEC.It is found that dipolar interactions increase Ω 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 ⊥ but in the oblate case they increase it.The rotational energy density (1/2)n(ρ, z)v 2 (ρ) is lower at larger radii because v ∝ 1/ρ 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], 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 ⊥ in Eq. ( 47) is replaced by its dipolar version as given in Eq. ( 27), the resulting prediction for Ω 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.

Vortex in a quasi-two-dimensional dipolar Bose-Einstein condensate
We now review the vortex solutions in the simpler context of the homogeneous quasi-two-dimensional dipolar condensate, following the work of Refs.[129,134].Figure 8 [129] plots these solutions, found by numerical solution of the quasi-two-dimensional dipolar GPE, as a function of ε dd .
The vortex profile has a non-trivial dependence on the polarization angle and ε dd .For α = 0, the vortex density is axisymmetric.For ε 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 ξ = / √ mµ [108].Since the length scaling applied in Figure 8 is the dipolar healing length, the vortex core structure shown in Figure 8(a) for ε dd = 0 remains, for the most part, similar to that for ε 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 ε 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 0 , and their wavelength is of order the roton wavelength ≈ 4ξ.
For α = 0, see Figures 8(b,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 dd > 0 and perpendicular for C dd < 0).These anisotropic ripples are related to the anisotropic mixing of the roton into the ground state [130].

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 0 , with the density expressed using the decomposition as in Eq. ( 43).Then, by noting that the dipolar potential Φ is a linear functional of density, the total dipolar potential can be written as, (48) i.e. a contribution Φ 0 , which is just a constant, from the uniform background n 0 and a spatially dependent contribution Φ[n v ] 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.
Consider a vortex in the quasi-two-dimensional dipolar condensate.
For α = 0 and away from the roton and phonon instabilities, the vortex profile is well-approximated by the non-dipolar ansatz given by Eqs. ( 42) and (43), where the healing length is taken to be the dipolar healing length ξ = / √ mµ.The dipolar potential generated by the density defect can be determined via the convolution result Φ ⊥ (ρ, t) = F −1 Ũ ⊥ dd (q)ñ ⊥ (q, t) .Due to the cylindrical symmetry of the α = 0 case, one can perform the Fourier transforms through Hankel transforms.
The Hankel transform of the twodimensional equivalent of Eq. ( 43) is, where K 0 (•) is a modified Bessel function of the second kind.Expanding the dipolar interaction potential Ũ ⊥ dd (q) given in Eq. ( 34) with α = 0 as a series in the condensate width parameter, σ, gives Figure 9. Schematic of the decomposition of a density featuring a vortex into a uniform density n 0 and a negative vortex density nv.
Then, to first order in σ = l z /ξ and third order in 1/ρ = ξ/ρ , the dipolar potential due to a vortex is [134,129], with constants A = − 9π/8 ≈ −1.88, B = (ln 2 − 1)A ≈ 0.577, and Φ 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 Eq. ( 51) that the dipolar potential generated by the vortex relaxes predominantly as 1/ρ 2 to the background value Φ 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/ρ 2 , and the leading contribution to Φ 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/ρ 2 dependence [164].The long-range contribution to Φ v can be interpreted as arising from effective anti-dipoles in the vortex core, see Figure 9.This non-local contribution is represented in Eq. ( 51) by terms linear in σ.In the limit σ → 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 ρ /ρ 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].
The first (third) order expansion in σ (1/ρ ) allows the leading long-range dipolar contribution to Ũ ⊥ dd (q) to be evaluated (∝ σ).For dipoles tilted away from the vertical (α = 0), the vortex core and its dipolar potential become anisotropic and an analytic treatment is challenging.Figure 10(a,b) shows an example numerical solution of the vortex density and dipolar potential for α = 0.The dipolar potential is indeed anisotropic about the vortex.Remarkably, the modification to the dipolar potential induced by the vortex mimics the dipoledipole interaction itself, with an angular dependence which resembles 1−3 cos 2 θ, 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 α = 0 case, the dipolar potential is found to decay at long-range as 1/ρ 2 to the background value, as seen in Figure 10(c).At short range, ρ < ∼ 10ξ, the dipolar potential is dominated by the core structure.
6 Vortex Pairs: Interactions and Dynamics

Interaction between vortices
Two vortices (or indeed two elements of the same three-dimensional vortex line [154,155,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 ρ 1 and ρ 2 , a distance d apart.The vortices have charge q 1 and q 2 , and individual velocity fields v 1 and v 2 , respectively.The net velocity field of the two vortices is v 1 +v 2 .The energy of the vortices can be estimated by integrating the total kinetic energy across the system, For simplicity, one can ignore the vortex core density and set n(ρ) = n.Assuming ξ d R then the (kinetic) energy of the pair is, 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 0 , then the number of particles in this system can be expressed as, Now consider the system without a vortex, but with the same number of particles.It has constant density n 0 = N/A and its energy is where ) into Eq.(55) gives, where a term, negligible in the limit A ξ 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], where E 2 (ρ 1 , ρ 2 ) is the energy of the 2-vortex system.This is plotted in Figure 12 for α = 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].
Figure 11.The dipolar interaction between two vortices may be interpreted as the interaction between two collections of antidipoles.In the absence of dipoles, E 12 increases with d for the VA pair and decreases for the VV pair.For d ξ, E 12 closely follows the logarithmic scaling of the hydrodynamic prediction, while for d < ∼ ξ the overlap of the cores causes a breakdown of the logarithmic behaviour.With dipoles, E 12 is significantly modified at short and intermediate length scales up to d ≈ 5ξ, 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 nonlocal interactions.
When the dipoles are tilted in the plane, E 12 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 ξ, ∆E 12 = E 12 (η) − E 12 (η = 0) approaches a sinusoidal dependence on η, analogous to the dipolar interaction itself.

Dynamics of vortex pairs
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 intervortex axis.This solitary wave has speed v = /md for well-separated vortices.
When the VA separation is small, d ∼ ξ, 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 α = 0 the speed of sound varies with angle [130], one can expect the pair speed to be anisotropic in space.
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, α = 0.However, for α = 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, corotation 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.
7 Generation of Vortices

Summary of vortex generation methods
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,166,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 Ticknor et al [130] examined this process in a quasi-two-dimensional dipolar condensate.For dipoles tilted into the plane (α > 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 52 Cr [172].However, this shift may be more significant in 168 Er and 164 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 chargeneutral 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 mag-netic field, the potential geometry, and the orientation of the dipoles, with anisotropic interactions significantly altering vortex nucleation.

Stationary solutions of rotating dipolar condensates in elliptical traps
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,181], with a potential of the form, 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 Ω ∼ 0.3ω ⊥ .Surprisingly, in non-dipolar BEC experiments the observed nucleation of vortices occurs at considerably larger rotation frequencies Ω ∼ 0.7ω ⊥ .Theoretical analysis based on the hydrodynamic equations reveals the important role of collective modes.Specifically, for Ω 0.7ω ⊥ 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 [182,183].Evidence for this comes from comparison between experiments [181,184] and full numerical simulations of the GPE [185,186,187,188].
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 Ω where Ω = |Ω| is the rotation frequency, the Hamiltonian in the rotating frame is given by, where the Hamiltonian in the absence of rotation is H 0 and the quantum mechanical angular momentum operator is L = −i (r × ∇).Applying this result for Ω = (0, 0, Ω), the dipolar GPE in the rotating frame is [189,190], where in the rotating frame the trapping potential V , given by Eq. ( 58), is stationary.The spatial coordinates r are those of the rotating frame, with the momentum coordinates expressed in the laboratory frame [189,190,191].
Using the Madelung transform, as per Section 3.3, leads to the following dipolar hydrodynamical equations in the rotating-frame where the quantum pressure is assumed to be small and is neglected, i.e. the Thomas-Fermi limit.
Stationary solutions of Eqs. ( 61) and ( 62) satisfy the equilibrium conditions, Assuming the irrotational (∇ × v = 0) velocity field ansatz [182] permits us to examine the rotating solutions in terms of the velocity field amplitude α v .A physical interpretation of velocity field amplitude, α v , can be gained from the continuity equation ( 61).Specifically, it can be written as [109] where D is the deformation of the BEC in the x − y plane, where . . .denotes the expectation value in the stationary state and κ x = R x /R z and κ y = R y /R z represent the aspect ratios of the BEC along x and y with respect to z.
Combining Eqs. ( 62) and ( 64) gives, where the effective trap frequencies ωx and ωy are given by, 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, where n cd = 15N/(8πR x R y R z ).This is an exact solution of the stationary dipolar hydrodynamic equations given in Eqs. ( 61) and ( 62) with velocity field (64), and it only remains to find the radii {R x , R y , R z } and velocity amplitude α v .The exact dipolar potential due to a parabolic density distribution of general ellipsoidal symmetry, with dipoles aligned in the zdirection, i.e. the generalized version of that given in Eq. ( 30), is derived in the Appendices of Refs.[92] and [117], and is given by, where the coefficients β ijk are, x + s) for integer-valued i, j and k.Thus, Eq. ( 67) can be rearranged to obtain an expression for the density profile [165,166,167], .
(73) By equating the x 2 , y 2 and z 2 terms in Eq. ( 70) and Eq. ( 73) three self-consistency conditions are found.These conditions define the size and shape of the condensate, where ζ = 1 − ε dd 1 − 9κxκy 2 β 002 .Furthermore, by inserting Eq. ( 73) into Eq.( 61), the stationary solutions are seen to satisfy the condition [165,166,167], Equation (77) gives the velocity field amplitude α v for a given ε dd , Ω and trap geometry.For ε dd = 0, α v is independent of the s-wave interaction strength g and the trap ratio γ.Dipolar interactions qualitatively alter this scenario with α v becoming dependent on ε dd and γ.The solutions to Eq. ( 77), as a function of Ω, have significantly different properties depending on whether the traps are circular ( = 0) or elliptical ( > 0) in the x − y plane.We restrict our analysis to Ω < ω ⊥ , since the static solutions are known to disappear for Ω ∼ ω ⊥ due to centre of mass instabilities [182].Below, the circular ( = 0) and elliptical ( > 0) cases are considered.
7.2.1 Circular trapping in the x − y plane: = 0 In Figure 15(a) [167] the solutions of Eq. ( 77) are plotted as a function of rotation frequency Ω for a spherically-symmetric trap, γ = 1 and = 0, for various values of ε dd .For a given value of ε dd the solutions have the same qualitative structure.Specifically, only one solution exists (α v = 0) up to some critical frequency.Two additional solutions (α v > 0 and α v < 0) bifurcate from this single solution at the critical rotation frequency, denoted as the bifurcation frequency Ω b .
As expected, when ε dd = 0 the results of Refs.[182,183]  Specifically, in the Thomas-Fermi regime, the surface excitation dispersion is given by (see page 183 of [109]) where R is the radius of the BEC.For a surface excitation with angular momentum l = q l R, in the absence of rotation, Eq. ( 78) reduces to ω l = √ lω ⊥ .Rotation shifts the mode frequency by −lΩ, see Eq. ( 59).Hence for a rotating BEC the quadrupole surface excitation (l = 2) frequency is ω 2 (Ω) = √ 2ω ⊥ − 2Ω [109].The bifurcation frequency, Ω b = ω ⊥ / √ 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 Ω ≥ ω ⊥ / √ 2 the two additional solutions arise from the excitation of the quadrupole mode.
Referring back to Eq. ( 78), the inclusion of non-local dipolar interactions implies that the force −∇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 ε dd .Indeed, the dependence on ε dd can be seen in Figure 15(a) where the introduction of dipolar interactions leads to a shift in Ω b .Specifically, for ε dd > 0 (ε dd < 0) the bifurcation frequency is increased (decreased).It is possible to evaluate Ω b analytically by realising that κ x = κ y = κ, for α v = 0.In this limit the aspect ratio κ is determined by the transcendental Eq. ( 25) [91,92].As α v → 0 + , the first order corrections to κ x and κ y with respect to κ from Eqs. ( 74) and ( 75) can be calculated and inserted into Eq.( 77).Solving for Ω (noting that Ω → Ω b as which is plotted in Figure 15(b) as a function of γ for various values of ε dd .When ε dd = 0 then Ω b = ω x / √ 2, which is independent of γ = ω z /ω x [182,183].However, for ε dd = 0 the value of γ for which Ω b reaches a minimum changes.Specifically, as ε dd is increased from −0.5 the minimum value of Ω b changes from occurring at values of γ where the trap shape is oblate (γ > 1) to ones where it is prolate (γ < 1).Fixing γ and increasing ε dd leads to a monotonic decrease in Ω b .Increasing ε dd can lead to a significant reduction in the bifurcation frequency, i.e. for ε dd = 0.99, Ω b is reduced by 20% compared to its non-dipolar value to ≈ 0.55ω ⊥ .It is tempting to consider the case where ε dd > 1 to induce shifts in Ω b to even lower values, however, in these regimes the premise of the calculation, i.e. the Thomas-Fermi approximation, may not be valid.

Elliptical trapping in the x − y plane:
> 0 Rotating elliptical traps have been created experimentally with lasers and magnetic fields [35,181].Following the experiment by Madison et al. [35], below we consider a trap with a weak ellipticity of = 0.025.Figure 16(a) [167] shows the solutions to Eq. ( 77) for various values of ε dd for γ = 1.As in the non-dipolar case [182,183], the solutions become heavily modified for > 0. There is an upper branch solution (α v > 0) which extends over the range 0 ≤ Ω ≤ ω ⊥ and a lower solution (α 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 Ω b , which in the limit = 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 α v = 0 solutions for Ω > 0. In the absence of dipolar interactions increasing the trap ellipticity results in an increase of the back-bending frequency Ω b .As shown in Figure 16(b) dipolar interactions, as in the case of = 0, reduce (increase) Ω b for ε dd > 0 (ε dd < 0).
Due to the anisotropy of the dipolar interactions increasing ε dd decreases both κ x and κ 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 ε dd is increased, for Ω > 0, α v increases and hence, see Eq. ( 65), the deformation (D) of the dipolar BEC increases.

Dynamical stability of stationary solutions
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 eq + δn and S = S eq + δS.By linearizing the dipolar hydrodynamic equations Eqs.(61,62), the dynamics of such perturbations can be described as [165,166,167], ∂ ∂t where v c = v − Ω × r and the integral operator K is defined as To investigate the stability of the BEC the eigenfunctions and eigenvalues of the operator in Eq. ( 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 [192], 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 Eq. ( 80), including K [93,94,96], result in polynomials which are of degree N or less.This enables Eq. ( 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,166,167,183].
Below, the stability of the rotating solutions are considered for > 0 [Figure 16(a)].For α v < 0 there are two static solutions for Ω > Ω b .The solution nearest the α v = 0 axis is always dynamically stable except when Ω → ω ⊥ , where the BEC is susceptible to a centre-of-mass instability [193].The other solution, in the α v < 0 half-plane, is always dynamically unstable and hence experimentally irrelevant.The dynamical stability of solutions in the upper half-plane  (α v > 0) is more interesting.In Figure 17 [167] the maximum positive real eigenvalues of the upperbranch solution is plotted as a function of Ω for a fixed trap geometry and various values of ε dd .To obtain these results a maximum polynomial perturbation of δn = x p y q z r with p + q + r ≤ 3 = N was considered.
In Figure 17(inset) an example of the unstable region (for the upper-branch solution) is shown, in the − Ω plane, for ε dd = 0 and γ = 1.The shaded crescents [183] denote the regimes of dynamical instability (Re(λ) > 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 [188] (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 [188].These numerical results showed that the narrow instability regions have negligible effect when ramping Ω at rates greater than dΩ/dt = 2 × 10 −4 ω 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 0.1 and the unstable regime can be quantified with N = 3 for the perturbations.Denoting the lower bound for the dynamical instability as Ω i , Figure 17 indicates that as ε dd is increased the lower bound for Ω i decreases.
As the size of the scalar matrix operator ( 80) is increased to N = 4, 5, . .., the higher lying modes develop real eigenvalues as Ω 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 [194] considered the case of superfluid 4 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(p) → E(p)−pv.There is, therefore, a critical velocity v c = min E(p)/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 4 He has a roton minimum in its spectrum, the excitations created at v c have a specific momentum p ≈ p r , where p r is the roton minimum, triggering an instability to the formation of a density wave with wavelength ∝ p −1 r .For a rotating system similar arguments can be made with the Galilean-shifted energy taking the form E → E − LΩ 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 Ω 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 ≤ ε dd ≤ 1) this is not surprising since roton minima in dipolar BECs do not occur for −0.5 ≤ ε dd ≤ 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.

Routes to instability and vortex lattice formation
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 Ω can be utilised to trigger dynamical instability.
For non-dipolar BECs this has been realised both experimentally [181,184] and numerically [185,186,187], with excellent agreement with the hydrodynamic predictions.
For dipolar BECs the static solutions and their instability depend on the rotation frequency Ω, the trap ellipticity , the trap ratio γ and the interaction parameter ε 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 α v of Eq. ( 77) as a function of Ω [Figure 18 a dynamical instability (dashed grey arrows), or (ii) back-bending of the stationary solution such that it ceases to exist (solid grey arrows).For α v > 0, the instability always arises from dynamical instabilities in the rotating frame stationary solution.For α v < 0 the instability always arises from the rotating frame stationary solution ceasing to exist.
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 Ω v defined in Eq. (37), where in the Thomas-Fermi limit (assuming cylindrical symmetry) Ω v is given by Eq. ( 47).For typical condensate parameters Ω v ∼ 0.3ω ⊥ .This result is inconsistent with experimental observations of a much higher threshold for vortex lattice formation, Ω ∼ 0.7ω ⊥ .This discrepancy can be understood by considering the mechanism for vortices to be seeded into the BEC.For example, consider the case where Ω is increased adiabatically from zero [Figure 18(a)].Then the pertinent stationary solution is in the α v > 0 half-plane, which becomes dynamically unstable when Ω ≈ 0.7ω ⊥ .This instability provides a mechanism for vortices to be seeded into the BEC and ultimately relax into a vortex lattice.Alternatively, for Ω > Ω 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 Ω is adiabatically increased from zero, for Ω < Ω v the global energy minimum is defined by the rotating frame stationary solutions.As such as Ω passes through Ω 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 Ω ≈ 0.7ω ⊥ .
Of course, the details of how a vortex or a vortex lattice forms once a hydrodynamical instability occurs is nontrivial, however the instability is the first step in the process [185,186,187].For example, consider a BEC undergoing an adiabatic introduction of Ω (from zero), see dashed grey arrow in Figure 18(a).At a critical rotation frequency, Ω i , the BEC becomes dynamically unstable.This leads to the exponential growth of surface ripples in the BEC [35,184,187].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 ε 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 [187,195].
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 ε dd as discussed in Section 5.3 based upon the results in Ref. [161].There it was assumed that the system was cylindrically-symmetric, however, if we assume a very weak ellipticity = 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 ε dd .Initially an oblate trap is considered with γ = 10, see Figure 19(a) which plots the instability frequencies Ω i (short dashed curve) and Ω b (long dashed red curve) as a function of the dipolar interactions ε dd .For adiabatic changes in Ω (vertical path) or ε 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 Ω i ≈ Ω b ≈ 0.75ω ⊥ .Also shown in Figure 19(a) is the rotation frequency Ω 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 Ω v , with Ω v ≈ 0.1ω ⊥ 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.
For a prolate trap with γ = 0.1, Figure 19(b) plots the instability frequencies Ω i (short dashed curve) and Ω b (long dashed red curve) and Ω v (solid curve).
The instability frequencies follow a similar trend to the oblate case.However, Ω v is drastically different, increasing significantly with ε dd [161].Depending on the strength of the dipolar interactions, two regimes are predicted.For ε dd 0.8 both Ω i and Ω b are larger than Ω 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 ε dd 0.8, both Ω i and Ω b are lower than Ω 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 [196,197,198].However, the true nature of the final state warrants further investigation.
Numerical results of the dipolar GPE in the quasitwo-dimensional regime [199,200] 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.

Generalisation of Thomas-Fermi analysis
In principle the analysis presented in Sections 7.2, 7.3 and 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.

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.3 and 7.4 for dipolar BECs can be modified to include synthetic magnetic fields.Let us outline the calculation: The starting point is, where the generalized velocity is ν = v−q * A * /m, with A * 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, ∂ ∂t (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.

Off-axis dipole orientation
In Sections 7.2, 7.3 and 7.4 the analysis was restricted to the case where the dipoles were aligned perpendicular to the plane of rotation.In the absence of rotation it is possible to find Thomas-Fermi solutions when the dipoles are not aligned along one trap axis [201].Although non-trivial, it is in principle possible to generalise the work of Sapina et al. [201] to include rotation, thereby providing a framework to quantify under what regimes it is favourable for a vortex to be seeded into a dipolar BEC when the dipoles are not perpendicular to the axis of rotation † †.

Vortex Lattices
To date there have been several examinations of the properties of vortex lattices in quasi-two-dimensional rotating dipolar BECs [144,202,203,204,205,206,207].Work by Cooper et al. [202,204,205] found, by numerical minimization of the interaction energy with the wavefunction constrained to states in the Lowest Landau Level (LLL), that the dipolar interaction could modify the symmetry of the vortex lattice.Specifically, they found new phases could emerge with square, stripe (rectangular) and bubble phases, in contrast to the conventional triangular lattice structure of nondipolar BECs.These new phases emerge in the regime a s −0.13C dd m/ 2 .This result coincided with the work of Zhang and Zhai [203] who also found that the triangular phase is not favoured when the contact interactions are attractive (a s < 0), with stronger dipolar interactions leading to a square then rectangular lattice.Additionally, recent work by Kishor Kumar et al. [207], using numerical solutions of the three-dimensional purely (a s = 0) dipolar GPE, found both triangular and square vortex lattice configurations, with both the strength of the dipolar interactions and the rotation frequency determining the symmetry of the final state.Each of these investigations assumed that the axis of rotation of the BEC was the same as the alignment direction of the dipoles.Numerical simulations of the dipolar GPE, carried out by Yi and Pu [144] did not find any evidence of this change of lattice structure for such a configuration.However, for dipoles aligned off axis, they did find evidence of a change in the symmetry of † † To be able to carry out such a calculation, along the lines presented in Sections 7.2, 7.3 and 7.4, one must consider the case where the dipole alignment is stationary in the rotating frame, i.e. rotating with the trap.the vortex lattice.This work also found evidence for ripply vortex lattices, where density modulations arise in the vicinity of a vortex.This can be understood in the context of the analysis presented in Section 5 where we deduced that as the BEC approaches the roton instability density ripples appear around a single vortex.This is consistent with numerical work presented by Jona-Lasinio et al. [124].
Below we present new results for the structure of a vortex lattice in a dipolar BEC in the quasi-twodimensional regime.We consider both the case where the dipole alignment is perpendicular to the condensate plane (in reference to Figure 4 α = 0), and generalise this to include configurations where the alignment has a component into the plane (in reference to Figure 4 α > 0).The treatment considered draws on the works of Cooper et al. [202,205], Zhang and Zhai [203] and [208,209].For clarity, we shall begin our analysis by considering vortex lattices for the case of a nondipolar BEC.Although this case has been covered rather extensively in the literature, it will be useful for the reader to provide a comprehensive, self-contained review here.We shall then build on the methodology presented to consider dipolar BECs.

Vortex lattice in a non-dipolar BEC
Assuming that a condensate is confined in a cylindrically-symmetric harmonic trap, the energy functional in the rotating frame is given by, where, (86) There are two distinct contributions to the total energy: a single-particle energy contribution, E 0 (the first term in Eq. ( 85)) and an interaction energy contribution, which in the absence of dipolar interactions is E vdW (the second term in Eq. ( 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 ( ω z 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 Eq. (31).
The resulting form for the energy functional, in the absence of dipolar interactions is, where, Consider the fast-rotating limit where Ω → ω ⊥ .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 singleparticle Hamiltonian in quasi-two dimensions tends to , up to a constant.The eigenfunctions of this Hamiltonian are well known: they are the Landau level orbitals u m ,n (x, y) with corresponding eigenenergies n = Ω(n + 1/2).The n quantum number labels the Landau level, and may take on any non-negative integer value.Each Landau level is infinitely degenerate since n does not depend on m (which also takes on non-negative integer values).
Although the Landau level orbitals do not necessarily form a complete basis for the full quasi-twodimensional Hamiltonian, in the limit of weak interactions they are a good approximate basis choice [196].
In searching for the ground state of the condensate it is assumed that it is adequately described by a superposition of n = 0 Landau level orbitals, i.e. the lowest Landau level (LLL) approximation.Using an unrestricted minimization this assumption implies that [52], In this limit the two-dimensional ground state wavefunction of the condensate can be written as, where the right-hand side follows from the explicit form of u m ,0 (x, y).The length scale l ⊥ = /mΩ characterises the radial extent of the condensate and is effectively equal to the transverse trap length since Ω → ω ⊥ .
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 ρ = xx + yŷ is mapped onto the complex number w = x + iy.Then the ground state wavefunction of Eq. ( 90) may be written as, where h is an analytic function of w.With this definition, the coefficients of the superposition c m 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 [209].
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 wellstudied function in complex analysis which will be useful here: the Jacobi theta function θ 1 (z, ζ).The roots of θ 1 (z, ζ) are able to describe any regular lattice up to a rotation by making an appropriate choice for the parameter ζ.This means that h may be expressed in terms of θ 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: b 1 and b 2 .The points of the lattice are obtained from these basis vectors by constructing all possible linear combinations of the form m 1 b 1 + m 2 b 2 where m 1 and m 2 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 b 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 b 1 is specified solely in terms of its magnitude, which is denoted by the parameter b 1 , since its direction is fixed along x.The second basis vector b 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 1 , τ and η, results in b 1 = b 1 x and b 2 = τ b 1 (cos ηx + sin ηŷ).
The other parameter that appears in Figure 20 is the area of the unit cell, given by v c = b 2 1 τ sin η.This parameter is redundant in specifying the lattice because it depends on the other parameters: b 1 , τ 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 q 1 and q 2 , we have 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 q 1 + m 2 q 2 where m 1 and m 2 are integers.Each one of these linear combinations is a reciprocal lattice vector, which is denoted in general by q.
From this, it is possible to describe the two-dimensional density as, where g(ρ) = q gq exp(iq • ρ) and, with q = m 1 q 1 + m 2 q 2 .In other words, the two-dimensional condensate density is the product of a Gaussian envelope and a function g(ρ) = q gq exp(iq • ρ) which is periodic on the vortex lattice.The Fourier coefficients of g(ρ) are rescaled compared to those in Eq. ( 93) so that gq = g q / v gv exp(−iv 2 χ 2 /4).The radial extent of the condensate cloud is quantified by the length-scale χ = (l −2 ⊥ − πv −1 c ) −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 ψ ⊥ (ρ) 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, Here dρ is a dimensionless analogue of the interaction energy contribution.Substituting the ansatz for the vortex lattice ground state from Eq. ( 92) gives, This quantity depends on τ and η implicitly through its dependence on the reciprocal lattice vectors.
Considering the limit of large vortex number where χ 2 q 2 1, it is possible to simplify this expression significantly.
In fact, one may assume that the is so sharply peaked at q = −v that it may be approximated by a Kronecker delta, δ q,−v .This collapses the double sum to a single sum over q, which is much simpler to evaluate, In order to minimise I(τ, η) as defined in Eq. ( 96) a numerical treatment is required, which means cutting off the sums over m 1 and m 2 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 I(τ, η) 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 I(τ, η) graphically, by generating a contour plot of I(τ, η) as a function of τ and η.The result is shown in Figure 21.In this plot, the dark purple shading corresponds to regions where I(τ, η) approaches its minimum value.The white areas correspond to regions where I(τ, η) 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 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.
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 τ = 1 and η = π/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.
Figure 22.The Bravais lattice may be classified as one of five types according to the geometry of the unit cell.The three pertinent Bravais lattice structures considered are triangular, square and rectangular lattices which are all special cases of the parallelogrammic lattice.For the triangular, square and rectangular lattice, we have given the corresponding (τ ,η) values in the so-called standard parametrisation.

Vortex lattice in a quantum ferrofluid: dipoles perpendicular to the plane of rotation
For the case of dipolar BECs the dipolar interaction potential, U dd (r), 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 where σ = l z /R ⊥ .The above has been obtained by considering the expansion of the Ũ ⊥ dd (q) in terms of the width of the BEC to first order, i.e.Eq. ( 50), and then calculating the dipolar potential to second order in ρ/R ⊥ .
The additional contribution to the energy functional arises due to the dipolar interactions, Adding this to the single-particle and contact interaction energy contributions gives the total energy in the rotating frame for a dipolar BEC, The results obtained in Section 8.1 for the singleparticle (E 0 ) and contact interaction (E vdW ) energy contributions are unchanged for the dipolar case.All that remains then, is to perform the calculation for E dd .
Rewriting the dipolar interaction energy contribution in reciprocal space by applying the Fourier convolution theorem leads to, To evaluate this in quasi-two dimensions, we substitute the Fourier transform of the quasi-two-dimensional condensate density, ñ(k) = e −k 2 z l 2 z /4 ñ⊥ (q), such that where Ũ ⊥ dd (q) is given by Eq. ( 34), with α = 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 Eq. ( 92) results in, ñ(q) = N q gq e −χ 2 |q −q| 2 /4 . ( This enters into the expression for the dipolar interaction energy, leading to, gq gv e −χ 2 (q 2 +v 2 )/4 A(v − q), (103) where, and I 0 (•) is the modified Bessel function of the first kind.At this point, the integral may be separated into two terms A 1 (q) + A 2 (q).The first term has a simple solution, Returning to the expression for the dipolar interaction energy given in Eq. ( 103), and substituting the simplified result for A 1 (q) gives where I(τ, η) is as defined in Eq. ( 95) and, W(τ, η) = q,v gq gv e −χ 2 (q 2 +v) 2 /4 A 2 (v − q). ( This expression shows that the dipolar energy has been separated into two distinct contributions: a local contribution which is proportional to I(τ, η) and a non-local contribution which is proportional to W(τ, η).In principle, it is now possible to calculate the dipolar interaction energy as a function of the lattice parameters, however the function for the nonlocal contribution, W(τ, η) 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, 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 Eq. ( 107) to a single sum in the limit of large vortex number, resulting in, where, where β = χ/l z and L n (•) is the n th Laguerre polynomial.
With the results of Section 8.1 it is now feasible, from a computational perspective, to numerically minimise the condensate energy, E(τ, η) = E 0 + E int (τ, η), 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, since the single-particle contribution does not depend on τ and η.Assuming C dd > 0 so that the factor outside the square brackets is positive the vortex lattice configuration is determined by the the minimization of 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 ε −1 dd .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.
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 ε −1 dd = −1.82 in Figure 23(a), a square lattice is favourable since the optimal values of τ and η at that point are 1 and π/2 respectively.Continuing in this way, three types of lattice emerge depending on the value of ε −1 dd : 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 ε −1 dd ≈ −2 indicates that the transition arises in the regime where the long-range contribution to the interactions [W(τ, η)] dominates, see Eq. (113).Comparing Figures 23(a), with l ⊥ /l z = 40, and (b), with l ⊥ /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.
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.

Phase
Figure 23(a) Figure 23 Table 1.Definition of the four regions in phase space shown in Figure 23.
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 ⊥ ) = 1.0191, l ⊥ /lz = 80 and χ/lz = 584.66.In each plot, the black solid (dashed) line represents the optimal value of τ (η) with reference to the scale on the left (right) vertical axis.By classifying the (τ ,η) parameters, four distinct regions in phase space are identified: a collapse phase, rectangular lattice (stripe) phase, square lattice phase, and triangular lattice phase.These regions are indicated by distinct shading.The red line with arrows on the left side specifies the region for which the interaction energy is less than zero.
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 ε −1 dd 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 ε −1 dd and the values of the length scales l ⊥ , l z and χ.In general, a triangular lattice geometry is favoured in regions where the local interaction contribution dominates, as was seen in the nondipolar case.However, in regions where the nonlocal interaction contribution becomes significant, the favoured lattice geometry changes from triangular to square or rectangular.Below a certain value of ε −1 dd (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 [203] and Cooper et al. [202,204].Zhang and Zhai also find that the lattice geometry undergoes a transition from triangular → square → rectangular → collapse as the value of ε −1 dd 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 [203] and Cooper et al. [202,204].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.

Vortex lattice in a quantum ferrofluid: dipoles not perpendicular to the plane of rotation
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 xdirection.If the density profile of the cloud is expressed in the form exp[−x 2 /l 2 x − y 2 /l 2 y ], then the aspect ratio is written as λ = l y /l x .For a cylindrically symmetric BEC, the width of the cloud along the x-and ydirections must be the same, which implies that λ = 1.For α > 0 [see Figure 4], since dipolar BECs elongate along the direction of polarisation, it is expected that l x > l y and hence 0 < λ < 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, b 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 noncylindrically symmetric dipolar BECs (α > 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, where R represents the standard two-dimensional rotation operator, χ = [(l x l y ) −1 − πv −1 c ] −1/2 and, From this starting point is possible to generalize the approach presented in Section 8.2 to the case where α > 0. The total energy can be broken into two components, The non-interacting component of the energy is given by, which in the limit Ω → ω ⊥ tends to a constant.Following the same procedure as in Section 8.2, using Eq. ( 92), the interaction energy is given by, where q  lxly/lz = 40 and χ/lz = 292.33.The optimal value of ζ is not included because Ẽ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.
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: (ε −1 dd = 0.9, α = π/2) and (ε −1 dd = 0.95, α = π/2) we find evidence for a transition.In the limit Ω → ω ⊥ (E 0 (λ) → N [ Ω + ω z /2]) Table 2 shows E int (τ, η, λ, ζ) 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 ε dd .It is also found that the variational parameter ζ is essentially irrelevant, since the minimum value of E int (τ, η, λ, ζ) 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.
Numerical studies [144], based on the dipolar GPE, suggest that for α > 0 the vortex lattice can undergo a phase transition from triangular structure to a non-triangular structure as ε 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 α > 0.

Vortex lattices in two-component dipolar Bose-Einstein condensates
The theoretical study of non-dipolar two-component BECs [209,210] has shown how interspecies interactions g 12 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 [209] were able to quantify various regimes through the parameter β = g 12 / √ g 1 g 2 , where g 12 is the interspecies interaction parameter and g 1( 2) is the intraspecies interaction parameter for component 1 (2) of the two-component BEC.This treatment and subsequent numerical analysis [210] shows that for β < 0 the two components overlap and a single triangular vortex lattice arises.For β > 0 the two components separate to form interlaced triangular vortex lattices.As β ∼ 1 the triangular lattice distorts to form square or rectangular arrays.This theoretical analysis is consistent with experiment [211].
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. [212] 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 [213]), 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 [214,215,216] and has been extended to consider how dipole alignment in the plane of rotation [216] and component-dependent optical lattices [217] influences the phase diagram.9 Summary and Outlook

Summary
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 longrange and anisotropic component to the interactions, making a significant departure from conventional swave 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 quasitwo-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 longrange and anisotropic contribution to the vortexvortex interaction.This can, for instance, lead to the suppression of the annihilation of vortex-antivortex pairs and induce the co-rotational dynamics of vortexvortex 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.

Dipolar Bose-Einstein condensates in toroidal traps
The experimental study of persistent superfluid flow in BECs confined to toroidal traps [218,219,220,221,222,223,224,225,226,227,228,229] has matured significantly over the last decade.As such, ring shaped BECs in toroidal traps have been the subject of many experimental and theoretical investigations [230,231,232,233,234,235] focusing on persistent currents [218,221,236], weak links [219,222], formation of matterwave patterns by rotating potentials [237], solitary waves [232,238], and the decay of the persistent current via phase slips [220,239,240].In these studies the transference of angular momentum from optical fields [218,222] or stirring with an optical potential [222,223] is used to generate persistent flow.
Within the context of this review we consider a persistent flow in a toroidal BEC as a giant vortex state.One might suppose that beyond studying the density profile and stability of a dipolar condensate within a toroidal trap [241,242] that dipolar interactions have a limited influence on the superflow properties in such a geometry.This assumption arises since to a close approximation the wavefunction can be considered to have the form ψ where θ is the azimuthal angle and q v is the vortex charge (charaterising the persistent flow) around the toroid, and the density is independent of q v .As such, when considering the energy difference between q v and q v + 1 the dipolar interactions play no role.
Despite this observation work has been carried out on the properties of dipolar BECs in toroidal traps focusing on the generation of persistent flows via the He-McKellar-Wilkens or Aharonov-Casher effect [243] and the properties of two-component dipolar BECs in concentrically coupled toroidal traps [244].In the former case it was shown that for atomic dipolar BECs, where the dipolar interaction is mediated via a magnetic dipole moment, that although it is theoretically possible to induce persistent flow, via the Aharonov-Casher effect [245,246,247], the strength of electric field required is prohibitive.In the case of polar molecules, with significant electric dipole moments, the He-McKellar-Wilkens effect [248,249,250] could ultimately be used to generate a persistent flow in a dipolar BEC in a toroidal geometry.For the case of a two-component dipolar BECs [244] in concentrically coupled toroidal traps, various vortex structures can arise depending on the strength of the dipolar interactions and the rotation frequency.The interesting vortex structures predicted include polygonal vortex clusters and vortex necklaces.

Fractional quantum Hall physics in dipolar Bose-Einstein condensates
In Section 8 we considered the LLL regime to investigate vortex lattice structures in dipolar BECs.Ultimately, in the limit Ω → ω ⊥ , a rotating BEC is predicted to make a quantum phase transition to a highly correlated, non-superfluid, fractional quantum Hall groundstate.This state emerges when the LLL meanfield vortex lattice melts, i.e. when the number of vortices (N v ) in the BEC is the same as or greater than the number of atoms (N).This regime occurs approximately when Ω/ω ⊥ ∼ 0.999.In the absence of dipolar interactions theoretical evidence for such a transition has come from exact two-dimensional groundstate calculations for a small number of bosons with a large angular momentum [196,251,252].For a review of quantum Hall physics in rotating BECs see Ref. [253].
The natural question which arises is: do dipolar interactions influence this phase transition and the properties of the highly correlated state?Work by Rezayi et al. [254] showed that at a filling factor of N/N v = 3/2, with N = 18, dipolar interactions support an incompressible fluid ground state which possesses non-Abelian statistics for the quasiparticle excitations.
Additionally, dipolar interactions in lattice systems have been shown to increase the gap between the ground state and the first excited state [255], for N/N v = 1/2.
It is expected that this increase in the gap will be maintained in the thermodynamic limit and hence relevant for experiments.Numerical simulations by Chung and Jolicoeur [256] showed that at N/N v = 1 the ground state is a Moore-Read paired state, as is the case for bosons with purely contact interactions.This state is destabilized when the contact interactions are small enough, i.e. dipolar interactions alone can not support this state.For N/N v = 1/3 a composite fermion sea emerges, where each boson is bound with three vortices.The robustness of fractional quantum Hall states in artificial gauge fields, in the presence of dipolar interactions, has been investigated by Graß et al. [257].

Dipolar fermions
There is considerable interest in the properties of dipolar Fermi gas systems.The partially attractive nature of the dipolar interaction in single-component dipolar Fermi gases opens up the possibility of BCS pairing resulting in superfluid states in three dimensions [22,23,24,28,31] and two dimensions [25,26,29] at sufficiently low temperatures.In these systems, the anisotropy of the superfluid order parameter causes a major change in properties as compared to the case of a two-component BCS superfluid, dominated by van der Waals interactions, where the superfluid order parameter is isotropic (swave).It is expected that the anisotropic gap will lead to significant differences in the properties of single and multiple vortex states in these systems.For example, Levinson et al. [29] have proposed a scheme to construct a topological p x + ip y superfluid phase in a quasi-two-dimensional single component dipolar Fermi gas which can support vortices which carry zero energy Majorana modes on their cores [258,259,260].However, to date, there has been very little research into the properties of vortices and vortex lattices of such states.Assuming that such a state can be experimentally achieved there is a significant opportunity to revisit much of what has been discussed in this review within the context of vortices and vortex lattices in the BCS state of a dipolar Fermi gas.
There have been extensive studies of the properties of rotating single-component fermionic quantum ferrofluids [261,262,263,264,265,266] away from the BCS superfluid transition.These studies have primarily focused on the emergence of fractional quantum Hall states in the Ω → ω ⊥ regime.Specifically, it has been shown [261,262] that for a filling fraction of ν = 2πl 2 ⊥ n 0 = 1/3 the many-body state is well described by the Laughlin wave function with a significant gap between the ground and the excited states.Further studies [262,263,265] have shown that as the filling fraction is reduced further (ν < 1/7) Wigner crystal [267] states may emerge.To our knowledge, these studies have all been carried out in the regime where the plane of rotation is perpendicular to the orientation of the dipoles.As such an interesting question to ask may be how do such states change when the dipole orientation has a component in the plane of rotation.

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 [268,269], but even this is destroyed when the temperature is raised through the Berezinskii-Kosterlitz-Thouless transition [270,271] 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 [272], vortices [273] and changes in the density profile due to the onset of superfluidity [274].The Berezinskii-Kosterlitz-Thouless transition in the presence of dipolar interactions has been studied using, for a homogeneous system, Monte Carlo methods [275], the mean-field Hartree-Fock-Bogoliubov-Popov model [276], the dipolar XY-model [277] 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 vortexantivortex pairs in a uniform two-dimensional nondipolar BEC (of two-dimensional background density n 0 ), the energy of a pair, separated by a distance d and calculated from hydrodynamical arguments, is approximately given by V (d) = 2π 2 n 0 ln(d/ξ)/m.The critical temperature associated with the transition is given by the relation 2π 2 n 0 /(mk B T ) = 4.This is calculated by determining the average distance between the pairs, where we have used the short-hand notation Γ = 2π 2 n 0 /(mk B T ).This result diverges at Γ = 4, signifying the Berenzinksii-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 vortexantivortex 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) → V (d, η) [see Figure 13(b)].The above generalisations can be incorporated into Eq.( 122) by considering the following analysis [278].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 ∝ (0, sin α, cos α).If the separation between the vortices is d = d(cos η, sin η, 0) then the interaction between the vortex and the anti-vortex can be written as For simplicity, we have assumed that the interaction between vortices arising from the absence of dipoles in their cores is ∝ 1/d 3 .Given the results in Section 5. where When λ = 0 the result obtained from Eq. ( 122) is regained, i.e. d 2 = ξ 2 (Γ − 2)/(Γ − 4), with the Berezinskii-Kosterlitz-Thouless transition arising from the divergence at Γ = 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 Eq. (125), diverges first at Γ = 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 [275,277] 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 [277] 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 finitesize 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.

Vortex lattices in the supersolid phase
A supersolid phase is characterised by the spontaneous formation of a periodic structure in a system which also supports superfluidity, i.e. both on-and off-diagonal long-range order in the single particle density matrix.For more than half a century considerable effort has been expended theorizing about the possibility of supersolidity in condensed matter systems [279,280,281,282,283].In particular, the investigation of a supersolid phase in 4 He [284,285,286] has been the primary focus of the more recent work.However, the most credible claim for the experimental observation of a supersolid phase in 4 He [287] has now been withdrawn [288].The relatively recent realization of dipolar cold atoms and molecules now provides an alternative venue for investigating supersolid phases.Apart from the strongly correlated dipolar systems [70,72,73,74,75,76] cited in Section 2, extended Bose-Hubbard lattice models are also thought to support supersolidity [289,290,291,292,293,294,295,296,297,298,299].The common ingredient in both sets of investigations is the presence of long-range interactions.There have also been various works focusing on how an artificial gauge field in a dipolar BEC influences the boundaries between Mott-insulator and supersolid regimes [300,301], and how staggered fluxes lead to supersolid phases with staggered vortex phases [302].However, there does not appear to have been much investigation into the structural properties of vortices [303] and vortex lattices in the supersolid state.

The Onsager vortex phase transition
By studying a two-dimensional point vortex model Onsager predicted that negative temperature states may be relevant for two-dimensional fluids [304].While intended as a model of two-dimensional fluids in general, Onsager noted that the model was potentially particularly relevant for two-dimensional superfluids, whose vortices have quantized circulation and uniform size.Simulations of the two-dimensional GPE have shown how it is possible to dynamically evolve to the negative temperature Onsager vortex state [305,306].Starting with a random configuration of vortices and antivortices one might expect the vortices and antivortices to evaporate from the BEC via pairwise annihilation and re-thermalization of the emitted sound, simply resulting in a BEC with an increased temperature.However, Simula et al. [305] showed that only some vortices annihilate, and the remaining vortices self-organise into two ordered clusters of like-sign circulation, which represent the Onsager vortices.This outcome was found to be the result of the evaporative heating of quantized vortices via vortex-antivortex pair annihilation, leaving the remaining vortices to re-thermalize to a state with higher energy per vortex.This process drives the vortex component of the superfluid to ever higher energies, leading to the Onsager vortex phase transition.The final Onsager vortex state emerges as a clustering of vortices and antivortices.The dynamical process which leads to this final state is non-trivial, but is underpinned by the interaction between the constituent vortices and antivortices in the system.Introducing dipolar interactions, as shown in Section 6, can significantly modify this interaction and as such may result in significant changes to the Onsager vortex phase transition.As such we suspect that dipolar interactions may offer the opportunity to significantly modify the Onsager vortex phase transition.

Figure 1 .
Figure 1.(a) For a prolate trapped condensate the net dipolar interaction is attractive.(b) For an oblate system the net dipolar interaction is repulsive.

Figure 2 .
Figure 2. Illustration of the dipole-dipole interaction.The magic angle, at which the dipole-dipole interaction reduces to zero, is indicated by the black dashed line.

Figure 3 .
Figure 3. Top: The aspect ratio, κ, (solid curves) of harmonically trapped, cylindrically-symmetric, dipolar condensates in the Thomas-Fermi regime.Each line corresponds to an equally spaced (on a logarithmic scale) trap aspect ratio, γ (γ = [0.1,10]).White, light grey and dark grey shading correspond to regimes of global, metastable and unstable solutions respectively.Reprinted figure with permission from [117].Copyright 2010 by the American Physical Society.Bottom: Stability diagram of the purely dipolar harmonically-trapped condensate (ground state), as a function of the trap aspect ratio λ ≡ γ = ωz/ω ⊥ and the dipolar interaction parameter is D = N mC dd /(4π 2 l ⊥ ).The shaded region denotes stability against collapse.The dark shaded regions indicate biconcave condensates.Reprinted figure with permission from [119].Copyright 2007 by the American Physical Society.

Figure 4 .
Figure 4. Schematic of the quasi-two-dimensional dipolar condensate, with strong harmonic trapping along z.The dipoles are taken to be polarized at angle α to the z-axis in the x − z plane.The condensate is assumed to follow the static ground harmonic oscillator state along z. Figure reproduced from Ref. [129] under a CC BY licence.

Figure 5 .
Figure 5.The stability diagram in ε dd − α space, for (a) g > 0 and (b) g < 0, of a homogeneous dipolar BEC in the quasitwo-dimensional regime (σ = 0.5).Shown are the regions of stability (white), phonon instability (pink) and roton instability (blue).The vertical dashed line indicates the magic angle αm.For α > π/2 the results are the mirror image of the presented region.Figure adapted from Ref. [129].

Figure 6 .
Figure 6.Schematic of a three-dimensional density (red isosurface plot) of a trapped non-dipolar condensate featuring a vortex line along the z-axis.The corresponding two-dimensional phase profile (grey scale plot at the base of the figure, with white corresponding to a condensate phase of 0 and black corresponding to a condensate phase of 2π) and central onedimensional density profile (solid black curve) are also depicted.

Figure 7 .
Figure 7. Top: Stability diagram of the trapped purely dipolar condensate, with dipoles polarized along the z-axis and featuring an axial vortex, as a function of the trap ratio (λ ≡ γ = ωz/ω ⊥ ) and dipolar interaction strength [D = N mC dd /(4π 2 l ⊥ )].Below the solid line the condensate is dynamically stable.Ripples about the vortex arise in the pink region.The inset shows an isosurface of the density for such a solution.Reprinted figures with permission from [120].Copyright 2009 by the American Physical Society.Bottom: Density (a-b) and phase (c) profiles of a trapped quasi-two-dimensional condensate with dipoles polarized along x.Reprinted figures with permission from [144].Copyright 2006 by the American Physical Society.

Figure 8 .
Figure 8. Vortex solutions in an infinite dipolar condensate, as a function of ε dd , in the quasi-two-dimensional regime (σ = 0.5).Along x (for y = 0) the normalised density profile (n/n 0 ) is shown on the left-hand side of the main plots.Along the right-hand side the normalised density profile is plotted along y (with x = 0).(a) Dipoles polarized along z (α = 0).(b) Dipoles polarized off-axis at α = π/4.(c) Dipoles polarized off-axis at α = π/2.In each of the main plots grey bands indicate the unstable regimes of ε dd .Insets: normalised density profile over an area (40ξ) 2 for indicated values of ε dd .Figure reproduced from Ref. [129] under a CC BY licence.

Figure 16 .
Figure 16.(a)Irrotational velocity field amplitude αv of the rotating frame stationary solutions as a function of the rotation frequency of the trap, Ω, with γ = 1 and = 0.025 and ε dd = −0.49,0, 0.5 and 0.99, with increasing ε dd denoted by the arrow.Insets illustrate the deformation of the condensate in the x − yplane, for αv > 0 and αv < 0. (b) Back-bending frequency Ω b versus ε dd for = 0.025 and γ = 0.5 (solid curve), 1.0 (long dashed curve) and 2.0 (short dashed curve).Reprinted figure with permission from[167].Copyright 2009 by the American Physical Society.

Figure 17 .
Figure 17.The maximum positive real eigenvalues of Eq. (80) (solid curves) for the upper-branch solutions of αv as a function of Ω, for = 0.025, γ = 1, N = 3, ε dd = −0.49,0, 0.5 and 0.99, with ε dd increasing in the direction of the arrow.The inset shows the region of dynamical instability in the − Ω plane for ε dd = 0.The narrow regions, around Ω/ω ⊥ < 0.6 with < 0.1 and Ω/ω ⊥ < 0.56 with < 0.25, have negligible effect and so only the main instability region is considered (bounded by the dashed line).Reprinted figure with permission from [167].Copyright 2009 by the American Physical Society.

Figure 19 .
Figure 19.The relation between the instability frequencies, Ω b (long dashed red curve) and Ω i (short dashed curve), and the critical rotation frequency for vorticity Ωv (solid curve) for (a) an oblate trap γ = 10 and (b) a prolate trap γ = 0.1.The instability frequencies are based on a trap with ellipticity = 0.025 while Ωv is obtained from Eq. (47) under the assumption of a 52 Cr BEC with 150, 000 atoms and scattering length as = 5.1nm in a circularly symmetric trap with ω ⊥ = 2π × 200Hz= 200rad/s.Reprinted figure with permission from [167].Copyright 2009 by the American Physical Society.

1 Figure 20 .
Figure 20.Illustration of the Bravais lattice basis vectors and the lattice parameters.

Figure 21 .
Figure 21.Contour plot of I(τ, η), a dimensionless analogue of the interaction energy.The white areas are regions where I(τ, η) is greater than the cut-off value of 0.6.The purple areas correspond to regions where I(τ, η) approaches its minimum value.Multiple local minima are visible, although they become difficult to see in the lower left-hand region of the plot.The three minima labelled by white numbers are mentioned in the text.

Figure 23 .
Figure 23.Plots showing the optimal values of τ and η in dipolar BECs with on-axis polarisation.Plot (a) assumes the length scales vc/(πl 2 ⊥ ) = 1.0191, l ⊥ /lz = 40 and χ/lz = 292.33(b) assumes the length scales vc/(πl 2⊥ ) = 1.0191, l ⊥ /lz = 80 and χ/lz = 584.66.In each plot, the black solid (dashed) line represents the optimal value of τ (η) with reference to the scale on the left (right) vertical axis.By classifying the (τ ,η) parameters, four distinct regions in phase space are identified: a collapse phase, rectangular lattice (stripe) phase, square lattice phase, and triangular lattice phase.These regions are indicated by distinct shading.The red line with arrows on the left side specifies the region for which the interaction energy is less than zero.