Dark solitons in atomic Bose-Einstein condensates: from theory to experiments

# Dark solitons in atomic Bose-Einstein condensates: from theory to experiments

## Abstract

This review paper presents an overview of the theoretical and experimental progress on the study of matter-wave dark solitons in atomic Bose-Einstein condensates. Upon introducing the general framework, we discuss the statics and dynamics of single and multiple matter-wave dark solitons in the quasi one-dimensional setting, in higher-dimensional settings, as well as in the dimensionality crossover regime. Special attention is paid to the connection between theoretical results, obtained by various analytical approaches, and relevant experimental observations.

Topical Review
###### pacs:
03.75.Lm, 03.75.Kk, 05.45.Yv
1

## 1 Introduction

A dark soliton is an envelope soliton having the form of a density dip with a phase jump across its density minimum. This localized nonlinear wave exists on top of a stable continuous wave (or extended finite-width) background. Dark solitons are the most fundamental nonlinear excitations of a universal model, the nonlinear Schrödinger (NLS) equation with a defocusing nonlinearity, and, as such, they have been studied in diverse branches of physics. Importantly, apart from a vast amount of literature devoted to relevant theoretical works, there exist many experimental results on dark solitons, including the observation of optical dark solitons, either as temporal pulses in optical fibers [1, 2], or as spatial structures in bulk media and waveguides [3, 4], the excitation of a non-propagating kink in a parametrically-driven shallow liquid [5], dark soliton standing waves in a discrete mechanical system [6], high-frequency dark solitons in thin magnetic films [7], dissipative dark solitons in a complex plasma [8], and so on.

Theoretical studies on dark solitons started as early as in 1971 [9] in the context of Bose-Einstein condensates (BECs). In particular, in Ref. [9], exact soliton solutions of the Gross-Pitaevskii (GP) equation (which is a variant of the NLS model) [10] were found and connected, in the small-amplitude limit, with the solitons of the Korteweg-de Vries (KdV) equation. Later, and shortly after the integration of the focusing NLS equation [11], the defocusing NLS equation was also shown [12] to be completely integrable by means of the Inverse Scattering Transform (IST) [13]; this way, single as well as multiple dark soliton solutions of arbitrary amplitudes were found analytically. The IST approach allowed for an understanding of the formation of dark solitons [14, 15, 16, 17, 18, 19], the interaction and collision between dark solitons [12, 20] (see also Refs. [21, 22, 23, 24, 25] and [26] for relevant theoretical and experimental studies, respectively), and paved the way for the development of perturbation methods for investigating their dynamics in the presence of perturbations [25, 27, 28, 29, 30, 31, 32]. From a physical standpoint, dark solitons were mainly studied in the field of nonlinear optics — from which the term “dark” was coined. The first theoretical work in this context, namely the prediction of dark solitons in nonlinear optical fibers at the normal dispersion regime [33], was subsequently followed by extensive studies of optical dark solitons [34, 35].

A new era for dark solitons started shortly after the realization of atomic BECs [36, 37, 38]; this achievement was awarded the Nobel prize in physics of 2001 [39, 40], and has been recognized as one of the most fundamental recent developments in quantum and atomic physics over the last decades (see, e.g., the books [41, 42] for reviews). In an effort to understand the properties of this exciting state of matter, there has been much interest in the macroscopic nonlinear excitations of BECs (see reviews in Refs. [43, 44]). In that regard, the so-called matter-wave dark solitons, were among the first purely nonlinear states that were experimentally observed in BECs [45, 46, 47, 48, 49].

The interest on matter-wave dark solitons is not surprising due to a series of reasons. First of all, for harmonically confined BECs, these structures are the nonlinear analogues of the excited states of a “prototype” quantum system [50, 51], namely the quantum harmonic oscillator [52]. On the other hand, the topological nature of matter-wave dark solitons (due to the phase jump at their density minimum) renders them a “degenerate”, one-dimensional (1D) analogue of vortices, which are of paramount importance in diverse branches of physics [53]. Additionally, and perhaps more importantly, matter-wave dark solitons are — similarly to vortices [54, 55, 56] — quite fundamental structures arising spontaneously upon crossing the BEC phase-transition [57, 58], with properties which may be used as diagnostic tools probing the rich physics of a purely quantum system (BEC) at the mesoscale [59]. Finally, as concerns applications, it has been proposed that the dark soliton position can be used to monitor the phase acquired in an atomic matter-wave interferometer in the nonlinear regime [60, 61] (see also relevant experiments of Refs. [62, 63] devoted to atom-chip interferometry of BECs).

The early matter-wave dark soliton experiments, as well as previous works on dark solitons in optics, inspired many theoretical efforts towards a better understanding of the stability, as well as the static and dynamical properties of matter-wave dark solitons. Thus, it is probably not surprising that a new series of experimental results from various groups have appeared [64, 65, 66, 67, 68, 69, 70, 71], while still other experiments — not directly related to dark solitons — reported observation of these structures [62, 63, 72]. These new, very recent, experimental results were obtained with an unprecedented control over the condensate and the solitons as compared to the earlier soliton experiments. Thus, these “new age” experiments were able not only to experimentally verify various theoretical predictions, but also to open new exciting possibilities. Given this emerging interest, and how new experiments in BEC physics inspire novel ideas — both in theory and in experiments — new exciting results are expected to appear.

The present paper aims to provide an overview of the theoretical and experimental progress on the study of dark solitons in atomic BECs. The fact that there are many similarities between optical and matter-wave dark solitons [73], while there exist excellent reviews on both types of dark solitons (see Ref. [34] for optical dark solitons and Ch. 4 in Ref. [43] for matter-wave dark solitons), provides some restrictions in the article: first, the space limitations of the article, will not allow for an all-inclusive presentation; in that regard, important entities — relevant to dark solitons — such as vortices [53, 74, 75] and vortex rings [76, 77] will only be discussed briefly. In fact, this review (which obviously entails a “personalized” perspective on dark solitons) will cover the basic theory emphasizing, in particular, to the connection between the theoretical results and experimental observations; this way, in most cases, theoretical discussion will be immediately followed by a presentation of pertinent experimental results. In that regard, it is also relevant to note that our theoretical approach will basically be based on the mean-field theory: as will be shown, the latter can be used as a basis of understanding of most effects and experimental findings related to matter-wave dark solitons; this way, thermal and quantum effects — which may be particularly relevant and important in certain cases — will only be briefly covered. Following the above limitations, the structure of the manuscript will be as follows.

Section 2 is devoted to the mean-field description of BECs. Particularly, we first present the GP equation and discuss its connection with the respective full quantum many-body problem. Next, we present the ground state of the condensate and discuss how its small-amplitude excitations can be studied by means of the Bogoliubov-de Gennes (BdG) equations. Lower-dimensional versions of the GP model, pertinent to highly anisotropic trapping potentials, are also discussed; this way, depending on the shape of the trap, we start from purely three-dimensional (3D) BECs and introduce elongated (alias “cigar-shaped”) BECs, quasi one-dimensional (1D) BECs and quasi two-dimensional (2D) (alias “disk-shaped”) ones, as well as discuss cases relevant to the dimensionality crossover regimes. The topics of strongly-interacting Bose gases, and their relevant mean-field description, are also briefly covered.

Section 3 provides the theoretical basis for the study of matter-wave dark solitons. Specifically, first we present the completely integrable 1D NLS equation, its basic properties and the dark soliton solutions. Relevant mathematical tools, such as Inverse Scattering Transform (IST), the renormalization of the integrals of motion of dark solitons and the small-amplitude approximation — leading to the connection of matter-wave dark solitons to Korteweg-de Vries (KdV) solitons — are discussed. Furthermore, the generation of matter-wave dark solitons by means of the phase-, density- and quantum-state-engineering methods are also presented. We also provide the multiple-dark soliton solutions of the NLS equation, and discuss their interactions and collisions.

Section 4 deals with matter-wave dark solitons in quasi-1D Bose gases. Particularly, we first discuss the adiabatic dynamics of dark solitons in the presence of the harmonic trap by means of different analytical techniques; these include the Hamiltonian and Lagrangian approaches of the perturbation theory, the Landau dynamics and the small-amplitude approximation approaches. Next, a connection between the stability, statics and dynamics of dark solitons is presented, relying on a study of the Bogoliubov spectrum of single- and multiple-dark solitons and the role of the pertinent anomalous modes. Non-adiabatic effects, namely emission of radiation of solitons in the form of sound waves, as well as rigorous results concerning the persistence and stability of matter-wave dark solitons, are also discussed.

Section 5 studies matter-wave dark solitons in higher-dimensional settings. Considering, at first, the case of purely 2D or 3D geometries, the transverse (alias “snaking”) instability of rectilinear dark solitons, and the concomitant soliton decay into vortex pairs or vortex rings, is presented. The theme of matter-wave dark solitons of radial symmetry, namely ring dark solitons and spherical shell solitons, is also covered. Furthermore, we present results concerning the stability of dark solitons in cigar-shaped (3D) BECs, and in BECs in the dimensionality crossover regime from 3D to 1D; in the latter experimentally relevant setting, both single- and multiple- dark soliton statics and dynamics are analyzed.

In Section 6, we discuss various experimentally relevant settings and parameter regimes for matter-wave dark solitons. In particular, we first present results concerning matter-wave dark solitons in multi-component (pseudo-spinor and spinor) BECs. Next, we discuss how matter-wave interference and the breakdown of BEC superfluidity are connected to the generation of matter-wave dark solitons. We continue by referring to matter-wave dark solitons in periodic potentials, namely optical lattices (OLs) and superlattices, and conclude this Section by discussing the statics and dynamics of dark solitons at finite temperatures.

Finally, in Section 7 we briefly summarize our conclusions and discuss future challenges.

## 2 Mean-field description of Bose-Einstein condensates

Bose-Einstein condensation of dilute atomic gases is an unambiguous manifestation of a macroscopic quantum state in a many-body system. As such, this phenomenon has triggered an enormous amount of experimental and theoretical work [41, 42]. Importantly, this field is intimately connected with branches of physics such as superfluidity, superconductivity, lasers, coherent optics, nonlinear optics, and physics of nonlinear waves. Many of the common elements between BEC physics and the above areas, and in particular optics, rely on the existence of macroscopic coherence in the many-body state of the system. From a theoretical standpoint, this can be understood by the fact that many effects related to BEC physics can be described by a mean-field model, namely the Gross-Pitaevskii (GP) equation [10]. The latter is a partial differential equation (PDE) of the NLS type, which plays a key role — among other fields — in nonlinear optics [35]. Thus, BEC physics is closely connected to nonlinear optics (and the physics of nonlinear waves), with vortices and solitons being perhaps the most prominent examples of common nonlinear structures arising in these areas [43, 44].

Below we will briefly discuss the theoretical background for the description of BECs. We emphasize, in particular, lowest-order mean-field theory, as this can be used as a basis to understand the nonlinear dynamics of matter-wave dark solitons. Interesting effects naturally arise beyond the GP mean-field, both due to thermal and quantum fluctuations. Such effects become particularly relevant in extremely tightly confining geometries, or when the Bose-Einstein condensation transition region is approached.

### 2.1 The Gross-Pitaevskii equation.

In order to describe theoretically the statics and dynamics of BECs a quantum many-body approach is required [41, 42] (see also Ref. [78] for a recent review on the many-body aspects of BECs). Particularly, a sufficiently dilute ultracold atomic gas, composed by interacting bosons of mass confined by an external potential , can be described by the many-body Hamiltonian; the latter can be expressed, in second quantization form, through the boson annihilation and creation field operators, and (which create and annihilate a particle at the position ) namely,

 ^H = ∫dr^Ψ†(r,t)^H0^Ψ(r,t) (1) + 12∫drdr′^Ψ†(r,t)^Ψ†(r′,t)V(r−r′)^Ψ(r′,t)^Ψ(r,t),

where is the single-particle operator and is the two-body interatomic potential. Apparently, the underlying full many-body problem is very difficult to be treated analytically (or even numerically) as increases and, thus, for convenience, a mean-field approach can be adopted. The mean-field approach is based on the separation of the condensate contribution from the boson field operator as follows [79]:

 ^Ψ(r,t)=⟨^Ψ(r,t)⟩+^Ψ′(r,t)=Ψ(r,t)+^Ψ′(r,t). (2)

In the above expression, the expectation value of the field operator, , is known as the macroscopic wave function of the condensate, while describes the non-condensate part, which accounts for quantum and thermal fluctuations. Considering the case of a dilute ultracold gas with binary collisions at low energy, characterized by the -wave scattering length , the interatomic potential can be replaced by an effective delta-function interaction potential, [41, 42], with the coupling constant given by . Under these assumptions, a nontrivial zeroth-order theory for the BEC wave function can be obtained by means of the Heisenberg evolution equation , upon replacing the field operator with the classical field , i.e., ignoring the quantum and thermal fluctuations described by . Such a consideration leads to the Gross-Pitaevskii (GP) equation [10], which has the form:

 iℏ∂tΨ(r,t)=[−ℏ22m∇2+Vext(r)+g|Ψ(r,t)|2]Ψ(r,t). (3)

In the above equation, is normalized to the number of atoms , namely,

 N=∫|Ψ(r,t)|2dr, (4)

and the nonlinearity (which is obviously introduced by interatomic interactions) is characterized by the -wave scattering length , which is or for repulsive or attractive interatomic interactions, respectively. Notice that Eq. (3) can be written in canonical form, (with denoting complex conjugate), where the dynamically conserved energy functional is given by

 E=∫dr[ℏ22m|∇Ψ|2+Vext|Ψ|2+12g|Ψ|4], (5)

with the three terms in the right-hand side representing, respectively, the kinetic energy, the potential energy and the interaction energy.

A time-independent version of the GP equation can be obtained upon expressing the BEC wave function as , where is the chemical potential. This way, Eq. (3) yields the following equation for the stationary state :

 [−ℏ22m∇2+Vext(r)+g|Ψ0|2(r)]Ψ0(r)=μΨ0(r). (6)

### 2.2 The mean-field approach vs. the many-body quantum mechanical problem.

Although the GP equation is known from the early 60’s [10], it was only recently shown that it can be derived rigorously from a self-consistent treatment of the respective many-body quantum mechanical problem [80]. In particular, in Ref. [80] — which dealt with the stationary GP Eq. (6) — it was proved that the GP energy functional describes correctly the energy and the particle density of a trapped Bose gas to the leading-order in the small parameter 2 , where is the average density of the gas. The above results were proved in the limit where the number of particles and the scattering length , such that the product stays constant. Importantly, although Ref. [80] referred to the full three-dimensional (3D) Bose gas, extensions of this work for lower-dimensional settings were also reported (see the review [81] and references therein).

The starting point of the analysis of Ref. [80] is the effective Hamiltonian of identical bosons, which can be expressed (in units so that ) as follows:

 H=N∑j=1[−∇2j+Vext(rj)]+∑i

where is a general interaction potential assumed to be spherically symmetric and decaying faster than at infinity. Then, assuming that the quantum-mechanical ground-state energy of the Hamiltonian (7) is (here is the number of particles and is the dimensionless two-body scattering length), the main theorem proved in Ref. [80] is the following. The GP energy is the dilute limit of the quantum-mechanical energy:

 ∀~a1>0:limn→∞1NEQM(N,~a1n)=EGP(1,~a1), (8)

where is the energy of a solution of the stationary GP Eq. (6) (in units such that ), and the convergence is uniform on bounded intervals of .

The above results (as well as the ones in Ref. [81]) were proved for stationary solutions of the GP equation, and, in particular, for the ground state solution. More recently, the time-dependent GP Eq. (3) was also analyzed within a similar asymptotic analysis in Ref. [82]. In this work, it was proved that the limit points of the -particle density matrices of (which is the solution of the -particle Schrödinger equation) satisfy asymptotically the GP equation (and the associated hierarchy of equations) with a coupling constant given by , where describes the interaction potential.

These rigorous results, as well as a large number of experimental results related to the physics of BECs, indicate that (under certain conditions) the GP equation is a good starting point for analyzing the statics and dynamics of BECs.

### 2.3 Ground state and excitations of the condensate.

Let us now consider a condensate confined in a harmonic external potential, namely,

 Vext(r)=12m(ω2xx2+ω2yy2+ω2zz2), (9)

where , , and are the (generally different) trap frequencies along the three directions. In this setting, and in the case of repulsive interatomic interactions () and sufficiently large number of atoms , Eq. (6) can be used to determine analytically the ground state of the system. In particular, in the asymptotic limit of (where is the harmonic oscillator length associated with the geometrical average of the trap frequencies), it is expected that the atoms are pushed towards the rims of the condensate, resulting in slow spatial variations of the density profile . Thus, the latter can be obtained as an algebraic solution stemming from Eq. (6) when neglecting the kinetic energy term — the so-called Thomas-Fermi (TF) limit [41, 42, 43]:

 n(r)=g−1[μ−Vext(r)], (10)

in the region where , and outside, and the value of being determined by the normalization condition [cf. Eq. (4)]. Notice that the TF approximation becomes increasingly accurate for large values of .

Small-amplitude excitations of the BEC can be studied upon linearizing Eq. (6) around the ground state. Particularly, we consider small perturbations of this state, i.e.,

 Ψ(r,t)=e−iμt/ℏ[Ψ0(r)+∑j(uj(r)e−iωjt+υ∗j(r)eiωjt)], (11)

where , are the components of the linear response of the BEC to the external perturbations that oscillate at frequencies [the latter are (generally complex) eigenfrequencies]. Substituting Eq. (11) into Eq. (6), and keeping only the linear terms in and , we obtain the so-called Bogoliubov-de Gennes (BdG) equations:

 [^H0−μ+2g|Ψ0|2(r)]uj(r)+gΨ20(r)υj(r)=ℏωjuj(r), [^H0−μ+2g|Ψ0|2(r)]υj(r)+gΨ∗20(r)uj(r)=−ℏωjυj(r), (12)

where is the single-particle Hamiltonian. Importantly, these equations can also be used, apart from the ground state, for any other stationary state (including, e.g., solitons) with the function being modified accordingly. In such a general context, the BdG equations provide the eigenfrequencies and the amplitudes and of the normal modes of the system. Note that due to the Hamiltonian nature of the system, if is an eigenfrequency of the Bogoliubov spectrum, so are , and . In the case of stable configurations with , the solution of BdG equations with frequency represent the same physical oscillation with the solution with frequency [42].

In the case of a homogeneous gas () characterized by a constant density , the amplitudes and in the BdG equations are plane waves, , of wave vector . Then, Eqs. (12) lead to the dispersion relation,

 (ℏω)2=(ℏ2k22m)(ℏ2k22m+2gn0). (13)

In the case of repulsive interatomic interactions (), Eq. (13) indicates that small-amplitude harmonic excitations of the stationary state

 Ψ=√n0exp(−iμt/ℏ), (14)

with , are always stable since for every . Thus, this state is not subject to the modulational instability (see, e.g., Ref. [83] and references therein). This fact is important, as the wave function of Eq. (14) can serve as a stable background (alias “pedestal”), on top of which strongly nonlinear localized excitations may be formed; such excitations may be, e.g., matter-wave dark solitons which are of particular interest in this work. Notice that the above mentioned small-amplitude harmonic excitations are in fact sound waves, characterized by the phonon dispersion relation [see Eq. (13) for small momenta ], where

 cs=√gn0/m, (15)

is the speed of sound. We should note in passing that in the case of attractive interatomic interactions (), the speed of sound becomes imaginary, which indicates that long wavelength perturbations grow or decay exponentially in time. Thus, the stationary state of Eq. (14) is subject to the modulational instability, which is responsible for the formation of matter-wave bright solitons [84, 85, 86] in attractive BECs (see also the reviews [43, 44, 83, 87] and references therein).

### 2.4 Lower-dimensional condensates and relevant mean-field models.

Let us consider again a condensate confined in the harmonic trap of Eq. (9). In this case, the trap frequencies set characteristic length scales for the spatial size of the condensate through the harmonic oscillator lengths (). Another important length scale, introduced by the effective mean-field nonlinearity, is the so-called healing length defined as (with being the maximum condensate density). The healing length, being the scale over which the BEC wave function “heals” over defects, sets the spatial widths of nonlinear excitations, such as matter-wave dark solitons.

Based on the above, as well as the form of the ground state [cf. Eq. (10)], it is clear that the shape of the BEC is controlled by the relative values of the trap frequencies. For example, if (i.e., for an isotropic trap), the BEC is almost spherical, while for (i.e., for an anisotropic trap) the BEC is “cigar shaped”. It is clear that such a cigar-shaped BEC (a) may be a purely 3D object, (b) acquire an almost 1D character (for strongly anisotropic traps with and ), or (c) being in the so-called dimensionality crossover regime from 3D to 1D. These regimes can be described by the dimensionless parameter [88],

 d=NΩaa⊥, (16)

where is the so-called “aspect ratio” of the trap. Particularly, if the dimensionality parameter is , the BEC locally retains its original 3D character (although it may have an elongated, quasi-1D shape) and its ground state can be described by the TF approximation in all directions. On the other hand, if , excited states along the transverse direction are not energetically accessible and the BEC is effectively 1D. Apparently, this regime is extremely useful for an analytical study of matter-wave dark solitons. Finally, if , the BEC is in the crossover regime between 1D and 3D, which is particularly relevant as recent matter-wave dark soliton experiments have been conducted in this regime [69, 71].

Let us now discuss in more detail lower-dimensional mean-field models describing cigar-shaped BECs. First, we consider the quasi-1D regime () characterized by an extremely tight transverse confinement. In this case, following Refs. [50, 89, 90], the BEC wave function is separated into transverse and longitudinal components, namely . Then, the transverse component is described by the Gaussian ground state of the transverse harmonic oscillator (and, thus, the transverse width of the condensate is set by the transverse harmonic oscillator length ), while the longitudinal wave function obeys the following effectively 1D GP equation:

 iℏ∂tψ(z,t)=[−ℏ22m∂2z+V(z)+g1D|ψ(z,t)|2]ψ(z,t), (17)

where the effective 1D coupling constant is given by and . Notice that in the case under consideration, if the additional condition is fulfilled, then the longitudinal condensate density can be described by the TF approximation — see Eq. (10) with now being the 1D chemical potential (and ) [88]. Following the terminology of Ref. [69], this regime will hereafter be referred to as the TF-1D regime.

Next, let us consider the effect of the deviation from one-dimensionality on the longitudinal condensate dynamics. In this case, the wave function can be factorized as before, but with the transverse component assumed to depend also on the longitudinal variable (and time ) [91, 92, 93, 94]. Physically speaking, it is expected that the transverse direction will no longer be occupied by the ground state, but would still be approximated by a Gaussian function with a width that can be treated as a variational parameter [92, 93, 94]. This way, it is possible to employ different variational approaches and derive the following NLS equation for the longitudinal wave function,

 iℏ∂ψ∂t=[−ℏ22m∂2∂z2+V(z)+f(n)]ψ. (18)

The nonlinearity function in Eq. (18) depends on the longitudinal density and may take different forms. Particularly, in Ref. [92] (where variational equations related to the minimization of the action functional were used) is found to be:

 f(n)=g2πa2⊥n√1+2an+ℏω⊥2(1√1+2an+√1+2an), (19)

and the respective NLS equation is known as the non-polynomial Schrödinger equation (NPSE). On the other hand, in Refs. [93, 94] (where variational equations related to the minimization of the transverse chemical potential were used) the result for is:

 f(n)=ℏω⊥√1+4an. (20)

Since, as explained above, the derivation of the mean-field models with the nonlinearity functions in Eqs. (19) and (20) is based on different approaches, these nonlinearity functions are quite different. Nevertheless, they can be “reconciled” in the weakly-interacting limit of : in this case, the width of the transverse wave function becomes and Eq. (18) — with either the nonlinearity function of Eq. (19) or that of Eq. (20) — is reduced to the 1D GP model of Eq. (17).

The above effective 1D models predict accurately ground state properties of quasi-1D condensates, such as the chemical potential, the axial density profile, the speed of sound, collective oscillations, and others. Importantly, these models are particularly useful in the dimensionality crossover regime, where they describe the axial dynamics of cigar-shaped BECs in a very good approximation to the 3D Gross-Pitaevskii (GP) equation (see, e.g., theoretical work related to matter-wave dark solitons in Ref. [95] and relevant experimental results in Ref. [69]).

On the other hand, extremely weak deviations from one-dimensionality can also be treated by means of a rather simple non-cubic nonlinearity that can be obtained by Taylor expanding , namely:

 f(n)=g1n−g2n2, (21)

where and depends on the form of . In this case, Eq. (18) becomes a cubic-quintic NLS (cqNLS) equation. This model was derived self-consistently in Ref. [91], where dynamics of matter-wave dark solitons in elongated BECs was considered; there, the coefficient was found to be equal to .

Here, it is worth mentioning that the quintic term in the cqNLS equation may have a different physical interpretation, namely to describe three-body interactions, regardless of the dimensionality of the system. In this case, the coefficients and in Eq. (21) are generally complex, with the imaginary parts describing inelastic two- and three-body collisions, respectively [96]. As concerns the rate of the three-body collision process, it is given by [41], where is the loss rate (which is of order of cms for various species of alkali atoms [97]). Accordingly, the decrease of the density is accounted for by the term in the time dependent GP equation, i.e., to the quintic term in the cqNLS equation.

It is also relevant to note that the NLS Eq. (18) has also been used as a mean-field model describing strongly-interacting 1D Bose gases and, particularly, the so-called Tonks-Girardeau gas of impenetrable bosons [98] (see also Refs. [99, 100] for recent experimental observations). In this case, the function takes the form [101],

 f(n)=π2ℏ22mn2, (22)

and, thus, Eq. (18) becomes a quintic NLS equation. Although the applicability of this equation has been criticized (as in certain regimes it fails to predict correctly the coherence properties of the strongly-interacting 1D Bose gases [102]), the corresponding hydrodynamic equations for the density and the phase arising from the quintic NLS equation under the Madelung transformation are well-documented in the context of the local density approximation [103]. Additionally, it should be noticed that the time-independent version of the quintic NLS equation has been rigorously derived from the many-body Schrödinger equation [104].

We finally mention that another lower-dimensional version of the fully 3D GP equation can be derived for “disk-shaped” (alias “pancake”) condensates confined in strongly anisotropic traps with and . In such a case, a procedure similar to the one used for the derivation of Eq. (17) leads to the following -dimensional NLS equation:

 iℏ∂tψ(x,y,t)=[−ℏ22m∇2⊥+V(r)+g2D|ψ(x,y,t)|2]ψ(x,y,t), (23)

where , , the effectively 2D coupling constant is given by , while the potential is given by . It should also be noticed that other effective 2D mean-field models (involving systems of coupled 2D equations [92] or 2D GP equations with generalized nonlinearities [92, 94]) have also been proposed for the study of the transverse dynamics of disk-shaped BECs.

The above models will be used below to investigate the static and dynamical properties of matter-wave dark solitons arising in the respective settings.

## 3 General background for the study of matter-wave dark solitons

### 3.1 NLS equation and dark soliton solutions.

We start by considering the case of a quasi-1D condensate described by Eq. (17). The latter, can be expressed in the following dimensionless form:

 i∂tψ(z,t)=[−12∂2z+V(z)+|ψ(z,t)|2]ψ(z,t), (24)

where the density , length, time and energy are respectively measured in units of , , and , while the potential is given by

 V(z)=12Ω2z2. (25)

In the case under consideration, the normalized trap strength (aspect ratio) is and, thus, as a first step in our analysis, the potential is ignored. 3 In such a case, the condensate is homogeneous and can be described by the completely integrable defocusing NLS equation [12] (see also the review [34]):

 i∂tψ(z,t)=[−12∂2z+|ψ(z,t)|2]ψ(z,t). (26)

This equation possesses an infinite number of conserved quantities (integrals of motion); the lowest-order ones are the number of particles , the momentum , and the energy , respectively given by:

 N = ∫−∞−∞|ψ|2dz, (27) P = i2∫−∞−∞(ψ∂zψ∗−ψ∗∂zψ)dz, (28) E = 12∫−∞−∞(|∂zψ|2+|ψ|4)dz. (29)

It is also noted that the NLS Eq. (26) can be obtained by the Euler-Lagrange equation =, where the Lagrangian density is given by:

 L=i2(ψ∂tψ∗−ψ∗∂tψ)−12(|∂zψ|2+|ψ|4). (30)

The simplest nontrivial solution of Eq. (26) is a plane wave of wave number and frequency , namely,

 ψ=√n0exp[i(kz−ωt+θo)],ω=12k2−μ, (31)

where the constant BEC density sets the chemical potential, i.e., and is an arbitrary constant phase. This solution, which is reduced to the stationary state of Eq. (14) for , is also modulationally stable as can be confirmed by a simple stability analysis (see, e.g., Refs. [34, 83]). For small densities, , the above plane wave satisfies the linear Schrödinger equation, , and the pertinent linear wave solutions of the NLS equation are characterized by the dispersion relation . Notice that if the system is characterized by a length , then the integrals of motion for the stationary solution in Eq. (31) take the values:

 N=2n0L,P=kn0L,E=12(k2−n0)n0L. (32)

The NLS equation admits nontrivial solutions, in the form of dark solitons, which can be regarded as strongly nonlinear excitations of the plane wave solution (31). In the most general case of a moving background [ in Eq. (31)], a single dark soliton solution may be expressed as [12],

 ψ(z,t)=√n0(Btanhζ+iA)exp[i(kz−ωt+θo)], (33)

where ; here, is the soliton center, is an arbitrary real constant representing the initial location of the dark soliton, is the relative velocity between the soliton and the background given by , the frequency is provided by the dispersion relation of the background plane wave, [cf. Eq. (31)] 4 , and, finally, the parameters and are connected through the equation . In some cases it is convenient to use one parameter instead of two and, thus, one may introduce

 A=sinϕ,B=cosϕ, (34)

where is the so-called “soliton phase angle” (). Notice that although the asymptotics of the dark soliton solution (33) coincide with the ones of Eq. (31), the plane waves at have different phases; as a result, there exists a nontrivial phase jump across the dark soliton, given by:

 Δϕ=2[tan−1(BA)−π2]=−2tan−1(AB). (35)

Note that, hereafter, we will consider the simpler case where the background of matter-wave dark solitons is at rest, i.e., ; then, the frequency actually plays the role of the normalized chemical potential, namely , which is determined by the number of atoms of the condensate.

The soliton phase angle describes also the darkness of the soliton, namely,

 |ψ|2=n0(1−cos2ϕsech2ζ). (36)

This way, the cases and correspond to the so-called black and gray solitons, respectively. The amplitude and velocity of the dark soliton are given (for ) by and , respectively; thus, the black soliton

 ψ=√n0tanh(√n0z)exp(−iμt), (37)

is characterized by a zero velocity, (and, thus, it is also called stationary kink), while the gray soliton moves with a finite velocity . Examples of the forms of a black and a gray soliton are illustrated in Fig. 1.

In the limiting case of a very shallow (small-amplitude) dark soliton with , the soliton velocity is close to the speed of sound which, in our units, is given by:

 cs=√n0. (38)

The speed of sound is, therefore, the maximum possible velocity of a dark soliton which, generally, always travels with a velocity less than the speed of sound. We finally note that the dark soliton solution (33) has two independent parameters (for ), one for the background, , and one for the soliton, , while there is also a freedom (translational invariance) in selecting the initial location of the dark soliton 5 .

In the case of a condensate confined in a harmonic trap [cf. Eq. (9)], the background of the dark soliton is, in fact, of finite extent, being the ground state of the BEC [which may be approximated by the Thomas-Fermi cloud, cf. Eq. (10)]. For example, in the quasi-1D setting of the 1D GP Eq. (24) with the harmonic potential in Eq. (25), the “composite” wave function (describing both the background and the soliton) can be approximated as , where is the TF background and is the dark soliton wave function of Eq. (33), which satisfies the 1D GP equation for .

### 3.2 Dark solitons and the Inverse Scattering Transform.

The single dark soliton solution of the NLS Eq. (26) presented in the previous Section, as well as multiple dark soliton solutions (see Sec. 3.6 below), can be derived by means of the Inverse Scattering Transform (IST) [12]. A basic step of this approach is the solution of the Zakharov-Shabat (ZS) eigenvalue problem, with eigenvalue , for the auxiliary two-component eigenfunction , namely,

 LU=(i∂zψ(z,0)ψ∗(z,0)−i∂z)(u1u2)=λU, (39)

with the boundary conditions , for , and , for . Here, is the amplitude of the background wave function and is a constant phase. Since the operator is self-adjoint, the ZS eigenvalue problem possesses real discrete eigenvalues , with magnitudes . Importantly, each real discrete eigenvalue corresponds to a dark soliton of depth and velocity . To make a connection to the dark soliton solutions of the NLS equation presented in the previous Section, we note that the dark soliton of Eq. (33) corresponds to a single eigenvalue .

Although the system of ZS Eqs. (39) is linear, its general solution for arbitrary initial condition is not available. Thus, various methods have been developed for the determination of the spectrum of the ZS problem, such as the so-called quasi-classical method [14, 15] (see also Ref. [19]), the variational approach [105], as well as other techniques that can be applied to the case of dark soliton trains [106, 107]. In any case, the generation of single- as well as multiple-dark solitons (see Sec. 3.6 below) can be studied in the framework of the IST method, and many useful results can be obtained. In that regard, first we note that a pair of dark solitons — corresponding to a discrete eigenvalue pair in the associated scattering problem — can always be generated by an arbitrary small dip on a background of constant density [14] (see also Ref. [15]). This means that the generation of dark solitons is a thresholdless process, contrary to the case of bright solitons which are created when the number of atoms exceeds a certain threshold [108]. In another example, as dark solitons are characterized by a phase jump across them, we may assume that they can be generated by an anti-symmetric initial wave function profile of the form,

 ψ(z,0)=√n0tanh(αz), (40)

characterized by a background density and a width (the ratio is assumed to be arbitrary). In such a case, the ZS eigenvalue problem (39) can be solved exactly [16, 17, 18] and the resulting eigenvalues of the discrete spectrum are given by and , where positive are defined as , , and is the largest integer such that . These results show that for arbitrary the initial wave function profile of Eq. (40) will always produce a black soliton [cf. Eq. (37)] at (corresponding to the first, zero eigenvalue) and additional pairs of symmetric gray solitons (corresponding to the even number of the secondary, nonzero eigenvalues), propagating to the left and to the right of the primary black soliton. Apparently, the total number of eigenvalues and, thus, the total number of solitons, is and depends on the ratio . Apart from the above example, dark soliton generation was systematically studied in Ref. [15] for a variety of initial conditions (such as box-like dark pulses, phase steps, and others). Notice that, generally, initial wave function profiles with odd symmetry will produce an odd number of dark solitons, while profiles with an even symmetry (as, e.g., in the study of Ref. [14]) produce pairs of dark solitons; this theoretical prediction was also confirmed in experiments with optical dark solitons [109]. Furthermore, the initial phase change across the wave function plays a key role in dark soliton formation, while the number of dark solitons that are formed can be changed by small variations of the phase.

### 3.3 Integrals of motion and basic properties of dark solitons.

Let us now proceed by considering the integrals of motion for dark solitons. Taking into regard that Eqs. (27)–(29) refer to both the background and the soliton, one may follow Refs. [27, 28, 110, 111], and renormalize the integrals of motion so as to extract the contribution of the background [see Eqs. (32)]. This way, the renormalized integrals of motion become finite and, when calculated for the dark soliton solution (33), provide the following results (for ). The number of particles of the dark soliton reads:

 Nds=∫−∞−∞(n0−|ψ|2)dz=2√n0B. (41)

The momentum of the dark soliton is given by,

 Pds = i2∫−∞−∞(ψ∂zψ∗−ψ∗∂zψ)dz−n0Δϕ (42) = Missing or unrecognized delimiter for \right = −2v(c2s−v2)1/2+2c2stan−1[(c2s−v2)1/2v],

where is given by Eq. (35) and is the speed of sound. Furthermore, the energy of the dark soliton is given by,

 Eds=12∫−∞−∞[|∂zψ|2+(|ψ|2−n0)2]dz=43(c2s−v2)3/2, (43)

while the renormalized Lagrangian density takes the form [25]:

 Lds=i2(ψ∂tψ∗−ψ∗∂tψ)(1−n0|ψ|2)−12[|∂zψ|2+(|ψ|2−n0)2]. (44)

The renormalized integrals of motion can now be used for a better understanding of basic features of dark solitons. To be more specific, one may differentiate the expressions (42) and (43) over the soliton velocity to obtain the result,

 ∂Eds∂Pds=v, (45)

which shows that the dark soliton effectively behaves like a classical particle, obeying a standard equation of classical mechanics. Furthermore, it is also possible to associate an effective mass to the dark soliton, according to the equation . This way, using Eq. (42), it can readily be found that

 mds=−4√n0B, (46)

which shows the dark soliton is characterized by a negative effective mass. The same result, but for almost black solitons () with sufficiently small soliton velocities (), can also be obtained using Eq. (43) [112]: in this case, the energy of the dark soliton can be approximated as or, equivalently,

 Eds=E0+12mdsv2, (47)

where , and the soliton’s effective mass is .

### 3.4 Small-amplitude approximation: shallow dark solitons as KdV solitons.

As mentioned above, the case of corresponds to a small-amplitude (shallow) dark soliton, which travels with a speed close to the speed of sound, i.e., . In this case, it is possible to apply the reductive perturbation method [113] and show that, in the small-amplitude limit, the NLS dark soliton can be described by an effective KdV equation (see, e.g., Ref. [114] for various applications of the KdV model). The basic idea of this, so-called, small-amplitude approximation can be understood in terms of the similarity between the KdV soliton and the shallow dark soliton’s density profile: indeed, the KdV equation for a field expressed as,

 ∂tu+6u∂zu+∂3zu=0, (48)

possesses a single soliton solution (see, e.g., Ref. [13]):

 u(z,t)=2κ2sech2[κ(z−4κ2t)] (49)

(with being an arbitrary constant), which shares the same functional form with the density profile of the shallow dark soliton of the NLS equation [see Eqs. (49) and (36)]. The reduction of the cubic NLS equation to the KdV equation was first presented in Ref. [9] and later the formal connection between several integrable evolution equations was investigated in detail [115]. Importantly, such a connection is still possible even in cases of strongly perturbed NLS models, a fact that triggered various studies on dark soliton dynamics in the presence of perturbations (see, e.g., Refs. [116, 117, 118] for studies in the context of optics, as well as the recent review [119] and references therein). Generally, the advantage of the small-amplitude approximation is that it may predict approximate analytical dark soliton solutions in models where exact analytical dark soliton solutions are not available, or can only be found in an implicit form [116].

Let us now consider a rather general case, and discuss small-amplitude dark solitons of the generalized NLS Eq. (18); in the absence of the potential (), this equation is expressed in dimensionless form as:

 i∂tψ=−12∂2zψ+f(n)ψ, (50)

where the units are the same to the ones used for Eq. (24). Then, we use the Madelung transformation (with and representing the BEC density and phase, respectively) to express Eq. (50) in the hydrodynamic form:

 ∂tφ+f(n)+12(∂zφ)2−12n−1/2∂2zn1/2=0, (51) ∂tn+∂z(n∂zφ)=0. (52)

The simplest solution of Eqs. (51)–(52) is and , where . Note that in the model of Eq. (19) one has , for the model of Eq. (20), , and so on. Next, assuming slow spatial and temporal variations, we define the slow variables

 Z=ϵ1/2(z−ct),T=ϵ3/2t, (53)

where is a formal small parameter () connected with the soliton amplitude. Additionally, we introduce asymptotic expansions for the density and phase:

 n = n0+ϵn1(Z,T)+ϵ2n2(Z,T)+⋯, (54) φ = −f0t+ϵ1/2φ1(Z,T)+ϵ3/2φ2(Z,T)+⋯. (55)

Then, substituting Eqs. (54)–(55) into Eqs. (51)–(52), and Taylor expanding the nonlinearity function as (where ), we obtain a hierarchy of equations. In particular, Eqs. (51)–(52) lead, respectively, at the order and , to the following linear system,

 −c∂Zφ1+f′0n1=0,n0∂2Zφ1−c∂Zn1=0. (56)

The compatibility condition of the above equations is the algebraic equation , which shows that the velocity in Eq. (53) is equal to the speed of sound, . Additionally, Eqs. (56) connect the phase and the density through the equation:

 ∂Zφ1=csn0n1. (57)

To the next order, viz. and , Eqs. (51) and (52), respectively, yield:

 ∂Tφ1−cs∂Zφ2+f′0n2+12f′′0n21+12(∂Zφ1)2−14n−10∂2Zn1=0, (58) ∂Tn1−cs∂Zn2+∂Z(n1∂Zφ1)+n0∂2Zφ2=0. (59)

The compatibility conditions of Eqs. (58)–(59) are the algebraic equation , along with a KdV equation [see Eq. (48)] for the unknown density :

 2cs∂Tn1+(3f′0+n0f′′0)n1∂Zn1−14∂3Zn1=0. (60)

Thus, the density of the shallow dark soliton can be expressed as a KdV soliton [see Eq. (49)]. In terms of the original time and space variables, is expressed as follows:

 n1(z,t)=−3κ22(3f′0+n0f′′0)sech2[ϵ1/2κ(z−vt)], (61)

where is (as before) an arbitrary parameter [assumed to be of order ], while is the soliton velocity; the latter, is given by

 v=cs−ϵκ22cs, (62)

and, clearly, . Apparently, Eq. (61) describes a small-amplitude dip [of order  — see Eq. (54)] on the background density of the condensate, with a phase that can be found using Eq. (57); in terms of the variables and , the result is:

 φ1(z,t)=−3κcs2n0(3f′0+n0f′′0)tanh[ϵ1/2κ(z−vt)]. (63)

The above expression shows that the density dip is accompanied by a