Dark solitons in atomic BoseEinstein condensates: from theory to experiments
Abstract
This review paper presents an overview of the theoretical and experimental progress on the study of matterwave dark solitons in atomic BoseEinstein condensates. Upon introducing the general framework, we discuss the statics and dynamics of single and multiple matterwave dark solitons in the quasi onedimensional setting, in higherdimensional 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.
type:
Topical Reviewpacs:
03.75.Lm, 03.75.Kk, 05.45.YvContents:
 1 Introduction
 2 Meanfield description of BoseEinstein condensates

3 General background for the study of matterwave dark solitons
 3.1 NLS equation and dark soliton solutions.
 3.2 Dark solitons and the Inverse Scattering Transform.
 3.3 Integrals of motion and basic properties of dark solitons.
 3.4 Smallamplitude approximation: shallow dark solitons as KdV solitons.
 3.5 On the generation of matterwave dark solitons
 3.6 Multiple dark solitons and dark soliton interactions
 4 Matterwave dark solitons in quasi1D Bose gases
 5 Matterwave dark solitons in higherdimensional settings
 6 Matterwave dark solitons in various settings and parameter regimes
 7 Conclusions and perspectives
1 Introduction
A dark soliton is an envelope soliton having the form of a density dip with a phase jump across its density minimum. This localized nonlinear wave exists on top of a stable continuous wave (or extended finitewidth) 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 nonpropagating kink in a parametricallydriven shallow liquid [5], dark soliton standing waves in a discrete mechanical system [6], highfrequency 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 BoseEinstein condensates (BECs). In particular, in Ref. [9], exact soliton solutions of the GrossPitaevskii (GP) equation (which is a variant of the NLS model) [10] were found and connected, in the smallamplitude 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 socalled matterwave dark solitons, were among the first purely nonlinear states that were experimentally observed in BECs [45, 46, 47, 48, 49].
The interest on matterwave 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 matterwave dark solitons (due to the phase jump at their density minimum) renders them a “degenerate”, onedimensional (1D) analogue of vortices, which are of paramount importance in diverse branches of physics [53]. Additionally, and perhaps more importantly, matterwave dark solitons are — similarly to vortices [54, 55, 56] — quite fundamental structures arising spontaneously upon crossing the BEC phasetransition [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 matterwave interferometer in the nonlinear regime [60, 61] (see also relevant experiments of Refs. [62, 63] devoted to atomchip interferometry of BECs).
The early matterwave 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 matterwave 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 matterwave 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 matterwave dark solitons), provides some restrictions in the article: first, the space limitations of the article, will not allow for an allinclusive 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 meanfield theory: as will be shown, the latter can be used as a basis of understanding of most effects and experimental findings related to matterwave 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 meanfield description of BECs. Particularly, we first present the GP equation and discuss its connection with the respective full quantum manybody problem. Next, we present the ground state of the condensate and discuss how its smallamplitude excitations can be studied by means of the Bogoliubovde Gennes (BdG) equations. Lowerdimensional 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 threedimensional (3D) BECs and introduce elongated (alias “cigarshaped”) BECs, quasi onedimensional (1D) BECs and quasi twodimensional (2D) (alias “diskshaped”) ones, as well as discuss cases relevant to the dimensionality crossover regimes. The topics of stronglyinteracting Bose gases, and their relevant meanfield description, are also briefly covered.
Section 3 provides the theoretical basis for the study of matterwave 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 smallamplitude approximation — leading to the connection of matterwave dark solitons to Kortewegde Vries (KdV) solitons — are discussed. Furthermore, the generation of matterwave dark solitons by means of the phase, density and quantumstateengineering 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 matterwave dark solitons in quasi1D 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 multipledark solitons and the role of the pertinent anomalous modes. Nonadiabatic effects, namely emission of radiation of solitons in the form of sound waves, as well as rigorous results concerning the persistence and stability of matterwave dark solitons, are also discussed.
Section 5 studies matterwave dark solitons in higherdimensional 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 matterwave 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 cigarshaped (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 matterwave dark solitons. In particular, we first present results concerning matterwave dark solitons in multicomponent (pseudospinor and spinor) BECs. Next, we discuss how matterwave interference and the breakdown of BEC superfluidity are connected to the generation of matterwave dark solitons. We continue by referring to matterwave dark solitons in periodic potentials, namely optical lattices (OLs) and superlattices, and conclude this Section by discussing the statics and dynamics of dark solitons at finite temperatures.
Finally, in Section 7 we briefly summarize our conclusions and discuss future challenges.
2 Meanfield description of BoseEinstein condensates
BoseEinstein condensation of dilute atomic gases is an unambiguous manifestation of a macroscopic quantum state in a manybody 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 manybody 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 GrossPitaevskii (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, lowestorder meanfield theory, as this can be used as a basis to understand the nonlinear dynamics of matterwave dark solitons. Interesting effects naturally arise beyond the GP meanfield, both due to thermal and quantum fluctuations. Such effects become particularly relevant in extremely tightly confining geometries, or when the BoseEinstein condensation transition region is approached.
2.1 The GrossPitaevskii 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 manybody aspects of BECs). Particularly, a sufficiently dilute ultracold atomic gas, composed by interacting bosons of mass confined by an external potential , can be described by the manybody Hamiltonian; the latter can be expressed, in second quantization form, through the boson annihilation and creation field operators, and (which create and annihilate a particle at the position ) namely,
(1)  
where is the singleparticle operator and is the twobody interatomic potential. Apparently, the underlying full manybody problem is very difficult to be treated analytically (or even numerically) as increases and, thus, for convenience, a meanfield approach can be adopted. The meanfield approach is based on the separation of the condensate contribution from the boson field operator as follows [79]:
(2) 
In the above expression, the expectation value of the field operator, , is known as the macroscopic wave function of the condensate, while describes the noncondensate part, which accounts for quantum and thermal fluctuations. Considering the case of a dilute ultracold gas with binary collisions at low energy, characterized by the wave scattering length , the interatomic potential can be replaced by an effective deltafunction interaction potential, [41, 42], with the coupling constant given by . Under these assumptions, a nontrivial zerothorder theory for the BEC wave function can be obtained by means of the Heisenberg evolution equation , upon replacing the field operator with the classical field , i.e., ignoring the quantum and thermal fluctuations described by . Such a consideration leads to the GrossPitaevskii (GP) equation [10], which has the form:
(3) 
In the above equation, is normalized to the number of atoms , namely,
(4) 
and the nonlinearity (which is obviously introduced by interatomic interactions) is characterized by the wave scattering length , which is or for repulsive or attractive interatomic interactions, respectively. Notice that Eq. (3) can be written in canonical form, (with denoting complex conjugate), where the dynamically conserved energy functional is given by
(5) 
with the three terms in the righthand side representing, respectively, the kinetic energy, the potential energy and the interaction energy.
A timeindependent version of the GP equation can be obtained upon expressing the BEC wave function as , where is the chemical potential. This way, Eq. (3) yields the following equation for the stationary state :
(6) 
2.2 The meanfield approach vs. the manybody 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 selfconsistent
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 leadingorder in the small parameter
The starting point of the analysis of Ref. [80] is the effective Hamiltonian of identical bosons, which can be expressed (in units so that ) as follows:
(7) 
where is a general interaction potential assumed to be spherically symmetric and decaying faster than at infinity. Then, assuming that the quantummechanical groundstate energy of the Hamiltonian (7) is (here is the number of particles and is the dimensionless twobody scattering length), the main theorem proved in Ref. [80] is the following. The GP energy is the dilute limit of the quantummechanical energy:
(8) 
where is the energy of a solution of the stationary GP Eq. (6) (in units such that ), and the convergence is uniform on bounded intervals of .
The above results (as well as the ones in Ref. [81]) were proved for stationary solutions of the GP equation, and, in particular, for the ground state solution. More recently, the timedependent GP Eq. (3) was also analyzed within a similar asymptotic analysis in Ref. [82]. In this work, it was proved that the limit points of the particle density matrices of (which is the solution of the particle Schrödinger equation) satisfy asymptotically the GP equation (and the associated hierarchy of equations) with a coupling constant given by , where describes the interaction potential.
These rigorous results, as well as a large number of experimental results related to the physics of BECs, indicate that (under certain conditions) the GP equation is a good starting point for analyzing the statics and dynamics of BECs.
2.3 Ground state and excitations of the condensate.
Let us now consider a condensate confined in a harmonic external potential, namely,
(9) 
where , , and are the (generally different) trap frequencies along the three directions. In this setting, and in the case of repulsive interatomic interactions () and sufficiently large number of atoms , Eq. (6) can be used to determine analytically the ground state of the system. In particular, in the asymptotic limit of (where is the harmonic oscillator length associated with the geometrical average of the trap frequencies), it is expected that the atoms are pushed towards the rims of the condensate, resulting in slow spatial variations of the density profile . Thus, the latter can be obtained as an algebraic solution stemming from Eq. (6) when neglecting the kinetic energy term — the socalled ThomasFermi (TF) limit [41, 42, 43]:
(10) 
in the region where , and outside, and the value of being determined by the normalization condition [cf. Eq. (4)]. Notice that the TF approximation becomes increasingly accurate for large values of .
Smallamplitude 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.,
(11) 
where , are the components of the linear response of the BEC to the external perturbations that oscillate at frequencies [the latter are (generally complex) eigenfrequencies]. Substituting Eq. (11) into Eq. (6), and keeping only the linear terms in and , we obtain the socalled Bogoliubovde Gennes (BdG) equations:
(12) 
where is the singleparticle Hamiltonian. Importantly, these equations can also be used, apart from the ground state, for any other stationary state (including, e.g., solitons) with the function being modified accordingly. In such a general context, the BdG equations provide the eigenfrequencies and the amplitudes and of the normal modes of the system. Note that due to the Hamiltonian nature of the system, if is an eigenfrequency of the Bogoliubov spectrum, so are , and . In the case of stable configurations with , the solution of BdG equations with frequency represent the same physical oscillation with the solution with frequency [42].
In the case of a homogeneous gas () characterized by a constant density , the amplitudes and in the BdG equations are plane waves, , of wave vector . Then, Eqs. (12) lead to the dispersion relation,
(13) 
In the case of repulsive interatomic interactions (), Eq. (13) indicates that smallamplitude harmonic excitations of the stationary state
(14) 
with , are always stable since for every . Thus, this state is not subject to the modulational instability (see, e.g., Ref. [83] and references therein). This fact is important, as the wave function of Eq. (14) can serve as a stable background (alias “pedestal”), on top of which strongly nonlinear localized excitations may be formed; such excitations may be, e.g., matterwave dark solitons which are of particular interest in this work. Notice that the above mentioned smallamplitude harmonic excitations are in fact sound waves, characterized by the phonon dispersion relation [see Eq. (13) for small momenta ], where
(15) 
is the speed of sound. We should note in passing that in the case of attractive interatomic interactions (), the speed of sound becomes imaginary, which indicates that long wavelength perturbations grow or decay exponentially in time. Thus, the stationary state of Eq. (14) is subject to the modulational instability, which is responsible for the formation of matterwave bright solitons [84, 85, 86] in attractive BECs (see also the reviews [43, 44, 83, 87] and references therein).
2.4 Lowerdimensional condensates and relevant meanfield models.
Let us consider again a condensate confined in the harmonic trap of Eq. (9). In this case, the trap frequencies set characteristic length scales for the spatial size of the condensate through the harmonic oscillator lengths (). Another important length scale, introduced by the effective meanfield nonlinearity, is the socalled healing length defined as (with being the maximum condensate density). The healing length, being the scale over which the BEC wave function “heals” over defects, sets the spatial widths of nonlinear excitations, such as matterwave dark solitons.
Based on the above, as well as the form of the ground state [cf. Eq. (10)], it is clear that the shape of the BEC is controlled by the relative values of the trap frequencies. For example, if (i.e., for an isotropic trap), the BEC is almost spherical, while for (i.e., for an anisotropic trap) the BEC is “cigar shaped”. It is clear that such a cigarshaped BEC (a) may be a purely 3D object, (b) acquire an almost 1D character (for strongly anisotropic traps with and ), or (c) being in the socalled dimensionality crossover regime from 3D to 1D. These regimes can be described by the dimensionless parameter [88],
(16) 
where is the socalled “aspect ratio” of the trap. Particularly, if the dimensionality parameter is , the BEC locally retains its original 3D character (although it may have an elongated, quasi1D shape) and its ground state can be described by the TF approximation in all directions. On the other hand, if , excited states along the transverse direction are not energetically accessible and the BEC is effectively 1D. Apparently, this regime is extremely useful for an analytical study of matterwave dark solitons. Finally, if , the BEC is in the crossover regime between 1D and 3D, which is particularly relevant as recent matterwave dark soliton experiments have been conducted in this regime [69, 71].
Let us now discuss in more detail lowerdimensional meanfield models describing cigarshaped BECs. First, we consider the quasi1D regime () characterized by an extremely tight transverse confinement. In this case, following Refs. [50, 89, 90], the BEC wave function is separated into transverse and longitudinal components, namely . Then, the transverse component is described by the Gaussian ground state of the transverse harmonic oscillator (and, thus, the transverse width of the condensate is set by the transverse harmonic oscillator length ), while the longitudinal wave function obeys the following effectively 1D GP equation:
(17) 
where the effective 1D coupling constant is given by and . Notice that in the case under consideration, if the additional condition is fulfilled, then the longitudinal condensate density can be described by the TF approximation — see Eq. (10) with now being the 1D chemical potential (and ) [88]. Following the terminology of Ref. [69], this regime will hereafter be referred to as the TF1D regime.
Next, let us consider the effect of the deviation from onedimensionality on the longitudinal condensate dynamics. In this case, the wave function can be factorized as before, but with the transverse component assumed to depend also on the longitudinal variable (and time ) [91, 92, 93, 94]. Physically speaking, it is expected that the transverse direction will no longer be occupied by the ground state, but would still be approximated by a Gaussian function with a width that can be treated as a variational parameter [92, 93, 94]. This way, it is possible to employ different variational approaches and derive the following NLS equation for the longitudinal wave function,
(18) 
The nonlinearity function in Eq. (18) depends on the longitudinal density and may take different forms. Particularly, in Ref. [92] (where variational equations related to the minimization of the action functional were used) is found to be:
(19) 
and the respective NLS equation is known as the nonpolynomial Schrödinger equation (NPSE). On the other hand, in Refs. [93, 94] (where variational equations related to the minimization of the transverse chemical potential were used) the result for is:
(20) 
Since, as explained above, the derivation of the meanfield 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 : in this case, the width of the transverse wave function becomes and Eq. (18) — with either the nonlinearity function of Eq. (19) or that of Eq. (20) — is reduced to the 1D GP model of Eq. (17).
The above effective 1D models predict accurately ground state properties of quasi1D 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 cigarshaped BECs in a very good approximation to the 3D GrossPitaevskii (GP) equation (see, e.g., theoretical work related to matterwave dark solitons in Ref. [95] and relevant experimental results in Ref. [69]).
On the other hand, extremely weak deviations from onedimensionality can also be treated by means of a rather simple noncubic nonlinearity that can be obtained by Taylor expanding , namely:
(21) 
where and depends on the form of . In this case, Eq. (18) becomes a cubicquintic NLS (cqNLS) equation. This model was derived selfconsistently in Ref. [91], where dynamics of matterwave dark solitons in elongated BECs was considered; there, the coefficient was found to be equal to .
Here, it is worth mentioning that the quintic term in the cqNLS equation may have a different physical interpretation, namely to describe threebody interactions, regardless of the dimensionality of the system. In this case, the coefficients and in Eq. (21) are generally complex, with the imaginary parts describing inelastic two and threebody collisions, respectively [96]. As concerns the rate of the threebody collision process, it is given by [41], where is the loss rate (which is of order of – cms for various species of alkali atoms [97]). Accordingly, the decrease of the density is accounted for by the term in the time dependent GP equation, i.e., to the quintic term in the cqNLS equation.
It is also relevant to note that the NLS Eq. (18) has also been used as a meanfield model describing stronglyinteracting 1D Bose gases and, particularly, the socalled TonksGirardeau gas of impenetrable bosons [98] (see also Refs. [99, 100] for recent experimental observations). In this case, the function takes the form [101],
(22) 
and, thus, Eq. (18) becomes a quintic NLS equation. Although the applicability of this equation has been criticized (as in certain regimes it fails to predict correctly the coherence properties of the stronglyinteracting 1D Bose gases [102]), the corresponding hydrodynamic equations for the density and the phase arising from the quintic NLS equation under the Madelung transformation are welldocumented in the context of the local density approximation [103]. Additionally, it should be noticed that the timeindependent version of the quintic NLS equation has been rigorously derived from the manybody Schrödinger equation [104].
We finally mention that another lowerdimensional version of the fully 3D GP equation can be derived for “diskshaped” (alias “pancake”) condensates confined in strongly anisotropic traps with and . In such a case, a procedure similar to the one used for the derivation of Eq. (17) leads to the following dimensional NLS equation:
(23) 
where , , the effectively 2D coupling constant is given by , while the potential is given by . It should also be noticed that other effective 2D meanfield 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 diskshaped BECs.
The above models will be used below to investigate the static and dynamical properties of matterwave dark solitons arising in the respective settings.
3 General background for the study of matterwave dark solitons
3.1 NLS equation and dark soliton solutions.
We start by considering the case of a quasi1D condensate described by Eq. (17). The latter, can be expressed in the following dimensionless form:
(24) 
where the density , length, time and energy are respectively measured in units of , , and , while the potential is given by
(25) 
In the case under consideration, the normalized trap strength (aspect ratio)
is and, thus, as a first step in our analysis, the potential is
ignored.
(26) 
This equation possesses an infinite number of conserved quantities (integrals of motion); the lowestorder ones are the number of particles , the momentum , and the energy , respectively given by:
(27)  
(28)  
(29) 
It is also noted that the NLS Eq. (26) can be obtained by the EulerLagrange equation =, where the Lagrangian density is given by:
(30) 
The simplest nontrivial solution of Eq. (26) is a plane wave of wave number and frequency , namely,
(31) 
where the constant BEC density sets the chemical potential, i.e., and is an arbitrary constant phase. This solution, which is reduced to the stationary state of Eq. (14) for , is also modulationally stable as can be confirmed by a simple stability analysis (see, e.g., Refs. [34, 83]). For small densities, , the above plane wave satisfies the linear Schrödinger equation, , and the pertinent linear wave solutions of the NLS equation are characterized by the dispersion relation . Notice that if the system is characterized by a length , then the integrals of motion for the stationary solution in Eq. (31) take the values:
(32) 
The NLS equation admits nontrivial solutions, in the form of dark solitons, which can be regarded as strongly nonlinear excitations of the plane wave solution (31). In the most general case of a moving background [ in Eq. (31)], a single dark soliton solution may be expressed as [12],
(33) 
where ; here,
is the soliton center,
is an arbitrary real constant representing the initial location
of the dark soliton,
is the relative velocity between the soliton and the background given
by , the frequency is provided by
the dispersion relation of the background plane wave,
[cf. Eq. (31)]
(34) 
where is the socalled “soliton phase angle” (). Notice that although the asymptotics of the dark soliton solution (33) coincide with the ones of Eq. (31), the plane waves at have different phases; as a result, there exists a nontrivial phase jump across the dark soliton, given by:
(35) 
Note that, hereafter, we will consider the simpler case where the background of matterwave dark solitons is at rest, i.e., ; then, the frequency actually plays the role of the normalized chemical potential, namely , which is determined by the number of atoms of the condensate.
The soliton phase angle describes also the darkness of the soliton, namely,
(36) 
This way, the cases and correspond to the socalled black and gray solitons, respectively. The amplitude and velocity of the dark soliton are given (for ) by and , respectively; thus, the black soliton
(37) 
is characterized by a zero velocity, (and, thus, it is also called stationary kink), while the gray soliton moves with a finite velocity . Examples of the forms of a black and a gray soliton are illustrated in Fig. 1.
In the limiting case of a very shallow (smallamplitude) dark soliton with , the soliton velocity is close to the speed of sound which, in our units, is given by:
(38) 
The speed of sound is, therefore, the maximum possible velocity of a dark soliton which,
generally, always travels with a velocity less than the speed of sound. We finally note
that the dark soliton solution (33) has two independent parameters (for ),
one for the background, , and one for the soliton, , while there is also
a freedom (translational invariance) in selecting the initial location of the dark
soliton
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 ThomasFermi cloud, cf. Eq. (10)]. For example, in the quasi1D setting of the 1D GP Eq. (24) with the harmonic potential in Eq. (25), the “composite” wave function (describing both the background and the soliton) can be approximated as , where is the TF background and is the dark soliton wave function of Eq. (33), which satisfies the 1D GP equation for .
3.2 Dark solitons and the Inverse Scattering Transform.
The single dark soliton solution of the NLS Eq. (26) presented in the previous Section, as well as multiple dark soliton solutions (see Sec. 3.6 below), can be derived by means of the Inverse Scattering Transform (IST) [12]. A basic step of this approach is the solution of the ZakharovShabat (ZS) eigenvalue problem, with eigenvalue , for the auxiliary twocomponent eigenfunction , namely,
(39) 
with the boundary conditions , for , and , for . Here, is the amplitude of the background wave function and is a constant phase. Since the operator is selfadjoint, the ZS eigenvalue problem possesses real discrete eigenvalues , with magnitudes . Importantly, each real discrete eigenvalue corresponds to a dark soliton of depth and velocity . To make a connection to the dark soliton solutions of the NLS equation presented in the previous Section, we note that the dark soliton of Eq. (33) corresponds to a single eigenvalue .
Although the system of ZS Eqs. (39) is linear, its general solution for arbitrary initial condition is not available. Thus, various methods have been developed for the determination of the spectrum of the ZS problem, such as the socalled quasiclassical 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 multipledark 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 antisymmetric initial wave function profile of the form,
(40) 
characterized by a background density and a width (the ratio is assumed to be arbitrary). In such a case, the ZS eigenvalue problem (39) can be solved exactly [16, 17, 18] and the resulting eigenvalues of the discrete spectrum are given by and , where positive are defined as , , and is the largest integer such that . These results show that for arbitrary the initial wave function profile of Eq. (40) will always produce a black soliton [cf. Eq. (37)] at (corresponding to the first, zero eigenvalue) and additional pairs of symmetric gray solitons (corresponding to the even number of the secondary, nonzero eigenvalues), propagating to the left and to the right of the primary black soliton. Apparently, the total number of eigenvalues and, thus, the total number of solitons, is and depends on the ratio . Apart from the above example, dark soliton generation was systematically studied in Ref. [15] for a variety of initial conditions (such as boxlike dark pulses, phase steps, and others). Notice that, generally, initial wave function profiles with odd symmetry will produce an odd number of dark solitons, while profiles with an even symmetry (as, e.g., in the study of Ref. [14]) produce pairs of dark solitons; this theoretical prediction was also confirmed in experiments with optical dark solitons [109]. Furthermore, the initial phase change across the wave function plays a key role in dark soliton formation, while the number of dark solitons that are formed can be changed by small variations of the phase.
3.3 Integrals of motion and basic properties of dark solitons.
Let us now proceed by considering the integrals of motion for dark solitons. Taking into regard that Eqs. (27)–(29) refer to both the background and the soliton, one may follow Refs. [27, 28, 110, 111], and renormalize the integrals of motion so as to extract the contribution of the background [see Eqs. (32)]. This way, the renormalized integrals of motion become finite and, when calculated for the dark soliton solution (33), provide the following results (for ). The number of particles of the dark soliton reads:
(41) 
The momentum of the dark soliton is given by,
(42)  
where is given by Eq. (35) and is the speed of sound. Furthermore, the energy of the dark soliton is given by,
(43) 
while the renormalized Lagrangian density takes the form [25]:
(44) 
The renormalized integrals of motion can now be used for a better understanding of basic features of dark solitons. To be more specific, one may differentiate the expressions (42) and (43) over the soliton velocity to obtain the result,
(45) 
which shows that the dark soliton effectively behaves like a classical particle, obeying a standard equation of classical mechanics. Furthermore, it is also possible to associate an effective mass to the dark soliton, according to the equation . This way, using Eq. (42), it can readily be found that
(46) 
which shows the dark soliton is characterized by a negative effective mass. The same result, but for almost black solitons () with sufficiently small soliton velocities (), can also be obtained using Eq. (43) [112]: in this case, the energy of the dark soliton can be approximated as or, equivalently,
(47) 
where , and the soliton’s effective mass is .
3.4 Smallamplitude approximation: shallow dark solitons as KdV solitons.
As mentioned above, the case of corresponds to a smallamplitude (shallow) dark soliton, which travels with a speed close to the speed of sound, i.e., . In this case, it is possible to apply the reductive perturbation method [113] and show that, in the smallamplitude 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, socalled, smallamplitude approximation can be understood in terms of the similarity between the KdV soliton and the shallow dark soliton’s density profile: indeed, the KdV equation for a field expressed as,
(48) 
possesses a single soliton solution (see, e.g., Ref. [13]):
(49) 
(with being an arbitrary constant), which shares the same functional form with the density profile of the shallow dark soliton of the NLS equation [see Eqs. (49) and (36)]. The reduction of the cubic NLS equation to the KdV equation was first presented in Ref. [9] and later the formal connection between several integrable evolution equations was investigated in detail [115]. Importantly, such a connection is still possible even in cases of strongly perturbed NLS models, a fact that triggered various studies on dark soliton dynamics in the presence of perturbations (see, e.g., Refs. [116, 117, 118] for studies in the context of optics, as well as the recent review [119] and references therein). Generally, the advantage of the smallamplitude 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 smallamplitude dark solitons of the generalized NLS Eq. (18); in the absence of the potential (), this equation is expressed in dimensionless form as:
(50) 
where the units are the same to the ones used for Eq. (24). Then, we use the Madelung transformation (with and representing the BEC density and phase, respectively) to express Eq. (50) in the hydrodynamic form:
(51)  
(52) 
The simplest solution of Eqs. (51)–(52) is and , where . Note that in the model of Eq. (19) one has , for the model of Eq. (20), , and so on. Next, assuming slow spatial and temporal variations, we define the slow variables
(53) 
where is a formal small parameter () connected with the soliton amplitude. Additionally, we introduce asymptotic expansions for the density and phase:
(54)  
(55) 
Then, substituting Eqs. (54)–(55) into Eqs. (51)–(52), and Taylor expanding the nonlinearity function as (where ), we obtain a hierarchy of equations. In particular, Eqs. (51)–(52) lead, respectively, at the order and , to the following linear system,
(56) 
The compatibility condition of the above equations is the algebraic equation , which shows that the velocity in Eq. (53) is equal to the speed of sound, . Additionally, Eqs. (56) connect the phase and the density through the equation:
(57) 
To the next order, viz. and , Eqs. (51) and (52), respectively, yield:
(58)  
(59) 
The compatibility conditions of Eqs. (58)–(59) are the algebraic equation , along with a KdV equation [see Eq. (48)] for the unknown density :
(60) 
Thus, the density of the shallow dark soliton can be expressed as a KdV soliton [see Eq. (49)]. In terms of the original time and space variables, is expressed as follows:
(61) 
where is (as before) an arbitrary parameter [assumed to be of order ], while is the soliton velocity; the latter, is given by
(62) 
and, clearly, . Apparently, Eq. (61) describes a smallamplitude dip [of order — see Eq. (54)] on the background density of the condensate, with a phase that can be found using Eq. (57); in terms of the variables and , the result is:
(63) 
The above expression shows that the density dip is accompanied by a