Vortex Structures in Ultra-Cold Atomic Gases

Abstract

In this chapter a basic introduction to the theory of vortices in ultra-cold (superfluid) atomic gases is given. The main focus will be on bosonic atomic gases, since these contain the same basic physics, but with simpler formulas. Towards the end of the chapter, the difference between bosonic and fermionic atomic gases is discussed. This discussion will allow the reader to make the conceptual step from bosonic to fermionic gases, while pinpointing the main differences and difficulties when working with fermionic gases rather than bosonic gases. The goal of this chapter is to provide a good and general starting point for researchers, or other interested parties, who wish to start exploring the physics of ultra-cold gases.

1. Introduction

Superfluidity was discovered independently by Kapitza [1] and by Allen and Misener [2] in 1938, who found that liquid helium II flows without resistance through capillaries. The theoretical framework to describe superfluids was set up in the next decades, and a central concept in this description is the macroscopic (complex-valued) wave function. A normal fluid, when described at the quantum mechanical level, requires a many-body wave function that depends on the coordinates of all the particles in the fluid. A superfluid, however, can be well described by a macroscopic wave function that depends on a single-position coordinate. This wave function has a straightforward ‘hydrodynamic’ interpretation: the density of the superfluid is given by the modulus square of the macroscopic wave function, and the velocity field $v(r\mathrm{,}t)$ is proportional to the gradient of the phase $S(r\mathrm{,}t)$ of the wave function,

$v(r\mathrm{,}t)=\frac{\hslash }{m}\nabla S(r\mathrm{,}t),$ E1

with m the mass of the atoms. The single-valuedness of the wave function, in conjunction with Eq. (1) leads to remarkable properties for rotating superfluids. Since the phase can wind only in multiples of 2π, circulation is quantized. Contrary to vortices in classical fluids, in a quantum fluid, vortices carry an integer number of circulation quanta—it is not possible to divide a singly-quantized vortex into smaller units. Upon rotation, a superfluid will not behave as a solid body, but will generate a number of singly-quantized vortices that together carry the total circulation of the rotating fluid. In a (neutral) superfluid, vortices with the same sign of the circulation repel due to kinetic energy of the super-flow. This leads to the formation of a lattice of singly quantized vortices, known as an Abrikosov lattice [3].

In helium-II, the quantization of vorticity and the presence of a lattice of singly-quantized vortices were observed by Vinen [4]. His seminal experiments initiated the study of vortices in liquid helium, which still represents a very active and intense field of research. A thorough review of this field can be found in reference [5]. Liquid helium at non-zero temperature is not a ‘pure’ superfluid in the sense that thermal and quantum fluctuations lead to a finite viscosity in the bulk liquid, and this complicates the study of superfluid vortices. However, two decades ago, a new superfluid system became experimentally available: ultra-cold atomic quantum gases. These consist of a cloud of typically 10⁵ atoms that are optically and/or magnetically trapped, and cooled down to nano-kelvin temperatures. The system is very controllable and versatile: temperature, atom number, trapping geometry and dimensionality, and interaction strength are experimentally tunable. Bosonic atoms were cooled down below the critical temperature for Bose-Einstein condensation (BEC) in 1995 (for a review, see reference [6]), and fermionic atoms were cooled down to a regime where the atoms form superfluid pairs in 2003 (reviewed in reference [7]). The tunable interaction strength allows to bring Fermi gases in the regime where they form Cooper pairs, or in the regime where they form strongly bound molecules that Bose condense, or in the crossover regime between those two limits.

The present chapter aims to familiarize researchers in other fields of vortex physics with the ongoing research on vortices specifically in quantum gases, especially Bose-Einstein condensates. First, we present an (non-exhaustive) overview on the experiments related to these superfluid vortices, and in the second part, we outline the theoretical description of vortices in condensates.

2. Review of experiments on vortices in quantum gases

In the first few years after the initial observation of Bose-Einstein condensation in atomic gases [8–10], it was not clear how to form vortices in these new superfluids. Early experiments with rotating traps failed to produce the vortices (see the contribution of Ketterle and co-workers in reference [6]), and the question was raised on whether the vortex core might be too small to image. The healing length (and core size) is estimated to be of the order of 100 nm. To make observations of quantum gases, the mainstream method [6] is to shine a laser on the cloud, and image the shadow of that cloud. That way, a map of the column density of the cloud can be made, and vortices with their core parallel to the line of sight should appear as holes in the cloud. In this imaging technique, the trap holding the quantum gas is switched off, so the cloud—initially a micron-sized object—can expand to a size large enough to provide a good resolution. The initial difficulties to observe vortices were surprising since theory predicts that when the trap is switched off the vortex core expands faster than the radius of the atomic cloud [11, 12], and this should facilitate the observation, unless the vortex line twists away from the line of sight. An interferometric technique was proposed [13, 14] to overcome possible difficulties with direct imaging. The phase coherence of the superfluid implies matter-wave interference: when two adjacent condensates are allowed to expand, then a series of parallel matter interference fringes appear in the region of overlap of the expanding clouds. The phase pattern of a cloud containing a vortex is different, and this leads to a ‘fork dislocation’ in the fringe pattern [13]. This idea of interferometry is used in many different branches of physics, e.g. plasma physics to measure the electron density [15, 16] and in the study of optical vortices [17, 18], where the ‘fork dislocation’ is also found. In both cases, this is done by making an interference pattern of the ‘probing beam’ (which probes your object) with a ‘reference beam’. However, note that the interference pattern of two overlapping condensates is a pure quantum mechanical effect [19], while the interference pattern in conventional interferometry is a pure classical interference effect.

The first successful attempt, in 1999, to create a topological excitation similar to a vortex was realized not by rotating the cloud, but by quantum state manipulation [20]. This requires a multi-component condensate, consisting of atoms in different internal states, typically different hyperfine levels, confined to the same-atom trap. In the experiment given in reference [20], initially the trap contains only a non-rotating one-component BEC. A microwave field is then used to couple two hyperfine states and coherently transform a given number of atoms from the first component into the second, so that a BEC can be grown in the second component. An additional laser beam can be used to imprint a phase profile in the second component in such a way that a single quantum of circulation is present in the second-component, ring-like BEC. The core of the vortex-like structure is filled by the atoms still in the (non-rotating) first component. It can be argued that due to the multi-component nature of the condensate, it should be interpreted as a spinor condensate and the resulting object should consequently be classified rather as a skyrmion than as a ‘true’ vortex. Indeed, in these cases, the phase can ‘untwist’ itself, and the excitation is not topologically protected as expected for a quantum vortex [21]. In later experiments using quantum-state manipulation, vortices were created by transferring orbital momentum directly from a laser beam [22]. In other related experiments, the (hyperfine) spin orientation of the atoms in a BEC is manipulated by adiabatically changing an external magnetic field in a position-dependent way to imprint a vortex [23]. In this way, vortices carrying multiple quanta of circulation were also created and their decay into singly quantized vortices was studied [24]. The latest development in vortex physics using internal state manipulation is the possibility to create synthetic magnetic fields [25], i.e. to quantum engineer an effective Lagrangian in which a vector potential term appears, not due to physical rotation or a physical magnetic field, but due to a position-dependent coupling between internal atomic states.

The second class of vortex experiments in quantum gases uses setups in which the condensate is physically rotated. Note that the circulation of a vortex and thus the condensate are quantized. This means that a minimal (critical) frequency is needed in order to have a vortex in a quantum gas. As soon as the critical frequency is reached, the rotational frequency of the condensate snaps to this critical frequency, and starts rotating as a whole. Although initially this method was unsuccessful, it is now the most widely used method to create and study vortices in superfluid quantum gases. The first observation of a ‘true’ vortex in quantum gases—distinct from the coreless vortex given in reference [20]—can be found in reference [26]. In that experiment a time-dependent optical potential (a ‘laser spoon’) is used in conjunction with a magnetic trap to stir the condensate, and this proved to be more effective than rotating only the magnetic trap. Below a critical frequency the condensate shakes (surface oscillation modes are present) but it does not rotate (no circulation). Above the critical frequency a singly quantized vortex is nucleated: it arises from a surface mode at the edge of the cloud and moves towards the centre where it remains as long as the cloud is kept in rotation. Using this technique, by rotating faster, an array of vortices was nucleated [27]. Also a ‘field-cooled’ version of this experiment was performed: rather than rotating a condensate, one cools a rotating normal cloud below the critical temperature [28]. When the temperature drops below the critical temperature, a condensate is nucleated in which the circulation ‘snaps’ to the nearest quantized value of circulation, with a corresponding number of vortices appearing. A large number of subsequent experiments have further refined the rotation technique for vortex creation—some noteworthy cases are the investigation of vortex nucleation [29], the development of a purely magnetic technique [30] and the use of a rotating optical lattice [31]. These experiments rely largely on the direct imaging technique, in which the vortex core is visible as a hole in the cloud, and the success of these experiments have dispelled the early doubts that the vortex cores might be too small to observe. Nevertheless, the matter-wave interferometric technique was also used to confirm the presence of the phase singularity, and to make sure that there is indeed circulation around the visible hole in the condensate [32, 33].

With techniques to produce vortices in quantum gases well-established, the path was open to study various aspects of vortex physics, such as the properties and dynamics of single vortex lines [34]. In the next three paragraphs we highlight three specific lines of experimental research that have been pursued very intensively: one that focuses on the properties of the vortex lattice, one that looks for spontaneous generation of vortices and anti-vortices, and finally one that focuses on quantum turbulence.

The equilibrium properties of vortex lattices were studied [35], as well as their formation dynamics [36, 37]. The response of the lattice to external perturbations such as shaking and removal of atoms was investigated, and revealed various unusual vortex configurations, such as giant vortices and vortex sheets. This also allowed to excite the ‘phonons’ or vibrations of the lattice, known as Tkatchenko modes [38]. Increasing the rotation frequency, attempts were made to reach the Landau-level regime, where the number of quantum vortices becomes comparable to the number of particles. Although this regime could not be reached, interesting physics regarding the lattice structure and equilibrium cloud shapes emerged for fast-rotated condensates [39].

By adapting the trapping potential, the atomic cloud can be squeezed into a pancake thin enough to reveal the physics of the two-dimensional (2D) Bose gas—in this case superfluidity is destroyed through the spontaneous appearance of vortices and anti-vortices above a critical temperature, the so-called Kosterlitz-Thouless transition, as given in reference [40]. The creation of vortex-anti-vortex pairs can also be induced in a three-dimensional (3D) cloud by moving an obstacle (such as a blue de-tuned laser beam) through the condensate [41]. The spontaneous generation of vortices and anti-vortices also occurs when a Bose gas is quenched through the phase transition. This quench leads to the nucleation of individual phase-coherent domains with different phases that upon merging lead to phase singularities at the interface of the merged domains. This phenomenon is called the Kibble-Zurek mechanism and with the advent of uniform (‘flat box’) trapping potentials [42] it has become possible to observe this transition in quantum gases [43].

By tilting the vortex array by a small angle around the axis of symmetry, the precession of the array was observed, and Kelvin modes (transverse oscillations of the vortex line) were seen [44]. The technique of tilting the axes of rotation would later prove useful in the studies on quantum turbulence. This is a regime of turbulent flow characterized by a tangle of vortices, and has been extensively studied in helium [45]. The difficulty in studying vortex tangles in liquid helium is that the core size is very small, so it is hard to directly image the vortices. In quantum gases, vortex cores can be readily imaged, and in addition new in situ [46] and tomographic [34] imaging techniques were developed. However, in quantum gases, the ratio between the condensate size and the vortex core is much smaller than in helium (typically it is of the order of 10–50 in BECs), limiting the size of the vortex tangle and inducing edge and finite-size effects. Quantum turbulence was achieved [47] by rotating the condensate around different axes simultaneously, or by rotating in conjunction with shaking. The possibility to flatten the atomic cloud has opened the way to also study quantum turbulence in two dimensions [48]. These techniques offer the potential to resolve many outstanding questions [49] regarding the transition from classical to quantum turbulence, the role of vortex reconnections and the dimensional crossover.

The above examples of recent vortex studies relate to Bose gases. Also in superfluid Fermi gases under rotation, vortex lattices have been observed [50]. The appearance of arrays of singly quantized vortices has been used as a hallmark of superfluidity in these systems. This in turn has been used to investigate how superfluidity breaks down when a population imbalance exists between the two states that combine to form the superfluid pairs [51]. Finally, to conclude this review section, we provide some references to more extensive reviews [52, 53] and a book [54] on vortices in quantum gases.

3. Building a theory of ultra-cold bosonic gases

Before we take a look at vortices, we will briefly discuss and derive the theory that we will use to describe ultra-cold bosonic gases. More elaborated discussions can be found in references [55–57]. The results of this section will be the Gross-Pitaevskii (GP) equation, a hydrodynamic description of the system and the energy of a general ultra-cold, dilute bosonic gas.

3.1. Choosing the right framework: classical or quantum?

In order to give a successful theoretical description of ultra-cold quantum gases, we first need to ask ourselves which framework to use. Is it sufficient to use the classical kinetic theory, or should we resort to quantum theory to get reliable answers? In order to answer this question, we have to look at the de Broglie wavelength

This wavelength is derived by comparing the kinetic energy of a free particle with its thermal energy, using the de Broglie hypothesis p = hk. This wavelength also emerges naturally when looking at ideal gases, see reference [58].

(or thermal wavelength):

${\lambda }_{\text{dB}}=\sqrt{\frac{h^2}{2\pi m{k}_{B}T}}\mathrm{.}$ E2

As soon as the de Broglie wavelength becomes comparable (or larger) than the typical distance between atoms, we need quantum mechanics. In this case the wavelength associated with a particle will become so large that the quantum mechanical wave functions will overlap.

For typical systems of ultra-cold quantum gases, we deal with clouds of around:

• 10⁵ atoms,

• with a scattering length in the order of nm,

• trapped in a magneto-optical trap with a typical length scale of μm,

• at temperatures of the order of nK.

The typical de Broglie wavelength of this system can now be estimated using the above, this results in:

${\lambda }_{\text{dB}}\propto {10}^{-5}\text{m}=10\text{μm}\mathrm{.}$ E3

The distance between the atoms can be estimated from the inverse particle density (n), which yields:

$d_{\text{atoms}}\propto \sqrt[3]{n^{-1}}\propto {10}^{-8}\text{m}=10\text{nm}\mathrm{.}$ E4

As we see the distance between the atoms is around a factor 1000 smaller than the de Broglie wavelength, which means that we should use quantum mechanics.

Since we are working with a gas of N atoms (many body system), we will start in second quantization [59, 60] since it is then easier to write down (and simplify) our Hamiltonian. The most general Hamiltonian for a system with two-body interactions is [61]:

${\displaystyle \widehat{H}}={\displaystyle \int }dr{{\displaystyle \widehat{\psi }}}^{†}(r)\left[-\frac{{\hslash }^2{\nabla }^2}{2m}+V_1(r)\right]{\displaystyle \widehat{\psi }}(r)+\frac{1}{2}{\displaystyle \int }dr{\displaystyle \int }d{\mathbf{r}}^{\prime }{{\displaystyle \widehat{\psi }}}^{†}(r){{\displaystyle \widehat{\psi }}}^{†}({\mathbf{r}}^{\prime })V_2(r-{\mathbf{r}}^{\prime }){\displaystyle \widehat{\psi }}({\mathbf{r}}^{\prime }){\displaystyle \widehat{\psi }}(r),$ E5

where $V_1(r)$ is a global potential (e.g. a trapping potential) and $V_2(r-{\mathbf{r}}^{\prime })$ describes the interactions.

Of course in reality three (or more) body interactions can be a possibility. Since ultra-cold atomic gases are usually dilute gases, interactions between more than two particles are negligible.

The operators ${\displaystyle \widehat{\psi }}(r)$ are bosonic operators, and the extension to fermionic atomic gases is discussed towards the end of the chapter. Determining the two potentials will fix the system, allowing the calculation of the dynamics of the system.

3.2. Simplifying by assumptions: ultra-cold dilute bosonic gases

For ultra-cold gases the main contribution to the partial wave expansion will come from the s-wave scattering. Since the range of the inter-atomic potential is much smaller than the inter-atomic distance (dilute gases), the potential can be approximated as a contact potential. Using partial wave analysis [62], the two-body potential can now be simply written as [63, 64]:

$V_2(r-{\mathbf{r}}^{\prime })=g\delta (r-{\mathbf{r}}^{\prime })\text{with}g=\frac{4\pi {\hslash }^2{a}_s}{m}\mathrm{,}$ E6

with a _s the (boson-boson

Also sometimes noted as a _BB . For now this does not matter, when you study fermions this will become more important.

) s-wave scattering length. Physically this s-wave scattering length can be interpreted as the effective radius a _s of our interacting particles if they would be approximated as hard spheres.

The larger the spheres (value of a _s ), the more interactions, and the larger the value of g.

A second approximation can be made if we assume that we have a Bose-Einstein condensate (BEC). The mechanism behind Bose-Einstein condensation can be easily seen when we look at the Bose-Einstein distribution:

$n(\epsilon ;\mu \mathrm{,}T)={\left[\exp \left(\frac{\epsilon -\mu }{k_{B}T}\right)-1\right]}^{-1}\text{with}\epsilon >0\text{and}\mu \le \mathrm{0,}$ E7

where the energy was chosen in such a way that the lowest energy level is situated in ε = 0. This distribution has a divergence in ε = μ and drops exponentially for large values of ε. This means that we can store an infinite amount of particles in the lowest energy level, while we can only have a finite amount in the higher energy levels. The total number of particles in our system of bosons is given by:

$N_{\text{tot}}=N_0+{\displaystyle \underset{N_{\text{exc}}}{{\displaystyle \underbrace{{\displaystyle {\displaystyle \underset{0}{\overset{\infty }{\int }}}}d\epsilon g(\epsilon )n(\epsilon ;\mu \mathrm{,}T)}}}}\mathrm{,}$ E8

where N ₀ is the number of particles in the ground state, $N_{\text{exc}}$ the number of particles in the excited state and $g(\epsilon )$ the density of states. When we start putting more and more particles in our system, we will see at a certain point that the number of particles in excited states will start saturating (the chemical potential becomes maximal). At the point where the number of excited states is saturated, every new particle will be placed in the ground state. If we keep increasing the number of particles, then at a certain point the number of particles in the ground state will become macroscopically large, dominating the many-body wave function. This means that for a BEC all particles are in the same quantum state and hence the many-particle wave function $\psi (r)$ can be approximated by a macroscopic wave function $\text{Ψ}(r)$ .

Notice that in order to achieve a BEC we want as many particles as possible, while for our approximation of the interaction potential we want a dilute gas. How can we combine this? The answer is to make the gas ultra-cold, this way the number of particles required to get a BEC becomes small enough to make the gas dilute. If our atomic gas is not diluted, higher order interactions will play a more important role. These higher order interactions will lead to particle being expelled from the atomic gas.

Combining both assumptions, it is now possible to simplify our Hamiltonian to make it equal to:

${\displaystyle \widehat{H}}=-\frac{{\hslash }^2{\nabla }^2}{2m}+V_1(r)+\frac{1}{2}g{\text{Ψ}}^{*}(r)\text{Ψ}(r).$ E9

In order to get our equation of motion for our system (9), we use the Heisenberg equation of motion, resulting in:

$i\hslash \frac{\partial \text{Ψ}(r\mathrm{,}t)}{\partial t}=\left(-\frac{{\hslash }^2{\nabla }^2}{2m}+V_1(r)+g{\left|\text{Ψ}(r\mathrm{,}t)\right|}^2\right)\text{Ψ}(r\mathrm{,}t).$ E10

The above equation (10) is known as the time dependent Gross-Pitaevskii equation (or the nonlinear Schrödinger equation). Note that we substituted the operator ${\displaystyle \widehat{\psi }}(r\mathrm{,}t)$ by the macroscopic field value $\text{Ψ}(r\mathrm{,}t)$ . Mathematically this substitution is not done since we go from an operator on Hilbert space to a function in space and time. However from a physical point of view this substitution makes sense. In order to see why the substitution of the operator is sensible, we look at the Penrose-Onsager criterion for BEC. The Penrose-Onsager criterion for BEC states that when the ground state becomes macroscopically occupied (which is the case), our one-particle density matrix should have an eigenvalue that becomes macroscopically large. This allows us to approximate the density matrix by its biggest contribution (the macroscopic wave function):

${\rho }_1(r;{\mathbf{r}}^{\prime })=⟨{{\displaystyle \widehat{\psi }}}^{†}(r){\displaystyle \widehat{\psi }}({\mathbf{r}}^{\prime })⟩\approx {\text{Ψ}}^{*}(r)\text{Ψ}({\mathbf{r}}^{\prime }).$ E11

Since all expectation values for our system can be rewritten in terms of the density matrix ${\rho }_1$ , we can replace all operators by the macroscopic field. With the two above Eqs. (9) and (10), we are ready to describe our condensate of bosons.

Note that if you want to minimize your energy (9) to find the ground state of your system with a given normalization, you explicitly have to use Lagrange multipliers to fix your normalization. This is because of the fact that our Hamiltonian and thus equations of motion are non-linear. For the ground state this means that you define the Hamiltonian:

${{\displaystyle \widehat{H}}}_{\mu }=-\frac{{\hslash }^2{\nabla }^2}{2m}-\mu +V_1(r)+\frac{1}{2}g{\text{Ψ}}^{*}(r)\text{Ψ}(r),$ E12

which now contains the chemical potential as a Lagrange multiplier for the number of particles. Doing a functional minimization then yields the differential equation:

$\mu {\text{Ψ}}_{\mu }(r\mathrm{,}t)=\left(-\frac{{\hslash }^2{\nabla }^2}{2m}+V_1(r)+g{\left|{\text{Ψ}}_{\mu }(r\mathrm{,}t)\right|}^2\right){\text{Ψ}}_{\mu }(r\mathrm{,}t).$ E13

So depending on whether you work with a fixed number of particles or a fixed chemical potential you use respectively, Eqs.(12) and (13) or Eqs.(9) and (10). The relation between both views can be easily established by comparison. For the energy we get the known Legendre transform from statistical mechanics:

${\displaystyle \widehat{H}}={{\displaystyle \widehat{H}}}_{\mu }+\mu {\displaystyle \widehat{N}}\mathrm{.}$ E14

For the macroscopic wave function, we have an additional phase rotation:

$\text{Ψ}(r\mathrm{,}t)=\exp \left(\frac{i\mu t}{\hslash }\right){\text{Ψ}}_{\mu }(r\mathrm{,}t).$ E15

With the above two relations it is possible to switch between fixed number of particles or fixed chemical potential if needed.

3.3. The hydrodynamic description

In order to simplify the calculations one usually works in the hydrodynamic description. This means that our complex field $\text{Ψ}(r\mathrm{,}t)$ will be written as the product of two real fields $n(r\mathrm{,}t)$ and $S(r\mathrm{,}t)$ as:

$\text{Ψ}(r\mathrm{,}t)=\sqrt{n(r\mathrm{,}t)}{e}^{iS(r\mathrm{,}t)}\mathrm{.}$ E16

The field $n(r\mathrm{,}t)$ is known as the superfluid density (field), whereas $S(r\mathrm{,}t)$ is known as the phase field.

If we now take the Gross-Pitaevskii equation, Eq. (10) and treat it canonically

Multiply (10) by the complex conjugate field ${\text{Ψ}}^{*}(r\mathrm{,}t)$ , and subtract the product ${\text{Ψ}}^{*}{\partial }_t\text{Ψ}$ by the complex conjugate of the whole $\text{Ψ}{\partial }_t{\text{Ψ}}^{*}$ . In other words, calculate ${\text{Ψ}}^{*}{\partial }_t\text{Ψ}-\text{Ψ}{\partial }_t{\text{Ψ}}^{*}$ .

we get, after simplification:

$\frac{\partial {\left|\text{Ψ}(r\mathrm{,}t)\right|}^2}{\partial t}+\nabla \cdot \left(\frac{\hslash }{2mi}\left[{\text{Ψ}}^{*}(r\mathrm{,}t)\nabla \text{Ψ}(r\mathrm{,}t)-\text{Ψ}(r\mathrm{,}t)\nabla {\text{Ψ}}^{*}(r\mathrm{,}t)\right]\right)=0.$ E17

This is a result known from standard quantum mechanics, where it is seen as a continuity equation since ${\left|\text{Ψ}(r\mathrm{,}t)\right|}^2=n(r\mathrm{,}t)$ is interpreted as a (probability) density. If (17) is a continuity equation, then the part with a divergence is given by a current $j(r\mathrm{,}t)=n(r\mathrm{,}t)v(r\mathrm{,}t)$ . Using the fact that the argument of the divergence is a current, we can define our velocity field as:

$v(r\mathrm{,}t)=\frac{\hslash }{2mi}\frac{{\text{Ψ}}^{*}(r\mathrm{,}t)\nabla \text{Ψ}(r\mathrm{,}t)-\text{Ψ}(r\mathrm{,}t)\nabla {\text{Ψ}}^{*}(r\mathrm{,}t)}{{\left|\text{Ψ}(r\mathrm{,}t)\right|}^2}\mathrm{.}$ E18

Of course we want to work in the hydrodynamic equations so we will substitute the hydrodynamic description (16) in the continuity equation, Eq. (17) and the definition for the velocity field, Eq. (18). For the continuity equation, Eq. (17) this yields (as expected) that:

$\frac{\partial n(r\mathrm{,}t)}{\partial t}+\nabla \cdot \left[n(r\mathrm{,}t)v(r\mathrm{,}t)\right]=\mathrm{0,}$ E19

which is known from classical hydrodynamics. For the velocity field (18), we find that

$v(r\mathrm{,}t)=\frac{\hslash }{m}\nabla S(r\mathrm{,}t),$ E20

so the velocity field is determined entirely by the phase field. Note that our velocity field is also irrotational:

$\nabla \times v(r\mathrm{,}t)=\nabla \times \left(\frac{\hslash }{m}\nabla S(r\mathrm{,}t)\right)=\mathrm{0,}$ E21

meaning that we will only have irrotational flow patterns.

Of course since we have two real fields $n(r\mathrm{,}t)$ and $S(r\mathrm{,}t)$ , we also need two equations of motion. To find these equations of motion, we need to substitute the hydrodynamic description (16) into the Gross-Pitaevskii equation, Eq. (10). After the substitution, the Gross-Pitaevskii equation, Eq. (10) splits up into real and imaginary parts which can be considered as two independent equations since both fields are real-valued. The real and imaginary parts after minor simplifications are given by:

$\hslash \frac{\partial \sqrt{n(r\mathrm{,}t)}}{\partial t}=-\frac{{\hslash }^2}{2m}\left[2\left(\nabla \sqrt{n(r\mathrm{,}t)}\right)\cdot \left(\nabla S(r\mathrm{,}t)\right)+\sqrt{n(r\mathrm{,}t)}\left({\nabla }^{2}S(r\mathrm{,}t)\right)\right]\mathrm{,}$ E22

$-m\frac{\partial }{\partial t}\left(\frac{\hslash }{m}S(r\mathrm{,}t)\right)=-\frac{{\hslash }^2}{2m}\frac{{\nabla }^2\left(\sqrt{n(r\mathrm{,}t)}\right)}{\sqrt{n(r\mathrm{,}t)}}+\frac{m}{2}{\left(\frac{\hslash }{m}\nabla S(r\mathrm{,}t)\right)}^2+V_1(r)+gn(r\mathrm{,}t),$ E23

where the first line (22) is the imaginary part and the second line (23) is the real part of the Gross-Pitaevskii equation, Eq. (10). The imaginary part (22) can be simplified and rewritten in terms of the velocity field and density field, this will yield the continuity equation, Eq. (19). The real part (23) can also be rewritten in terms of the velocity and density fields, this is most easily done by taking the gradient of the equation (23). After simplifying the gradient of the real part (23) we get:

$m\frac{\partial v(r\mathrm{,}t)}{\partial t}+m\left(v(r\mathrm{,}t)\cdot \nabla \right)v(r\mathrm{,}t)+\nabla \cdot \left(-\frac{{\hslash }^2}{2m}\frac{{\nabla }^2\left(\sqrt{n(r\mathrm{,}t)}\right)}{\sqrt{n(r\mathrm{,}t)}}+V_1(r)+gn(r\mathrm{,}t)\right)=\mathrm{0,}$ E24

which is the same as the Euler equation (or frictionless Navier-Stokes equation) from hydrodynamics. In other words, in a BEC we have an irrotational, frictionless flow. By looking at the Euler equation (24), we can identify three kinds of pressure that will occur in our system:

Mechanical pressure: The first pressure we can recognize is the one particle potential $V_1(r)$ . This pressure is called the mechanical pressure, and this is the pressure used to keep the condensate together.

Interaction pressure: The second pressure is given by the interaction term $gn(r\mathrm{,}t)$ and is known as the interaction pressure. This pressure is a consequence of the collisions between the atoms, leading to a pressure that is proportional to the density.

Quantum pressure: The remaining term

$-\frac{{\hslash }^2}{2m}\frac{{\nabla }^2\left(\sqrt{n(r\mathrm{,}t)}\right)}{\sqrt{n(r\mathrm{,}t)}}$ E25

is known as the quantum pressure. This pressure is a consequence of quantum mechanics where your kinetic energy is proportional to the second derivative of the wave function (the curvature). Note that this pressure is a pure quantum mechanical effect, indeed if we take the limit h → 0 this term becomes equal to zero.

3.4. Summary

So for our ultra-cold bosonic gases forming a BEC, we get one macroscopic wave function $\text{Ψ}(r\mathrm{,}t)$ which is subject to the Gross-Pitaevskii equation, Eq. (10). This complex wave function $\text{Ψ}(r\mathrm{,}t)$ can be rewritten in terms of two real valued fields (16): a density field $n(r\mathrm{,}t)$ and a phase field $S(r\mathrm{,}t)$ . The description in terms of the real fields is called the hydrodynamic description, where one usually uses the velocity field (20) instead of the phase field. The relations for the density field and velocity field are given by the continuity equation, Eq. (19) and the Euler equation, Eq. (24), where the quantum mechanics enters through the quantum pressure. From the hydrodynamic description we concluded that a BEC typically describes a frictionless, irrotational flow.

4. Basic properties of a vortex in a superfluid

Now that we know the basic physics for a BEC, we can look at vortices in ultra-cold gases. Note that the results discussed here will be general for both vortices in superfluid Bose gases as well as in superfluid Fermi gases.

This is because of the fact that fermions have to pair up and become bosons in order to become superfluid. So the effective superfluid behaves like a BEC.

In this section we will not look at the structure of the vortex (this is done in the next section), but we will rather define what we mean with a vortex and what this implies for our vortex flow. By thoroughly studying the definition of a vortex we can already derive a few basic properties of a vortex. The rest of the section is devoted to two analytically solvable systems of a BEC and on how vortices are detected in experiments.

4.1. Definition of a (quantum) vortex

In order to study a vortex, we first have to explain what we mean when we talk about a vortex. Just as in classical fluid dynamics, we start from the definition of the circulation

In fluid dynamics the circulation is usually denoted by Γ.

κ, defined as:

$\kappa ={\displaystyle {\oint }_C}ds\cdot v(r\mathrm{,}t),$ E26

where C is a closed contour containing the vortex core. For the circulation we then have three possibilities: counter-clockwise rotation (κ > 0), clockwise rotation (κ < 0) and no rotation κ = 0. The value of the circulation κ should be the same for any contour C that encloses the vortex core.

This is because of the fact that we assume that the different fields (density field, velocity field,…) are meromorphic functions and hence the singularities are localized. According to Morera’s theorem we can freely deform the contour (as long as we do not cross a singularity) and still have the same contour integral.

Given the velocity field $v(r\mathrm{,}t)$ , it is always possible to define a vorticity field as:

$\mathbf{\omega }(r\mathrm{,}t)=\nabla \times v(r\mathrm{,}t),$ E27

which is a quantity telling us how the fluid rotates locally. The vorticity follows the vorticity equation which can be derived by taking the curl of the Euler equation, Eq. (24). Applying Stokes’ theorem in Eq. (26), we can relate the circulation κ in terms of the vorticity (24), this results in:

$\kappa ={\displaystyle {\int }_S}dS\cdot \omega (r\mathrm{,}t),$ E28

making the circulation equal to the flux of vorticity. Note that we can use Stokes’ theorem as long as S is a surface without any holes.

Using the above two formulas (26) and (28) we can calculate the circulation of a BEC. To calculate the circulation using the velocity field, we will use the definition of the velocity field in terms of the phase field (20). The circulation then becomes:

$\kappa =\frac{\hslash }{m}{\displaystyle {\oint }_C}\nabla S(r\mathrm{,}t)=\frac{\hslash }{m}\left(S(r_e\mathrm{,}t)-S(r_b\mathrm{,}t)\right)\mathrm{,}$ E29

where the gradient theorem was used. In (29) r _b is the point where the closed loop begins and r _e is the point where the loop ends, although the physical positions are the same the points r _b and r _e can be different. An example of two points with a different coordinate but the same physical location is given for example when we describe the phase in polar coordinates $S(\rho \mathrm{,}\varphi \mathrm{,}z\mathrm{,}t)$ . Making one complete circle in the angular coordinate then yields for the start and end points that:

$S(r_b\mathrm{,}t)=S(\rho \mathrm{,}\varphi \mathrm{,}z\mathrm{,}t)\text{and}S(r_e\mathrm{,}t)=S(\rho \mathrm{,}\varphi +2\pi \mathrm{,}z\mathrm{,}t).$ E30

Since we expect to find the same physics again when we make one round, we also expect our quantum mechanical wave function to be the same up to a global phase that is a multiple of 2π. The fact that we expect to find the same physics again tells us that the phase difference in (29) has to be given by:

$S(r_e\mathrm{,}t)-S(r_b\mathrm{,}t)=n2\pi \text{with}n\in \mathbb{Z}\mathrm{.}$ E31

Using the above restriction (31) on the phase difference the equation for the circulation (29) then yields:

$\kappa =n\frac{h}{m}\text{with}n\in \mathbb{Z}\mathrm{,}$ E32

in other words the circulation is quantized in multiples of h / m.

The circulation can now also be calculated by using the formula containing the vorticity field (28), yielding

$\kappa ={\displaystyle {\int }_S}dS\cdot \omega (r\mathrm{,}t)={\displaystyle {\int }_S}dS\cdot \left[\nabla \times \left(\frac{\hslash }{m}\nabla S(r\mathrm{,}t)\right)\right]=\mathrm{0,}$ E33

where we got zero since the curl of a gradient field is always zero. The fact that our two results for the circulation Eqs. (32) and (33) are in contradiction means we did something (mathematically) wrong. Since Eq. (26) is the definition and the restriction Eq. (31) is exact our result, Eq. (32) should be the right one. The only assumption that we made to go from the definition of the circulation Eq. (26) to the formula with the vorticity field, Eq. (28) is our condensate (and thus the surface we integrate over) has no holes. The fact that the results are different means that it was wrong to assume that the condensate has no holes. Or in other words, having a vortex (so a non-vanishing circulation: n ≠ 0 in Eq. (32)) will create a hole in the condensate. We can even go further and say that it is impossible for a vortex to begin and end in two different points inside of the BEC, and this can be illustrated by contradiction. Suppose we have a vortex that ends (or begins) inside of the BEC, then we can choose a surface S belonging to the contour C in such a way that the surface does not intersect with the vortex, an example is given in Figure 1 . If we can choose our surface as illustrated in Figure 1 , then Stokes’ theorem should hold, leading to a contradiction when calculating the circulation using the two formulas (26) and (28). What we can conclude from this is that vortices extend through the entire BEC,

The exception to this is vortex rings which are closed rings with a non-vanishing circulation. But since a ring has no beginning or ending this is actually no exception to the rule.

with end and begin points at the end of the condensate.

Figure 1.

In blue a cylindrical vortex that ends/begins inside the BEC. In red the contour C and its surface S is chosen in such a way that the surface does not intersect with the vortex core.

From the definition of the circulation we can now make three general conclusions about vortices in ultra-cold gases:

• Vortices will make a hole in the BEC.

• The circulation is quantized with quanta h / m.

• Vortices cannot have a beginning or ending within the BEC.

The second property is also the reason why vortices in ultra-cold gases are called quantum vortices, and throughout this chapter, we will however drop the prefix quantum. The last property also tells us that vortices have to be formed at the edges of the BEC, since they cannot end the start inside the BEC. This can be seen in experiments where a condensate that is removing its vortices just pushes them to the outside of the condensate, where they leave.

In this chapter we will look at the simplest kind of vortices: vortices which have a cylindrically symmetric. This means that we will be looking at straight vortices with through the centre the (rotational) symmetry axis, for which we will choose the z axis. If our macroscopic wave function $\text{Ψ}(r\mathrm{,}t)$ has a cylindrical symmetry, it is possible to rewrite the wave function as:

$\text{Ψ}(r\mathrm{,}t)=f(\rho \mathrm{,}z\mathrm{,}t)e^{in\varphi }\mathrm{,}$ E34

where $f(\rho \mathrm{,}z\mathrm{,}t)$ is a real-valued structure function. Using the symmetric form of the wave function, Eq. (34), we can calculate the velocity field using Eq. (20), this yields:

$v(r\mathrm{,}t)=\frac{\hslash }{m}\left(\frac{1}{\rho }\frac{\partial (n\varphi )}{\partial \varphi }{e}_{\varphi }\right)=n\frac{\hslash }{m\rho }{e}_{\varphi }\mathrm{.}$ E35

What we see is that the velocity field lies along the ϕ direction, which is indeed a spinning motion in the counter-clockwise direction

Since we are working in a right-handed coordinate system $(\rho \mathrm{,}\varphi \mathrm{,}z)$ , the direction of rotation along ϕ is given by the right-hand rule. This results in a counter-clockwise direction for the velocity field of the condensate when viewed from above (along the −e _z direction).

for positive values of n (vortices). If n is a negative number, the singularity is called an anti-vortex, and we have a clockwise rotation of the condensate. We also see that the velocity field drops as $\propto 1/\rho$ as we move away from the vortex core. Using the velocity field Eq. (35) we can calculate the angular momentum, this becomes equal to:

$L_z(r\mathrm{,}t)=r\times mv(r\mathrm{,}t)=m\rho v(r\mathrm{,}t)=n\hslash \mathrm{.}$ E36

This tells us that the angular momentum is quantized, every time that we make a vortex the condensate gets an angular momentum kick of the size nℏ. Finally we can also look at the circulation (26) of the wave function $\text{Ψ}(r\mathrm{,}t)$ , this yields:

$\kappa ={\displaystyle {\oint }_C}ds\cdot v(r\mathrm{,}t)={\displaystyle {\displaystyle \underset{0}{\overset{2\pi }{\int }}}}\rho d\varphi v(r\mathrm{,}t)={\displaystyle {\displaystyle \underset{0}{\overset{2\pi }{\int }}}}\rho d\varphi n\frac{\hslash }{m\rho }{e}_{\varphi }=n\frac{h}{m}\mathrm{.}$ E37

Comparing the above result (37) with the circulation of a BEC (32), we can indeed see that the wave function, Eq. (34) can describe a cylindrically symmetric vortex. The number n in the phase factor of the wave function, Eq. (34) is the number of circulation quanta h / m the vortex will have.

With now the macroscopic wave function (34) as a starting point, we can now start our study of vortices in a BEC. In order to determine the vortex structure, we should determine the functional dependence of $f(\rho \mathrm{,}z\mathrm{,}t)$ either analytically (by using a variational model) or numerically, this is done in the next section. To conclude this section we will look at two analytic cases: the homogeneous condensate, and the condensate near a hard wall. The first analytic case will tell us how the condensate behaves far away from the vortex, and this asymptotic behaviour will come in handy both analytically and numerically. The second analytic case tells us how the condensate heals if you make a hole in it, which is the basis of a commonly used variational model.

4.2. A stationary homogeneous BEC

The first system that can be solved analytically is the homogeneous BEC. For this system we expect that the density is given by:

${\left|\text{Ψ}(r)\right|}^2=n_{\infty }\mathrm{,}$ E38

with $n_{\infty }$ the particle density of the BEC. If we want to work with a fixed number of N particles, we need to determine the chemical potential. Substituting the wave function, Eq. (38) into the Gross-Pitaevskii equation for a fixed number of particles, Eq. (13) then yields:

$g{n}_{\infty }\text{Ψ}(r)=\mu \text{Ψ}(r).$ E39

Solving the above equation then yields for the chemical potential:

$\mu =g{n}_{\infty }\mathrm{.}$ E40

The energy of the condensate can be found by taking the expectation value of the Hamiltonian (9), resulting in:

$E_N[\text{Ψ}]=\frac{1}{2}g{\displaystyle {\int }_V}{d}^{3}r{\left|\text{Ψ}(r)\right|}^4=\frac{1}{2}g{n}_{\infty }^{2}V=\frac{1}{2}g\frac{N^2}{V}=\frac{1}{2}\mu N\mathrm{,}$ E41

with V the volume of the BEC and N the total number of particles in the BEC. The above formula for the energy (41) tells us that a BEC likes to spread out and maximize its volume (this decreases the energy). This means that in experiments we will need an external potential $V_1(r)$ to keep the BEC together. Usually this is a harmonic trapping potential of the form

$V_1(r)=\frac{1}{2}m{\omega }_{\perp }^2{\rho }^2+\frac{1}{2}m{\omega }_{\parallel }^2{z}^2\mathrm{,}$ E42

which follows from the magneto-optical trap. This harmonic trapping potential leads to cigar-shaped condensates. Another type of confinement potential is the hard wall potential, this trapping potential is realized by using a hollow laser [65], this leads to a potential of the form:

$V_1(r)=\text{Θ}(R-\rho ).$ E43

In this chapter, we will mainly focus on the latter hard-wall potential (43) because this yields the simplest calculations and still contains all of the physics. For the harmonic-trapping potential, it is also possible to make simple approximations in order to get easier calculations. The most common approximation made when working with a harmonic trapping potential is the Thomas-Fermi approximation. In the Thomas-Fermi approximation, it is assumed that the kinetic energy is negligible, and this is the case if the number of particles is large enough. In typical experiments this is usually the case.

4.3. A stationary BEC near a hard wall and the healing length

The next analytic case that will be considered is if we put a hard wall edge in the BEC. In other words we will punch a hole in the BEC and see how it heals. This analytic case is important because it will give us the length scale over which the condensate heals, also called the healing length. This healing length will be useful to use as a typical length scale of the system in order to make numerical equations dimensionless. To introduce the healing length, one usually looks at the Gross-Pitaevskii equation, Eq. (10) and fills in the length scale for which the kinetic energy term becomes equal to the interaction term. When substituting the second derivative in Eq. (10) by $1/{\xi }^2$ this then yields the equation:

$\frac{{\hslash }^2}{2m{\xi }^2}=g{n}_{\infty }\mathrm{,}$ E44

for the healing length ξ this solves to:

$\xi =\frac{\hslash }{\sqrt{2m{n}_{\infty }g}}=\frac{1}{\sqrt{8\pi n_{\infty }{a}_s}}\mathrm{.}$ E45

Given the healing length, we can now again look at a free BEC (in x > 0) in which we place a hard wall potential at x = 0. Filling in the Gross-Pitaevskii equation, Eq. (10) for a 1D BEC now yields:

$-\frac{{\hslash }^2}{2m}\frac{d^2\text{Ψ}(x)}{d{x}^2}+g{\left|\text{Ψ}(x)\right|}^2\text{Ψ}(x)=\mu \text{Ψ}(x).$ E46

In the limit x → ∞ we get our homogeneous BEC back, meaning that the chemical potential is given by Eq. (40). If we now define the structure function f(x) as

$\text{Ψ}(x)=\sqrt{n_{\infty }}f(x)$ E47

and we fill in the chemical potential, then the Gross-Pitaevskii equation can be written as:

${\xi }^2\frac{d^{2}f(x)}{d{x}^2}-\left[f{(x)}^2-1\right]f(x)=\mathrm{0,}$ E48

where we divided the Gross-Pitaevskii equation by gn _∞ in order to get the healing length. The above non-linear differential equation is solved by substituting a hyperbolic tangent. The solution for the condensate wave function with a hard wall potential is then given by:

$\text{Ψ}(x)=\sqrt{n_{\infty }}\tanh \left(\frac{x}{\sqrt{2}\xi }\right)\mathrm{.}$ E49

What we see is that the condensate ‘heals’ from a hole over a typical length given by the healing length ξ, and this is also illustrated in Figure 2 .

Figure 2.

The healing of a BEC near a hard wall over a typical length of ξ.

This solution for a hard wall edge of the BEC will play an important role when we will be studying vortices. As said at the beginning of this section, a vortex will make a hole in the condensate from which the condensate also heals. It is to be expected that when the condensate heals from a hole of a vortex, it happens in an analogous way as when it heals from a hole due to a hard wall. This is because of the fact that both holes follow the same condensate physics, with the same boundary condition: $\text{Ψ}(x=0)=0$ , where x = 0 is either the centre of the vortex or the location of the hard wall. The difference for the vortex is that there is an extra rotation of the condensate, this rotation tends to ‘smear out’ the hole of the vortex, resulting in a larger length over which the condensate heals from a vortex: $\alpha \xi$ (with $\alpha \ge 1$ ). Thanks to the analogous behaviour as a vortex, the hard wall solution will seem to be a very good variational model to study vortices in ultra-cold gases, even for fermionic gases [66].

4.4. Detection of vortices in superfluids

In essence, every experiment on atomic clouds is done by scattering probe laser light off it. So essentially experimenters will take a ‘picture’ of the condensate by using a laser in either a destructive or non-destructive manner. This ‘picture’ that is obtained is in essence the plot of the superfluid density field $n(r\mathrm{,}t)$ of the hydrodynamic description given in Eq. (16) at a given time. A phase can also be measured by the way the condensate interacts with either the laser or another condensate.

This means that the main source of information is given by visual information. For the detection of vortices this is problematic since the healing length is of the order of 10 nm (if we fill in the quantities of the third section). To get an idea of how the healing lengths compare to the scattering lengths for several elements that are used to make a BEC, we have given a few values in Table 1 for typical condensate atoms.

Condensate atom	Scattering length as	Healing length ξ	as / ξ
Cesium (^133Cs)	300	3.646 × 10^2	1.214
Rubidium (^87Rb)	100	6.31 × 10^2	6.31
Sodium (^23Na)	50	8.92 × 10^2	17.84
Lithium (^7Li)	30	1.15 × 10^3	38.3
Hydrogen (^1H)	1	6.31 × 10^3	6310

Table 1.

The approximate scattering length a _s and approximate healing length ξ (44) for a few typical atoms that undergo BEC [55].

Both lengths are given in units of a ₀ = 10⁻⁸cm (or 0.1 nm), for a typical particle density equal to n _∞ = 10¹⁵cc⁻¹ (or 10²¹m⁻³). Note that although the healing length is very small, it is always considerably larger than the s-wave scattering length.

From Table 1 , we see that the healing length is indeed a very small value, which makes vortices hard to detect (and even measure their structure). The common method to detect vortices in a superfluid is to turn off the confining potential. As we know a BEC likes to maximize its volume, so turning off the confinement will lead to an expanding cloud of atoms (also called ‘Ballistic expansion’ by experimenters). Since the BEC expands, the cores of the vortices will start expanding as well, making the vortex cores easier to detect. Typical images that we get from this kind of experiments are given in Figure 3 . As we see the vortices become visible after expansion. The two things that we can see in an experiment are how many vortices we have and how they are ordered.

Figure 3.

Four experimental outcomes for vortices, depending on how fast the condensate is stirred. From left to right (a–d), the stirring frequency increases, leading to the formation of more vortices in the BEC. The diameter of the cloud in figure (d) is 1 mm after the ballistic expansion, which is a magnification of 20 times. Figures taken from Ref. [67].

During ballistic expansion, we can also extract extra information by using interference of matter waves (condensates). When we make two BECs overlap they will give rise to a matter interference pattern. This means that if we take a reference BEC (without any vortices) and make this interfere with a BEC that contains a vortex, the effect of the vortex can be seen in the interference pattern. The vortex shows up as a fork dislocation, an example of a BEC with one vortex overlapping with a BEC with no vortices is given in Figure 4 . As seen in the figure, the presence of the vortex leads to an irregularity in the interference pattern that is known as a fork dislocation. The fork dislocation was also already known in optical systems where optical vortices are also a field of study [69].

Figure 4.

Two BECs in a trap, when the trap is turned off the condensates will expand and eventually overlap leading to an interference pattern. The top figures show two condensates without a vortex, in the bottom figures the left condensate has a vortex. In the interference pattern with the vortex, the fork dislocation becomes apparent. Figure taken from Ref. [68].

The additional information that we get with respect to the simple ballistic expansion is the fact that we can now also see the number of vortex quanta (or circulation quanta) that are present in the vortex. The more circulation quanta a vortex has, the more fringes will be involved in the fork dislocation. The interference pattern of two condensates can also tell us something about the sign of the circulation. A vortex (velocity field gives a counter-clockwise rotation) will have a fork dislocation which is upside down compared to the fork dislocation of an anti-vortex (spinning clockwise).

5. Vortex formation in Bose-Einstein condensates

With two basic solutions for ultra-cold gases and the general properties of a vortex from the previous section we can now continue. In this section we will first look at how vortices appear in a condensate. Afterwards, we will look at the structure of a single vortex and two simple variational models.

5.1. Energy of a rotating BEC

Since vortices classically appear when we start stirring our system, a rotating condensate might seem a good setting for a vortex system. In atomic gases this rotating system is achieved by stirring the condensate using a (blue de-tuned) laser which acts like a spoon. We will look at a system that rotates with a rotational velocity equal to Ω. Looking back at the previous section it only makes sense that we should be looking for vortices in a rotating condensate. The only way to get a rotating BEC (and hence a non-zero circulation in Eq. (26)) is by having a vortex in the BEC.

A condensate with the same amount of vortices and anti-vortices does not need to be rotating since these circulations cancel each other out.

From the third section the energy of a BEC is already known and given by the expectation value of Eq. (9). This leads to three different energy contributions given by:

$E_{\text{kin}}[\text{Ψ}(r)]=\frac{-{\hslash }^2}{2m}{\displaystyle {\int }_V}{d}^{3}r{\text{Ψ}}^{*}(r){\nabla }^2\text{Ψ}(r)=\frac{{\hslash }^2}{2m}{\displaystyle {\int }_V}{d}^{3}r{\left|\nabla \text{Ψ}(r)\right|}^2\mathrm{,}$ E50

$E_{\text{trap}}[\text{Ψ}(r)]={\displaystyle {\int }_V}{d}^{3}r{\text{Ψ}}^{*}(r)V_1(r)\text{Ψ}(r),$ E51

$E_{\text{int}}[\text{Ψ}(r)]=\frac{1}{2}g{\displaystyle \int }{d}^{3}r{\left|\text{Ψ}(r)\right|}^4\mathrm{,}$ E52

where the second equality in Eq. (50) follows from partial integration and the fact that the boundary term drops out.

We either will look at an infinite condensate or a condensate with a hard-wall edge. Both cases yield a zero valued boundary term

The time dependency was explicitly left out since we will be looking at stationary solutions for the vortex in the remainder of the chapter.

Since for vortices we will be working in the hydrodynamic description, we will substitute the hydrodynamic form of the wave function in Eq. (16) in the different energy contributions in Eqs. (50)–(52), this yields:

$E_{\text{kin}}[\text{Ψ}(r)]=\frac{{\hslash }^2}{2m}{\displaystyle {\int }_V}{d}^{3}r{\left(\nabla \sqrt{n(r)}\right)}^2+\frac{m}{2}{\displaystyle {\int }_V}{d}^{3}r\left[n(r)v^2(r)\right]\mathrm{,}$ E53

$E_{\text{trap}}[\text{Ψ}(r)]={\displaystyle {\int }_V}{d}^{3}rn(r)V_1(r),$ E54

$E_{\text{int}}[\text{Ψ}(r)]=\frac{1}{2}g{\displaystyle \int }{d}^{3}r{\left(n(r)\right)}^2\mathrm{,}$ E55

The only missing ingredient is the Hamiltonian in a rotating frame of reference, this Hamiltonian is given by:

${{\displaystyle \widehat{H}}}^{\prime }={\displaystyle \widehat{H}}-\mathbf{\Omega }\cdot {\displaystyle \widehat{\mathbf{L}}}\mathrm{.}$ E56

This additional term $-\mathbf{\Omega }\cdot {\displaystyle \widehat{\mathbf{L}}}$ is nothing else but the rotational energy (also known from classical mechanics). What we see is that in order to form a vortex the profit in energy when rotating $\mathbf{\text{Ω}}\cdot {\displaystyle \widehat{\mathbf{L}}}$ should be larger than the energy cost of making the vortex. In other words, a vortex will be formed if in Eq. (56) we have that $⟨{{\displaystyle \widehat{H}}}^{\prime }⟩$ with a vortex is smaller than $⟨{\displaystyle \widehat{H}}⟩$ without a vortex. In order to give an estimate of the minimal required rotational velocity Ω needed to form a vortex we will need a model for the vortex core. To get a fairly rough, but close estimate we introduce the simplest possible model for a vortex: the cylindrical vortex hole. This model will allow us to easily make a few statements about vortex stability and interactions.

5.2. Simplest model for a vortex

The simplest model for a vortex that we can imagine is to represent the vortex by a hollow cylindrical hole. The radius of the vortex hole is chosen to be equal to the healing length

You can change this by some effective length R _cylinder, this will not change the calculations or results

and the height of the vortex given by H. As we know the velocity field of the vortex should be given by Eq. (35), only inside the vortex core there is no condensate and we are free to choose our velocity field. In order to have a correct circulation (even within the vortex core) we assume the velocity field to take a constant value. The simple vortex model is depicted in Figure 5 .

Figure 5.

The simplest possible vortex model: a hollow cylindrical vortex core with radius ξ in a cylindrical condensate with height H and radius R, shown in figure (a). Figure (b) shows the corresponding vortex velocity field.

In formulas this means that the simple vortex model is given by:

$n(r)=n_{\infty }\text{Θ}(\rho -\xi )\text{Θ}(R-\rho ),$ E57

$v(r)=n\frac{\hslash }{m\rho }{e}_{\varphi }\text{Θ}(\rho -\xi )\text{Θ}(R-\rho )+n\frac{\hslash }{m\xi }{e}_{\varphi }\text{Θ}(\xi -\rho ).$ E58

With this model (and the above formulas) we can now calculate the different energies as given in Eqs. (41–44) for the rotating condensate.

5.3. Formation of one and more vortices

Since only the energy difference compared to the homogeneous (vortex free) state is relevant, we will only calculate the energy differences between the homogeneous and vortex state. We will look at a condensate confined by an infinite potential well, so $V_1(r)=0$ within the condensate. This means that we have two possible energy contributions, one for the velocity field and one for the interactions.

If we have a vortex in our condensate, then we will get a velocity field of the form (58), which will give an extra kinetic energy. The difference kinetic energy (41) then becomes

The first term drops out since we have a homogeneous condensate and thus $\nabla \left(\sqrt{n(r)}\right)=0$ .

:

$\text{Δ}{E}_{\text{kin}}={\displaystyle {\int }_V}{d}^{3}r\left(\frac{m}{2}n(r)v^2(r)\right)=2\pi \frac{1}{2}m{n}_{\infty }H{\displaystyle {\displaystyle \underset{\xi }{\overset{R}{\int }}}}\rho d\rho {\left(n\frac{\hslash }{m\rho }\right)}^2=n^2\pi n_{\infty }\frac{{\hslash }^2}{m}H\ln \left(\frac{R}{\xi }\right)\mathrm{.}$ E59

We also have a contribution in energy of the interaction energy, namely the energy that is needed to punch a hole in the BEC. This needed energy is given by

$\text{Δ}{E}_{\text{int}}=\frac{1}{2}g{\displaystyle {\int }_V}{d}^{3}r{\left(n(r)\right)}^2=\frac{1}{2}g{V}_{\text{hole}}\mathrm{,}$ E60

with

$V_{\text{hole}}=\pi {\xi }^{2}H\mathrm{,}$ E61

the volume of the hole that has been made in the condensate. If we now assume (for simplicity) that

This is not a bad assumption since in most condensates we have that R / ξ ≈ 30.

$R>>\xi$ (the condensate is much larger than the volume of the hole), then the energy contribution of Eq. (60) will be negligible. The total cost of energy for making a vortex is thus given by:

$\text{Δ}E=\text{Δ}{E}_{\text{kin}}=n^2\pi n_{\infty }\frac{{\hslash }^2}{m}H\ln \left(\frac{R}{\xi }\right)\mathrm{.}$ E62

In order to be advantageous, the energy cost of making a vortex from Eq. (62) should be smaller as the energy reduction of making a vortex rotate. In order to calculate the energy reduction of rotating the condensate, we need to calculate the second term of Eq. (56), for this we need the expected angular momentum $⟨{\displaystyle \widehat{\mathit{L}}}⟩$ . Within our simple model the expected angular momentum $⟨{\displaystyle \widehat{\mathit{L}}}⟩$ can be easily calculated:

$⟨{\displaystyle \widehat{\mathit{L}}}⟩=⟨{\displaystyle \widehat{r}}\times mv(r)⟩={\displaystyle {\int }_V}{d}^{3}r\left[\left(\rho e_{\rho }\right)\times \left(mn(r)\frac{n\hslash }{m\rho }{e}_{\varphi }\right)\right]=n\pi \hslash n_{\infty }H\left(R^2-{\xi }^2\right)e_z\mathrm{.}$ E63

With Eq. (63) the energy reduction of rotating the condensate is given by:

$\text{Δ}{E}_{\text{rot}}=⟨\mathbf{\Omega }\cdot {\displaystyle \widehat{\mathit{L}}}⟩=n\pi \hslash n_{\infty }H\left(R^2-{\xi }^2\right)e_z\mathrm{,}$ E64

where we assumed that Ω is parallel to the vortex core. The fact that Ω is parallel to the vortex core is also needed because otherwise we can have circulation in the BEC without having a vortex, in the previous section we showed this was not possible.

A vortex will be formed in the BEC if it is energetically favourable; this means that the cost of making a vortex as in Eq. (62) is smaller than (or equal to) the energy profit or rotating the condensate as in Eq. (64). Setting the energy cost (62) equal to the energy profit (64) we can now calculate the critical frequency:

${\text{Ω}}_{\text{cr}}=n\frac{\hslash }{m}\frac{\ln (R/\xi )}{R^2-{\xi }^2}\mathrm{.}$ E65

The critical frequency yields the minimal rotational frequency that is required in order to get a vortex. The critical frequency is plotted as a function of the condensate size in Figure 6 .

Figure 6.

The critical frequency in Eq. (65) as a function of R / ξ.

A first thing that can be noted is that the critical frequency given in Eq. (65) is directly proportional to the number of circulation quanta present in the vortex. This means that a n = 1 vortex (one circulation quantum) is the first possible vortex that can be formed. If we now rotate twice as fast we can in principle, a vortex can be formed with n = 2 (two circulation quanta) or we can make two n = 1 vortices. Which one of the two that will be formed depends on the relative energy of both vortex configurations. In order to see which configuration will be beneficial, it is worthwhile to look at how the energy cost (62) and energy reduction (64) depend on the number of circulation quanta n. For the energy cost (62) we see that $\text{Δ}{E}_{\text{kin}}\propto n^2$ , while the energy reduction (64) goes as $\text{Δ}{E}_{\text{rot}}\propto n$ . This means that as your number of circulation quanta rises, the energy cost will rise faster than the energy reduction. Indeed if we have two n = 1 vortices, we see that the energy cost will rise with a factor two, and the energy reduction will also rise with a factor two. While if we have a n = 2 vortex the energy cost will rise with a factor four, while the energy reduction will still remain a factor two. This means that for the condensate it is more advantageous to make more vortices as you rotate faster, rather than making one big vortex. This is also seen in experiments (for example Figure 3 ), due to interactions the vortices will form lattices (Abrikosov lattices) in the BEC. Note that for both cases (many vortices of one big vortex) the critical frequency (65) will be the same.

A second aspect that can be noticed is that from our simple derivation it follows that two vortices will repel each other, while a vortex and anti-vortex will attract each other. The fact that two vortices repel each other can be seen from the fact that if we try to push them in the same place, they will behave like one vortex with the double amount of circulation quanta. As we have seen one vortex with the double amount of circulation quanta is energetically less advantageous than two separate vortices. This means that two vortices will repel each other in order to avoid an energetically less advantageous situation. For a vortex-anti-vortex pair this is different, these will attract in order to reduce the energy. If we have a vortex and an anti-vortex the total circulation equals zero, meaning that the condensate as a whole is not rotating. The fact that the condensate is not rotating means that there is no energy reduction (64), but only an energy cost (62) for making the holes. By pushing the vortex and anti-vortex together, their velocity fields will cancel each other out, leading to an annihilation of both vortices. Interactions between vortex pairs and vortex-anti-vortex pairs will be studied in more detail in the next section.

The results about vortex formation in a BEC can be summarized as:

• Energetically is more favourable to only make vortices with one circulation quantum (n = 1).

• Two vortices will repel each other.

• A vortex and anti-vortex will attract each other.

Since for a condensate it is advantageous to only make n = 1 vortices, we will assume for the remainder of the chapter that n = 1. So we will only look at vortices with one circulation quantum.

5.4. The exact structure of a single vortex

For our derivation of the vortex formation we assumed that our vortex core was a cylindrical hole. As it will turn out this is not a bad approximation as long as $R>>\xi$ , which was exactly what we assumed in our previous derivation for the critical frequency. Of course we should always check whether our assumptions are good, in order to do this we will determine the exact structure of a n = 1 vortex.

To determine the vortex structure we start from the Gross-Pitaevskii equation, Eq. (13), where we substitute the Laplacian in cylindrical coordinates $(\rho \mathrm{,}\varphi \mathrm{,}z)$ :

$-\frac{{\hslash }^2}{2m}\left[\frac{1}{\rho }\frac{\partial }{\partial \rho }\left(\frac{\partial \text{Ψ}(r)}{\partial \rho }\right)+\frac{1}{{\rho }^2}\frac{{\partial }^2\text{Ψ}(r)}{\partial {\varphi }^2}+\frac{{\partial }^2\text{Ψ}(r)}{\partial z^2}\right]+V_1(r)\text{Ψ}(r)+g{\left|\text{Ψ}(r)\right|}^2\text{Ψ}(r)=\mu \text{Ψ}(r).$ E66

In the Gross-Pitaevskii equation, Eq. (66) we can now substitute the general cylindrical symmetric wave function (34). After eliminating the phase factor $e^{i\varphi }$ we find that our equation of motion (66) becomes equal to:

$-\frac{{\hslash }^2}{2m}\left[\frac{1}{\rho }\frac{\partial }{\partial \rho }\left(\rho \frac{\partial f(\rho \mathrm{,}z)}{\partial \rho }\right)+\frac{{\partial }^{2}f(\rho \mathrm{,}z)}{\partial z^2}\right]+\left[\frac{{\hslash }^2}{2m{\rho }^2}+V_1(\rho \mathrm{,}z)\right]f(\rho \mathrm{,}z)+g{f}^3(\rho \mathrm{,}z)=\mu f(\rho \mathrm{,}z).$ E67

The extra term next to the potential is the rotational energy due to the vortex flow, note that this term fill force $f(\rho \mathrm{,}z)$ to zero as we approach the vortex core in $\rho =0$ .

The Gross-Piteavskii equation, Eq. (67) will now be solved for the uniform BEC, so $V_1(\rho \mathrm{,}z)=0$ . Since there is no more z dependency in the equation of motion (67) for the chosen potential, the structure of the vortex will also no longer depend on z. For a uniform BEC the Gross-Pitaevskii equation then becomes:

$-\frac{{\hslash }^2}{2m}\frac{1}{\rho }\frac{\partial }{\partial \rho }\left(\rho \frac{df(\rho )}{d\rho }\right)+\frac{{\hslash }^2}{2m{\rho }^2}f(\rho )+g{f}^3(\rho )=\mu f(\rho ).$ E68

We are now left with an ordinary non-linear differential equation (68) that needs to be solved. This means that the solution will be most likely given by a numerical solution. In order to easily solve the equation of motion (68) numerically, it should be made dimensionless. If we look at the limit $\rho \to \infty$ (far away from the vortex core) we find that the solution to the differential equation is given by:

$f(\rho \to \infty )=\sqrt{\frac{\mu }{g}}=\sqrt{n_{\infty }}\mathrm{,}$ E69

where the last equality follows from the solution of a homogeneous BEC in Section 4. We now have two relevant quantities to make the differential equation dimensionless, namely the value at infinity f ₀ and the healing length ξ:

$f_0=\sqrt{\frac{\mu }{g}}=\sqrt{n_{\infty }}\text{and}\xi =\frac{1}{\sqrt{8\pi n_{\infty }{a}_s}}=\frac{1}{f_0}\frac{\hslash }{\sqrt{2gm}}\mathrm{.}$ E70

Using the asymptotic value and healing length (70), we can define a dimensionless radial coordinate x and a wave function $\chi (x)$ as:

$x=\frac{\rho }{\xi }\text{and}\chi (x)=\frac{f(x\xi )}{f_0}\mathrm{.}$ E71

Substituting our dimensionless quantities (70) and variables (71) into the differential equation (68) yields:

$-\frac{1}{x}\frac{d}{dx}\left(x\frac{d\chi (x)}{dx}\right)+\frac{\chi (x)}{x^2}+{\chi }^3(x)-\chi (x)=0.$ E72

This remaining differential equation, Eq. (72) can now be solved. In order to solve the differential equation, Eq. (72) for the vortex structure, we impose the boundary conditions:

$\chi (x=0)=0\text{and}\chi (x\to \infty )=1.$ E73

The last boundary condition can also be replaced by

${\frac{d\chi (x)}{dx}|}_{x\to \infty }=\mathrm{0,}$ E74

depending on which numerical algorithm is used. In the latter case one should check whether the stable value equals one. In order to solve the equation of motion, Eq. (72) together with the boundary conditions (73) the shooting method can be used. The shooting method is a fast and popular way to solve non-linear ODEs with boundary conditions. Using the shooting method one starts in one of the boundary conditions and then ‘shoots’ out paths at different slopes. In order for a solution to become accepted, it should end within a given tolerance near the other boundary condition. At first the possible paths for $\chi (x)$ are shot at various large slopes, the range of slopes under which the paths are shot out reduces in each step. The reduction of the range of slopes is done by choosing the slopes the two closest paths near the other boundary condition. The great advantage of the shooting method is that it reduces a problem with boundary conditions to a problem with initial conditions, which is more easily solved. The shooting method is demonstrated in Figure 7 for our equation of motion, Eq. (58). The big problem with the shooting method is that for unstable (non-linear) problems (as our vortex structure for example) the algorithm is also fairly unstable. A trick to overcome this instability is to start at a low x value (e.g. at x = 3 in Figure 7 ) for comparing with the boundary condition. At lower values of x, the value of $\chi (x)$ is smaller and hence also the value of ${\chi }^3(x)$ , leading to a more stable program. Each time the interval of slopes becomes narrower, the x-value at which the solutions are compared to the boundary value increases since the solution then becomes more stable.

Figure 7.

The shooting method solution of the ODE (72) with boundary conditions (73). The eventual vortex structure is plotted as a thick red dashed line.

As we see from the above discussion solving a non-linear ODE is quite involved and unstable. A more stable method to find the vortex structure is to numerically minimize the energy functional. In polar coordinates the total energy of the BEC (expectation value of Eq. (9)) is given by:

$E[f(\rho )]={\displaystyle {\int }_V}{d}^{3}r\left[\frac{{\hslash }^2}{2m}{\left(\frac{df(\rho )}{d\rho }\right)}^2+\frac{{\hslash }^2}{2m}\frac{f^2(\rho )}{{\rho }^2}+\frac{g}{2}{f}^4(\rho )\right]\mathrm{.}$ E75

Since we have no ϕ of z dependency in the energy (75), these can be integrated out yielding an energy per height H of the condensate:

$\frac{E[f(\rho )]}{H}={\displaystyle {\displaystyle \underset{0}{\overset{R}{\int }}}}2\pi \rho d\rho \left[\frac{{\hslash }^2}{2m}{\left(\frac{df(\rho )}{d\rho }\right)}^2+\frac{{\hslash }^2}{2m}\frac{f^2(\rho )}{{\rho }^2}+\frac{g}{2}{f}^4(\rho )\right]\mathrm{.}$ E76

One thing that can already be noted is the fact that for large values of ρ we have that f(ρ) becomes constant, this will yield a logarithmic divergence of the integral. In order to get sensible results we can either renormalize the integral by subtracting the ρ → ∞ behaviour, or (as is the case in Eq. (76)) keeping the integration interval finite. In order to simplify the numerics, the energy of the uniform BEC is subtracted from the energy of the cylindrical vortex (76). For a homogeneous BEC we know from Section 4 the energy of a homogeneous BEC (41), per unit length this becomes:

$\frac{E_0[\text{Ψ}]}{H}=\frac{1}{2}g{n}_{\infty }^2\frac{V}{H}=\frac{1}{2}g{n}_{\infty }^2\pi R^2=\frac{1}{2}g{f}_0^4\pi R^2\mathrm{,}$ E77

with V the volume of the vortex core. Since the number of particles should be preserved when we go from a uniform BEC to a BEC containing a vortex, the number of particles per unit length N / H should also be preserved. The number of particles per unit length (per cross section of the cylindrical condensate) can be written as ν, related to n _∞ as:

$n_{\infty }=\frac{\nu }{\pi R^2}\mathrm{.}$ E78

Since in the system with the vortex the number of particles per unit length ν should be the same as in the uniform case, we can now calculate ν using the vortex solution:

$\nu =\frac{1}{H}{\displaystyle {\int }_V}{d}^{3}r{f}^2(\rho )={\displaystyle {\displaystyle \underset{0}{\overset{R}{\int }}}}2\pi \rho d\rho f^2(\rho )=\pi R^2{f}_0^2-{\displaystyle {\displaystyle \underset{0}{\overset{R}{\int }}}}2\pi \rho d\rho \left(f_0^2-f^2(\rho )\right)\mathrm{,}$ E79

where we added and subtracted the uniform solution in the last step. Using both Eqs. (78) and (79), the energy of a homogeneous BEC per unit length (77) can be rewritten as:

$\frac{E_0[\text{Ψ}]}{H}=\frac{1}{2}g{f}_0^4\pi R^2-g{f}_0^2{\displaystyle {\displaystyle \underset{0}{\overset{R}{\int }}}}2\pi \rho d\rho \left(f_0^2-f^2(\rho )\right)\mathrm{,}$ E80

where the term

$\frac{g}{\pi R^2}{\left({\displaystyle {\displaystyle \underset{0}{\overset{R}{\int }}}}2\pi \rho d\rho \left(f_0^2-f^2(\rho )\right)\right)}^2$ E81

was dropped. The term that was dropped is of the order

$\frac{f_0^{4}g{\xi }^4}{R^2}\mathrm{,}$ E82

Since we assume that $R>>\xi$ , this term is negligible. If we now subtract the energy of the homogeneous BEC (80) from our total energy (76) we get:

$\frac{E[f(\rho )]}{H}={\displaystyle {\displaystyle \underset{0}{\overset{R}{\int }}}}2\pi \rho d\rho \left[\frac{{\hslash }^2}{2m}{\left(\frac{df(\rho )}{d\rho }\right)}^2+\frac{{\hslash }^2}{2m}\frac{f^2(\rho )}{{\rho }^2}+\frac{g}{2}{\left(f_0^2-f^2(\rho )\right)}^2\right]\mathrm{.}$ E83

The advantage of writing the energy as Eq. (83) is that the last term now vanishes as ρ → ∞ together with the first term. The form (83) is a convenient form for numerical calculations of the vortex energy. By using our dimensionless quantities (70) and variables (71), the energy (83) can be made dimensionless yielding:

$\frac{E[\chi (x)]}{H}=\frac{\pi {\hslash }^2}{m}{n}_{\infty }{\displaystyle {\displaystyle \underset{0}{\overset{R/\xi }{\int }}}}xdx\left[{\left(\frac{d\chi (x)}{dx}\right)}^2+\frac{{\chi }^2(x)}{x^2}+\frac{1}{2}{\left(1-{\chi }^2(x)\right)}^2\right]\mathrm{.}$ E84

Note that the upper boundary of the energy integral equals $R/\xi$ , since $R>>\xi$ this upper boundary is usually (with variational models) approximated by ∞. If we now integrate to a very large value of R / ξ then the second term in (84) will have a dominant contribution. As a first approximation for the energy of the vortex we get:

$\frac{E[\chi (x)]}{H}\approx \frac{\pi {\hslash }^2}{m}{n}_{\infty }\ln \left(\frac{R}{\xi }\right)\mathrm{.}$ E85

Note that this is the energy of the cylindrical vortex (62). To get a numerically even nicer form for the energy (83) the first (asymptotic) approximation (85) can be pulled out of the integral. After rewriting the numerically most optimal form is given by:

$\begin{array}{l}\frac{E[\chi (x)]}{H}=\frac{\pi {\hslash }^2}{m}{n}_{\infty }\ln \left(\frac{R}{\xi }\right)\\ +\frac{\pi {\hslash }^2}{m}{n}_{\infty }{\displaystyle {\displaystyle \underset{0}{\overset{R/\xi }{\int }}}}xdx\left[{\left(\frac{d(1-\chi (x))}{dx}\right)}^2-\frac{1-{\chi }^2(x)}{x^2}+\frac{1}{2}{\left(1-{\chi }^2(x)\right)}^2\right]\mathrm{.}\end{array}$ E86

The main advantage of Eq. (86) is that the integral integrates over $(1-\chi (x))$ which drops to zero as x → ∞, this makes Eq. (86) numerically more stable than Eq. (84). This numerical stability is because of the fact that the logarithmic divergence was analytically removed in Eq. (86), allowing for an integration up to ∞. We can substitute the vortex exact vortex solution (shown in Figure 7 ) in the energy, the dependency as a function of R / ξ is shown in Figure 8 . By making an asymptotic fit of the form $\ln (a\times x)$ the asymptotic expansion for the exact solution was found:

Figure 8.

The energy in Eq. (86) as a function of R / ξ for the exact solution. The dashed line represents the asymptotic expansion (87).

$\frac{E[{\chi }_{\text{exact}}]}{H}=\frac{\pi \hslash }{m}{n}_{\infty }\ln \left(1.466\frac{R}{\xi }\right)\mathrm{.}$ E87

As we see in Figure 8 , the asymptotic solution of Eq. (87) becomes quite accurate quite fast (already around $R/\xi =2$ ). The fact that the asymptotic solution is accurate so fast supports the claim at the beginning of the section that the cylindrical vortex already gives quite a good estimate. In Figure 8 , we in fact see that for condensates that are large enough and do not contain too many vortices, the cylindrical solution will give good results.

5.5. Variational model for the vortex core profile

In the previous we solved for the vortex structure exactly. Once we start looking at more complicated situations this will become quite a time-consuming activity. In order to reduce the numerical work and allow for more analytical work, it is possible to use a variational model. A good variational model should be subjected to the boundary conditions (73). Three possible candidates are given by:

${\chi }_1(x)=\frac{x}{\sqrt{{\alpha }_1+x^2}}\mathrm{,}$ E88

${\chi }_2(x)=\tanh \left(\frac{x}{\sqrt{2}{\alpha }_2}\right)\mathrm{,}$ E89

${\chi }_3(x)=1-\exp \left(-\frac{x}{{\alpha }_3}\right)\mathrm{.}$ E90

In order to make the calculations as easy as possible, it is wise to choose a variational function consisting of analytical functions, so either Eq. (89) or Eq. (90). Remember from Section 4 that a condensate heals from a hole like a hyperbolic tangent (49), this might be a possible argument to choose (89) as the variational vortex function. In order to determine the values of $\left\{{\alpha }_1\mathrm{,}{\alpha }_2\mathrm{,}{\alpha }_3\right\}$ the variational models are substituted into the energy functional (86) and this energy is minimized with respect to $\left\{{\alpha }_1\mathrm{,}{\alpha }_2\mathrm{,}{\alpha }_3\right\}$ . The values for that are found for the variational parameters are given by:

${\alpha }_1=\mathrm{2,}$ E91

${\alpha }_2=\sqrt{\frac{3}{4\ln (2)-1}}\approx \mathrm{1.3,}$ E92

${\alpha }_3=\frac{12}{\sqrt{89}}\approx \mathrm{1.27.}$ E93

Given these variational parameters we can now see how they compare to the exact solution. The three variational solutions for Eqs.(88)–(90) are plotted with the exact solution in Figure 9 .

Figure 9.

The exact solution for the vortex structure together with the three variational solutions (88)–(90).

For the remainder of this chapter we will work with either the cylindrical vortex core (57) or the hyperbolic tangent (89). This choice of variational model is made because for the hyperbolic tangent it is already proven that it gives good results for the vortex core structure (88). For larger values of R / ξ, the hyperbolic tangent reduces to the cylindrical vortex, as it should since the exact solution does the same, in order to show that the hyperbolic tangent indeed reduces to the cylindrical case.

6. Several basic vortex configurations

Now that we know how reliable our variational models are and their range of applicability, we can look at three different vortex configurations:

• A vortex-antivortex system

• An off-centre vortex

• A vortex-vortex system

This will give us an idea of how stable vortices are in a BEC and how they interact. The interactions between vortices happen because the vortices feel each other’s velocity field, which leads to a Magnus force. For the first part, the easiest variational model (cylindrical hole) is used. Once the forces are known, we can also look at the energy (and thus potential) of the system in a variational manner. For this system, the variational solutions using the hyperbolic tangent are presented.

Note that although these methods might seem specific to vortices and hydrodynamics, there are a lot of other applications. An example is given by the application of the Magnus force to model the interaction of the solar wind plasma with the upper layers of planetary atmospheres. In the absence of particle-particle interactions between both plasma media (solar wind and atmosphere), the kinetic momentum of the solar wind will be transferred to the planetary particles through wave-particle interactions. These wave-particle interactions will produce vortex motion, which will cause these Magnus force interactions [70, 71].

6.1. A vortex in a velocity field: Magnus force

The interaction of vortices is because they feel each other’s velocity field. If we have a rotating system (vortex) that moves with a relative velocity with respect to a fluid, it will feel a Magnus force. A nice qualitative picture of the Magnus force is given in Figure 10 . From the image on the right side, it is intuitively clear how the Magnus force will emerge. The rotating cylinder will bend the indecent velocity field. In order to bend the velocity field, an acceleration is needed, which results in a force (second law of Newton). This force will both act on the fluid as well as on the cylinder.

Figure 10.

Left: a flow (velocity field) around a rotating (circular) body together with the resulting Magnus force (red), right: the simple qualitative picture.

In order to calculate the force on a vortex, we will use the variational model of the cylindrical vortex core. In this case the force can even be analytically calculated. The starting point for the Magnus force is the Bernoulli equation

$p(r)+\frac{1}{2}n(r)v^2(r)+\varphi (r)=\text{const}\mathrm{,}$ E94

where $\varphi (r)$ is a potential term. The Bernoulli equation (94) can be derived from the Euler equation ((24)) by integration. Note that since our vortex flow is irrotational the Bernouilli equation, Eq. (94) is valid everywhere in the BEC and not only along the same stream line. For simplicity we also look in a homogeneous BEC (no potentials) so that our potential term $\varphi (r)$ in the Bernoulli equation drops out. Taking our reference values far away from the vortex core (at infinity), the Bernoulli equation becomes:

$p(r)+\frac{1}{2}n(r)v^2(r)=p_{\infty }+\frac{1}{2}{n}_{\infty }{v}_{\infty }^2\mathrm{.}$ E95

Note that $v_{i}nfty$ is equal to the translational velocity $v_{\text{trans}}$ of the vortex. Since the velocity field $v(r)$ and density field $n(r)$ of the vortex are known the Bernoulli equation, Eq. (95) gives us a prescription for the pressure. Given the pressure we can calculate the resulting force as:

$F=-{\displaystyle {\oint }_S}dSp(r)=-{\displaystyle {\int }_V}{d}^{3}r\left(\nabla p(r)\right)\mathrm{,}$ E96

where the Gauss integration theorem was used. To calculate the Magnus force we now fill in the pressure using the Bernoulli equation, Eq. (95). The density field of our cylindrical vortex is known (57), the velocity field is given by [72]:

$v(r)=\left(2v_{\text{trans}}\sin (\varphi )+\frac{\kappa }{2\pi \xi }\right)e_{\varphi }\mathrm{,}$ E97

which is a combination of a doublet flow with a vortex flow in a uniform velocity field. For the Bernoulli equation, Eq. (95) we need the square of the velocity, leading to

$v^2(r)=4v_{\text{trans}}^2{{\displaystyle \sin }}^2(\varphi )+{\left(\frac{\kappa }{2\pi \xi }\right)}^2+2\frac{\kappa }{\pi \xi }{v}_{\text{trans}}\sin (\varphi ).$ E98

In order to calculate the Magnus force (96) we will use the first formula, expanding this in the three parts of the cylinder this gives:

$F=-{\displaystyle {\oint }_S}dSp(r)=-\left({\displaystyle {\int }_{\text{Upper plane}}}+{\displaystyle {\int }_{\text{Lower plane}}}+{\displaystyle {\int }_{\text{Mantle}}}\right)dSp(r).$ E99

The first two integrals in Eq. (99) will be equal to zero since the upper and lower parts of the cylindrical vortex are not inside of the BEC, yielding:

$F=-{\displaystyle {\int }_{\text{Mantle}}}dSp(r)=H{\displaystyle {\displaystyle \underset{0}{\overset{2\pi }{\int }}}}R\cos (\varphi )d\varphi e_x+H{\displaystyle {\displaystyle \underset{0}{\overset{2\pi }{\int }}}}R\sin (\varphi )d\varphi e_y\mathrm{,}$ E100

the e _z component is zero since this comes from the lower and upper planes which are zero. Note that constant terms (p _∞,…) do not contribute since they gets averaged out in the integration. Substituting the pressure $p(r)$ using the Bernoulli equation (95) the Magnus force becomes (100):

$F=-{\displaystyle {\oint }_S}dS\left[-\frac{1}{2}{n}_{\infty }\left(4v_{\text{trans}}^2{{\displaystyle \sin }}^2(\varphi )+\frac{2\kappa }{\pi \xi }{v}_{\text{trans}}\sin (\varphi )\right)\right]=H\xi n_{\infty }\frac{\kappa }{\xi }{v}_{\text{trans}}{e}_y\mathrm{.}$ E101

The Magnus force (per unit of length) can now be rewritten as

Note that in aerodynamics this equation is also called the Kutta-Joukowski theorem, which gives a good estimate for the lift of the vortex flow.

$\frac{F_{\text{Magnus}}}{H}=n_{\infty }\kappa v_{\text{trans}}{e}_y=n_{\infty }\kappa \times v_{\text{trans}}\mathrm{,}$ E102

where $\kappa =\kappa e_z$ and $v_{\text{trans}}=v_{\text{trans}}{e}_x$ . Note that this means that the z-axis should point along the vortex core. Given the above formula for the Magnus force, we can now look at our three different vortex configurations.

6.2. Force between a vortex-anti-vortex

The first structure that we will study is a vortex and anti-vortex that can feel each other’s velocity field. This is the simplest system that we can have since for this system the condensate does not have to be rotating (total circulation is zero). We look at a vortex-anti-vortex pair where the two vortices are separated by a distance d, a top-view sketch of this system is given in Figure 11 .

Figure 11.

A top-view sketch of the vortex-anti-vortex system.

In order to determine the force, we need to fill in the formula of the Magnus force (102). What we need to know is the circulation and the translational velocity at which the vortex is moving. The circulation is known since we only look at vortices with one vortex quantum. For the translational velocity we say that the vortices are moving in each other’s velocity field (35), this allows us to estimate $v_{\text{trans}}$ by the velocity of the other vortex. Using the above the Magnus force (102) on a vortex can be calculated as:

$\frac{F_{\text{vortex}}}{H}=n_{\infty }\left(\kappa e_z\right)\times \left(\frac{\kappa }{2\pi d}\right){{\mathit{e}}^{}}_{\varphi }=-\frac{n_{\infty }{\kappa }^2}{2\pi d}{e}_{\rho }\mathrm{.}$ E103

As we see from the force on both vortices in the vortex-anti-vortex pair (103), the force points along $-e_{\rho }$ . The fact that the force points along $-e_{\rho }$ means that both the vortex and anti-vortex are being pulled towards the centre of the condensate. Once both vortices meet in the middle, they will annihilate each other. This again verifies the fact that vortices and anti-vortices will attract each other, this result was derived in the previous section based on energy statements. With Eq. (103) we now have an idea on the factors influencing the force between the vortices.

6.3. Force on an off-centre vortex

Given the previous vortex case, we can now look at a vortex that is displaced from the origin with a distance d. A top-view sketch of the system is given in Figure 12 .

Figure 12.

A top-view sketch of our single vortex system.

Now the translational velocity of the off-centre vortex is given by the rotational velocity of the trapping potential $v=\mathbf{\Omega }\times r$ . Substituting the translational velocity into the formula for the Magnus force then yields that:

$\frac{F_{\text{vortex}}}{H}=n_{\infty }\mathit{\kappa }\times \left(\mathbf{\text{Ω}}\times r\right)=n_{\infty }\kappa \text{Ω}\rho \left[e_z\times (e_z\times e_{\rho })\right]=-n_{\infty }\kappa \text{Ω}\rho e_{\rho }\mathrm{.}$ E104

As we see from Eq. (104) the off-centre vortex is being pushed towards the centre of the vortex due to the rotational velocity.

Now there is still one problem with our setup ( Figure 12 ) that comes from the velocity field (35). Since we are working in a cylindrical trap with a hard wall potential, our velocity at the edge of the condensate cannot have a radial component (flow in or out the condensate). In order to get the correct velocity field, we should in principle solve the equations of motion again, given the new boundary conditions. However, since the velocity equation is a linear equation in this case we can apply the method of images. In order to get the right velocity field, that follows the boundary conditions, an anti-vortex is added at a distance b > R from the centre of the trap. The full situation with the image vortex is shown in Figure 13 .

Figure 13.

A top-view sketch of our single vortex system with an anti-vortex image.

To get a correct estimate of the force on the displaced vortex, we need to determine the position b of the image anti-vortex. To do this, we write down the total velocity field given by:

$v_{\text{tot}}(r)=\frac{\kappa }{2\pi }\frac{1}{{\rho }^{\prime }}{{\mathit{e}}^{\prime }}_{\varphi }-\frac{\kappa }{2\pi }\frac{1}{{{\rho }^{″}}^{}}{{\mathit{e}}^{″}}_{\varphi }\mathrm{.}$ E105

Here, the notation ${\rho }^{\prime }$ and ${{\rho }^{″}}^{}$ are the distances to the vortex and anti-vortex, respectively. The vectors ${{\mathit{e}}^{\prime }}_{\varphi }$ and ${{\mathit{e}}^{″}}_{\varphi }$ are the basis vectors $e_{\varphi }$ with the vortex and anti-vortex, respectively, in the origin of the polar coordinate system. Rewriting the above in the chosen system of coordinates $(x\mathrm{,}y\mathrm{,}z)$ in Figure 13 yields:

$\begin{array}{l}{v}_{\text{tot}}(x\mathrm{,}y)=\frac{\kappa }{2\pi }\frac{1}{\sqrt{{(x-d)}^2+y^2}}\left[\frac{-y}{\sqrt{{(x-d)}^2+y^2}}{e}_x+\frac{x-d}{\sqrt{{(x-d)}^2+y^2}}{e}_y\right]\\ -\frac{\kappa }{2\pi }\frac{1}{\sqrt{{(x-b)}^2+y^2}}[\frac{-y}{\sqrt{{(x-b)}^2+y^2}}{e}_x+\frac{x-b}{\sqrt{{(x-b)}^2+y^2}}{e}_y.]\end{array}$ E106

Our Boundary condition (no radial flow) is given by

$v_{\text{tot}}(x\mathrm{,}y)\cdot e_{\rho }=0.$ E107

Writing out the left-hand side of the boundary condition Eq. (73) using Eq. (106) then yields:

$\begin{array}{l}{v}_{\text{tot}}(x\mathrm{,}y)\cdot e_r=\frac{\kappa }{2\pi }\frac{1}{{(x-d)}^2+y^2}\left[\frac{-xy}{\sqrt{x^2+y^2}}+\frac{(x-d)y}{\sqrt{x^2+y^2}}\right]\\ -\frac{\kappa }{2\pi }\frac{1}{{(x-b)}^2+y^2}\left[\frac{-xy}{\sqrt{x^2+y^2}}+\frac{(x-b)y}{\sqrt{x^2+y^2}}\right]\mathrm{.}\end{array}$ E108

For the boundary condition Eq. (108) should equal zero, using the fact that on the boundary $x^2+y^2=R^2$ yields that:

$0=d{b}^2-\left[R^2+d^2\right]b+d{R}^2\mathrm{.}$ E109

Solving the second order equation in b then yields two possible solutions:

$b_{+}=d\mathrm{,}$ E110

$b_{-}=\frac{R^2}{d}\mathrm{.}$ E111

The first solution is an anti-vortex in the same place as the vortex, hence an annihilation of the vortex. The second solution is a vortex outside of the condensate (since $d\le R$ ). The resulting velocity field then becomes:

$\begin{array}{l}{v}_{\text{tot}}(x\mathrm{,}y;d\mathrm{,}R)=\frac{\kappa }{2\pi R}\frac{1}{{(x/R-d/R)}^2+{(y/R)}^2}\left[-\frac{y}{R}{e}_x+\left(\frac{x}{R}-\frac{d}{R}\right)e_y\right]\\ -\frac{\kappa }{2\pi R}\frac{1}{{(x/R-R/d)}^2+{(y/R)}^2}\left[-\frac{y}{R}{e}_x+\left(\frac{x}{R}-\frac{R}{d}\right)e_y\right]\mathrm{.}\end{array}$ E112

The size of the velocity field $v_{\text{tot}}=\sqrt{v_{\text{tot}}\cdot v_{\text{tot}}}$ is given by:

$v_{\text{tot}}(x\mathrm{,}y;d\mathrm{,}R)=\frac{\kappa }{2\pi R}\sqrt{\frac{{(d/R-R/d)}^2}{[(x/R-d/R{)}^2+{(y/R)}^2][(x/R-R/d{)}^2+{(y/R)}^2]}}\mathrm{.}$ E113

Note that in the limit d → 0 we get our single vortex case back, and in the case d → R the vortex and anti-vortex annihilate yielding a zero-velocity field.

Using the above result, we see that in order to get a correct velocity field, we should add an anti-vortex. The addition of the anti-vortex means that our Magnus force gets an extra contribution, pulling the vortex out of the condensate. Since the distance between the vortex and anti-vortex is given by:

$a=\frac{R^2}{d}-d=\frac{R^2-d^2}{d}\mathrm{,}$ E114

the extra term in the Magnus force becomes equal to:

$\frac{F_{\text{Vortex}}}{L}=\frac{n_{\infty }{\kappa }^2}{2\pi a}{e}_x=\frac{n_{\infty }{\kappa }^2}{2\pi }\frac{d}{R^2-d^2}{e}_x\mathrm{.}$ E115

Combining both forces on the vortex (rotating condensate and anti-vortex), the total force on the vortex becomes equal to:

$\frac{F_{\text{total}}}{L}=\left(\frac{n_{\infty }{\kappa }^2}{2\pi }\frac{d}{R^2-d^2}-n_{\infty }\kappa \text{Ω}d\right)e_x=n_{\infty }\kappa d\left(\frac{\kappa }{2\pi }\frac{1}{R^2-d^2}-\text{Ω}\right)e_x\mathrm{.}$ E116

In the total force we see a competition between the two forces, by looking for the positions of the vortex d where the force equals zero, we can find the equilibrium positions. Setting Eq. (116) to zero then yields two solutions:

$d_1=\mathrm{0,}$ E117

$d_2=R\sqrt{1-\frac{\kappa }{2\pi R^2\text{Ω}}}\mathrm{.}$ E118

The first solution yields a stable position at the centre of the vortex, while the second solution gives an equilibrium outside of the centre. Note that the second solution only exists if the condensate is rotating fast enough. The second solution is thus some kind of maximum (barrier) needed to keep the vortex inside of the condensate when it is rotating. The minimal velocity needed to keep the vortex inside is given by:

${\text{Ω}}_{\text{min}}=\frac{\hslash }{m{R}^2}\mathrm{.}$ E119

Note that this frequency (119) is almost the same as the critical frequency (65). The big difference with the critical frequency is that there we have a length scale (the healing length) that also plays a role. The reason for this difference in frequency is because of the fact that when the Magnus force (102) is calculated, the typical length scale of the vortex (healing length) gets eliminated. This means that the obtained results are excellent as long as the vortices are far enough away from each other. When the vortices do get close to each other, for example if the off centre vortex reaches the edge, we need a more precise (variational or exact) calculation that also takes the vortex structure into account.

Using the total force on the vortex (116) a potential can also be defined as:

$\frac{U}{L}=-{\displaystyle {\displaystyle \underset{0}{\overset{d}{\int }}}}\left(\frac{F_{\text{total}}}{L}\right)d(d)e_{\rho }=n_{\infty }\kappa \left(\frac{\kappa }{4\pi }\ln \left(\frac{R^2-d^2}{R^2}\right)+\text{Ω}\frac{d^2}{2}\right)\mathrm{.}$ E120

This potential (120) is also plotted in Figure 14 . As we see in Figure 14 , as long as the condensate spins fast enough, a vortex will always be pushed towards the centre.

Figure 14.

The resulting potential for the off-centre vortex. The black dashed line and dots yield the second stationary position of the off-centre vortex (118).

6.4. Vortex-vortex

The last situation we can look at is the vortex-vortex interaction. As we know from the previous section two vortices will repel each other, meaning that now the centre of the trap will not yield a stable solution. As we know from the off-centre vortex, we want our boundary condition to hold Eq. (107), which means that each vortex should have an anti-vortex. A sketch of the two-vortex system is given in Figure 15 .

Figure 15.

A top-view sketch of the two vortex system with the corresponding anti-vortices.

As earlier, we can now write down the Magnus force of the three different vortices on one of the vortices in the condensate. Since the system is symmetric, it does not matter which vortex is chosen for the calculation; we chose the right one, yielding:

$\frac{F_{\text{vortex}}}{L}=\frac{n_{\infty }{\kappa }^2}{2\pi }\left(\frac{d}{R^2-d^2}+\frac{1}{2d}-\frac{d}{R^2+d^2}\right)e_x=\frac{n_{\infty }{\kappa }^2}{2\pi }\frac{R^4+3d^4}{2d(R^4-d^4)}{e}_x\mathrm{.}$ E121

As we see from Eq. (121) the vortices push each other out of the condensate. Since the vortices give a total circulation to the condensate, we should also have a rotation of the superfluid.

Note that as before this rotation is needed to keep the vortices inside of the condensate.

Adding the attractive rotational force then yields for the total force on the (right) vortex that:

$\frac{F_{\text{total}}}{L}=n_{\infty }\kappa \left(\frac{\kappa }{2\pi }\frac{R^4+3d^4}{2d(R^4-d^4)}-\text{Ω}d\right)e_x\mathrm{.}$ E122

Again we can look for stable positions of the vortex pair by solving the equation

$0=\frac{\kappa }{2\pi }\frac{R^4+3d^4}{2d(R^4-d^4)}-\text{Ω}\mathrm{.}$ E123

This equation (123) is analytically solvable by substituting the variable y = d ² and solving the remaining equation of third degree by Cardano’s method. This will however lead to a large (non-insightful) equation. For simplicity we assume that $d<<R$ (vortices close to the centre), this then yields for the equilibrium positions that:

$\frac{\kappa }{4\pi }=\text{Ω}{d}^2\Rightarrow d=\sqrt{\frac{\kappa }{4\pi \text{Ω}}}=\sqrt{\frac{\hslash }{2m\text{Ω}}}\mathrm{.}$ E124

The equilibrium position, both exact and the above approximation (124), is plotted in Figure 16 . As we see, the faster the condensate rotates, the more the vortices are pushed together. In the limit $d<<R$ a very high rotational frequency is needed in order to get accurate results. The limit $d<<R$ corresponds to ignoring the two image vortices, and these get pushed to infinity as d decreases. In principle we can again calculate a critical frequency as before, but this would require the full formula and some tedious calculating work that is left for the reader. A nice estimate of the critical frequency can also be given by the plot in Figure 16 by looking at the frequency for which the solid curve stops, this yields:

Figure 16.

The equilibrium position of a vortex-vortex pair.

$\frac{R^2\text{Ω}}{\kappa }\approx \mathrm{0.351.}$ E125

Notice that this is slightly higher than the critical frequency we would expect from the discussion in Section 5, where we would have:

$\frac{R^2\text{Ω}}{\kappa }=\frac{1}{\pi }\approx \mathrm{0.318.}$ E126

This overestimation is again an artefact of the fact that the vortex core structure is ignored.

Again it is possible to define a potential which yields:

$\frac{U}{L}=-{\displaystyle {\displaystyle \underset{R/2}{\overset{d}{\int }}}}\left(\frac{F_{\text{total}}}{L}\right)d(d)e_x=n_{\infty }\kappa \left[\frac{\kappa }{4\pi }\ln \left(\frac{R/2}{d}\frac{R^4-d^4}{R^4-{(R/2)}^4}\right)+\frac{\text{Ω}}{2}\left(d^2-{\left(\frac{R}{2}\right)}^2\right)\right]\mathrm{.}$ E127

This potential is plotted for two different rotational frequencies in Figure 17 . In Figure 17 , we see that the vortices indeed repel each other but stay in the trap if the trapping frequency is high enough.

Figure 17.

The potential for a vortex-vortex pair for two different frequencies. The solid curves yield the full solution, the dashed curves are the limit d ≪ R.

6.5. Variational calculations for the energy functional

From the previous sections we now know four variational models for the vortex structure $\chi (r)$ given by Eqs. (57), (88)–(90). Given one of these variational, we can write the variational macroscopic wave function for a vortex or anti-vortex as:

$\text{Ψ}(r)=\sqrt{n_{\infty }}\chi (r)e^{iS(r)}\mathrm{.}$ E128

The above macroscopic wave function (128) is not limited to one-vortex states. We can also construct a structure function $\chi (r)$ for a multi-vortex state. Assume we have a total of $N_v$ vortices, then the variational structure function is written as

The reason that we have to multiply the vortex solutions for the structure rather than add them has several reasons. First of all the particle density is given by the macroscopic wave function (and thus of the structure function), if we have a sum here this will yield extra terms. Second the Gross-Pitaevskii equation, Eq. (10) is a non-linear equation, which means that a sum of two solutions is not a solution of the differential equation. By multiplying the structure functions we create one state (solution) with multiple vortices.

$\chi (r)={\displaystyle {\displaystyle \prod _{j=1}^{N_v}}}{\chi }_j(r),$ E129

where ${\chi }_j(r)$ are the one-vortex structure functions. Usually one takes the hyperbolic tangent (89) for the one-vortex structure, resulting for Eq. (129) in:

$\chi (r)={\displaystyle {\displaystyle \prod _{j=1}^{N_v}}}\tanh \left(\frac{|r-r_j|}{\alpha \sqrt{2}\xi }\right)\text{with}\alpha =\sqrt{\frac{3}{4\ln (2)-1}}\approx \mathrm{1.3,}$ E130

with $r_j$ the positions for the vortices. The velocity field for one vortex in Eq. (35) can be rewritten in terms of a circulation vector $\mathbf{\kappa }=\kappa e_z$ and a position vector $r=\rho e_{\rho }$ . Note that the z-axis is oriented parallel to the rotational velocity Ω. Using the circulation vector κ and position vector r the velocity field (27) for a vortex or anti-vortex becomes:

$v(r)=\frac{\pm \hslash }{m\rho }{e}_{\varphi }=\frac{\kappa }{2\pi \rho }{e}_{\varphi }=\frac{\kappa }{2\pi \rho }(e_z\times e_{\rho })=\frac{\mathbf{\kappa }\times r}{2\pi r^2}\mathrm{.}$ E131

where the sign of κ depends on whether you have a vortex κ > 0 or anti-vortex κ < 0. If we look at a multi-vortex structure with $N_v$ vortices, then the velocity field is given by

In the case of velocity fields we take the sum. This also has two reasons. First of all we know that each vortex hole gives rise to such a velocity field (35). Second in the case of the cylindrical hole we have an irrotational and incompressible flow for which the velocity follows a linear differential equation. So for the cylindrical vortex hole we can add velocity solutions thanks to linearity. For other vortex structures this variational model is again approximate, but appropriate for the chosen variational model for the density (129).

:

$v(r)={\displaystyle {\displaystyle \sum _{j=1}^{N_v}}}\frac{\mathbf{\kappa }\times (r-r_j)}{2\pi |r-r_j{|}^2}\mathrm{.}$ E132

Since we are going to look at stability of the system, we will again add the rotational term according to Eq. (56). For the off-centre vortex and vortex-vortex structures, this will be needed. For the vortex-anti-vortex this will not be needed since the condensate will not be rotating then.

Finally since we will be working in polar coordinates, our structure functions $\chi (r)$ and velocity field should be converted to polar coordinates. Since we will be looking at one- and two-vortex systems we can make the calculations simpler by choosing the x-axis of the Cartesian coordinate system such that all of the vortices lie on the x-axis. For the structure function with the hyperbolic tangent given in Eq. (89) we get in polar coordinates for a vortex in $(x\mathrm{,}y)=(d\mathrm{,0})$ :

$\begin{array}{c}\tanh \left(\frac{|r-d{e}_x|}{\alpha \sqrt{2}\xi }\right)=\tanh \left(\frac{\sqrt{{(\rho \cos (\varphi )-d)}^2+{(\rho \sin (\varphi ))}^2}}{\alpha \sqrt{2}\xi }\right)\\ =\tanh \left(\frac{\sqrt{{\rho }^2-2\rho d\cos (\varphi )+d^2}}{\alpha \sqrt{2}\xi }\right)\end{array}$ E133

Of course in order to do a variational calculation we need an energy functional, this is given by:

$E[n(r),v(r)]={\displaystyle {\int }_V}{d}^{3}r\left(\frac{1}{2}mn(r)v^2(r)+n(r)V_{\text{trap}}(r)+g{n}^2(r)\right)\mathrm{.}$ E134

In this chapter we look at a homogeneous condensate ( $V_{\text{trap}}(r)$ ). For the variational calculations, we also assumed that the vortices are small compared to the size of the condensate ( $R>>\xi$ ), if that is the case then our interaction term is constant and can be neglected. Finally we need to add the rotational energy, doing this results in the following form for the rotational energy:

$E[n(r),v(r)]={\displaystyle {\int }_V}{d}^{3}r\left(\frac{1}{2}mn(r)v^2(r)\right)-\mathbf{\text{Ω}}\cdot {\displaystyle {\int }_V}dr\left(r\times (mn(r)v(r))\right)\mathrm{.}$ E135

Choosing the rotational velocity along the z-axis is possible to do the z integration resulting in:

$\frac{E[n(r),v(r)]}{mL}=\frac{1}{2}{\displaystyle {\int }_S}{d}^{2}r\left(n(r)v^2(r)\right)-\text{Ω}{\displaystyle {\int }_S}{d}^{2}r\left|r\right|n(r)|v(r)|.$ E136

The above energy can be used for further variational calculations; the calculations themselves are left to the reader.

6.6. Variational results for the energy of the systems

Given the variational energy (136) it is possible to look at our three different vortex cases again. To conclude the section the different energies for the three-vortex systems above will be calculated variationally. For the different calculations a trap of the size R = 30 was chosen with a healing length of ξ = 1. As a variational model, the hyperbolic tangent was chosen.

For the vortex pair the energy as a function of the distance between the vortices is given in Figure 18 . As we see the vortex-anti-vortex pair attracts as well in the variational model so our results from the simple calculations on the Magnus force still hold. For the off-centre vortex and vortex pair, a rotating condensate is needed. This means that the two integrals of Eq. (136) should be integrated. The first integral yields the energy of the stationary condensate, the second integral yields the energy reduction of rotating the condensate. The total energy is then given by a competition between both energies.

Figure 18.

The variational energy for a vortex-anti-vortex pair the solid line gives the full solution, whereas the dashed line yields the asymptotic solution.

The potentials for the off-centre vortex and vortex-vortex systems are given in Figure 19 . As seen on Figure 19 , a repulsive and attractive energy are in competition. Depending on the values of the rotational frequency we will have different situations. The situation for the off-centre vortex is given for a few values of the rotational frequency Ω in Figure 20 . What we can notice in Figure 20 is that for a single vortex we find a non-centred vortex in the stable situation. This is because of the fact that the variational velocity field of the cylindrical model was used for the hyperbolic tangent, yielding small deviations. For the vortex-vortex case the same can be done, resulting in the variational energies as given in Figure 21 . Here we again see that as the rotational velocity increases, the vortices get pushed more towards the centre of the trap, while if the rotational velocity decreases below a certain value the vortices get pushed out of the trap.

Figure 19.

The kinetic (a) and rotational (b) energy for the off-center vortex.

Figure 20.

The variational energy for the off-centre vortex for several values of Ω. The dashed lines trace the minima and maxima of the potential.

Figure 21.

The variational energy for the vortex-vortex for several values of Ω. The dashed lines trace the minima and maxima of the potential.

7. The Berezinskii-Kosterlitz-Thouless phase transition

To conclude the discussion about vortices in BECs we will look at a vortex induced phase transition, the Berezinskii-Kosterlitz-Thouless (BKT) phase transition (Nobel Prize in Physics 2016).

7.1. Vortex entropy

Vortices carry entropy. The entropy of a given state can be computed from the free energy,

$S={\frac{\partial F}{\partial T}|}_{V\mathrm{,}N}\mathrm{,}$ E137

but in this simplified presentation we will use some qualitative arguments to estimate the entropy. The ‘ size’ of a vortex is given by its healing length ξ, and the number of ways that an object of this size can be placed in a condensate of radius R is estimated by ${\xi }^2/R^2$ , where we consider either a two-dimensional condensate, or we enforce the vortex lines to be parallel. An illustration of a possible microstate—a way of placing vortices—is shown in Figure 22 . With this simple argument we can estimate the entropy of one or more vortices and anti-vortices, again by counting the number Ω of microstates, and using Boltzmann’s formula $S=k_B\log \text{Ω}\mathrm{.}$ For example, the entropy of a single vortex becomes:

Figure 22.

A schematic view of a possible microstate contributing to the vortex positional entropy.

Figure 23.

A qualitative picture of the BEC-BCS crossover, figure adapted from Ref. [83].

$S_{\text{1vortex}}=k_B\ln ({\text{Ω}}_{\text{1vortex}})\approx k_B\ln \left[{\left(\frac{R}{\xi }\right)}^2\right]=2k_B\ln \left(\frac{R}{\xi }\right)\mathrm{.}$ E138

When two vortices are present, or a vortex and an anti-vortex are present, we consider them to be independent. From the previous section, it is clear that this is an approximation that fails when they are close together. However, the number of microstates where this happens is small, and we can still assume the entropy is extensive in the number of vortices:

$S_{\text{2vortices}}=2S_{\text{1vortex}}=4k_B\ln \left(\frac{R}{\xi }\right)\mathrm{.}$ E139

When information (or order) is added to the system, the entropy must decrease. For example, we might know that the vortex and anti-vortex are bound, and must be in adjacent patches of Figure 22 . Then the number of microstates that are compatible with this constraint is again proportional to ${\xi }^2/R^2$ , so that up to a constant, negligible for $R\gg \xi$ ,

$S_{\text{vort}\exp \text{air}}\approx 2k_B\ln \left(\frac{R}{\xi }\right)=\frac{1}{2}{S}_{\text{2vortices}}\mathrm{.}$ E140

7.2. BKT temperature

In a three-dimensional Bose gas, superfluidity vanishes above the critical temperature for Bose-Einstein condensation. In the normal state, phase correlations vanish exponentially fast as a function of distance, and phase coherence is absent. In the BEC regime, the phase is locked by long-range off-diagonal order [73], and the phase correlation function is a constant as a function of distance. However, this is not the case for a two-dimensional Bose gas. Due to quantum fluctuations, there is no phase transition to the Bose-Einstein condensation (more accurately: the critical temperature is zero). Nevertheless, this does not mean that phase coherence is lost altogether: the phase correlation function gets a power-law decay rather than the exponentially fast decay in the normal state [74]. This regime is called ‘quasi-condensation’ and allows for superfluidity.

Obviously, the mechanism whereby superfluidity is destroyed must be different. Berezinskii [75] and independently Kosterlitz and Thouless [76, 77] proposed that phase coherence gets destroyed by the unbinding of vortex-anti-vortex pairs, leading to a sudden proliferation of vortices and anti-vortices that scramble the phase coherence. Here, the results for the interaction energy between vortices and anti-vortices calculated above will be of use. The binding energy of a vortex-anti-vortex pair can be defined as the energy it takes to move a vortex and an anti-vortex that form a pair, away from each other to infinity.

Heuristically, the Berezinskii-Kosterlitz-Thouless (BKT) transition can be thought of as a competition between the internal energy and the entropy for vortices. The realized state is that which minimizes the free energy

$F=U-TS\mathrm{.}$ E141

When the temperature is low T ≈ 0, the internal energy (here, the binding energy of the vortex-anti-vortex pair) dominates. In this case, the pairs remain bound and no uncompensated phase singularities perturb the coherence. Indeed, as can be seen from Figure 18 , the internal energy is minimal when the vortex and anti-vortex are co-located. The phase windings then cancel out, and U = 0, F = 0. When the temperature is increased, the entropy term TS can become large enough to dominate the internal energy, leading to a transition to the more disordered state. Focusing on the vortex-anti-vortex pair, the disordered state (the state with highest entropy) is that in which the vortex and the anti-vortex are unbound and can be placed freely anywhere in the available volume. In that disordered state, we will also take U to be equal to twice the kinetic energy as in Eq. (62)—the number of microstates where the vortex is close to the anti-vortex is negligible in comparison to the number of microstates when they are well separated, so we can neglect the interaction. This allows to approximate the internal energy of the disordered state by

$U=2\pi n_{\infty }\frac{{\hslash }^2}{m}H\ln \left(\frac{R}{\xi }\right)\mathrm{.}$ E142

Comparing the free energy of the state with bound pairs (F = 0) with the free energy of the disordered state yields the critical temperature where superfluidity vanishes:

$0=U-T_{C}S\mathrm{.}$ E143

Using our earlier results (142) for U and (139) for S, the critical temperature is found to be

$T_C=\frac{2\pi n_{\infty }\frac{{\hslash }^2}{m}H\ln (R/\xi )}{4k_B\ln (R/\xi )}=\frac{\pi n_{\infty }{\hslash }^{2}H}{2m{k}_B}\mathrm{.}$ E144

In the units used for (143), we can see that the right-hand side indeed describes a temperature, and that it is determined by the superfluid density n _∞. Above this temperature, it is advantageous for the free energy to spontaneously create and unbind vortex-anti-vortex pairs that scramble the phase and turn the gas into the normal state. As mentioned in the introduction, the BKT transition was observed experimentally [40]. The observed critical phase space density differs from the back of the envelope estimate given here by a factor of about four, mostly due to the fact that the experimental condensate does not have a uniform density, and also due to finite size effects.

8. Concluding section

With this chapter an introduction to vortices in ultra-cold gases was given. The main focus of the discussion was aimed towards atomic gases of bosons. However the derived (qualitative) results are also applicable to fermionic quantum gases. In this concluding section the main results of the chapter are summarized together with a first look at a condensate of fermionic atoms.

8.1. Conclusions

Concerning vortices in superfluid atomic gases, the following conclusions can be made:

• A vortex will always make a hole that extends throughout the entire condensate.

• The circulation of a vortex is quantized in multiples of h / m see Eq. (32).

• Vortices in a BEC will always carry one quantum of circulation, this is energetically more favourable.

• For vortices a minimal rotational frequency (critical frequency) is required as given in Eq. (65).

• Two vortices will repel each other.

• A vortex and anti-vortex will attract each other.

• Rotating a BEC faster will push vortices to the centre of the trap more, for a single vortex.

• A single vortex in a trap will be most stable in the centre of the trap, for a single vortex this means that it is stable in the centre of the trap.

• For a 2D quasi-condensate superfluidity is destroyed by the unbinding of vortex-anti-vortex pairs that destroy the overall phase coherence. This happens at a certain critical temperature T _C .

8.2. One step further: fermionic gases

This chapter was restricted to superfluidity in Bose gases since these already (qualitatively) describe the behaviour of vortices in general superfluids very well. If we are working with fermionic gases a few things change, these are briefly discussed in this section. More elaborate discussions about fermionic quantum gases can be found in references [78–81].

First of all when we have fermions, we have to take the ‘spin’. With spin we label the different quantum states of our atom, in fermionic atomic gases these will be two different hyperfine states of the considered atom. These spins will be called ‘up’ and ‘down’ as an analogy with the spin known from standard quantum mechanics. In order to fix the number of particles in both spin states, we need two chemical potentials labelled μ _↑ and μ _↓. If we now say that we have an ultra-cold gas and thus have s-wave scattering, the Pauli exclusion principle tells us that only opposite spins can interact. Since the s-scattering wave function is symmetric, the spin wave function should be the anti-symmetric singlet state. Applying this to our Hamiltonian we get:

${\displaystyle \widehat{H}}={\displaystyle \int }{d}^{3}r({\displaystyle {\displaystyle \sum _{\sigma =\uparrow \mathrm{,}\downarrow }}}{{\displaystyle \widehat{\psi }}}_{\sigma }^{†}(r\mathrm{,}t)\left[-\frac{{\hslash }^2{\nabla }^2}{2m}-{\mu }_{\sigma }\right]{{\displaystyle \widehat{\psi }}}_{\sigma }(r\mathrm{,}t)+\frac{1}{2}{\displaystyle {\displaystyle \sum _{\sigma =\uparrow \mathrm{,}\downarrow }}}g{{\displaystyle \widehat{\psi }}}_{\sigma }^{†}(r\mathrm{,}t){{\displaystyle \widehat{\psi }}}_{\neg \sigma }^{†}(r^{\prime }\mathrm{,}t){{\displaystyle \widehat{\psi }}}_{\neg \sigma }(r^{\prime }\mathrm{,}t){{\displaystyle \widehat{\psi }}}_{\sigma }(r\mathrm{,}t)).$ E145

Given the above Hamiltonian (144) atomic Fermi gases can be studied. In order to simplify calculations, one usually defines the average chemical potential μ and polarization ζ as:

$\mu =\frac{{\mu }_{\uparrow }+{\mu }_{\downarrow }}{2}\mathrm{,}$ E146

$\zeta =\frac{{\mu }_{\uparrow }-{\mu }_{\downarrow }}{2}\mathrm{.}$ E147

Of course the question we can ask now is ‘what is the criterion for superfluidity now’? In Section 3, this criterion was given by the Penrose-Onsager criterion which states that the one-body density matrix should have a macroscopic eigenvalue, indicating that there was one state being occupied macroscopically. For fermions, the criterion for bosons was expanded by C.N. Yang, hence named the Penrose-Onsager-Yang criterion. Yang stated that for fermions the second order density matrix should have a macroscopically large eigenvalue, indicating a condensation of fermion pairs. These fermion pairs are also referred to as Cooper pairs, named after L.N. Cooper who first proposed this for electrons [82]. Once this is assumed it is possible to rewrite the entire system in terms of a macroscopic two-body wave function $\text{Ψ}(y\mathrm{,}x)$ . This macroscopic wave function is usually rewritten by assuming the following product form:

$\text{Ψ}(y\mathrm{,}x)=\psi (R)\chi (\mathit{\rho }),$ E148

where R is the centre of mass coordinate and ρ is the relative position. So in the end $\psi (R)$ will describe superfluidity of the Cooper pairs and $\chi (\mathit{\rho })$ will describe the interactions between the atoms in a different spin state.

The study of superfluid fermionic atoms opens up a few extra parameters to play with in experiments, which can help explore the physics of ultra-cold gases. The two most prominent are:

Spin-imbalance: A first parameter to play with is the number of ‘spin-up’ and ‘spin-down’ atoms. If the numbers are unequal (or if the polarisation ζ in Eq. (146) is non-zero), then you will have excess atoms when the condensation occurs. These excess atoms will need to go somewhere. Usually excess atoms are pushed towards the edge of the condensate, or placed inside holes in the condensate, for example vortices. When the number of ‘spin-up’ and ‘spin-down’ atoms is unequal one also talks about pair frustration since the unbound (frustrated atoms) will still try to pair up with other paired atoms.

BEC-BCS crossover: A second parameter that can be played with is the interaction strength a _s , this can also be done with bosons but is more interesting to do with fermions. With fermions it is possible to change the sign of the interaction from negative (repulsive) to positive (attractive). By doing this it is possible to go from a gas of interacting fermions (Cooper pairs or BCS), to a gas of interacting bosons (BEC) by making the interactions so strong that the particles are quasi bosons. The point at which the interaction changes in sign is also referred to as unitatiry. A qualitative picture of the transition form a repulsive interaction (BCS) to an attractive interaction (BEC) is given in Figure 23 .

Acknowledgments

The authors gratefully acknowledge the extensive, weekly discussions with S.N. Klimin, G. Lombardi, W. Van Alphen and J.P.A. Devreese. Financial support was provided by the Fund for Scientific Research Flanders, through FWO projects: G042915N (Superfluidity and super-conductivity in multi-component quantum condensates) and G011512N (Quantum turbulence in atomic and solid state Bose-Einstein condensates). One of us (N.V.) acknowledges a Ph.D. fellowship of the University of Antwerp (UAntwerp).