{'fraction_non_alphanumeric': 0.0536435639070443, 'fraction_numerical': 0.02228576615831518, 'mean_word_length': 4.491727683439788, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dark solitons in atomic Bose-Einstein condensates: from theory to experiments

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.


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 Kortewegde 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 matterwave 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 multipledark 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 smallamplitude 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.

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 meanfield 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.

The Gross-Pitaevskii equation.
In order to describe theoretically the statics and dynamics of BECs a quantum manybody 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 N interacting bosons of mass m confined by an external potential V ext (r), 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,Ψ(r, t) andΨ † (r, t) (which create and annihilate a particle at the position r) namely, H = drΨ † (r, t)Ĥ 0Ψ (r, t) whereĤ 0 = −( 2 /2m)∇ 2 + V ext (r) is the single-particle operator and V (r − r ′ ) is the two-body interatomic potential. Apparently, the underlying full many-body problem is very difficult to be treated analytically (or even numerically) as N 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]:Ψ In the above expression, the expectation value of the field operator, Ψ (r, t) ≡ Ψ(r, t), is known as the macroscopic wave function of the condensate, whileΨ ′ (r, t) 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 s-wave scattering length a, the interatomic potential can be replaced by an effective delta-function interaction potential, V (r ′ − r) = gδ(r ′ − r) [41,42], with the coupling constant g given by g = 4π 2 a/m. Under these assumptions, a nontrivial zeroth-order theory for the BEC wave function can be obtained by means of the Heisenberg evolution equation i (∂Ψ/∂t) = [Ψ,Ĥ], upon replacing the field operatorΨ with the classical field Ψ, i.e., ignoring the quantum and thermal fluctuations described byΨ ′ (r ′ , t). Such a consideration leads to the Gross-Pitaevskii (GP) equation [10], which has the form: In the above equation, Ψ(r, t) is normalized to the number of atoms N, namely, and the nonlinearity (which is obviously introduced by interatomic interactions) is characterized by the s-wave scattering length a, which is a > 0 or a < 0 for repulsive or attractive interatomic interactions, respectively. Notice that Eq. (3) can be written in canonical form, i ∂ t Ψ = δE/δΨ * (with * denoting complex conjugate), where the dynamically conserved energy functional E is given by 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 Ψ(r, t) = Ψ 0 (r) exp(−iµt/ ), where µ = ∂E/∂N is the chemical potential. This way, Eq. (3) yields the following equation for the stationary state Ψ 0 : 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 manybody 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 parametern|a| 3 † , wheren is the average density of the gas. The above results were proved in the limit where the number of particles N → ∞ and the scattering length a → 0, such that the product Na stays constant. Importantly, although Ref. [80] referred to the full three-dimensional (3D) Bose gas, extensions of this work for lowerdimensional 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 N identical bosons, which can be expressed (in units so that = 2m = 1) as follows: where v(|r|) is a general interaction potential assumed to be spherically symmetric and decaying faster than |r| −3 at infinity. Then, assuming that the quantum-mechanical ground-state energy of the Hamiltonian (7) is E QM (N,ã) (here N 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 quantummechanical energy: where E GP (N,ã) is the energy of a solution of the stationary GP Eq. (6) (in units such that = 2m = 1), and the convergence is uniform on bounded intervals ofã 1 . 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 k-particle density matrices of Ψ N,t (which is the solution of the N-particle Schrödinger equation) satisfy asymptotically the GP equation (and the associated hierarchy of equations) with a coupling constant given by v(x)dx, where v(x) 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.

Ground state and excitations of the condensate.
Let us now consider a condensate confined in a harmonic external potential, namely, where ω x , ω y , and ω z are the (generally different) trap frequencies along the three directions. In this setting, and in the case of repulsive interatomic interactions (a > 0) and sufficiently large number of atoms N, Eq. (6) can be used to determine analytically the ground state of the system. In particular, in the asymptotic limit of Na/a ho ≫ 1 (where a ho = /(mω ho ) is the harmonic oscillator length associated with the geometrical average ω ho = (ω x ω y ω z ) 1/3 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 n(r) ≡ |Ψ 0 (r)| 2 . Thus, the latter can be obtained as an algebraic solution stemming from Eq. (6) when neglecting the kinetic energy termthe so-called Thomas-Fermi (TF) limit [41][42][43]: in the region where µ > V ext (r), and n = 0 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., where u j , υ j are the components of the linear response of the BEC to the external perturbations that oscillate at frequencies ±ω j [the latter are (generally complex) eigenfrequencies]. Substituting Eq. (11) into Eq. (6), and keeping only the linear terms in u j and υ j , we obtain the so-called Bogoliubov-de Gennes (BdG) equations: whereĤ 0 ≡ −( 2 /2m)∇ 2 + V ext (r) 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 Ψ 0 being modified accordingly. In such a general context, the BdG equations provide the eigenfrequencies ω ≡ ω r + iω i and the amplitudes u j and υ j 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 ω i = 0, 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 (V ext (r) = 0) characterized by a constant density n 0 = |Ψ 0 | 2 , the amplitudes u j and υ j in the BdG equations are plane waves, ∼ exp(ik·r), of wave vector k. Then, Eqs. (12) lead to the dispersion relation, In the case of repulsive interatomic interactions (g > 0), Eq. (13) indicates that smallamplitude harmonic excitations of the stationary state with µ = n 0 , are always stable since ω i = 0 for every k. 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 ω = |k|c s [see Eq. (13) for small momenta k], where is the speed of sound. We should note in passing that in the case of attractive interatomic interactions (g < 0), 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).

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 a j ≡ ( /mω j ) 1/2 (j ∈ {x, y, z}). Another important length scale, introduced by the effective mean-field nonlinearity, is the so-called healing length defined as ξ = (8πn 0 a) −1/2 (with n 0 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 ω x = ω y ≡ ω ⊥ ≈ ω z (i.e., for an isotropic trap), the BEC is almost spherical, while for ω z < ω ⊥ (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 ω z ≪ ω ⊥ 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], where Ω = ω z /ω ⊥ is the so-called "aspect ratio" of the trap. Particularly, if the dimensionality parameter is d ≫ 1, 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 d ≪ 1, 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 d ≈ 1, 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 (d ≪ 1) 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 Ψ(r, t) = Φ(r; t)ψ(z, t). Then, the transverse component Φ(r; t) 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 a ⊥ ), while the longitudinal wave function ψ(z, t) obeys the following effectively 1D GP equation: where the effective 1D coupling constant is given by g 1D = g/2πa 2 ⊥ = 2a ω ⊥ and V (z) = (1/2)mω 2 z z 2 . Notice that in the case under consideration, if the additional condition [(N/ √ Ω)(a/a ⊥ )] 1/3 ≫ 1 is fulfilled, then the longitudinal condensate density n(z, t) ≡ |ψ(z, t)| 2 can be described by the TF approximation -see Eq. (10) with µ now being the 1D chemical potential (and g → g 1D ) [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 z (and time t) [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 w = w(z, t) 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, The nonlinearity function f (n) in Eq. (18) depends on the longitudinal density n(z, t) and may take different forms. Particularly, in Ref. [92] (where variational equations related to the minimization of the action functional were used) f (n) is found to be: 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 f (n) is: 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 weaklyinteracting limit of an ≪ 1: in this case, the width of the transverse wave function becomes w = a ⊥ 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 f (n), namely: where g 1 = g 1D and g 2 depends on the form of f (n). 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 g 2 was found to be equal to g 2 = 24 ln(4/3)a 2 ω ⊥ . 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 g 1D and g 2 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 (dn/dt) = −Ln 3 [41], where L is the loss rate (which is of order of 10 −27 -10 −30 cm 6 s −1 for various species of alkali atoms [97]). Accordingly, the decrease of the density is accounted for by the term −(L/2)|ψ| 4 ψ 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 f (n) takes the form [101], 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 n and the phase ϕ arising from the quintic NLS equation under the Madelung transformation ψ = √ n exp(iϕ) 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 ω ⊥ ≪ ω z and µ ≪ ω z . In such a case, a procedure similar to the one used for the derivation of Eq. (17) leads to the following (2 + 1)dimensional NLS equation: where r 2 = x 2 + y 2 , ∇ 2 ⊥ = ∂ 2 x + ∂ 2 y , the effectively 2D coupling constant is given by g 2D = g/ √ 2πa z = 2 √ 2πaa z ω z , while the potential is given by V (r) = (1/2)mω 2 ⊥ r 2 . 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.

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: where the density |ψ| 2 , length, time and energy are respectively measured in units of 2a, a ⊥ , ω −1 ⊥ and ω ⊥ , while the potential V (z) is given by In the case under consideration, the normalized trap strength (aspect ratio) is Ω ≪ 1 and, thus, as a first step in our analysis, the potential V (z) is ignored. † In such a † Note that in the limit of z → ±∞ this approximation always breaks down.
case, the condensate is homogeneous and can be described by the completely integrable defocusing NLS equation [12] (see also the review [34]): This equation possesses an infinite number of conserved quantities (integrals of motion); the lowest-order ones are the number of particles N, the momentum P , and the energy E, respectively given by: It is also noted that the NLS Eq. (26) can be obtained by the Euler-Lagrange equation where the Lagrangian density L is given by: The simplest nontrivial solution of Eq. (26) is a plane wave of wave number k and frequency ω, namely, where the constant BEC density n 0 sets the chemical potential, i.e., n 0 = µ and θ o is an arbitrary constant phase. This solution, which is reduced to the stationary state of Eq. (14) for k = 0, is also modulationally stable as can be confirmed by a simple stability analysis (see, e.g., Refs. [34,83]). For small densities, n 0 ≪ 1, the above plane wave satisfies the linear Schrödinger equation, i∂ t ψ + 1 2 ∂ 2 z ψ = 0, and the pertinent linear wave solutions of the NLS equation are characterized by the dispersion relation ω = 1 2 k 2 . Notice that if the system is characterized by a length L, then the integrals of motion for the stationary solution in Eq. (31) take the values: 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 [k = 0 in Eq. (31)], a single dark soliton solution may be expressed as [12], where ζ ≡ √ n 0 B [z − z 0 (t)]; here, z 0 (t) = vt + z o is the soliton center, z o is an arbitrary real constant representing the initial location of the dark soliton, v is the relative velocity between the soliton and the background given by v = A √ n 0 + k, the frequency ω is provided by the dispersion relation of the background plane wave, ω = (1/2)k 2 + n 0  [cf. Eq. (31)] † , and, finally, the parameters A and B are connected through the equation A 2 + B 2 = 1. In some cases it is convenient to use one parameter instead of two and, thus, one may introduce where φ is the so-called "soliton phase angle" (|φ| < π/2). Notice that although the asymptotics of the dark soliton solution (33) coincide with the ones of Eq. (31), the plane waves at z → ±∞ have different phases; as a result, there exists a nontrivial phase jump ∆φ across the dark soliton, given by: Note that, hereafter, we will consider the simpler case where the background of matterwave dark solitons is at rest, i.e., k = 0; then, the frequency ω actually plays the role of the normalized chemical potential, namely ω = µ = n 0 , which is determined by the number of atoms of the condensate. The soliton phase angle describes also the darkness of the soliton, namely, † Here, this dispersion relation implies that ω > k 2 and, thus, the allowable region in the (k, ω) plane for dark solitons is located above the parabola ω = 1 2 k 2 corresponding to the linear waves.
This way, the cases φ = 0 and 0 < φ < π/2 correspond to the so-called black and gray solitons, respectively. The amplitude and velocity of the dark soliton are given (for k = 0) by √ n 0 cos φ and √ n 0 sin φ, respectively; thus, the black soliton is characterized by a zero velocity, v = 0 (and, thus, it is also called stationary kink), while the gray soliton moves with a finite velocity v = 0. 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 cos φ ≪ 1, the soliton velocity is close to the speed of sound which, in our units, is given by: 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 k = 0), one for the background, n 0 , and one for the soliton, φ, while there is also a freedom (translational invariance) in selecting the initial location of the dark soliton z o † .
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 ψ(z, t) = Φ(z) exp(−iµt)ψ ds (z, t), where Φ(z) is the TF background and ψ ds (z, t) is the dark soliton wave function of Eq. (33), which satisfies the 1D GP equation for V (z) = 0.

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 U = (u 1 , u 2 ) T , namely, with the boundary conditions ψ(z, 0) → √ n 0 , for z → +∞, and ψ(z, 0) → √ n 0 exp(iθ), for z → −∞. Here, √ n 0 is the amplitude of the background wave function and θ is a constant phase. Since the operator L is self-adjoint, the ZS eigenvalue problem possesses real discrete eigenvalues λ j , with magnitudes |λ j | < √ n 0 . Importantly, each real discrete eigenvalue λ j = √ n 0 sin φ j corresponds to a dark soliton of depth √ n 0 cos φ j † Recall that the underlying model, namely the completely integrable NLS equation, has infinitely many symmetries, including translational and Galilean invariance. and velocity √ n 0 sin φ j . 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 λ = √ n 0 sin φ.
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, characterized by a background density n 0 and a width α −1 (the ratio √ n 0 /α 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 λ 1 = 0 and λ 2j = −λ 2j+1 = n 0 − µ 2 j , where positive µ j are defined as µ j = √ n 0 − jα, j = 1, 2, · · · , N 0 , and N 0 is the largest integer such that N 0 < √ n 0 /α. These results show that for arbitrary √ n 0 /α the initial wave function profile of Eq. (40) will always produce a black soliton [cf. Eq. (37)] at z = 0 (corresponding to the first, zero eigenvalue) and additional N 0 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 2N 0 + 1 and depends on the ratio √ n 0 /α. 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.

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 k = 0). The number of particles N ds of the dark soliton reads: The momentum P ds of the dark soliton is given by, where ∆φ is given by Eq. (35) and c s = √ n 0 is the speed of sound. Furthermore, the energy E ds of the dark soliton is given by, while the renormalized Lagrangian density takes the form [25]: 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 v ≡ A √ n 0 to obtain the result, 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 m ds = ∂P ds /∂v. This way, using Eq. (42), it can readily be found that which shows the dark soliton is characterized by a negative effective mass. The same result, but for almost black solitons (B ≈ 1) with sufficiently small soliton velocities (v 2 ≪ c 2 s ), can also be obtained using Eq. (43) [112]: in this case, the energy of the dark soliton can be approximated as E ds ≈ (4/3)c 3 s − 2c s v 2 or, equivalently, where E 0 ≡ 4 3 c 3 s , and the soliton's effective mass is m ds = −4 √ n 0 .
As mentioned above, the case of B 2 = cos 2 φ ≪ 1 corresponds to a small-amplitude (shallow) dark soliton, which travels with a speed v close to the speed of sound, i.e., v ≈ c s . 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 u(z, t) expressed as, possesses a single soliton solution (see, e.g., Ref. [13]): (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 (V (z) = 0), this equation is expressed in dimensionless form as: where the units are the same to the ones used for Eq. (24). Then, we use the Madelung transformation ψ(z, t) = n(z, t) exp[iϕ(z, t)] (with n ≡ |ψ| 2 and ϕ representing the BEC density and phase, respectively) to express Eq. (50) in the hydrodynamic form: The simplest solution of Eqs.
Note that in the model of Eq. (19) one has f 0 = 2+3n 0 2 √ 1+n 0 , for the model of Eq. (20), f 0 = √ 1 + 2n 0 , and so on. Next, assuming slow spatial and temporal variations, we define the slow variables where ǫ is a formal small parameter (0 < ǫ ≪ 1) connected with the soliton amplitude. Additionally, we introduce asymptotic expansions for the density and phase: Then, substituting Eqs. (54)-(55) into Eqs. (51)- (52), and Taylor expanding the , we obtain a hierarchy of equations. In particular, Eqs. (51)-(52) lead, respectively, at the order O(ǫ) and O(ǫ 3/2 ), to the following linear system, The compatibility condition of the above equations is the algebraic equation c 2 = f ′ 0 n 0 , which shows that the velocity c in Eq. (53) is equal to the speed of sound, c ≡ c s . Additionally, Eqs. (56) connect the phase ϕ 1 and the density n 1 through the equation: To the next order, viz. O(ǫ 2 ) and O(ǫ 5/2 ), Eqs. (51) and (52), respectively, yield: The compatibility conditions of Eqs.
Thus, the density n 1 of the shallow dark soliton can be expressed as a KdV soliton [see Eq. (49)]. In terms of the original time and space variables, n 1 is expressed as follows: where κ is (as before) an arbitrary parameter [assumed to be of order O(1)], while v is the soliton velocity; the latter, is given by and, clearly, v c s . Apparently, Eq. (61) describes a small-amplitude dip [of order O(ǫ) -see Eq. (54)] on the background density of the condensate, with a phase ϕ 1 that can be found using Eq. (57); in terms of the variables z and t, the result is: The above expression shows that the density dip is accompanied by a tanh-shaped phase jump. Thus, the wave function characterized by the density n 1 in Eq. (61) and the phase ϕ 1 in Eq. (63) is an approximate shallow dark soliton solution of the GP Eq. (50), obeying the effective KdV Eq. (60).
Notice that the above analysis applies for f (n) = n (i.e., for the cubic NLS model), as well as for all forms of the nonlinearity function in Eqs. (19)- (22). Furthermore, variants of the reductive perturbation method have also been applied for the study of matter-wave dark solitons in higher-dimensional settings [120,121], multi-component condensates [122,123] (see also Sec. 6.1) and combinations thereof [124].

On the generation of matter-wave dark solitons
Matter-wave dark solitons can be created in experiments by means of various methods, namely the phase-imprinting, density-engineering, quantum-state engineering (which is a combination of phase-imprinting and density engineering), the matter-wave interference method and by dragging an obstacle sufficiently fast through a condensate. In connection to Sec. 3.2 -and following the historical evolution of the subject -here we will discuss the phase-imprinting, density-engineering and quantum-state engineering methods (the remaining two methods will be presented in Secs. 6.2 and 6.3 below).
3.5.1. The phase-imprinting method. The earlier results of Sec. 3.2, as well as more recent theoretical studies in the BEC context [125,126] (see also Ref. [127]), paved the way for the generation of matter-wave dark solitons by means of the phase-imprinting method. This technique was used in the earlier [45,46,49] -but also in recent [67,68] -matter-wave dark soliton experiments. The phase-imprinting method involves a manipulation of the BEC phase, without changing the BEC density, which can be implemented experimentally by illuminating part of the condensate by a short offresonance laser beam (i.e., a laser beam with a frequency far from the relevant atomic resonant frequency -see details in the review [128]). This procedure can be described in the framework of Eq. (17), by considering a time-dependent potential of the form is the laser pulse envelope and φ(z) is the imprinted phase, given by [129], where ∆φ is the phase gradient, while the width W of the potential edge sets the steepness of the phase gradient at z * . Note that since experimentally relevant values correspond to a 10-20% absorption width of the phase step, an empirical factor b = 0.45 is also introduced in Eq. (64) [129]. From a theoretical standpoint, phase-imprinting can be studied (in the absence of the trapping potential) in the framework of the IST method, upon considering an initial wave function of the form ψ(z, 0) = exp[iφ(z)]; here the imprinted phase φ(z) is assumed to increase from left to right and approach constants as z → ±∞ [126] [as, e.g., in Eq. (64)]. The pertinent ZS eigenvalue problem can be solved by mapping Eqs. (39) to a damped driven pendulum problem. This way, a formula for the number of both the even and the odd number of generated dark solitons, traveling in both directions, can be derived analytically.
In some experiments (see, e.g., Ref. [45]), the generation of the "dominant" dark soliton is followed by the generation of a secondary wave packet traveling in the opposite direction with a velocity near the speed of sound. This effect can also be understood in the framework of IST: small perturbations of the dark soliton produce shallow "satellite" dark solitons moving with velocities v c s [18].
3.5.2. The density-engineering method. The density-engineering method involves a direct manipulation of the BEC density, without changing the BEC phase, such that local reductions of the density are created which eventually evolve into dark solitons. This technique was used in the Harvard experiments [47,65], where a compressed pulse of slow light was used to create a defect on the condensate density. This defect induced the formation of shock waves that shed dark solitons (or other higher-dimensional topological structures, such as vortex rings [65]). Notice that the use of a compressed pulse of slow light is not really necessary or beneficial in order to create dark solitons by means of the density engineering method: in fact, a local reduction of the BEC density can also be created by modifying the (harmonic) trapping potential with an additional barrier potential, which may be induced by an optical dipole potential or a far-detuned laser beam; this barrier can then be switched off non-adiabatically (while the harmonic trap is kept on), creating the desired local reduction of the density [129]. This technique was employed in a recent experiment [72], where such a dipole beam was used in different setups to induce merging and splitting rubidium condensates; depending on the parameters, this process leads to the formation of dark soliton trains, or a high density bulge and dispersive shock waves.
As in the case of phase-imprinting, the density-engineering technique can be studied by means of the IST method (in the absence of the trapping potential). In fact, earlier works [14,15] (see also Ref. [107]) have already addressed the problem of dark soliton generation induced by initial change of the density: for example, in the case of a box-like initial condition, namely ψ(z, 0) = √ n 0 for |z| > z 0 and ψ(z, 0) = √ n 1 for |z| < z 0 (with n 1 < n 0 ), the ZS spectral problem admits an explicit solution, as it can be solved exactly on the intervals |z| < z 0 and |z| > z 0 . In particular, it can be shown that there appear two discrete eigenvalues and, thus, two small-amplitude dark solitons are generated.
3.5.3. The quantum-state engineering method. A combination of the phase-imprinting and density-engineering methods is also possible, leading to the so-called quantum-state engineering technique [129,130]. This method, which involves manipulation of both the BEC density and phase, has been used in experiments at JILA [48] and Hamburg [67] with a two-component 87 Rb BEC (see Sec. 6.1 below): in the one component, a so-called "filled" dark soliton was created, with the hole in this component being filled by the other component. Depending on the trap geometry, the created filled dark soliton was found to be either unstable or stable. Particularly, in the JILA experiment [48], the dark soliton evolved in a quasi-spherical trap (after the filling from the other component was selectively removed) and, due to the onset of the so-called snaking instability, the soliton was found to decay into vortex rings (see Sec. 5.1 below). On the other hand, in the Hamburg experiment [67], the filled dark soliton in the one component was allowed to evolve (in the presence of the other component) in an elongated cigar-shaped trap; this way, a so-called dark-bright soliton pair was created (see Sec. 6.1), which was found to be stable, performing slow oscillations in the trap as predicted in theory [131].
Notice that a similar two-component engineering technique was also used for the creation of vortices [132], while earlier experimental results from the JILA group could be interpreted as a formation of a stack of filled dark solitons in a single BEC [133].
3.6. Multiple dark solitons and dark soliton interactions 3.6.1. The two-soliton state and dark soliton collisions. Apart from the single dark soliton solution, the NLS Eq. (26) possesses exact analytical multiple dark soliton solutions, which can be found by means of the IST [12,20] (see also Refs. [21,22]). Such solutions describe the elastic collision between dark solitons as, in the asymptotic limit of t → ±∞, the multiple-soliton solution can be expressed as a linear superposition of individual single-soliton solutions, which remain unaffected by the collision apart from a collision-induced phase-shift. To be more specific, let us consider the two-soliton wave function ψ = ψ(z, t), which can be asymptotically expressed as: where z ± 1,2 denote the position of each individual soliton (in the above expressions, the parameters A j and B j (j = 1, 2), with A 2 j + B 2 j = 1, characterize the velocity and depth of the soliton j). Apparently the shape and the parameters of each soliton are preserved, while the phase-shift of each soliton is given by: Note that if the soliton velocities are equal, i.e., then the phase-shift is equal for both solitons and is given by ∆z = (2B) −1 ln(1+B 2 /A 2 ). Equations (67)- (68) show that the spatial shift of each soliton trajectory is in the same direction as the velocity of each individual soliton and, thus, the dark solitons always repel each other. Here it should be mentioned, however, that this important result (as well as the collision dynamics near the collision point) can better be understood upon studying the explicit form of the two-soliton wave function rather than its asymptotic limit considered above. To do so, we consider again the case of a two-soliton solution, assuming for simplicity that the two solitons are moving with equal velocities (i.e., In such a case, the two-soliton wave function is given by [21,22]: where while n min = n 0 − n 0 B 2 = n 0 A 2 is the minimum density (i.e., the density at the center of each soliton). The density profile of the two-soliton solution in Eq. (69) is sketched in the top panel of Fig. 2.
To study analytically the interaction and collision between dark solitons, we follow the approach of Ref. [71] and find, at first, the trajectory of the soliton coordinate z 0 as a function of time: using the auxiliary equation ∂ z |ψ| 2 = 0 † [where the density |ψ| 2 is determined by Eq. (69)], the following result is obtained: Then, Eq. (72) determines the distance 2z * 0 between the two solitons at the point of their closest proximity, i.e., the collision point corresponding to t = 0: This equation [which holds for n min /n 0 = ν 2 ≤ 1/4, otherwise Eq. (73) provides a complex (unphysical) value for z * 0 ] shows that z * 0 = 0 for n min /n 0 = A 2 = 1/4. Thus, it † Recall that the dark soliton coordinate z 0 is the location of the minimum density (see Fig. 2).
is clear that there exists a critical value of the soliton velocity, namely v c = 1 2 √ n 0 ≡ 1 2 c s , which defines two types of dark solitons, exhibiting different behavior during their collision: "low-speed" solitons, with v < v c , which are reflected by each other, and "highspeed" solitons, with v > v c , which are transmitted through each other. In fact, as shown in the bottom panels of Fig. 2, the density profile of the low-speed (high-speed) twosoliton state exhibits two separate minima (a single non-zero minimum) at the collision point, namely n(z * 0 , 0) = 0 (n(z * 0 , 0) = 0) † . In other words, low-speed solitons are, in fact, well-separated solitons, which can always be characterized by two individual density minima -even at the collision point -while high-speed solitons completely overlap at the collision point. According to the nomenclature of Ref. [134], the collision between slow-speed (high-speed) solitons is called "black collision" ("gray collision"), since the dark solitons become black (remain gray) at t = 0. Notice that the case of gray collision can effectively be described -in the small-amplitude approximationby the collision dynamics of the KdV equation [134].
3.6.2. The repulsive interaction between slow dark solitons. Let us now investigate in more detail the case of well-separated solitons, which are always reflected by each other, with their interaction resembling the one of hard-sphere-like particles. In particular, we consider the limiting case of extremely slow solitons, i.e., n 0 /n min = A 2 ≪ 1 4 , for which the soliton separation is large for every time (i.e., z * 0 ≫ 0); in this case, the second term in the right-hand side of Eq. (72) is much smaller than the first one and can be ignored. This way, the soliton coordinate is expressed as: The above equation yields the soliton velocities: which, in the limit t → 0, become dz 0 /dt = 0. Thus, as the dark solitons approach each other, their depth (velocity) is increased (decreased), and become black at the collision point (t = 0), while remaining at some distance away from each other. Afterwards, the dark solitons are reflected by each other and continue their motion in opposite directions, with their velocities approaching the asymptotic values dz 0 /dt = ± √ n 0 A for t → ±∞ Next, differentiating Eq. (74) twice with respect to time, and using Eq. (72) (without the second term which is negligible for well-separated solitons), one may derive an equation of motion for the soliton coordinate in the form d 2 z 0 /dt 2 = −∂V int (z 0 )/∂z 0 , where the interaction potential V int (z 0 ) is given by: It is clear that V int is a repulsive potential, indicating that the dark solitons repel each other. If the separation between the dark solitons is sufficiently large (i.e., 2z 0 ≫ 1) then the hyperbolic sinh function in Eq. (76) can be approximated by its exponential asymptote, and the potential in Eq. (76) can be simplified as: The latter expression can also be derived by means of a Lagrangian approach [25]. Importantly, although the above result refers to a symmetric two-soliton collision, the results of Ref. [71] show that it is possible to use the repulsive potential (76) in the cases of non-symmetric collisions -using an "average depth" of the two solitonsand multiple dark solitons -with each soliton interacting with its neighbors (see also relevant discussion in Sec. 5.4).

Experiments on multiple dark solitons.
Multiple dark solitons were first created in a 23 Na BEC in the NIST experiment [46] by the phase-imprinting method (see Sec. 3.5.1), while the interaction and collision between two dark solitons in a 87 Rb BEC was first studied in the Hannover experiment of Ref. [49]. Nevertheless, in this early experiment the outcome of the collision was not sufficiently clear due to the presence of dissipation caused by the interaction of the condensate with the thermal cloud. In the more recent Hamburg experiment [68], the phase-imprinting method was also used to create two dark solitons in a 87 Rb BEC with slightly different depths. These solitons propagated to opposite sides of the condensate, reflected near the edges of the BEC, and subsequently underwent a single "gray" collision near the center of the trap. In addition, in the recent Heidelberg experiment [69] two dark solitons were created in a 87 Rb BEC by the so-called interference method (see Sec. 6.2 below). The solitons observed in this experiment, which were "well-separated" ones, propagated to opposite directions, reflected and then underwent multiple genuine elastic "black" collisions, from which the solitons emerged essentially unscathed. Notice that the experimentally observed dynamics of the oscillating and interacting dark-soliton pair of Ref. [69], as well as the one of multiple dark solitons in another Heidelberg experiment [71], was in a very good agreement with theoretical predictions based on the effective particle-like picture for dark solitons (see Sec. 4.2 and Sec. 5.4 below) and the interaction potential of Eq. (76).
3.6.4. Stationary dark solitons in the trap. At this point, it is relevant to briefly discuss the case where multiple dark solitons are considered in a trapped condensate. In this case, both the single dark soliton and all other multiple dark soliton states can be obtained in a stationary form from the non-interacting (linear) limit of Eq. (24), i.e., in the absence of the nonlinear term. In this case, Eq. (24) is reduced to a linear Schrödinger equation for a confined single-particle state. For the harmonic potential of Eq. (9), this Schrödinger equation describes the quantum harmonic oscillator, characterized by discrete energies and corresponding localized eigenmodes in the form of Hermite-Gauss polynomials [52]. As shown in Refs. [50,51], all these eigenmodes exist also in the fully nonlinear problem, and describe an analytical continuation of the above mentioned linear modes to a set of nonlinear stationary states. Additionally, analytical and numerical results of the recent work [135] suggest that in the case of a harmonic trapping potential there are no solutions of the 1D GP Eq. (24) without a linear counterpart. This actually means that interatomic interactions (i.e., the effective mean-field nonlinearity in the GP model) transforms all higher-order stationary modes into a sequence of stationary dark solitons confined in the harmonic trap [50,51]; note that as concerns its structure, this chain of, say n, stationary dark solitons shares the same spatial profile with the linear eigenmode of quantum number n. From a physical point of view, multiple stationary dark soliton states exist due to the fact that the repulsion between dark solitons is counter-balanced by the restoring force induced by the trapping potential.

General comments.
We consider again the quasi-1D setup of Eq. (24), but now incorporating the external potential V (z). In this setting, the dynamics of matter-wave dark solitons can be studied analytically by means of various perturbation methods, assuming that the trapping potential V (z) is smooth and slowly-varying on the soliton scale. This means that in the case, e.g., of the conventional harmonic trap [cf. Eq. (25)], the normalized trap strength is taken to be Ω ∼ ǫ, where ǫ ≪ 1 is a formal small (perturbation) parameter.
In such a case, Eq. (24) can be expressed as a perturbed NLS equation, namely, Then, according to the perturbation theory for solitons [136], one may assume that a perturbed soliton solution of Eq. (78) can be expressed in the following general form, Here, ψ s (z, t) has the functional form of the dark soliton solution (33), but with the soliton parameters depending on time, and ψ r is the radiation -in the form of sound waves -emitted by the soliton. Generally, the latter term is strong only for sufficiently strong perturbations (see, e.g., Refs. [137,138], as well as Ref. [73] and discussion in Sec. 4.4). Thus, the simplest possible approximation for a study of matter-wave dark solitons in a trap corresponds to the so-called adiabatic approximation of the perturbation theory for solitons [136], namely ψ(z, t) ≈ ψ s (z, t). In any case, the study of matter-wave dark solitons in a trap should take into regard that the trap changes the boundary conditions for the wave function, and BEC density, namely n → 0 [instead of n → n 0 in the homogeneous case -see, e.g., Eq. (36)] as z → ±∞. From a physical viewpoint, and based on the particle-like nature of dark solitons (see Sec. 3.3), one should expect that dark solitons could be reflected from the trapping potential; apparently, such a mechanism should then result in an oscillatory motion of dark solitons in the trap.
There exist many theoretical works devoted to the oscillations of dark solitons in trapped BECs. The earlier works on this subject reported that solitons oscillate in a condensate confined in a harmonic trap of strength Ω, and provided estimates for the oscillation frequency. In particular, in Ref. [139] soliton oscillations were observed in simulations and a soliton's equation of motion was presented without derivation; in the same work, it was stated that the solitons oscillate with frequency Ω (rather than the correct result which is Ω/ √ 2 -see below). The same result was derived in Ref. [140], considering the dipole mode of the condensate supporting the dark soliton. Other works [141][142][143] also considered oscillations of dark solitons in trapped BECs. An analytical description of the dark-soliton motion, and the correct result for the soliton oscillation frequency, Ω/ √ 2, were first presented in Ref. [144] by means of a multipletime-scale boundary-layer theory (this approach is commonly used for vortices [53]). The same result was obtained in Refs. [112,145] by solving the BdG equations (for almost black solitons performing small-amplitude oscillations around the trap center -see Sec. 4.3 below), using a time-independent version of the boundary-layer theory. Furthermore, in Ref. [112] a kinetic-equation approach was used to describe dissipative dynamics of the dark soliton due to the interaction of the BEC with the thermal cloud.
Matter-wave dark soliton dynamics in trapped BECs was also analyzed in other works by means of different techniques that were originally developed for optical dark solitons [34]. In particular, in Ref. [146] the problem was analyzed by means of the adiabatic perturbation theory for dark solitons devised in Ref. [28], in Ref. [147] by means of the small-amplitude approximation (see Sec. 3.4), while in Ref. [148] by means of the perturbation theory of Ref. [29]. Later, in Refs. [149,150] the so-called "Landau dynamics" approach was developed, based on the use of the renormalized soliton energy [cf. Eq. (43)], along with a local density approximation. Models relevant to the dynamics of matter-wave dark solitons in 1D strongly-interacting Bose gases, were also considered and analyzed by means of the small-amplitude approximation [151,152] (see also work for dark solitons in this setting in Refs. [153][154][155][156][157]). In other works, a Lagrangian approach for matter-wave dark solitons was presented [158] (see also Ref. [159]), and an asymptotic multi-scale perturbation method was used to describe dark soliton oscillations and the inhomogeneity-induced emission of radiation [160]. Recently, the motion of dark solitons was rigorously analyzed in Ref. [161] (where a wider class of traps was considered), while in Ref. [162] the same problem was studied in the framework of a generalized NLS model.
Finally, as far as experiments are concerned, the oscillations of dark solitons were only recently observed in the Hamburg [67,68] and Heidelberg [69,71] experiments. In these works, the experimentally determined soliton oscillation frequencies were found to deviate from the theoretically predicted value Ω/ √ 2. This deviation was explained in Refs. [67,68] by the anharmonicity of the trap, while in Refs. [69,71] by the dimensionality of the system and the soliton interactions (see also Sec. 5.4 below).

Adiabatic dynamics of matter-wave dark solitons
4.2.1. The perturbed NLS equation. The adiabatic dynamics of dark matter-wave solitons may be studied analytically by means of the Hamiltonian [27,28] or the Lagrangian [25] approach of the perturbation theory for dark solitons, which were originally developed for the case of a constant background. These approaches were later modified (see Ref. [146] for the Hamiltonian approach and Refs. [158,159] for the Lagrangian approach) to take into regard that, in the context of BECs, the background is inhomogeneous due to the presence of the external potential. The basic steps of these perturbation methods are: (a) determine the background wave function carrying the dark soliton, (b) derive from Eq. Here, we will present the first two steps of the above approach and, in the following two subsections, we will describe the adiabatic soliton dynamics in the framework of the Hamiltonian and Lagrangian approaches.
We consider again Eq. (78) and seek the background wave function in the form, where µ is the normalized chemical potential, θ o is an arbitrary phase, while the unknown real function Φ(z) satisfies the following equation, Then, we seek for a dark soliton solution of Eq. (78) on top of the inhomogeneous background satisfying Eq. (81), namely, ψ = Φ(z) exp(−iµt + iθ o )ψ s (z, t), where the unknown wave function ψ s (z, t) represents a dark soliton. This way, employing Eq. (81), the following evolution equation for the dark soliton wave function is readily obtained: It is clear that if the trapping potential V (z) is smooth and slowly-varying on the soliton scale, then the right-hand-side, and also part of the nonlinear terms of Eq. (82), can be treated as a perturbation. To obtain this perturbation in an explicit form, we use the TF approximation to express the background wave function as Φ(z) = 1 − V (z) [see Eq. (10) for g = 1 and µ = 1 † ] and approximate the logarithmic derivative of Φ as This way, Eq. (82) leads to the following perturbed NLS equation, where the perturbation Q(ψ s ) is given by: Then, considering an adiabatic evolution of the dark soliton, we assume that in the presence of the perturbation the dark soliton parameters become slowly-varying unknown functions of t [27,28,146]. Thus, the soliton phase angle becomes φ → φ(t) and, as a result, the soliton coordinate becomes In the latter expression, the dark soliton center z 0 (t) is connected to the soliton phase angle through the following equation: The evolution of the soliton phase angle can be found by means of the evolution of the renormalized soliton energy. In particular, employing Eq. (43) (for µ = 1), it is readily found that dE ds /dt = −4 cos 2 φ sin φ(dφ/dt). On the other hand, using Eq. (84) and its complex conjugate, it can be found that the evolution of the renormalized soliton energy is given by Then, the above expressions for dE ds /dt yield the evolution of φ, namely [28]: Next, we Taylor expand the potential V (z) around the soliton center z 0 , and assume that the dark soliton is moving in the vicinity of the trap center, i.e., µ ≡ 1 ≫ V , which means that the last two terms in the right-hand side of Eq. (85) can be neglected. This way, one may further simplify the expression for the perturbation in Eq. (85) which, when inserted into Eq. (88), yield the following result: To this end, combining Eq. (89) with Eq. (87), we obtain the following equation of motion for nearly stationary (black) solitons with cos φ ≈ 1, The above result indicates that the dark soliton center can be regarded as a Newtonian particle: Eq.
or, in physical units, with ω z / √ 2. An example of an oscillating matter-wave dark soliton is shown in Fig. 3. At this point, it is relevant to follow the considerations of Ref. [112] (see also [144,145]) and estimate the energy E ds of this almost dark soliton in the trap. Taking into regard that in the case of a homogeneous BEC this energy is given by Eq. (47), one may use a local density approximation and use in Eq. (43) the local speed of sound, c(z) = n 0 (z) [79] (here, n 0 (z) is the density of the ground state of the BEC), rather than the constant value c s = √ n 0 [cf. Eq. (38)]. Then, in the TF limit, the density is expressed as n 0 (z) = µ − 1 2 Ω 2 z 2 = c 2 s − 1 2 Ω 2 z 2 and, thus, one may follow the lines used for the derivation of Eq. (47) (for sufficiently slow solitons and weak trap strengths) and obtain the result: where E 0 ≡ 4 3 c 3 s and m ds = −4 √ n 0 as in Eq. (47). The above equation shows that the incorporation of the harmonic trap results in a decrease of the energy of the dark soliton by the potential energy term 1 4 |m ds |Ω 2 z 2 . Moreover, the ratio of the soliton mass over this potential energy is given by (Ω 2 z 2 /4) −1 , which is exactly two times the ratio of the atomic mass (which is equal to m = 1 in our units) over the external potential, namely (Ω 2 z 2 /2) −1 . This is another interpretation of the result that the effective mass of the dark soliton center is M eff = 2.
Here, A and B are unknown slowly-varying functions of time (with A 2 + B 2 = 1) representing, respectively, the velocity and amplitude of the dark soliton (which become time-dependent due to the presence of the perturbation), while ζ ≡ B(t) [z − z 0 (t)], where z 0 (t) is the dark soliton center. Note that in the unperturbed case, dz 0 /dt ≡ A, but in the perturbed case under consideration, this simple relationship may not be valid (see below). Next, the evolution of the unknown soliton parameters α j (t) (which is a generic name for z 0 (t) and A(t)) are obtained via the Euler-Lagrange equations [25,158]: whereα j ≡ dα j /dt and L ds = +∞ −∞ dzL ds {ψ s } represents the averaged Lagrangian of the dark soliton of the unperturbed NLS equation (namely for Q(ψ s ) = 0), with the Lagrangian density L ds being given by Eq. (44) (for n 0 = 1). The averaged Lagrangian can readily be obtained by substituting the ansatz (93) into Eq. (44): Therefore, substituting Eqs. (95) and (85) into Eq. (94), it is straightforward to derive evolution equations for the soliton parameters. For completeness, we will follow Ref. [158] and present the final result taking also into account the last two terms in the right-hand side of Eq. (85) -which were omitted in the previous subsection -so as to describe the motion of shallower solitons as well. This way, and employing a Taylor expansion of the potential around the soliton center (as in the previous subsection), we obtain the following evolution equations for z 0 (t) and A(t): Equations (96) This result is rigorously proved by means of the Landau dynamics approach that will be discussed below.

Landau dynamics of dark solitons.
The oscillations of dark solitons of arbitrary amplitudes in a trap can also be studied by means of the so-called Landau dynamics approach devised in Refs. [144,150]. This approach, which further highlights the particlelike nature of the matter-wave dark solitons, relies on a clear physical picture: when a dark soliton moves in a weakly inhomogeneous background, its local energy stays constant. Hence, one may employ the local density approximation, and rewrite the energy conservation law of Eq. (43) as c 2 (z 0 ) − v 2 = (3E ds /4) 2/3 , where c 2 (z 0 ) is the local speed of sound evaluated at the dark soliton center z 0 . Then, in the TF limit, one has c 2 (z 0 ) = c 2 s − 1 2 Ω 2 z 2 0 (as before), and taking into regard that the soliton velocity is v = dz 0 /dt, the following equation for the energy of the dark soliton is readily obtained: whereẼ ds = c 2 s − (3E ds /4) 2/3 and the effective mass of the dark soliton center is again found to be M eff = 2. It is readily observed that Eq. (98) can be reduced to Eq. (90) and, thus, it leads to the oscillation frequency of Eq. (91). Nevertheless, the result obtained in the framework of the Landau dynamics approach is more general, as it actually refers to dark solitons of arbitrary amplitudes. Moreover, the same approach can be used also in the case of more general models, including, e.g., the cases of non-harmonic traps and/or more general nonlinearity models, such as the physically relevant ones described by Eqs. (19)-(22) [150]. Nevertheless, it should be noted that in such more general cases the problem can be treated analytically for almost black solitons (v ≪ c), performing small-amplitude oscillations. In this case, the conservation law E ds (c(z 0 ), v) = E 0 can be Taylor expanded around z 0 = 0 and v = 0, leading to expressions for the soliton's effective soliton mass and oscillation frequency [150]. 4.2.5. The small-amplitude approximation. Next, we discuss the adiabatic dynamics of small-amplitude dark solitons in trapped 1D Bose gases. In this case, one may formally reduce the more general GP model of Eq. (50) (including the potential term V ψ) to a KdV equation with variable coefficients -see Ref. [44] for details and Ref. [163] for applications of this KdV model. The main result of such an analysis is that the density and the phase of the approximate shallow dark soliton solution of Eq. (50) have, to the leading-order of approximation, the functional form of their counterparts in Eqs. (61) and (63), but with the soliton parameter κ depending on a slow variable, say Z (see earlier work for a calculation of κ(Z) in Refs. [164,165]). This way, approximate analytical shallow soliton solutions have been found in various works [134,147,151,152] for different forms of the nonlinearity function f (n). Nevertheless, there are some subtle issues concerning the validity of this approximation, as discussed in Refs. [148,150,160], which is, strictly speaking, valid away from the turning points (where the soliton velocity vanishes). On the other hand, numerical results (see, e.g., Ref. [152]), illustrate that the range of validity of the above results is, in fact, wider than what may be expected based on the limitations of this approach.
The small-amplitude approximation, along with a local-density approximation, has also been used to estimate the shallow soliton's oscillation frequency: the shallow soliton's velocity v, which in the homogeneous problem was found to be close to the speed of sound, namely v ≈ c s = f ′ 0 n 0 [see Eq. (62)], can be approximated in the inhomogeneous system as follows: In some cases, Eq. (99) can be used for the derivation of physically relevant results. For example, following Ref. [150], we assume that f (n) = n α , where α = 1 for weaklyinteracting BECs, or α = 2 for strongly-interacting Tonks gases [101]. Then, in the TF limit, n 0 (Z) = [µ − V (Z)] 1/α , and Eq. (99) is reduced to the form dZ/ µ − V (Z) = √ αdt. The latter is integrated and yields [for V (Z) = (1/2)Ω 2 Z 2 ] the soliton trajectory: where L = √ 2µ/Ω is the length of the TF cloud. Equation (100) predicts that the shallow dark soliton will perform oscillations approximately in the entire spatial region occupied by the gas, with an oscillation frequency which takes the following values (in physical units): for α = 1, i.e., for quasi-1D BECs described by the cubic GP equation, ω osc = ω z / √ 2, while for α = 2, i.e., for the Tonks gas described by a quintic NLS equation, ω osc = ω z . Note that the latter result was first obtained via a many-body calculation [166], and later was derived by means of the KdV approximation [151].  91), can also be derived by means of a BdG analysis as was first demonstrated in Ref. [145] (see also results in Refs. [167][168][169]). Such an analysis can be done in the TF limit for the stationary dark soliton state, namely the black soliton ψ 0 [see Eq. (37)] located at the trap center, i.e., z 0 = 0, which is actually the first excited state of the condensate. Then, following the discussion in Sec. 2.3, the excitation spectrum can be found as follows: using the ansatz ψ(z, t) = [ψ 0 (z) + u(z)e −iωt + υ * (z)e iωt ] e −iµt [where ω = ω r + iω i is a (generally complex) eigenfrequency and (u, v) are perturbation eigenmodes], we derive from Eq. (24) the following BdG equations: whereĤ = −(1/2)∂ 2 z + (1/2)Ω 2 z 2 is the single particle operator. A typical example showing the initial configuration, i.e., the condensate and the stationary dark soliton (which can be found, e.g., by a Newton-Raphson method), as well as the corresponding spectral plane (ω r , ω i ), are shown in Fig. 4.
The BdG analysis reveals that all the eigenfrequencies of the spectrum are real, which indicates that the stationary dark soliton is dynamically stable. The four smallest magnitude eigenfrequency pairs † and their corresponding eigenmodes have the following physical significance (see, e.g., Ref. [42]). First, there exists a zero eigenfrequency, located at the origin of the spectral plane (ω r , ω i ), which reflects the phase invariance of the 1D GP equation. The respective eigenfunction is the so-called Goldstone mode and does not result in any physical excitation (oscillation) of the system. The solutions with eigenfrequencies ω r = ±Ω correspond to the so-called dipole mode (or Kohn mode), which is relevant to the motion of the center of mass of the system; note that as the system is harmonically confined, the center of mass oscillates with the frequency of the harmonic trap [170]. The solutions with eigenfrequencies having the next larger magnitude eigenfrequencies correspond to the quadrupole mode, with the location of the eigenfrequencies at ω r = ± √ 3Ω, being particular to the onedimensionality of the system [88]. Note that the excitation of the quadrupole mode (induced, e.g., by a time-modulation of the trap strength) results in a breathing behavior of the condensate, with its width oscillating with the above-mentioned frequency.
Of particular interest are the solutions with eigenfrequencies ω r = ω A ≡ Ω/ √ 2, which correspond to the so-called anomalous mode. This mode appears in the Bogoliubov analysis only when topological excitations of the condensate are involved, namely dark solitons or vortices [74]. A characteristic property of the anomalous mode is that the integral of the norm × energy product, (|u| 2 − |v| 2 )ωdz (in our units), is negative rather than positive as is the case for all the positive frequency modes associated with the ground state of the system [42]. Note that, from a mathematical viewpoint, the anomalous mode possesses a topological property of the so-called negative Krein signature [171], namely K ≡ sign{ (|u| 2 − |v| 2 )ωdz} < 0 (for positive eigenfrequencies ω). Practically, this means that the anomalous mode becomes structurally unstable (i.e., it becomes complex) upon collision with other eigenvalues, as is the case when dissipation is present [172]. In our case, finite temperature automatically implies the presence of dissipation which, as discussed in more detail in Sec. 6.5 below, may be described -in the simplest possible approach -by including a phenomenological temperature-induced damping term in the GP model.
In order to better clarify the above, and discuss in more detail the stability of the excitation corresponding to the anomalous mode, namely of the dark soliton, we note the following. At temperatures T → 0 (as is the case under consideration), the negative energy of the dark soliton does not imply any instability (e.g., a decay process) and, thus, the soliton is dynamically stable. Nevertheless, at finite temperatures, i.e., in the presence of a thermal cloud, the above mentioned properties of the anomalous mode indicate that the soliton will become unstable: in this case, the presence of the temperature-dependent damping results in the decay of the soliton (see discussion in Refs. [112,145] as well as in Sec. 6.5 below). From a physical point of view, the decay mechanism resembles the one of the low-energy excitations of trapped BECs [173] and originates from the scattering of thermal particles on the dark soliton. Thus, according to these arguments, matter-wave dark solitons can be regarded as thermodynamically unstable excitations as, in the presence of the temperature-induced dissipation, the system will be driven towards configurations with lower energy; in other words, the dark solitons will decay to the ground state. This scenario is also often referred to as energetic instability [174].
As mentioned above, the eigenfrequency of the anomalous mode is equal to the oscillation frequency of a dark soliton around the center of the harmonic trap in the TF limit. On the other hand, the eigenfunctions u A and υ A of the anomalous mode, shown in Fig. 5, are localized within the notch of the dark soliton [145,169], and their sum can be approximated as u A + υ A ∝ sech 2 ( √ n 0 z). Notice that in the case of a uniform condensate (i.e., in the absence of the trap), there exists a translational mode with zero frequency which has the same functional form, namely ∂ z ψ 0 , due to the translational invariance of the dark soliton solution. When the trap is present, however, this symmetry is broken, which suggests that the anomalous mode can be regarded as the "ghost" of the broken translational invariance of the dark soliton solution. We finally mention that the direct connection of the anomalous mode to the oscillation of the dark soliton can be better explained by the fact that an excitation of the stationary black soliton ψ 0 by the anomalous mode results in a displacement of the soliton from the trap center. In other words, the GP Eq. (24) with the initial condition ψ(z; t = 0) = ψ 0 (z) + u A (z) + υ * A (z), will naturally lead to dark soliton oscillations studied in Sec. 4.2.

4.3.2.
The multi-soliton state. The BdG analysis can also be performed in the more general case of multi-soliton states [175], which may be found in a stationary form (as explained in Sec. 3.6). In this case, starting from the non-interacting limit, it can be found that the Bogoliubov spectrum of the n-th excited state consists of one zero eigenvalue (corresponding to the Goldstone mode), n double eigenvalues (accounted for by the presence of the harmonic trap), and infinitely many simple ones. Then, in the nonlinear regime, one of the eigenvalues of each double pair becomes an anomalous mode of the system (characterized by a negative valued integral of the norm × energy product) and, thus, the number of anomalous modes in the excitation spectrum equals to the number of dark solitons [175]. This is in agreement with the fact that the number of eigenvalues with negative Krein signature equals to the number of the nodes of the stationary state [176]. Notice that in the framework of the 1D GP Eq. (24) -i.e., in the TF-1D regime -the first anomalous mode coincides with the oscillation frequency ω osc = Ω/ √ 2 of the single dark soliton. An example of a condensate with a stationary two-dark soliton state, as well as the pertinent spectral plane, are shown in Fig. 6.
The physical significance of the n-anomalous modes has been discussed in Refs. [71,175]: for example, in the case of a two-dark soliton state, the smallest of the anomalous modes corresponds to an in-phase oscillation (i.e., when the two dark solitons oscillate together without changing their relative spatial separation), the largest anomalous mode corresponds to an out-of-phase oscillation (i.e., when the two dark solitons move to opposite directions with the same velocity and undergo head-on collision), and so on. An example of the in-phase and the out-of-phase oscillation of the two-dark soliton state shown in the left panel of Fig. 6 (when the solitons are properly displaced from their equilibrium positions) is illustrated in Fig. 7. Here it should be noticed that since the starting point of the above considerations is the non-interacting limit, a similar analysis can also be performed in the case of other mean-field models with non-cubic nonlinearity, as the ones describing cigar-shaped BECs in the dimensionality crossover regime from 3D to 1D (see Sec. 2.4 and Sec. 5.4 below).
Finally, as concerns the stability of nonlinear modes (see relevant investigations in Refs. [51,71,177,178]), all higher-order nonlinear modes are unstable near the noninteracting limit [71,178] (but can be stabilized by using anharmonic traps [178]); nevertheless, the instability ceases to exist sufficiently deep inside the nonlinear regime (i.e., for sufficiently large BECs, with large N) [71].

Radiation effects: inhomogeneity-induced sound emission by the soliton.
As indicated in Sec. 4.1, the dark soliton experiences a background density gradient in the presence of the perturbation R(ψ) in Eq. (78) and, thus, it continuously emits energy in the form of sound waves. Here, we will study this effect in more detail, considering the case of R(ψ) = V (z)ψ with V (z) being a harmonic potential. A particular feature of this setup, is that the emitted sound energy remains confined within the spatial region of the trap and, hence, the soliton continuously re-interacts with the emitted sound waves [137]; in fact, this process is such that, on average, the dark soliton reabsorbs the radiation it emits. Thus, in the case of harmonic traps, an investigation of the inhomogeneity-induced sound emission, as well as an estimation of the rate of emission of energy, is relevant for short timescales, i.e., 0 ≤ t ≤ T osc ≡ 2π/ω osc [where ω osc is given by Eq. (91), with Ω ≪ 1]. On the other hand, if dark solitons evolve in the presence of non-harmonic potentials (such as localized barriers [73,138,146,159,179], disordered potentials [180,181], anharmonic traps [182], or other "properly designed" potentialssee below), the problem may be easier -at least in terms of a numerical investigation: in fact, as is explained below, it is possible to consider suitable setups that either damp off the emitted sound density or cause the emitted sound to dephase.
Various such setups were proposed and analyzed in the past; the most prominent example is, perhaps, a tight inner "dimple" trap, confining a dark soliton, located within a much weaker outer harmonic potential (such a configuration can be realized by focusing an off-resonant laser beam within a harmonic trap) [137]. In this case, if the depth of the dimple trap is sufficiently small, the sound waves can escape (to the outer trap), while the soliton can remain confined in this region. In this limit, sound energy is damped off for short enough timescales, until it bounces off the weaker outer trap and thus becomes forced to re-interact with the soliton. An alternative setup considered in Refs. [183,184] consists of a harmonic trap perturbed by an optical lattice potential (see Sec. 6.4). In this case, the optical lattice can confine a soliton within a single or a few lattice sites, with the sound (again for short enough times) escaping to neighboring sites. Although in this case the sound re-interacts with the soliton on faster timescales than in the case of the dimple trap mentioned above, the presence of the lattice dephases the emitted sound waves, and hence accelerates the soliton decay.
The radiation-induced dissipation of matter-wave dark solitons in harmonic traps was studied analytically in Ref. [160] by means of an asymptotic multi-scale expansion method. Particularly, assuming that R(ψ) = ǫ 2 z 2 ψ (with ǫ being a formal small parameter defined by the aspect ratio Ω), the following results were obtained. In the limit ǫ → 0, the dark soliton evolves adiabatically so that the dark soliton center z 0 (t) = vt + z o → s(T )/ǫ, i.e., it becomes a function of the slow timescale T = ǫt, while the soliton velocity is given by v(T ) =ṡ ≡ ds/dT . The adiabatic dynamics is followed by generation of sound waves, which can be taken into regard as per Eq. (79). In fact, in the decomposition of the wave function ψ into an inner and an outer asymptotic scale, the leading-order radiative effects are taken into account when the complex phase θ o [see Eq. (80)] depends also on T = ǫt, i.e., θ o ≡ θ(T ), and the first-order corrections to the dark soliton (33) grow linearly in z. Neglecting reflections from the trap, the extended dynamical equation for the position s(T ) of the dark soliton (33) takes the form: The left-hand-side of Eq.  103) is not applicable. Next, letting E = 1 2 (ṡ 2 + s 2 ) be the energy of the harmonic oscillator, one may employ Eq. (103) to calculate the rate of change of E due to the leading-order radiative effects appearing in the right-hand side of Eq. (103). The result is: and shows that due to the energy pumping (104), the amplitude of the harmonic oscillator increases in time. In the limit s 2 +ṡ 2 → 0, Eqs. (103) and (104) can be simplified. First, the energy of the dark soliton oscillations accelerates by the squared lawĖ = ǫṡ 2 /2, which was confirmed in numerical simulations in the setup of Ref. [137]. Second, the nonlinear equation (103) is linearized as follows: Equation (105) includes an anti-damping term accounting for the emission of radiation, indicating that the center point (0, 0) becomes an unstable spiral point on the plane (s,ṡ). Apparently, the leading-order solution reads s(T ) = s 0 e ǫT /4 cos(T + δ 0 ); thus, the amplitude of oscillations of a dark soliton increases while its own amplitude decreases. The above results of the asymptotic analysis were confirmed by the numerical simulations of Ref. [160], but also by numerical findings reported in other works: the radiation-induced effects were also observed for dark solitons oscillating between two Gaussian humps [138], or for dark solitons that are parametrically driven by a pair of two periodically-modulated Gaussian barriers, oscillating in anti-phase at a frequency close to the soliton frequency [185]. It is worth noticing that the mechanism proposed in Ref. [185] pumps energy into the dark soliton, which may compensate the inhomogeneity-induced emission of radiation, as well as the damping due to the presence of the thermal cloud [112] (see Sec. 6.5 below).

Persistence and stability of dark solitons.
As was highlighted in this Section, there exist many alternative approaches for the study of the statics and dynamics of matter-wave dark solitons in the quasi-1D setting. Nevertheless, rigorous results concerning the persistence and stability of dark solitons in a generalized GP-like model [cf. Eq. (18)] were obtained only recently [161,186]. Particularly, in Ref. [186], the existence and stability of a black soliton of Eq. (50) were studied in the absence of the potential term, while in Ref. [161] a more general model, incorporating the potential term, was considered. More specifically, the model used in Ref. [161] was of the following form, where ε is a formal small parameter setting the strength of the potential. The results obtained in Ref. [161] for bounded and exponentially decaying potentials (as, e.g., ones corresponding to red-detuned laser beams -see, e.g., the experiment of Ref. [187]) can be summarized as follows.
Let us consider that, in the absence of the potential, Eq. (106) admits a black soliton solution of the form ψ(z, t) = q(z − s) exp[−if (n 0 )t + iθ] (here, s is the soliton center and θ an arbitrary constant phase), with boundary conditions q 0 → ± √ n 0 as z → ±∞. Then, the above solution persists in the presence of the perturbation induced by the potential term in Eq. (106) provided that the function possesses at least one single root, say s 0 . Then, the stability of the dark soliton solution depends on the sign of the first derivative of the function in Eq. (107), evaluated at s 0 : an instability occurs, with one imaginary eigenfrequency pair for εM ′′ (s 0 ) < 0, and with exactly one complex eigenfrequency quartet for εM ′′ (s 0 ) > 0. In fact, this instability is dictated by the translational eigenvalue, which bifurcates from the origin as soon as the perturbation is present. For εM ′′ (s 0 ) < 0, the relevant eigenfrequency pair moves along the imaginary axis, leading to an instability associated with exponential growth of a perturbation along the relevant eigendirection. On the other hand, for εM ′′ (s 0 ) > 0, although the eigenfrequency should move along the real axis, it can not do so because the latter is filled with continuous spectrum; thus, since the translation mode and the eigenmodes of the continuous spectrum have opposite Krein signature, the collision of the eigenfrequency of the translational mode with the continuous spectrum results in a complex eigenvalue quartet, signalling the presence of an oscillatory instability. The relevant eigenfrequencies can be determined by a quadratic characteristic equation which, in the case of the cubic GP model (106) with f (n) = n and n 0 = 1, takes the form [161], and the eigenvalues λ are related to the eigenfrequencies ω through λ 2 = −ω 2 . For sufficiently small ε > 0, this equation has only one real root λ(ε) > 0 for M ′′ (s 0 ) < 0 and two complex-conjugate roots, with Re{λ(ε)} > 0 for M ′′ (s 0 ) > 0. It is interesting to observe that if the characteristic equation (108) is formally applied to the cubic GP model (106) with f (n) = n, n 0 = 1 and V (z) = z 2 , one obtains M ′′ (s 0 ) = 2 +∞ −∞ sech 2 (z)dz = 4 and, thus, Eq. (108) takes the form λ 2 − (ε/2)λ + ε = O(ε 2 ). Using appropriate rescalings, it can easily be shown that the latter characteristic equation can be derived from Eq. (105) of Sec. 4.4. Although the validity of the radiative boundary conditions for V (z) = z 2 cannot be verified by the analysis of Ref. [161], the above observation leads to the following conjecture [161]: in the most general GP model [cf. Eq.(106)], the two complex-conjugate eigenvalues with positive real part for M ′′ (s 0 ) > 0 result from the following Newton's particle equation of motion for the soliton center s(t): where M(s) is the effective potential implied by Eq. (107), while the constants µ 0 and λ 0 represent, respectively, the soliton's mass and anti-damping -as per the discussion of Sec. 4.4. The validity of Eq. (109), as well as the other theoretical predictions presented in this Section, were tested against numerical simulations in Ref. [161] for small decaying potentials, and the agreement between the analytical and numerical results was found to be very good. Notice that although the above results of Ref. [161] can only be rigorously applied to the case of small, bounded and exponentially decaying potentials, the basic qualitative features may formally persist for other types of external potentials as well. A pertinent example is the work of Ref. [188], where the the persistence and stability of matter-wave black solitons were studied in a condensate characterized by a periodic, piecewise-constant scattering length † : as shown in Ref. [188], a formal application of the predictions of Ref. [161] concerning the persistence and stability of dark solitons in this setting, was found to be in very good agreement with relevant numerical findings. Nevertheless, an analysis similar to the one presented in Ref. [161], but for other types of potentials (such as confining and periodic ones) is still missing.

Matter-wave dark solitons in higher-dimensional settings
Quasi-1D matter-wave dark solitons may naturally exist in higher-dimensional settings. For example, in the experimentally relevant case of a cigar-shaped trap, the actual dimensionality of the BEC density is 3D rather than 1D, despite of the fact that the BEC can be treated as an effectively 1D object using the NLS Eq.  Fig. 8. Furthermore, quasi-1D dark solitons may also exist in disk-shaped BECs (see Sec. 2.4), which are described by the (2 + 1)-dimensional GP Eq. (23). The latter, can be expressed in the following dimensionless form, where ∇ 2 = ∂ 2 x + ∂ 2 y , the density |ψ| 2 , length, time and energy are respectively measured in units of 2 √ 2πaa z , a z , ω −1 z and ω z , while the potential V (r) is given by with the aspect ratio being Ω = ω ⊥ /ω z ≪ 1. In this case, the soliton of Eq. (33), with the variable z being replaced by x, is an exact analytical solution of Eq. (110) for V (z) = 0. This "rectilinear" soliton has the form of a dark "stripe" on top of a 2D TF cloud, and the BEC wave function can be expressed (similarly to the 1D case) as ψ = ψ TF (r) exp(−iµt)ψ ds (x, t). It is also natural to consider the full (3 + 1)-dimensional version of Eq. (110), with ∇ 2 = ∂ 2 x + ∂ 2 y + ∂ 2 z , where quasi-1D dark soliton stripes exist as well. In this case, the potential and the 3D TF cloud are modified according to the relative values of the confining frequencies in the three directions. Examples of the densities of a disk-shaped BEC and a spherical BEC carrying a rectilinear dark soliton are shown, respectively, in the middle and right panels of Fig. 8.
Apart from the quasi-1D dark solitons, purely 2D dark solitons have also been predicted to occur in theory (but they have not been observed so far in experiments). Such dark soliton solutions of the GP Eq. (110), which have been derived in the framework of the small-amplitude approximation (see Sec. 3.4), may have the form of lumps satisfying an effective Kadomtsev-Petviashvili (KP) equation [120] or dromions satisfying an effective Davey-Stewartson system [121]; quasi-1D and 2D dark solitons of † BECs with spatially-varying coupling constant g, so-called collisionally inhomogeneous condensates [189], have attracted much attention, as they provide a variety of interesting phenomena [190][191][192][193][194][195][196][197][198]. the dromion type, have also been predicted to occur in disk-shaped multi-component condensates [124].

Snaking instability of rectilinear dark solitons.
5.1.1. Basic phenomenology and results. An important issue arising in higherdimensional settings is the stability of dark solitons which, for simplicity, will be studied at first in the (2 + 1)-dimensional geometry (relevant to disk-shaped BECs) and in the absence of the potential V (r). The stability of the 1D dark soliton stripe (lying, say, along the x-direction) in such a 2D setting was first studied in Ref. [199] (see also Refs. [200,201]). In this work, it was shown that the soliton is prone to transverse modulational instability, i.e., it is unstable against long-wavelength transverse periodic perturbations ∼ cos(Qy), where Q is the wave number of the perturbation. In particular, the instability band is defined by Q < Q cr , where the critical value of the perturbation wave number is given by (for µ = 1): The above expression is a result of a linear stability analysis, which indicates that the amplitude of the rectilinear dark soliton will grow exponentially in the transverse direction. Nonlinear regimes of this instability were also studied analytically by means of asymptotic expansion techniques (see, e.g., Ref. [202] and references therein). This instability was extensively studied in the context of nonlinear optics, both theoretically [34] and experimentally [203,204], and was found to be responsible for a possible decay of a plane dark soliton into a chain of vortices of opposite topological charges (vortex-antivortex pairs). Particularly, when the transverse modulational instability sets in, a plane black soliton undergoes a transverse "snake" deformation (hence the name "snaking instability") [34,202], causing the nodal plane to decay into vortex pairs. On the other hand, unstable gray solitons may not decay into vortices, but rather perform long lived oscillations accompanied by emission of radiation in the form of sound waves. The basic phenomenology and results described above, persist in the case of other higher-dimensional setups as well. For example, in Fig. 9, we show the onset of the snaking instability of a rectilinear dark soliton on top of a cigar-shaped BEC (see also Sec. 5.3 below), confined in a trap of the form V (z, r) = (1/2)(Ω 2 r r 2 + Ω 2 z z 2 ) (with r 2 = x 2 + y 2 ). In such higher-dimensional settings, the unstable soliton collapses into more stable vortex structures, namely vortex rings. From the viewpoint of experimental observations, the snaking instability and the decay of matter-wave dark solitons into vortex rings was first observed in a JILA experiment [48] with a two-component 87 Rb BEC (see Sec. 6.1). In particular, a quasi-1D dark soliton created in one component (see Sec. 3.5.3) evolved in a quasi-spherical trap and, thus, the onset of the snaking instability caused the soliton to decay into vortex rings -as predicted in theory [167].

Avoiding the snaking instability.
As is known from the context of nonlinear optics [34], the snaking instability can be avoided by using finite-sized background optical beams (see, e.g., relevant experimental results in Ref. [3]). Thus, one should expect that the suppression of the snaking instability may also be possible in the case of a trapped (disk-shaped) condensate, which also constitutes a background of a finite extent. Indeed, in such a case, a simple criterion for the suppression of the snaking instability can be found by means of scale competition arguments [205]. In particular, if the characteristic length scale of the condensate L BEC ≡ 2 √ 2/Ω (i.e., the TF diameter for µ = 1) is below the critical length L c ≡ 2π/Q cr stemming from Eq. (112), then the snaking instability will not manifest itself. Considering the case of a black (stationary) soliton with sin φ = 0, Eq. (112) yields Q cr = 1 and, thus, L c = 2π; in such a case, the above scale competition argument, L BEC < L c , leads to the prediction that a use of a sufficiently strong trap, such that Ω > Ω c ≡ √ 2/π ≃ 0.45, can suppress the snaking instability. The above criterion was tested against direct numerical simulations [205], and it was found that the critical value of the trap strength is less than the theoretically predicted, namely Ω c ≈ 0.31. This discrepancy can be understood by the fact that for small BECs (i.e., for tight traps) the presence of the dark soliton significantly modifies the maximum density which is less than µ by a "rescaling" factor f , found to be f ≈ 0.5.
On the other hand, it was recently predicted [206] that stable 3D stationary dark solitons may exist in dipolar condensates (for this type of BECs see, e.g., the recent review [207] and references therein). In particular, the special feature of dipolar condensates, namely the dipole-dipole interaction, together with the use of a sufficiently deep optical lattice in the soliton's nodal plane, allows for the existence of dark solitons of arbitrarily large transversal sizes, which are not prone to the snaking instability. In this case, the underlying reason for the suppression of the snaking instability is that the dipole-dipole interaction is long-range (it decays like r −3 , where r is the interparticle distance), which means that the respective GP equation incorporates a nonlocal nonlinear term. Generally, such a nonlocal response may arrest collapse and stabilize solitons in higher-dimensions, as was shown in the context of optics (see, e.g., Ref. [208], as well as Ref. [209] for a relevant recent work on dark solitons).
We also note that a more "exotic" dark soliton configuration in the 2D setting, which is not subject to the snaking instability, was presented in Ref. [210] (see also Ref. [211]). This configuration, which refers to a two-component BEC (see Sec. 6.1), consists of a "cross" formed by the intersection of two rectilinear domain walls, with the wave functions of the same species filling each pair of opposite quadrants having a π phase difference. This way, a dark soliton configuration is formed, which was found to be stable for long times -and even in the presence of rotation of the trap -in a large parametric region.

5.2.
Matter-wave dark solitons of radial symmetry 5.2.1. Ring dark solitons (RDS) and spherical shell solitons (SSS). A special class of dark solitons in higher-dimensional settings consists of dark solitons exhibiting a radial symmetry, that can be realized by wrapping the nodal plane around on itself. Such structures were first introduced in the context of nonlinear optics [212], with the motivation being that these dark solitons may not be prone to the snaking instability: indeed, if a dark stripe is bent so as to form a dark ring of length L < 2π/Q cr (in the (2 + 1)-dimensional geometry), then the snaking instability will be suppressed. Such ring dark solitons (RDSs) were studied in nonlinear optics both theoretically [213,214] and experimentally [215][216][217], while later were also predicted to occur in BECs [218]. In this context, and in the 2D setting (i.e., in a disk-shaped BEC), the RDS has the form of an annular "trough". On the other hand, in a 3D setting (i.e., in a spherical BEC), the radially symmetric dark soliton is called spherical dark soliton [158], or spherical shell soliton (SSS) -according to the nomenclature of Ref. [219] -and has the form of a nodal spherical "shell". Examples of the 3D densities of a disk-shaped and a spherical BEC, carrying a RDS and a SSS, are shown in Fig. 10.
As was originally proposed in Ref. [218], RDS may be generated in BECs by means of phase-engineering techniques (i.e., by a proper phase-imprinting method -see Sec. 3.5.1 below), similar to the ones used for the generation of optical RDS [215][216][217]. Another technique that has been proposed for the generation of RDS in BEC is the matter-wave interference method (see Sec. 6.2): if the condensate is initially trapped in a narrow cylindrical box-like potential, and then is allowed to coherently expand in the presence of a wider cylindrical impenetrable hard-wall potential, it is reflected from the boundary, and the self-interference pattern has the form of a sequence of non-stationary concentric RDS [220][221][222] (see also relevant work in Ref. [223]).

Dynamics and stability of RDS and SSS.
From a mathematical standpoint, matter-wave dark solitons of radial symmetry can be considered as quasi-1D objects and, accordingly, be analyzed by means of a quasi-1D GP equation. In particular, either RDS or SSS can be described by Eq. (110), with the Laplacian being in the form, with D = 1, 2, 3. In this setup, the simplest case of D = 1 reduces Eq. (110) to the 1D GP Eq. (24) describing a quasi-1D BEC (here, r ≡ z). The higher-dimensional setups correspond to the cases of D = 2 and D = 3: in the former case, Eq. (110) describes a disk-shaped BEC in the (x, y) plane (with r given by r = x 2 + y 2 ), while in the latter case Eq. (110) describes a spherical BEC (with r given by r = x 2 + y 2 + z 2 ).
In such a quasi-1D setup, the Hamiltonian or the Lagrangian perturbation theory for dark solitons (see Sec. 4.2) may also be applied for the study of the dynamics of RDS and SSS. In particular, Eq. (110), with the Laplacian of Eq. (113), can be treated as a perturbed 1D NLS equation [similar to Eq. (78)] provided that both the potential term and the term ∼ r −1 can be considered as small perturbations. The latter assumption is physically relevant for radially symmetric solitons of large radius r 0 . Then, approximating the functional form of the RDS or SSS as [cf. Eq. (86)], where r 0 (t) is the soliton radius, it can be found [218] (see also Ref. [158]) that r 0 is governed by the following Newtonian equation of motion: Here, the effective potential is given by and V (r 0 ) = (1/2)Ω 2 r 2 0 is the trapping potential evaluated at the soliton radius r 0 . In the 1D limit of D = 1, the last term in the right-hand side of Eq. (116) vanishes and Eq. (115) is reduced to Eq. (90). In the higher-dimensional cases of D = 2 or D = 3, the equation of the soliton motion (115) is clearly nonlinear (even for nearly black RDS or SSS) due to the presence of the repulsive curvature-induced logarithmic potential.
Equation (115) predicts the existence of both oscillating (grey) and stationary (black) RDS or SSS: the former, perform oscillations on top of the TF cloud, changing their radii between a minimum and a maximum value, while the latter correspond to the minimum of the effective potential Eq. (116) (such stationary states do not exist in the case of a uniform ground state -as in the context of nonlinear optics [212]). As was shown in Ref. [218] (see also Ref. [158]), such a particle-like approach can describe quite effectively the above generic scenarios and the RDS dynamics up to a certain time (before the development of instabilities -see below). Furthermore, in Refs. [224,225] it was shown that the dynamics of small-amplitude RDS, as well as the collisions between them, can be described in the framework of an effective cylindrical KdV equation [114] (see also Refs. [212][213][214] for similar findings in optics). On the other hand, numerical simulations of Ref. [218] revealed that RDS are generally unstable, as they either decay to radiation (the small-amplitude ones) or are subject to the snaking instability (the moderate-and large-amplitude ones). Interestingly, as shown in Fig. 11, the snaking instability of the RDS results in the formation of vortex-anti-vortex pairs in multiples of four, which are initially set along a ring, forming a so-called vortex necklace. Eventually, this pattern relaxes to a set of four pairs located on a ring, which oscillates in the radial direction between the same limits which confined the oscillations of the original RDS; simultaneously, the pairs perform oscillatory motion along the ring [218].
Matter-wave dark solitons of radial symmetry were also analyzed by means of other approaches. For example, in Ref. [219] (see also Ch. 7 in Ref. [43]) RDS and SSS were Figure 11. (color online) Contour plots of the density of a disk-shaped condensate carrying a stationary RDS, which develops the snaking instability. The initial condition used for the integration of the GP Eq. (110) is ψ(r, 0) = µ − (1/2)Ω 2 r 2 tanh(r − r 0 (0)), with µ = 1, Ω = 0.035 and initial soliton radius r 0 (0) = ( √ 2Ω) −1 = 20.2 (for a discussion concerning the value of r 0 (0) for stationary RDS see Ref. [158]). The top left (top right) panels show, respectively, the initial condition (t = 0) and a snapshot at t = 40, while the bottom left (bottom right) right panels show, respectively, the onset of the snaking instability (t = 80) and the formation of the vortex necklace (t = 100). considered as nonlinear Bessel functions, namely solutions of the equation, resulting from Eq. (110) (with V (r) = 0) when the ansatz ψ = q(r) exp(−iµt + iSφ) is introduced (in the latter expression, S is the topological charge of a central vortex). In this setting (i.e., in the absence of the trap), the solutions of Eq. (117) include, apart from the ground state, singly-and multiply-charged vortices, as well as infinitely many RDS; the nodes of the latter correspond to the nodes of the nonlinear Bessel function governed by Eq. (117). On the other hand, if an external harmonic potential is incorporated in Eq. (117), then it is possible to find infinite branches of nonlinear bound states, with each branch stemming from the respective mode of the underlying linear problem (the radially-symmetric 2D quantum harmonic oscillator) [226]. A stability analysis performed in Refs. [219,226] also revealed that the radiallysymmetric dark solitons are typically unstable but, in agreement to the findings of Ref. [218], their life-times may be considerably long. Additionally, as shown in a recent work [227], the life-time of RDS may be extended employing the so-called [228] Feshbach Resonance Management (FRM) technique, which is based on the use of external fields to periodically modulate in time the s-wave scattering length [229][230][231][232][233]. In any case, the theoretical investigations indicate that RDS and SSS have a good chance to be experimentally observed. In fact, structures similar to stationary SSS have already been observed as transients in the Harvard experiment of Ref. [65] (see also Ref. [234] were SSS are predicted to occur as a result of collisions of vortex rings).

Stability of dark solitons in cigar-shaped condensates.
The transverse (in)stability of dark solitons confined in a purely 3D setting, namely in a cylindrical trap of the form V (z, r) = (1/2)m[ω 2 z z 2 + ω 2 ⊥ (x 2 + y 2 )], was first studied in Ref. [145]. In that work, it was shown that dynamical stability of a black soliton (stationary kink), say ψ 0 , with a nodal plane perpendicular to the axis of the cylindrical trap (see left panel of Fig. 8), requires a strong radial confinement [as in the (2 + 1)dimensional case]. Particularly, it was shown that the instability can be suppressed if the transverse (radial) condensate component is not in the TF regime, which is guaranteed as long as ω ⊥ > µ (where µ is the 3D chemical potential). This criterion can physically be understood as follows: if the condensate is confined in a highly anisotropic (cigarshaped) trap, then the energy of the lowest possible radial excitation, ω ⊥ , must exceed the kink-related kinetic energy, K 0 = −(1/2)ψ 0 ∂ 2 z ψ 0 ∼ µ, so that the latter can not be transferred to the BEC's unstable transverse modes by the inter-atomic interaction. Moreover, a systematic study in Ref. [145] revealed a criterion of dynamical stability for the black soliton in terms of the ratio ω ⊥ /ω z , namely, The critical value γ c was calculated for various values of ω ⊥ /ω z and it was found that a characteristic value of γ c , pertinent to the limiting case of ω ⊥ ≫ ω z , is γ c ≈ 2.4. In a more quantitative picture, a detailed study of the BdG equations in Ref. [145] (see also relevant work in Refs. [167,168]) revealed the emergence of complex eigenvalues in the excitation spectrum and their connection to oscillatory dynamical instabilities, including the snaking instability. Additionally, in Ref. [235] it was found that the emergence of complex eigenvalues is directly connected to bifurcations of the rectilinear black solitons to other stationary states that may exist in BECs confined in cylindrical traps. In particular, an investigation of the dependence of the excitation energy on the dimensionality parameter d [cf. Eq. (16)] led to the following results: for sufficiently low excitation energies, the black soliton may bifurcate to a solitonic vortex or to an axisymmetric vortex ring (see also Ch. 7 in [43] and references therein), with the corresponding bifurcation points occurring at a low and a higher value of d. It was also found that the emergence of the first (second) complex eigenvalue in the Bogoliubov spectrum coincides with the above mentioned bifurcation points. Therefore, according to the above results, it can be concluded that the emergence of the complex eigenvalues in the excitation spectrum (a) denote the onset of dynamical instabilities of black solitons and (b) indicate the excitation of lower energy topological states. Since these states are energetically preferable, the onset of the dynamical instability will result in decay of the "high-energy" black soliton to these "low-energy" states carrying vorticity.
On the other hand, the stability of moving (gray) solitons was analyzed in Ref. [91]. According to this work, and following the arguments of Ref. [145], a criterion for not being in the radial TF regime (which is required for dynamical stability of the gray soliton) is w R ∼ ξ, where R is the radial size of the BEC and w is the soliton width (ξ is the healing length). Recalling that the soliton width w and velocity v depend on the soliton phase angle as w ∼ 1/ cos φ and v ∼ sin φ, it is clear that as the soliton is moving towards the boundaries of the BEC, its width (velocity) is increased (decreased). Thus, the instability border R ∼ w for the gray soliton is reached for larger values of the parameter γ [see Eq. (118)] than the ones pertaining to the black soliton. In the case of strongly anisotropic traps, the critical value of the chemical potential required for the dynamical stability of the gray soliton is proportional to the soliton velocity. In other words, the stability domain of gray solitons is wider than the one of black solitons. In fact, the shallower the soliton it is, it becomes more stable, similarly to the case of homogeneous systems: see Eq. (112), which indicates that the instability band vanishes for shallow solitons with cos φ → 0.
Numerical simulations of Ref. [91] have also revealed that for γ > 10, a phaseimprinted dark soliton -with a π-phase jump -always decays in a cigar-shaped BEC (on a time scale of order of ω −1 ⊥ ), while for γ 5 it transforms into a dark soliton characterized by a flat notch region and r-independent velocity. Here, it is relevant to mention that in the recent Technion experiment [70], an interesting nonlinear excitation that evolves periodically between a dark soliton and a vortex-ring was observed in a 87 Rb BEC, for γ = 4.95 (see also relevant theoretical work in Refs. [236,237]).

5.4.
Matter-wave dark solitons in the dimensionality crossover from 3D to 1D 5.4.1. The single-soliton state. As discussed in Sec. 2.4, if the dimensionality parameter [cf. Eq. (16)] takes values d ≈ 1, then a cigar-shaped condensate is in the so-called dimensionality crossover regime from 3D to 1D. However, in such an experimentally relevant regime † , exact analytical dark soliton solutions (of arbitrary amplitudes) of the pertinent effectively 1D mean-field models (see Sec. 2.4) are not available. As a result, the analytical techniques presented in Sec. 4.2 cannot be applied for the study of matterwave dark soliton dynamics in this regime. Nevertheless, the results of Secs. 4.2 and 5.3 indicate that the evolution of dark solitons should be similar to the one pertaining to the TF-1D regime, while the solitons would not be prone to the snaking instability, as in the case of the purely higher-dimensional setups.
The statics and dynamics of matter-wave dark solitons in the crossover regime between 1D and 3D were studied in Ref. [95]. Particularly, numerical simulations of the 3D GP equation revealed that matter-wave dark solitons are indeed dynamically stable, and perform harmonic oscillations in the harmonic trap. Importantly, in the † Note that the recent Heidelberg experiments [69,71] were conducted in this regime.
case of small-amplitude oscillations, the oscillation frequency ω osc resulting from the 3D GP equation was found to be equal to the eigenfrequency ω A of the anomalous mode of the effectively-1D NPSE model [cf. Eqs. (18) and (19)]; note that the latter, can be expressed in the dimensionless form: where units are the same to the ones used for Eq. (24). The above finding leads, in turn, to the following conclusion: an equation of motion for the center z 0 of a dark soliton in a condensate in the dimensionality crossover regime can be expressed as follows: As shown in Ref. [95], the soliton oscillation frequency is a decreasing function of the dimensionality parameter d, taking values ranging from ω osc = Ω (corresponding to the non-interacting limit of d → 0), to ω osc = Ω/ √ 2 [corresponding to the TF-1D regime of d ≪ 1 -cf. Eq. (91)], with Ω being the normalized strength of the harmonic trap. In any case, the oscillation frequency is up-shifted from its TF-1D value and, as a result, the effective trapping potential felt by the dark soliton during its motion [see Eq. (90)] will effectively become steeper. In that regard, it is relevant to mention that substantial shifts of the oscillation frequency (which may be of order of 10%) have been predicted in Ref. [95] and later confirmed in the Heidelberg experiment of Ref. [69] (see discussion in Sec. 5.4.3). It should also be noticed that the soliton oscillation frequency is also a decreasing function of the soliton amplitude (or, in other words, of the oscillation amplitude), contrary to the result corresponding to the TF-1D regime † [69,71].
Results similar to the ones obtained in the framework of the NPSE model [95] can also be obtained in the framework of the generalized NLS Eq. (18) with the nonlinearity function of Eq. (20) (see, e.g., the analysis of Ref. [71]). Furthermore, we should mention that the oscillations of dark solitons were also analyzed in the framework of a GP model with generalized nonlinearities, and specific results in the physically relevant case of a cubicquintic nonlinearity (modeling two-and three-body interactions -see Sec. 2.4) were presented [162]. The same model was also studied in Ref. [238], were various stationary states, including dark solitons, were found and analyzed in detail.

5.4.2.
The multiple-soliton state. Apart from the single-dark soliton state, the case of a multiple-dark soliton state can also be analyzed in the dimensionality crossover regime using, as a guideline, the methodology exposed in Secs. 3.6, 4.2 and 4.3. In particular, in the case of well-separated and symmetrically interacting dark solitons in a harmonic trap, one may follow the analysis of Refs. [69,71] and analyze this problem by adopting a simple physical picture: each soliton in the multi-soliton state follows the evolution of the single-soliton state, i.e., it oscillates in the trap with an oscillation † Recall that the results of Sec. 4.2 indicate that the soliton oscillation frequency does not depend on the soliton amplitude in the TF-1D regime. frequency ω osc equal to the eigenfrequency ω 1 of the first anomalous mode (determined by a BdG analysis of the pertinent 1D mean-field models of Sec. 2.4), and interacts with the neighboring solitons via the effective repulsive potential of Eq. (76).
In order to further elaborate on the above, let us consider -as an examplethe two-dark soliton state of the NPSE model of Eq. (119). This state (which, in the linear limit, corresponds to the second excited state of the quantum harmonic oscillator) can be obtained as a nonlinear stationary state of the system by means of, e.g., a Newton-Raphson method. The pertinent configuration, has the form of two overlapping dark solitons, placed at z i = ±2.185 (i = 1, 2) with a fixed relative distance δz 0 = 4.37. The corresponding Bogoliubov excitation spectrum, namely the spectral plane (ω r , ω i ) of the eigenfrequencies ω ≡ ω r + iω i , is shown in Fig. 12: as it can clearly be observed, among the lowest eigenfrequencies (such as the ones at ω = 0, ω d = Ω and ω q ≈ √ 3Ω, corresponding to the Goldstone, dipole and quadrupole mode, respectively), there exist two anomalous modes with eigenfrequencies ω 1 = 0.756Ω = 0.0378 and ω 2 = 2.094Ω = 0.1047.
According to the analysis of Sec. 4.3, small displacements of the dark solitons from their equilibrium points lead to the in-phase and out-of-phase oscillatory motion of the dark soliton pair (see Fig. 13), with the respective oscillation frequencies being equal to the eigenfrequencies ω 1 and ω 2 of the two anomalous modes. Importantly, the value of the eigenfrequency of the first anomalous mode, ω 1 = 0.0378, is quite close to the oscillation frequency of a single dark soliton, ω osc = 0.0375, in the same setup (i.e., with the same parameter values), with the percentage difference being ≈ 1.6%. This generic example suggests that, generally, the dynamics of the two dark soliton state can be described by an effective Lagrangian for the two solitons, namely L eff = T − V ; here, T and V are the kinetic and potential energies, respectively, depending on the soliton centers, z i (i = 1, 2), and soliton velocities,ż i ≡ dz i /dt, as follows, where V int (z 2 − z 1 ) (with z 2 − z 1 ≡ 2z 0 ) is the repulsive potential of Eq. (76). Then, the evolution of the soliton centers can readily be determined by the Euler-Lagrange equations d(∂ż i L eff )/dt − ∂ z i L eff = 0. The latter, may be simplified upon using the approximate form of the repulsive potential [cf. Eq. (77)], thus leading to the following equations of motion: where we have assumed well-separated, almost black solitons [i.e., in Eq. (77) we have set B ≈ 1]. Apparently, the above analysis can readily be generalized for multiple solitons [71], with each one interacting with its neighbors. Importantly, if ω osc in the system of Eqs. (122)-(123) was considered to be unknown, then it would be possible to be directly obtained in the form of the characteristic frequencies of the normal modes of this system (see details in Ref. [71]); these characteristic frequencies coincide to the ones determined via the BdG analysis. This results justifies a posteriori the considered decomposition of the principal physical mechanisms (oscillations and interactions of solitons) characterizing the system. We should also note that apart from the case of small-amplitude oscillations of two well separated, almost black solitons, the more general case of the dynamics of n-interacting dark solitons (which may also perform large amplitude oscillations -see Sec. 5.4.3 below) is possible using the full set of the above mentioned Euler-Lagrange equations; the latter, lead to the following n-coupled equations of motion [71]: The oscillation frequency of the dark solitons, ω osc , as a function of their initial relative distance δz. For δz → δz 0 (corresponding to the stationary state), we obtain ω osc → ω 2 , while for δz ≫ δz 0 , we obtain ω osc → 2ω 1 .
small -then ω osc differs from (in fact, it is quite smaller than) ω 2 : considering, e.g., that δz = 8 (corresponding to initial locations of the soliton centers z = ±4 -see Fig. 13), the out-of-phase oscillation of the two solitons is characterized by a frequency ω osc = 1.74Ω < ω 2 = 2.094Ω. A qualitative explanation for this difference is the following: As discussed above, the evolution of two initially overlapping dark solitons can be effectively described by the equations of motion (122)- (123), in the presence of the repulsive potential which depends exponentially on their relative distance. If the relative distance between the two solitons is not significantly different than the one pertaining to the corresponding stationary state, i.e., δz ≈ δz 0 , the effective repulsive force is strong and, as a result, their motion is strongly affected by their coupling. On the other hand, if their initial separation becomes larger (as, e.g., in the case of the example under consideration, with δz = 8), the repulsive force becomes exponentially small and, as a result, the motion of each individual soliton is not significantly affected by the presence of the other.
In the bottom panel of Fig. 13 we show the oscillation frequency of the dark soliton -obtained by the direct integration of the NPSE model of Eq. (119) -as a function of the initial relative distance between the two solitons. It is clear that the oscillation frequency, which takes values in the interval δz ≥ δz 0 (the value corresponding to the stationary state), exhibits two different asymptotic regimes: when the initial soliton separation is small enough, δz → δz 0 (i.e., for strong coupling between the two solitons), the oscillation frequency tends to the eigenfrequency ω 2 of the largest anomalous mode; on the other hand, when the initial soliton separation is large enough, δz ≫ δz 0 (i.e., when the solitons are actually decoupled), the oscillation frequency tends to 2ω 1 ; the latter value can be explained by the fact that the period of oscillation of each individual soliton is the half of the one that would correspond to a single soliton oscillation.
As concerns relevant experiments, an oscillating and interacting dark soliton pair in a 87 Rb BEC, in the dimensionality crossover regime between 1D and 3D, was experimentally observed in the Heidelberg experiment of Ref. [69]. In this experiment, large amplitude oscillations were induced by the method of matter-wave interference (see Sec. 6.2 below). The dependence of the oscillation frequency on the distance between the two solitons (see the left panel of Fig. 13) was found and compared with experimental data: the agreement between the theoretical predictions (based on a study of the NPSE model) and the experimentally observed oscillation frequencies was excellent. In accordance to the analysis of this Section, considerable upshifts -up to 16% -of the soliton oscillation frequency from the value of Ω/ √ 2 were observed in the study of Ref. [69], and they were quantitatively attributed to the dimensionality of the system and the soliton interactions.
6. Matter-wave dark solitons in various settings and parameter regimes 6.1. Matter-wave dark solitons in multi-component condensates.
Multi-component ultracold atomic gases and BECs may be composed by two or more atomic gases, which may have the form of mixtures of: (a) two different spin states of the same atom species (so-called pseudo-spinor condensates) [239][240][241]; (b) different Zeeman sub-levels of the same hyperfine level (so-called spinor condensates) [242]; (c) different atom species (so called heteronuclear mixtures) [243]; (d) degenerate bosonfermion clouds [244]; (e) purely degenerate fermion clouds [245] (see also [41][42][43] for reviews and references therein). Such multi-component systems support various types of matter-wave soliton complexes, with the type of soliton in one species being the same or different to the type of soliton in the other species. Here, of particular interest are the so-called vector solitons with the one component being a dark soliton, which have mainly been studied in the context of two-component and spinor condensates.
6.1.1. Dark solitons in two-component condensates. Generally speaking, a mixture of N purely bosonic components can be described in the framework of mean-field theory by a system of N coupled GP equations, which can be expressed in the following dimensionless form: g n,k |ψ k | 2 ψ n − κ n,k ψ k + ∆ n,k ψ n . (125) Here, ψ n is the wave function of the n-th component (n = 1, . . . , N ), V n (r) is the trapping potential confining the n-th component, ∆ n,k is the chemical potential difference between components n and k, the nonlinearity coefficients g n,k = g k,n characterize inter-atomic collisions, while the linear coupling coefficients κ n,k = κ k,n account for spin state inter-conversion, usually induced by a spin-flipping resonant electromagnetic wave (see, e.g., Ref. [246]). In some works (see, e.g., Refs. [247][248][249]) fermionic mixtures are also described in the framework of the mean-field theory, with the self-interacting nonlinear terms being replaced by g n,n |ψ n | 4/3 ψ n . Notice that in the GP Eqs. (125) both the energy E and the total number of atoms, |ψ k | 2 dr, are conserved; furthermore, in the absence of linear inter-conversions (κ n,k = 0), the number of atoms of each component N k is separately conserved.
Let us consider the case of two bosonic species (N = 2), and assume that the system is homogeneous (V n = 0). If, additionally, there is no spin-state inter-conversion (κ n,k = 0) and chemical potential difference ((∆µ n,k ) = 0), then the binary mixture is immiscible provided that the following immiscibility condition holds [250], ∆ ≡ (g 12 g 21 − g 11 g 22 )/g 2 11 > 0, (126) where the, so-called, miscibility parameter ∆ takes in practice the values ∆ ≈ 9 × 10 −4 or ∆ ≈ 0.036 for a mixture of two spin states of a 87 Rb BEC [239,240] (see also Ref. [262]) or a 23 Na BEC [241], respectively. The condition (126) indicates that if the mutual repulsion between species is stronger than the repulsion between atoms of the same species then the two species do not mix. In such a case, the two species tend to separate by filling two different spatial regions, thus forming a "ball" and "shell" configuration (see, e.g., Ref. [240] for relevant experimental results). This way, the ground state of the system -i.e., the state minimizing the energy -may take the form of domain-wall solutions of the GP Eqs. (125) [251][252][253][254][255][256]. In accordance to the experimental observations, these solutions represent configurations of the following form: in the Thomas-Fermi limit (where kinetic energy is negligible), one species occupies the region around the trap center, and it is separated by two domain-walls from side domains occupied by the other species; on the other hand, kinetic energy favors a configuration where a single domain-wall at the trap center separates two domains occupied by the different species [255]. The dynamics of phase-separation of two-component BECs has been studied in various works both theoretically (see, e.g., Refs. [257][258][259] and also Refs. [260,261] for proposed applications) and experimentally [262][263][264]. Importantly, magnetic-field Feshbach resonances can be used to controllably change the the inter-species [263] or the intra-species [264] scattering length, and thus controllably change the (im)miscibility between the two species [264]. Apart from domain-walls, a trapped two-component quasi-1D BEC supports vector solitons, with the one component being a dark soliton; in such a case, typically, the other component may be a dark soliton [122,124,258,[265][266][267][268][269][270][271] or a bright soliton [131,258,269,[272][273][274]. Additionally, apart from such dark-dark and dark-bright solitons, dark-antidark † solitons have also been predicted to exist in BECs, either in a stationary form [267] or as a dynamical entity resulting from instabilities [258,270]. Below we will focus on the most generic vector matter-wave solitons, namely the darkdark and dark-bright ones (note that the latter have also been observed experimentally [67]), presenting results corresponding to the simplest possible setup, which solely includes the traditional time-independent harmonic trap. Notice that in the absence of the trap, vector solitons of the above mentioned types have been extensively studied in the context of nonlinear optics: there, multi-component solitons occur when fields of one frequency, or one polarization, become coupled to fields of other frequencies, or other polarizations (see, e.g., the review [34] and references therein). Mathematically speaking, the existence (and stability) of multi-component optical solitons (and also matter-wave solitons in the miscible case) can be understood by the fact that the relevant coupled NLS equations rely on the so-called Manakov system [275]: the latter, has the form of a vector NLS equation, namely, where u(z, t) = (u 1 (z, t), u 2 (z, t), · · · , u n (z, t)) is a n-component vector. This system is known to be completely integrable [276][277][278] (in fact, it can be integrated by extending the IST method that has been used to integrate the scalar NLS equation [11,12]) and admits such vector N-soliton solutions [279][280][281].
As shown in Refs. [265,266], the dynamics and interaction of dark-dark solitons in a two-component quasi-1D BEC can be studied by means of a variational approach; in the case of equal chemical potentials, the latter is based on the use of the following ansatz for the single-component soliton wave-functions, where 2z 0 (t) denotes the relative distance between the two solitons, and the ∓ signs correspond, respectively, to a kink-anti-kink state (where the solitons' phase fronts are facing each other) and a kink-kink state (where the solitons' phase fronts are in the same direction). Both the miscible and immiscible cases where studied in Refs. [265,266] and the main results of the analysis can be summarized as follows. In the miscible case (g 11 = g 22 = g 12 ), and for the kink-anti-kink state, the trajectories in the (z 0 ,ż 0 ) phaseplane are either periodic surrounding the center (0, 0) [indicating the formation of a † An anti-dark soliton is actually a dark soliton with reverse-sign amplitude, i.e., it has the form of a hump (instead of a dip) on top of the background density (see, e.g., Refs. [117,118] In such a case, the soliton is accelerated after each reflection and eventually escapes from the domain. As mentioned above, dark-bright matter-wave solitons in a quasi-1D binary BEC are also possible. Particularly, in the miscible case (with all nonlinearity coefficients g n,k being normalized to unity), the wave functions ψ d (z, t) and ψ b (z, t) of the dark and bright soliton components may be expressed in the following form [131], Here, µ d = µ and µ b = µ + ∆ are the chemical potentials of the dark and bright components, φ is the dark soliton's phase angle, z 0 denotes the solitons' center, is the inverse width of the bright soliton, and θ b = (κ tan φ)x + [κ 2 (1 − tan 2 φ)/2 − ∆]t is the bright soliton's phase. According to the analysis of Ref. [131], if the external trapping potentials V d and V b for the dark and bright solitons are slowly varying on the soliton scale κ −1 , then the dynamics of the dark-bright soliton can be described by the effective particle approach of Sec. 4.2. Particularly, assuming that the solitons are sufficiently slow, a multiple-time-scale boundary-layer theory -similar to the one used in Ref. [144] -leads to the following equation of motion for the soliton center, where (132) is reduced to Eq. (90) (recall that the latter predicts dark soliton oscillations with a frequency Ω/ √ 2), while the motion of the vector soliton becomes more sensitive to the presence of the bright component as N b is increased. For example, in the case of equal harmonic traps of strength Ω, such that V b = V d ≪ µ (i.e., for soliton motion near the trap center), the oscillation frequency of the dark-bright soliton resulting from Eq. (132) reads, It is clear that Eq. (133) shows that the oscillation frequency is down-shifted as compared to the characteristic value of Ω/ √ 2, i.e., the dark-bright pair executes slower oscillations, as the bright component is enhanced.
The predictions of Ref. [131] can directly be compared to the findings of a Hamburg experiment [67], where long-lived dark-bright matter-wave solitons were observed in a two-component quasi-1D 87 Rb BEC. In particular, using the phase-imprinting method, a dark soliton was created in one spin state of the BEC and the density dip was filled by atoms, forming the bright soliton, in another spin state of the BEC (note that the number of atoms N b of the bright soliton was ≈ 10% of the total number of atoms). The created dark-bright soliton was then observed to perform slow oscillations with a frequency 0.24Ω, which is quite smaller than the frequency of the corresponding single dark soliton in the same setting. Moreover, due to the initial state preparation, an extra dark soliton was generated, which was allowed to interact with the co-existing dark-bright soliton; it was observed that this individual dark soliton was reflected off the slower dark-bright one, with the process resembling a hard-wall reflection.
6.1.2. Dark solitons in spinor condensates. The spin degree of freedom of spinor BECs gives rise to important new phenomena (including, among others, the formation of spin domains [242], spin textures [282] and vortices [283], as well as spin oscillations [284]), which are not present in other types of BECs. Generally, a spinor condensate formed by atoms with spin F can be described in the framework of mean-field theory by a (2F + 1)-component macroscopic wave function; accordingly, a spinor F = 1 condensate is characterized by a vector order parameter, with the three components corresponding to the values of the vertical spin projection, m F = −1, 0, +1. In a quasi-1D setting, the pertinent system of coupled GP equations for the wave functions ψ ±1,0 (z, t) can be expressed in the following dimensionless form (see, e.g., Refs. [123,285,286]), Here, H ≡ −(1/2)∂ 2 x + (1/2)Ω 2 z 2 + n, with Ω being the normalized trap strength and n = |ψ −1 | 2 + |ψ 0 | 2 + |ψ +1 | 2 the total density, while δ ≡ (a 2 − a 0 )/(a 0 + 2a 2 ) where a 0 and a 2 are the s-wave scattering lengths in the symmetric channels with total spin of the colliding atoms F = 0 and F = 2, respectively. Actually, the parameter δ represents the ratio of the strengths of the spin-dependent and spin-independent interatomic interactions, and may take negative or positive values for ferromagnetic or anti-ferromagnetic (alias polar) spinor BECs, respectively. Typically, in the relevant cases of 87 Rb and 23 Na atoms with F = 1, δ = −4.66 × 10 −3 [287] and δ = +3.14 × 10 −2 [288]; nevertheless, the above values may in principle be modified by employing the so-called confinement-induced Feshbach resonance [289].
In the limiting case of δ = 0 (and in the absence of the potential), the system of Eqs. (134)-(135) is reduced to the completely integrable Manakov system. On the other hand, as shown in Ref. [290], another completely integrable version of Eqs. (134)- (135) corresponds to the case δ = 1 (i.e., for interatomic and anti-ferromagnetic interactions of equal magnitude): in this case, the resulting matrix NLS equation with non-vanishing boundary conditions is completely integrable by means of the IST method [291] and admits exact analytical vector N-dark soliton solutions (i.e., single-and multiple-vector dark solitons of the dark-dark-dark type in terms of the m F = −1, 0, +1 spinor components) [292] (see also Ref. [293]). The one-dark soliton state of this system can be classified as (a) ferromagnetic (i.e., with nonzero total spin), which has domainwall shaped wave functions, and (b) polar (i.e., with zero total spin), characterized by the familiar hole soliton profile. Note that the collisions of two solitons give rise to interesting spin-dependent phenomena, such as spin-mixing or spin-transfer [292].
In the physically relevant case of small δ, mixed dark-bright solitons of the darkdark-bright or bright-bright-dark type (again in terms of the m F = −1, 0, +1 spinor components) were also predicted to occur in anti-feromagnetic spinor F = 1 BECs [123]. In the small-ampliude limit (and in the absence of the trap), these solitons were found to obey the completely integrable Yajima-Oikawa system [294], by means of which it was found that the functional form of the dark and bright components is similar to the one in Eqs. (130)- (131). Numerical simulations in Ref. [123] demonstrated that, for small-amplitudes, such dark-bright solitons feature genuine soliton behavior (i.e., they propagate undistorted and undergo quasi-elastic collisions), while for moderate and large amplitudes (and also for large values of δ) they can exist as long-lived objects as well. Furthermore, for sufficiently small number of atoms of the bright soliton, the bright component(s) are guided by the dark one(s), and the vector soliton performs harmonic oscillations; the oscillation frequency is different for small-and moderate-amplitude solitons, and it is respectively given by: In the above expression, Ω osc is given by Eq. (133), while the constants α 0,1 and ǫ 0,1 (with ǫ 0,1 ≪ α 0,1 ) depend on the normalized number of atoms of the bright component. It is clear that the characteristic oscillation frequency of the dark soliton (Ω/ √ 2), as well as the result of Ref. [131], is modified by the spin-dependent interactions that are present in the case of a spinor F = 1 BEC.

Matter-wave interference and dark solitons.
Matter-wave interference experiments (see, e.g., the seminal work of Ref. [295]) are known to demonstrate, apart from self-interference, the interference between two BECs confined in a trap, divided into separate parts by means of a barrier potential induced by a laser beam. In particular, the BECs are left to expand and overlap forming interference fringes, similar to the ones known in optics. Much interest has been drawn to a better understanding of this fundamental phenomenon, especially as concerns the coherence properties of the interfering BECs. In that regard, it is worth mentioning that the (incorrect) assumption that "when the interfering BECs have fixed atom numbers, there can be no phase", was resolved -shortly after the experimental realization of BECs [36][37][38] -in Ref. [296]. On the other hand, since most of the relevant experimental findings can be quantitatively reproduced in the framework of the GP mean-field theory [297], we will proceed by adopting this approach in order to discuss the connection between dark solitons and matter-wave interference.
An interesting variation of the interference process, which is naturally attributed to the inherent nonlinearity of BECs due to interatomic interactions, is thatunder certain conditions -the collision of two initially separated condensates can lead to the creation of dark solitons. This "nonlinear interference" effect was first observed in simulations [139], and was subsequently analyzed theoretically [298]. Other studies, basically relying on the self-interference of BECs, have also been proposed as well [220][221][222][223]299]. Importantly, relevant recent experiments employing this, socalled, matter-wave interference method have already been reported, demonstrating the generation of vortices [300] (see also theoretical work in Refs. [223,301,302]) and dark solitons [69][70][71][72] (see also the experiment of Ref. [62]).
To get a deeper insight into the physics of the matter-wave interference process, let us follow the arguments of Ref. [298] and consider the interference between two separated quasi-1D BECs colliding in the presence of a harmonic trap. There exist two different regimes characterizing this process, namely a linear and a nonlinear one, depending on the competition between the kinetic and the interaction energies. In the linear regime, the kinetic energy of the condensates exceeds the nonlinear interaction energy of the atoms. In this case, and at any time t, the total wave function of the system can be well approximated by a linear superposition of the wave functions that each individual condensate would have at t. The two initially well-separated BECs interfere at the trap center, produce a linear interference pattern, and then separate again regaining their initial shape. The fringe spacing l of the interference pattern is determined by the k vector that each individual condensate would have if performing a dipole oscillation alone in the trap and, at the time of maximum overlap, l = π/D( /2mω z ) (here D is the initial distance between the condensates). It is clear that the higher the kinetic energy is, the higher the number of fringes and the smaller the fringe spacing is. Approximating the individual wave functions in the TF limit, it can be found that the kinetic energy (estimated from the curvature of a cos 2 interference pattern) exceeds the peak nonlinear energy (at the center of the fringes) when the initial distance between the two BECs exceeds critical distance, namely D > D c , where D c is given by [298]: Here, N is the number of atoms, a the s-wave scattering length, ω z the longitudinal trap frequency, and m the atomic mass. If the above condition is not fulfilled, namely D < D c , then the system enters in the nonlinear regime. In the latter, the interference pattern consists of stable fringes with a phase jump of order of π across them, which can naturally be identified as genuine dark solitons. Notice that in the nonlinear regime, the initially individual condensates instead of reforming as separate objects, they form a combined condensate undergoing a quadrupole oscillation.
The results of Heidelberg experiments [69,71] can be compared directly to the above theoretical predictions. In these experiments, dark solitons were created by releasing a 87 Rb BEC from a double-well trap into a harmonic trap in the dimensionality crossover regime from 1D to 3D. For the parameters used, the initial distance of the individual condensates was approximately five times smaller than the critical distance and the trap frequencies were ramped, with ramping times chosen so as to minimize the excitation of the quadrupole mode. It was shown that, in accordance to the observations of Ref. [62], the number of created solitons is even for a zero phase-difference between the two initially separated condensates, while it is odd for a phase-difference close to π. If the phasedifference is exactly equal to π, a standing (black) dark soliton in the middle of the trap is always created. Notice that the total number of the created solitons depends on the momentum of the merging condensates, which may be controlled by varying, e.g., the distance between the condensate fragments, the number of atoms, or the aspect ratio of the trap [69,71,72].

BEC superfluidity and dark solitons
A flow past an obstacle is known to be one of the most fundamental contexts for studying superfluidity. Particularly, according to the Landau criterion for superfluidity [303], a superfluid flow past an obstacle, is stable (unstable) for group velocities smaller (larger) than the speed of sound. Actually, breakdown of superfluidity is caused by the opening of channels for emission of excitations in the fluid, whose formation manifests itself as an effective dissipation. In the BEC context, early experiments from the MIT group [187,304] demonstrated the onset of dissipation induced by the motion of an obstacle (in the form of a strongly repulsive dipole beam). From a theoretical standpoint, the problem can be studied by using a NLS (or a GP) equation that includes a localized external potential of the form V (r − vt) (with V (r) → 0 as |r| → ∞) accounting for the presence of the obstacle moving with velocity v; this potential may be naturally superimposed to the usual trapping potential confining the condensate. In relevant earlier studies [305], where the NLS equation as a model of superflow was used, vortex formation induced by the superfluid flow around an obstacle was predicted.
The lower-dimensional setting, namely the 1D flow of a repulsive NLS fluid in the presence of an obstacle, was also studied [306,307] (see also Ref. [308]). Specifically, in Ref. [306] it was shown that below an obstacle-dependent critical velocity, there exists a steady dissipationless flow solution, which disappears at the critical velocity by merging with an unstable solution in a saddle-node bifurcation. This unstable solution represents the transition state for emission of dark solitons, which are repeatedly generated above the critical velocity. In fact, the onset of dissipation corresponds to nonstationary flow with a wake asymptotically extending upstream to infinity, and downstream periodic emission of dark solitons [307]. Note that in both Refs. [306] and [307] the critical velocity was found to be smaller than the speed of sound, a result that may be explained by the fact that, in the region of the potential, the local fluid velocity can reach values higher than the local sound velocity (critical velocity values smaller than the speed of sound were also observed in Refs. [187,304]).
The above results paved the way for a better understanding of the BEC flow past an obstacle and inspired further investigations [309][310][311]. Importantly, in an recent experiment [66], the BEC flow induced by a broad, penetrable barrier (in the form of a laser beam) swept through an elongated 87 Rb condensate was systematically studied: it was demonstrated that at slow barrier speeds the flow is stable, at intermediate speeds becomes unstable and dark soliton generation is observed, while at faster speeds, remarkably, soliton formation completely ceases. Both repulsive and attractive barriers were used in the experiment of Ref. [66], and were found to lead to dark soliton formation; additionally, it was also found that the critical velocity for the breakdown of the BEC superfluidity and soliton generation was smaller than the speed of sound. Note that in a recent work [312], velocity regimes similar to the ones found in Ref. [66] were analytically predicted by using a hydrodynamic approach.
As shown theoretically [306,307,[309][310][311] and demonstrated experimentally [66], dark solitons (and vortices) are formed if the size of the "hypersonic" obstacle is of the order of, or greater than, the characteristic healing length of the condensate. On the other hand, if the size of the obstacle is much smaller than the healing length, the main loss channel, which opens at supersonic velocities of the obstacle, corresponds to the Cerenkov emission of Bogoliubov's excitations [313]. Notice that in the case of large hypersonic obstacles, two dispersive shock waves, which start propagating from the front and the rear parts of the obstacle, are formed. Far from the obstacle, the shock front gradually transforms into a linear "ship wave" located outside the Mach cone [314][315][316][317], whereas the rear zone of the shock is converted into a "fan" of oblique dark solitons located inside the Mach cone [318][319][320] (see also relevant experimental results in Refs. [314,321]). An important result reported in Ref. [322] is that although such dark solitons are unstable in higher-dimensional settings with respect to transverse perturbations (see Sec. 5.1), the instability becomes convective -rather than being absolute -for sufficiently large flow velocities and, thus, dark solitons are effectively stable in the region around the obstacle.
The flow of a multi-component BEC past an obstacle was also studied, and the cases of a two-component [270,271,323] and a spinor F = 1 condensate [324] were analyzed. It is interesting to note that, as shown in Ref. [270] in the case of a two-component BEC, the existence of two different speeds of sound provides the possibility for three dynamical regimes: when both components are subcritical, nucleation of coherent structures does not occur; when both components are supercritical they both form dark solitons in 1D and vortices or rotating vortex dipoles in 2D; in the intermediate regime, the nucleation of a dark-anti-dark soliton in 1D or a vortex-lump configuration in 2D is observed. Furthermore, as shown in Ref. [271], dark solitons can be convectively stabilized in the 2D setting at sufficiently high values of the obstacle velocity, similarly to the case of one-component BECs [322].

Matter-wave dark solitons in optical lattices.
Bose-Einstein condensates loaded into periodic optical potentials, so-called optical lattices (OLs), have attracted much attention as they demonstrate rich physical properties and nonlinear dynamics (see, e.g., Refs. [43,83,[325][326][327][328] for reviews). Optical lattices are generated by a pair of laser beams forming a standing wave which induces a periodic potential; thus, for a BEC confined in an optical lattice, the trapping potential in the GP model can be regarded as a superposition of a harmonic (magnetic or optical) trap and a periodic potential. Particularly, in a quasi-1D setting (generalization to higher-dimensional cases is straightforward) -cf. Eq. (78) -the trap takes the following dimensionless form (see, e.g., Ref. [44]): Here, Ω and V 0 denote, respectively, the harmonic trap and OL strengths, L ≡ π/k = (λ/2) sin(φ/2) is the periodicity of the lattice, with λ being the common wavelength of the two interfering laser beams, and φ the angle between them. In some cases (as, e.g., in the experiments of Refs. [329,330]), the harmonic potential is very weak as compared to the optical lattice and, thus, it can be ignored. Then, the stationary states of the pertinent GP equation -including solely the OL potential -can be found in the form of infinitely extended waves, with the periodicity of the OL, known as nonlinear Bloch waves (see, e.g., Ch. 6 in Ref. [43] and references therein). In the same case (i.e., in the absence of the harmonic potential), if the OL is very deep (compared to the chemical potential), the strongly spatially localized wave functions at the lattice sites are approximated by Wannier functions (see, e.g., Ref. [331]) and the tight-binding approximation can be applied; then, the continuous GP equation is reduced to the discrete NLS (DNLS) equation (see, e.g., Refs. [44,83,325] and Ch. 13 in Ref. [43], as well as Refs. [332,333] for reviews for the DNLS model). Dark solitons, which may naturally exist in all of the above settings and combinations thereof, have been studied both in combined harmonic and OL potentials [183,184,334,335] and in optical lattices (in the absence of the harmonic trap). In the latter case, various studies have been performed in the frameworks of the continuous GP equation, as well as its tight-binding approximation counterpart [336][337][338][339][340]. Notice that matter-wave dark solitons have also been studied in double-periodic optical superlattices † [340,341], while there exists a vast amount of work concerning dark solitons in periodic media arising in various contexts, such as nonlinear optics [343][344][345][346][347][348], solid-state physics [349] and the theory of nonlinear waves [350].
6.4.1. Dark solitons in combined harmonic and OL potentials. The stability of matterwave dark solitons in the combined harmonic and OL potential, was first studied in Ref. [334] by a means of a BdG analysis that was performed in the framework of both † Such a potential has the form V (z) = V 1 cos(k 1 z)+V 2 cos(k 2 z), where k 1 and k 2 > k 1 are the primary and secondary lattice wavenumbers, and V 1 and V 2 are the associated sublattice amplitudes [342].
the continuous quasi-1D GP equation and its DNLS counterpart. It was found that in the discrete model stationary dark solitons located at the minimum of the harmonic trap are, generally, subject to a weak oscillatory instability, which manifests itself as a shift of the soliton from its initial location, accompanied by quasi-periodic oscillations.
On the other hand, in the continuous GP model, dark solitons may be stable, with the (in)stability determined by the period and amplitude of the OL. In any case, the dark solitons are robust and if the oscillatory instability is present, it sets in at large times. The dynamics of dark solitons in the combined harmonic and OL potential can be studied upon distinguishing physically relevant cases, depending on the competition of the characteristic spatial scales of the problem [184]. Particularly, assuming that the harmonic trap varies slowly on the soliton scale, i.e., w = 1/ cos φ ∼ ξ ≪ Ω −1 (where w is the soliton width for chemical potential µ = 1, φ is the soliton phase angle, and ξ the healing length), the following three cases can readily be identified: (a) the case of a long-period OL, with L ≫ ξ, (b) the case of a short-period OL, with L ≪ ξ, and (c) the intermediate case, with L ∼ ξ. Then, if the OL strength is sufficiently small, the soliton dynamics in cases (a) and (b) can be treated in the framework of the adiabatic approximation (see Sec 4.2). Particularly, as shown in Ref. [184], case (a) can be studied by means of the Hamiltonian approach of the perturbation theory, and case (b) by means of a multi-scale expansion method (treating k −1 as a small parameter); this way, it can be shown that in both cases the dark soliton behaves as an effective classical particle, performing harmonic oscillations in the presence of the trap of Eq. (138). The oscillation frequency, which is different from its characteristic value Ω/ √ 2, is modified by the presence of the lattice according to the equations: for cases (a) and (b), respectively. As concerns the more interesting case (c) (see Refs. [183,184]), it can be shown that if the dark soliton is initially placed quite close to the bottom of a well of the OL potential, it remains there for a rather long time; eventually, however, it escapes due to the radiation-loss mechanism, and then performs large-amplitude oscillations in the condensate. Furthermore, if the harmonic trap is weak enough, the soliton eventually decays. In fact, as discussed in Ref. [183], the OL causes a dynamical instability (because the dark soliton has to "traverse" the potential humps caused by the lattice) resulting in a faster decay of the soliton than if it was evolving in the presence of the harmonic trap only: the presence of the lattice dephases the sound waves emitted by the soliton, hence reducing the effectiveness of the soliton to get stabilized by reabsorbing the sound waves (see discussion in Sec. 4.4). Nevertheless, according to the observations of Ref. [184] that the soliton may remain stationary for a relatively long time, in Ref. [335] (see also Ref. [184]) it was proposed that a timedependent OLs may either (i) capture a moving dark soliton or (ii) capture and drag a stationary soliton, bringing it to a pre-selected final destination. Notice that the transfer mechanism is robust as long as adiabaticity of the process is ensured (i.e., for sufficiently small speeds of the moving OL).
6.4.2. Dark solitons in optical lattices and superlattices. As mentioned above, dark solitons in OLs and superlattices have also been studied, in the absence of the harmonic trap, in the frameworks of the continuous and discrete NLS equations. Particularly, in Ref. [336], a DNLS model was derived in the tight-binding approximation -i.e., for a single-isolated band in the Floquet-Bloch spectrum -which was used to study matterwave dark solitons. Later, in Ref. [337], a continuous coupled-mode model was used to study the existence and stability of, so-called, dark lattice solitons, while a more general analysis was presented in Ref. [338]; in that work, a continuous GP model with periodic potential was shown to support stable stationary dark solitons for both attractive and repulsive interatomic interactions, which were found numerically [338] (see also relevant results in [339]).
In a more recent work [340], where both regular optical lattices and superlattices were considered, it was shown that each type of nonlinear Bloch wave can serve as a stable background supporting dark solitons. This way, different families of dark solitons, originating within the bands of the Floquet-Bloch spectrum, were found and their dynamical properties were analyzed. In particular, considering the continuous analogue of the Peierls-Nabarro potential (see, e.g., Ref. [351]) in discrete lattices, it was shown that the mobility and interaction properties of the dark solitons can be effectively controlled by changing the structure of the optical superlattice; moreover, following the ideas of Ref. [335], time-dependent superlattices were also shown to control the static and dynamical properties of matter-wave solitons [341].
Here we should point out that all the above mentioned studies on the dynamics of matter-wave dark solitons in optical lattices were carried out in the framework of the GP mean-field theory. Nevertheless, it is worth emphasizing that the GP equation is inadequate for dealing with several important aspects of ultracold bosons in optical lattices, such as the superfluid-to-Mott insulator phase transition (see, e.g., Refs. [352,353]) or, more generally, strong correlation effects (see, e.g., the review [328]). Thus, more recently, studies on the quantum dynamics of dark solitons have started to appear. In that regard, it is relevant to mention that matter-wave dark solitons were studied in the context of the Bose-Hubbard model [354,355], and it was found that dark soliton collisions become inelastic, in strong contrast to the predictions of mean-field theory. A conclusion of the above works [354,355] is that the lifetime and collision properties of matter-wave dark solitons in optical lattices may provide a clear signature of quantum effects. Additionally, in another recent work [356], dark solitons were studied in the vicinity of the superfluid-to-Mott insulator transition; particularly, in this work [356], antisymmetric eigenstates corresponding to standing solitons, as well as propagating solitons created by phase-imprinting, were presented and the soliton characteristics were found to depend on quantum fluctuations.
From the viewpoint of experiments, trains of stationary dark solitons were observed in a 87 Rb condensate confined in a 3D harmonic trap and a 1D OL [64]. The underlying mechanism for the formation of such structures were multiple Bragg reflections caused by displacing the harmonic trap and, thus, setting the BEC into motion. Due to the dimensionality of the system, the solitons were found to be subject to the snaking instability, giving rise to the subsequent formation of vortex rings (see Sec. 5.1), similarly to the observations of the pertinent JILA experiment (but without the OL) [48].
So far, we have considered the stability and dynamics of matter-wave dark solitons at zero temperature, T = 0. Nevertheless, as experiments are obviously performed at finite temperatures, it is relevant to consider the dissipative instability of dark solitons induced by the thermal excitations that naturally occur. This problem was first addressed in Ref. [112], where a kinetic-equation approach, together with a study of the Bogoliubovde Gennes (BdG) equations, was used. In this work, it was found that the dark soliton center obeys an equation of motion which includes an anti-damping term [similarly to Eq. (105)], which is nonzero (zero) for finite (zero)-temperature. The behavior of the solutions of this equation of motion incorporating the anti-damping term can be used to explain -at least qualitatively -the soliton dynamics observed in experiments: solitons either decay fast (for high temperatures) [45,46,49] or perform oscillations (for low temperatures) [67][68][69]71] of growing amplitude and eventually decay, so that the system finally relaxes to its ground state (see also discussion below for the role of the anti-damping term).
Dark soliton dynamics in BECs at finite temperatures was also studied in other works by means of different approaches. In particular, in Ref. [185], the problem was treated in the framework of a mean-field model, namely the so-called dissipative GP equation (see below); this equation incorporates a damping term, first introduced phenomenologically [357] and later justified from a microscopic perspective (see, e.g., the review [358]). On the other hand, in Refs. [359,360] the same problem was studied numerically, using coupled Gross-Pitaevskii and quantum Boltzmann equations, which include the mean field coupling and particle exchange between the condensate and the thermal cloud. Furthermore, in Refs. [58,361], finite-temperature dynamics of dark solitons was studied by means of the so-called stochastic GP equation (see, e.g., Ref. [358]), while in Ref. [362] quantum effects on dark solitons were additionally studied in the framework of the truncated Wigner approximation (see, e.g., Refs. [363,364] for this approach); we also note that in the recent work [365], the dissipative dynamics of a dark soliton at temperatures T , lower than the chemical potential µ of the background Bose liquid, was studied. In the work [361], it was shown that for sufficiently low temperatures and certain parameter regimes, averaged dark soliton trajectories obtained by the stochastic GP equation, are in a very good agreement with results obtained by the dissipative GP model. Thus, the results of Ref. [361] indicate that the use of the dissipative GP equation in studies of dark solitons in finite-temperature BECs (a) can reasonably be justified from a microscopic perspective, and (b) allows for an analytical description of the problem, by employing techniques exposed in Sec. 4.2, provided that the dynamics of the thermal cloud does not play the dominant role. bifurcation diagrams being almost identical. A similar situation occurs in the case of multiple-dark solitons as well: as shown in Ref. [366], eigenvalues derived by coupled equations of motions [similar to the ones in Eqs. (122)-(123)] for two-or three-dark solitons, were again found to be almost identical to the ones of the anomalous modes of the system.

Conclusions and perspectives
We have presented the recent progress on the study of dark solitons in atomic BECs, including analytical, numerical and experimental results. In fact, although the main body of this work was basically devoted to the theoretical aspects of this topic, we have tried to connect the theoretical results to pertinent experimental observations. In that regard, we have particularly tried to highlight the close connection between theory and experiments and the reasonable agreement between the two.
Matter-wave dark solitons were predicted to occur in BECs as early as 1971 [9], but were observed in experiments only 28 years later, in 1999 [45]. Although, till then, dark solitons had already a relatively long history in the context of nonlinear optics (where they were first observed in experiments on 1987 [1] and studied extensively in theory during the next years [34]), one can readily realize an emerging interest in them: during the last decade, there have been more than ten experiments on dark solitons [45][46][47][48][49][64][65][66][67][68][69][70][71], and half of them have been conducted very recently [66][67][68][69][70][71], with an unprecedented control over both the condensate and the solitons. As the experimental developments continuously inspire -and, at the same time, are guided by -a huge number of relevant theoretical works, one may expect that the interest in matter-wave dark solitons will still be growing in the near future.
Although there has been a tremendous progress on our understanding of matterwave dark solitons in atomic BECs over the last years, many important issues remain to be addressed or studied in a more systematic way. A relevant list is appended below.
• Beyond Mean-Field. Matter-wave dark solitons, being fundamental nonlinear macroscopic excitations of BECs, play an important role at probing the properties of the condensates at the mesoscale (see, e.g., discussion in Ref. [59]). In that regard, a quite interesting research direction is the study of these nonlinear structures, both in theory and in experiments, in various settings and regimes where thermal and quantum effects are important. In fact, mean-field theory can only account for averaged results (e.g., soliton decay times [361]), whereas recent experiments [67,69,71] indicate shot-to-shot variations that could be accounted for by stochastic approaches. There exist various experimentally relevant settings -such as the ones where the number of atoms is small, particularly those in optical lattices or at very low temperatures -which enhance the importance of quantum fluctuations; therefore, the latter should be appropriately included. One interesting question concerns, for example, the issue of the filling of the dark soliton due to averaging based on thermal or quantum fluctuations [169,367,368] and its relation to the optical counterparts [372].
• Ultracold Fermi gases. We finally note that, so far, matter-wave dark solitons have mainly been studied in the context of ultra-cold Bose gases. Nevertheless, recent progress in the area of ultra-cold Fermi gases (see, e.g., Refs. [373,374] for recent reviews), suggest that (similarly to vortices) dark solitons may be relevant in this context as well. In fact, pertinent theoretical studies have already started to appear [375][376][377], but there is still much work to be done towards this direction, both in theory and in experiments.