# Multivalley dark solitons in multicomponent Bose-Einstein condensates with
repulsive interactions

Yan-Hong Qin1,3 Li-Chen Zhao1,2,3<EMAIL_ADDRESS>Zeng-Qiang Yang4
Liming Ling5<EMAIL_ADDRESS>1School of Physics, Northwest University,
Xi’an 710127, China  2NSFC-SPTP Peng Huanwu Center for Fundamental Theory,
Xi’an 710127, China  3Shaanxi Key Laboratory for Theoretical Physics
Frontiers, Xi’an 710127, China  4Department of Physics, School of Arts and
Sciences, Shaanxi University of Science and Technology, Xi’an 710021, China
5School of Mathematics, South China University of Technology, Guangzhou
510640, China

###### Abstract

We obtain multivalley dark soliton solutions with asymmetric or symmetric
profiles in multicomponent repulsive Bose-Einstein condensates by developing
the Darboux transformation method. We demonstrate that the width-dependent
parameters of solitons significantly affect the velocity ranges and phase jump
regions of multivalley dark solitons, in sharp contrast to scalar dark
solitons. For double-valley dark solitons, we find that the phase jump is in
the range $[0,2\pi]$, which is quite different from that of the usual single-
valley dark soliton. Based on our results, we argue that the phase jump of an
$n$-valley dark soliton could be in the range $[0,n\pi]$, supported by our
analysis extending up to five-component condensates. The interaction between a
double-valley dark soliton and a single-valley dark soliton is further
investigated, and we reveal a striking collision process in which the double-
valley dark soliton is transformed into a breather after colliding with the
single-valley dark soliton. Our analyses suggest that this breather transition
exists widely in the collision processes involving multivalley dark solitons.
The possibilities for observing these multivalley dark solitons in related
Bose-Einstein condensates experiments are discussed.

###### pacs:

03.75. Lm, 03.75. Kk, 05.45.Yv, 02.30.Ik

##  I Introduction

Multicomponent Bose-Einstein condensates (BECs) provide a good platform for
the investigation of vector solitons both theoretically and experimentally
BEC1 ; BEC2 ; BEC3  due to the abundance of intra- and interatomic
interactions. Various vector solitons have been investigated in multicomponent
BECs with attractive or repulsive interactions BEC1 ; BEC2 ; BS1 ; BS2 ; BB ;
DD1 ; DD2 ; DD3 ; D-antiD ; spinzhao ; DB1 ; DB2 ; DB3 ; DB4 ; DDB ; DBB . The
major theme of research on attractive BECs is bright solitons BEC1 ; BS1 ; BS2
; BB , while studies on dark solitons (i.e., bright-dark solitons) are
considerably hampered by their background modulation instability MI . This
characteristic makes it difficult to observe dark solitons experimentally in
attractive BECs.

In contrast, many more dark vector solitons have been experimentally observed
in multicomponent repulsive BECs, such as dark-dark solitons DD1 ; DD2 ; DD3 ,
dark-bright solitons DB1 ; DB2 ; DB3 ; DB4 , dark-antidark solitons D-antiD ,
dark-dark-bright solitons and dark-bright-bright solitons DDB . Very recently,
experimental observations of the collisions of bright-dark-bright solitons
were realized in three-component BECs with repulsive interactions DBB .
Nevertheless, the dark solitons in the abovementioned works refer mainly to
single-valley dark solitons (SVDSs). Therefore, we aim to look for multivalley
dark soliton (MVDS) solutions in repulsive BECs.

In this work, we present the exact MVDS solutions in multicomponent BECs with
repulsive interactions by further developing the Darboux transformation (DT)
method. The explicit soliton solutions admit an MVDS in one of the components
and multihump bright solitons in the other components. In particular, the
soliton width-dependent parameters have a considerable impact on both the
velocity range and the phase jump of the MVDS. The phase jump of a double-
valley (triple-valley) dark soliton can vary in the range of $[0,2\pi]$
($[0,3\pi]$). These characteristics of MVDS are distinct from those of the
well-known scalar dark soliton, for which the width that depends on both the
velocity and the phase jump can be varied in the range $[0,\pi]$TVPzhao ; DS1
; DS2 . Furthermore, we explore the collision dynamics of MVDSs. The
interaction between two MVDSs reflects the density profile variations only
after a collision. Interestingly, one MVDS can transition to a breather after
colliding with an SVDS; this can occur because the mixture of the effective
energies of solitons in the three components emerges during the collision
process. This breather transition occurs extensively in collision processes
involving MVDSs. These findings provide an important supplement for recent
reports on nondegenerate vector solitons BB2 ; Stalin1 ; Stalin2 ; Stalin3 .
We expect that more abundant MVDSs could exist in coupled BECs comprising more
components and that the phase jump of an $n$-valley dark soliton could be in
the range of $[0,n\pi]$.

The remainder of this paper is organized as follows. In Sec.ii@, we introduce
the theoretical model and present the double-valley dark soliton (DVDS)
solutions in three-component repulsive BECs; we further show the density
profiles and analyze the phase features of DVDSs. In Sec.iii@, we investigate
the collision dynamics of DVDSs and report the striking state transition
dynamics when they collide with SVDSs. In Sec.iv@, we extend our study to
four-component repulsive BECs, where triple-valley dark solitons (TVDSs) can
be obtained. Finally, the conclusions and discussion are presented in<EMAIL_ADDRESS>
##  II Double-valley dark solitons in three-component repulsive condensates

###  II.1 Physical model and double-valley dark soliton solutions

We note that it is difficult to obtain dark solitons with more than a single
valley in $N(N\leqslant 2)$-component nonlinear systems based on the existing
methods, such as the inverse scattering method DD0 ; asymptotic1 , DT method
Mat ; DT1 , Hirota’s method Hirota1 ; Hirota2 , Kadomtsev-Petviashvili
reduction method KP1 ; KP2 , and Bäcklund transformation BT . Therefore, we
first attempt to find MVDSs in three-component BECs with repulsive
interactions. At sufficiently low temperatures and in the framework of the
mean-field approach, a three-component repulsive BEC (elongated along the x
direction) is well described by the three-component Manakov model DBB :

$\displaystyle\rm{i}q_{j,t}+\frac{1}{2}q_{j,xx}-(|q_{1}|^{2}+|q_{2}|^{2}+|q_{3}|^{2})q_{j}=0,$
(1)

where $q_{j}(x,t)(j=1,2,3)$ represents the mean-field wave functions of three-
component repulsive BECs. This model can also be used to describe the
evolution of light in defocusing nonlinear optical fibers DS1 ; DS2 . A recent
experiment on bright-dark-bright solitons in repulsive BECs strongly supports
the applicability of the above integrable three-component Manakov model to
three-component BECs with repulsive interactions DBB . For their attractive
counterpart, single-hump bright soliton solutions Hirota1 ; Lakshman ; Mat ;
Dok  and even multihump bright soliton solutions have been obtained in a
multicomponent Manakov system by the DT method BB2 , Hirota method Stalin1 ;
Stalin2 ; Stalin3  and other methods PCS1 ; MDS2 . Next, we systematically
seek the MVDS solutions in the multicomponent repulsive Manakov model.

The DT method is an effective and convenient way to derive localized wave
solutions lingDS ; DT2 ; DT3 ; DT4 ; DT5 ; DT6 . Recently, it was reported
that SVDSs can be obtained through the DT method for multicomponent repulsive
Manakov systems lingDS . In this paper, we further develop the DT method
lingDS  to derive MVDS solutions in combination with the multifold DT for
deriving nondegenerate bright solitons BB2 . We find that MVDSs can be derived
by performing a multifold DT with some special constraint conditions on the
eigenfunctions of the Lax pair. For example, one DVDS can be obtained by
performing a two-fold DT with the spectral parameters written as
$\lambda_{j}=\frac{1}{2}(\xi_{j}+\frac{1}{\xi_{j}})$ and adding some special
constraint conditions to the eigenfunctions. The complex parameter
$\xi_{j}=-v_{1}+\textrm{i}w_{j}(j=1,2)$ is introduced to simplify the soliton
solution and facilitate the analysis of physical meaning of each parameter;
the real part determines the soliton’s velocity, while the imaginary part is
called a soliton width-dependent parameter. The detailed derivation of exact
DVDS solutions is given in Appendix A. For the DVDS solutions given in Eq.
(A), one DVDS is admitted in the first component, and one double-hump bright
soliton is allowed in the other two components (DBBS). By simplifying Eq. (A),
the exact general soliton solutions can be expressed as follows:

$\displaystyle q_{1}=\frac{N_{1}}{M_{1}}e^{-\textrm{i}t}$ (2a) $\displaystyle
q_{2}=-\textrm{i}2w_{1}\sqrt{1\\!-\\!v_{1}^{2}\\!-\\!w_{1}^{2}}\frac{\alpha_{1}}{\xi_{1}}\frac{N_{2}}{M_{1}}e^{\textrm{i}\left[v_{1}x-\frac{1}{2}(2+v_{1}^{2}-w_{1}^{2})t\right]}$
(2b) $\displaystyle
q_{3}=-\textrm{i}2w_{2}\sqrt{1\\!-\\!v_{1}^{2}\\!-\\!w_{2}^{2}}\frac{\beta_{2}}{\xi_{2}}\frac{N_{3}}{M_{1}}e^{\textrm{i}\left[v_{1}x-\frac{1}{2}(2+v_{1}^{2}-w_{2}^{2})t\right]}$
(2c)

with

$\displaystyle N_{1}$
$\displaystyle\\!=\\!\left[\frac{\xi_{2}^{*}}{\xi_{2}}\alpha_{1}^{2}e^{\kappa_{1}-\kappa_{2}}+\frac{\xi_{1}^{*}}{\xi_{1}}\beta_{2}^{2}e^{\kappa_{2}-\kappa_{1}}+\alpha_{1}^{2}\beta_{2}^{2}e^{\kappa_{1}+\kappa_{2}}\right]$
$\displaystyle\times(w_{1}+w_{2})^{2}+\frac{\xi_{1}^{*}\xi_{2}^{*}}{\xi_{1}\xi_{2}}(w_{1}-w_{2})^{2}e^{-\kappa_{1}-\kappa_{2}},$
$\displaystyle N_{2}$
$\displaystyle\\!=\\!(w_{1}+w_{2})\left[\beta_{2}^{2}(w_{1}+w_{2})e^{\kappa_{2}}\\!+\\!(w_{1}-w_{2})e^{-\kappa_{2}}\right],$
$\displaystyle N_{3}$
$\displaystyle\\!=\\!(w_{1}+w_{2})\left[\beta_{2}^{2}(w_{1}+w_{2})e^{\kappa_{1}}\\!+\\!(w_{2}-w_{1})e^{-\kappa_{1}}\right],$
$\displaystyle M_{1}$
$\displaystyle\\!=\\!\left[\alpha_{1}^{2}\beta_{2}^{2}e^{\kappa_{1}+\kappa_{2}}\\!+\\!\alpha_{1}^{2}e^{\kappa_{1}-\kappa_{2}}\\!+\\!\beta_{2}^{2}e^{\kappa_{2}-\kappa_{1}}\right]\\!\\!(w_{1}+w_{2})^{2}$
$\displaystyle+(w_{1}-w_{2})^{2}e^{-\kappa_{1}-\kappa_{2}},$
$\displaystyle\kappa_{1}$
$\displaystyle=w_{1}(x-v_{1}t),\xi_{1}=-v_{1}+\textrm{i}w_{1},$
$\displaystyle\kappa_{2}$
$\displaystyle=w_{2}(x-v_{1}t),\xi_{2}=-v_{1}+\textrm{i}w_{2},$

where $\alpha_{1},$ $\beta_{2}$, $v_{1}$, $w_{1}$, and $w_{2}$ are arbitrary
real constants. Among them, the parameters $w_{1}$ and $w_{2}$ are two soliton
width-dependent parameters ($w_{1}\neq w_{2}$). The parameter $v_{1}$ is the
soliton velocity, which must satisfy the constraint condition
$v_{1}^{2}+w_{m}^{2}<1$, where $w_{m}$ is the larger of the two width-
dependent parameters, meaning that the soliton width-dependent parameters
affect the velocity range. In contrast, the width of a scalar dark soliton
depends on the velocity, which cannot exceed the speed of sound BEC2 ; TVPzhao
; DS1 ; DS2 . The parameters $\alpha_{1}$ and $\beta_{2}$ are two free
parameters related to the relative valley values and the center positions of
the two valleys; therefore, these parameters nontrivially contribute to the
soliton profiles.

###  II.2 Density profiles of double-valley dark solitons

Figure 1: The density profiles of DVDSs. (a) shows an asymmetric DVDSs in
component $q_{1}$. (b) and (c) show the density profiles of asymmetric bright
solitons in components $q_{2}$ and $q_{3}$, respectively. (d-f) show the
corresponding profiles for symmetric solitons. The parameters are
$\alpha_{1}=0.5,\beta_{2}=0.2$ for asymmetric solitons (a-c) and
$\alpha_{1}=\beta_{2}=1/\sqrt{5}$ for symmetric solitons (d-f). The other
parameters are chosen as $v_{1}=0.1$, $w_{1}=0.2$, and $w_{2}=0.3$.

The solution expressed in Eqs. (2) is generally asymmetric but becomes
symmetric when the two free parameters satisfy the condition
$\alpha_{1}=\beta_{2}=\sqrt{(w_{2}-w_{1})/(w_{1}+w_{2})}$. Examples of
asymmetric and symmetric DBBS solutions are shown in Fig. 1 with solid lines
and dashed lines, respectively. The plots from left to right correspond to
component $q_{1}$ (red line), component $q_{2}$ (green line) and component
$q_{3}$ (blue line). For the asymmetric case depicted in Figs. 1(a)-(c), the
parameters are $v_{1}=0.1$, $w_{1}=0.2$, $v_{3}=0.3$, $\alpha_{1}=0.5$, and
$\beta_{2}=0.2$. We can see that a DVDS emerges in component $q_{1}$,
manifesting two spatially localized density “dips” on a uniform background. An
asymmetric double-hump bright soliton with one node is exhibited in component
$q_{2}$, and an asymmetric single-hump bright soliton without a node presents
in component $q_{3}$ QM ; these two bright solitons are quite different from
the fundamental dark-bright-bright soliton observed in DDB ; DBB . As shown in
Fig. 1(a), the two valleys of the asymmetric DVDS are not equal in either the
width or the valley values, in sharp contrast to dark solitons reported in
previous studies BEC1 ; BEC2 ; BD ; DD1 ; DD2 ; DD3 ; D-antiD ; spinzhao ; DB1
; DB2 ; DB3 ; DB4 ; DDB ; DBB ; lingDS ; Hirota1 ; Hirota2 ; KP1 ; KP2 ; DS1 ;
DS2 . This characteristic is not observed for the stationary symmetric MVDS
solution PCS2 , which was reported to exist in focusing photorefractive media
(similar to attractive multicomponent BECs). The bright solitons that present
in components $q_{2}$ and $q_{3}$ are similar to those obtained in a two-
component attractive (focusing) system BB2 ; Stalin1 ; Stalin2 ; Stalin3 .

Figure 2: The asymmetry degree vs moving velocity for DVDSs. The red solid
line and blue dashed line correspond to the asymmetric and symmetric DVDSs,
respectively. The asymmetry degree increases with increasing moving velocity
for asymmetric DVDSs. However, that for symmetric DVDSs remains symmetric with
varying velocity. The other parameters are the same as in Fig. 1.

By setting the parameters
$\alpha_{1}=\beta_{2}=\frac{\sqrt{w_{2}-w_{1}}}{\sqrt{w_{2}+w_{1}}}$ and the
other parameters the same as in the asymmetric case, the solution given in
Eqs. (2) admits symmetric soliton profiles, as shown in Figs. 1(d)-(f) (dashed
lines). These profiles are similar to those reported in previous studies MDS2
; PCS2 , but the total density can differ from the “sech-squared” form.
Moreover, the DVDS is expressed by exponential functions, which is distinctly
different from the expressions with associated Legendre polynomials PCS2 .
Moreover, the above DVDS solution is more general than those given in PCS2 ;
MDS2  since it admits more free physical parameters. The soliton velocity of
the solution in Eqs. (2) changes in the range of
$\big{[}0,\sqrt{1-w_{m}^{2}}\big{)}$. In particular, the interaction between
the solitons can be investigated analytically with further iterations by
applying the developing DT method. When the soliton width-dependent parameters
$w_{1}$ and $w_{2}$ are close to each other, both components $q_{2}$ and
$q_{3}$ can show a symmetric (or an asymmetric) double-hump bright soliton.
The characteristics of bright solitons are similar to the results in a two-
component attractive (focusing) system BB2 ; PCS1 ; Stalin1 ; Stalin2 ;
Stalin3 . In fact, nondegenerate two-component bright soliton solutions can be
reduced from the above three-component solution with the aid of the close
relations between soliton solutions in attractive and repulsive systems MDS2 .

It is well known that the valley depth of an SVDS tends to decrease with
increasing moving velocity, and the symmetric properties do not change. A
similar characteristic is admitted by symmetric DVDSs. However, we find that
this characteristic no longer holds for asymmetric DVDSs, whose asymmetry can
significantly change with variation in the moving velocity, and the valley
depth also tends to become shallower with increasing moving velocity. The
asymmetry degree is defined as

$\displaystyle R(v_{1})=\frac{I_{1}-I_{2}}{I_{1}+I_{2}},$ (3)

where $I_{1}$ and $I_{2}$ are the values of the left and right valleys in
space, respectively. We exhibit the asymmetry degree vs the moving velocity in
Fig. 2 while keeping the other parameters the same as in Fig. 1. The red solid
curve (blue dashed line) corresponds to the asymmetric (symmetric) DVDS. The
asymmetry degree crucially depends on the velocity for asymmetric DVDSs,
whereas $R(v_{1})$ is constant for symmetric DVDSs. This is also a significant
difference between symmetric and asymmetric DVDSs.

###  II.3 Phase properties of double-valley dark solitons

SVDSs are well known to admit one phase jump across the dip, which varies in
the range $[0,\pi]$ BEC2 ; BEC1 ; DS1 ; DS2 ; TVPzhao . Here, we characterize
the phase jump properties of DVDSs. Generally, the phase value can be obtained
by calculating the argument $\phi_{ds}(x)=\rm{Arg}(q_{ds})$ by dropping the
momentum phase of the background. It should be emphasized that the phase
distribution in space cannot be directly obtained by calculating only the
argument based on the real and imaginary parts of the wave function.
Importantly, the phase jump direction must be judged by the phase gradient
flow
$\partial_{x}\arctan\left[\frac{\rm{Im}(q_{ds})}{\rm{Re}(q_{ds})}\right]$.
Namely, positive (negative) phase gradient flow corresponds to a rising
(falling) phase jump. Here, we define the phase jump as
$\triangle\phi_{ds}=\phi_{ds}(x\rightarrow-\infty)-\phi_{ds}(x\rightarrow+\infty)$.
As an example, we plot the corresponding phase distribution of Fig. 1(a) in
Fig. 3(a). Remarkably, the DVDS is characterized by two phase jumps through
two valleys, in contrast to the SVDS with one phase jump BEC1 ; BEC2 ; BD ;
DD1 ; DD2 ; DD3 ; D-antiD ; spinzhao ; DB1 ; DB2 ; DB3 ; DB4 ; DDB ; DBB ;
lingDS ; Hirota1 ; Hirota2 ; KP1 ; KP2 ; DS1 ; DS2 ; TVPzhao . Moreover, the
phase jump value across a DVDS is larger than $\pi$ in this case, which is
also distinctive from the phase jump of a SVDS (which cannot exceed $\pi$)
BEC1 ; BEC2 ; BD ; DD1 ; DD2 ; DD3 ; D-antiD ; spinzhao ; DB1 ; DB2 ; DB3 ;
DB4 ; DDB ; DBB ; lingDS ; Hirota1 ; Hirota2 ; KP1 ; KP2 ; DS1 ; DS2 ; TVPzhao
. We note that SVDSs can admit a phase jump greater than $\pi$ in a saturable
self-defocusing material darker1 ; darker2  but only one phase jump across the
soliton. The phase distribution of a symmetric DVDS is similar to that of an
asymmetric DVDS, with only a small difference in the spatial position.
Additionally, a bright soliton with one node in component $q_{2}$ always
maintains a $\pi$ phase jump, where the abrupt phase change occurs at the
node, while a bright soliton without a node has no phase jump in component
$q_{3}$.

Figure 3: The phase properties of a DVDS. (a) The phase distribution of an
asymmetric DVDS corresponding to the density profile in Fig. 1(a). We can see
that the DVDS admits two phase jumps across two valleys. (b) The phase jump vs
moving velocity for DVDSs. The red solid line shows the curve for the
asymmetric DVDS, while the blue dashed line shows the curve for the symmetric
DVDS. The other parameters for both cases are identical to those in Fig. 1.
The phase jump decreases gradually with increasing velocity. (c) The phase
jump (in $\pi$ unit) vs the soliton width-dependent parameter and velocity.
The parameter $w_{1}$ is fixed at $0.2$. This demonstrates that the soliton
width-dependent parameters significantly affect the phase jump and velocity
region of DVDSs. The other parameters are $\alpha_{1}=0.5$ and
$\beta_{2}=0.2$.

Furthermore, we find that the phase jump of a DVDS changes with variation in
the moving velocity. As an example, we display the variation in the phase jump
of a DVDS with the velocity in Fig. 3(b), where the red solid line and the
blue dashed line represent asymmetric and symmetric DVDSs, respectively. The
parameters are the same as in Fig. 1 except for the velocity $v_{1}$. The
phase jump value decreases dramatically with increasing velocity. For this
case, the velocity range is $|v_{1}|\in[0,\sqrt{91}/10)$, and the phase jump
is within the range $[0.32\pi,2\pi]$. Under the limit of zero velocity, i.e.,
for a stationary DVDS, the phase jump takes a value of $2\pi$. Fig. 3(b)
suggests that the change in the phase jump with velocity for a symmetric DVDS
is identical to that for an asymmetric DVDS, and the parameters $\alpha_{1}$
and $\beta_{2}$ do not affect the phase jump region. However, further analysis
indicates that the soliton width-dependent parameters essentially determine
the phase jump range. For example, we draw the variation in the phase jump
value with the soliton width-dependent parameter $w_{2}$ in Fig. 3(c), where
$w_{1}$ is fixed to be $0.2$ herein. With varying $w_{2}$, we can calculate
the velocity region under the constraint condition $w_{j}^{2}+v_{1}^{2}<1$.
Namely, when $w_{2}\in(0,0.2)$, $|v_{1}|\in[0,2\sqrt{6}/5)$; additionally,
when $w_{2}\in(0.2,1)$, $|v_{1}|\in[0,\sqrt{1-w_{2}^{2}})$. Fig. 3(c) clearly
shows the variation in the phase jump with an increase in the width-dependent
parameter $w_{2}$, always taking a value of $2\pi$ for a static DVDS. With an
increase in the width-dependent parameter $w_{2}$, the phase jump region
decreases since the phase jump value increases at the maximum velocity. With
$w_{1}=0.2$, the phase jump range is $[0.12\pi,2\pi]$. For smaller soliton
width-dependent parameters $w_{1}$ and $w_{2}$, the phase jump value tends to
be zero when the soliton’s velocity approaches the limit of the maximum
velocity, meaning that DVDSs admit the largest phase jump range of $[0,2\pi]$.

We emphasize that the soliton width-dependent parameters also affect the
velocity region; as a result, the maximum speed of the DVDS can be much
smaller than the speed of sound for DVDSs when choosing larger width-dependent
parameters (see Fig. 3(c)). This is a notable characteristic for MVDSs that is
absent for usual scalar SVDSs DS1 ; DS2 . For usual scalar SVDSs, the maximum
velocity is the speed of sound. For a bright-dark soliton in a two-component
BEC with attractive interactions, the width and velocity are completely
independent of each other BD , and there is no limit on the moving velocity.
However, for dark-bright solitons in repulsive condensates, a change in the
soliton width also causes the velocity range to vary BD . Our detailed
analysis indicates that the relation between the soliton width and velocity
satisfies an inequality for MVDSs and dark-bright solitons BD , but this
relation becomes an equality for usual scalar SVDSs DS1 ; DS2 . This
characteristic of dark-bright solitons has not been taken seriously in
previous work DB2 . Nevertheless, we expect that the soliton width plays an
important role in the motion of dark-bright solitons in harmonic traps, in
contrast to the motion of scalar dark solitons in harmonic traps Busch .
Therefore, the effects of the soliton width should be considered when
discussing the motion of MVDSs in external potentials.

Very recently, three-component vector solitons and their collisions were
experimentally observed in BECs with repulsive interactions DBB . Motivated by
these results and the developed density and phase modulation techniques in
ultracold atomic systems DB1 ; DDB ; DBB ; DDE ; BBC , we analytically
investigate the collision dynamics of DVDSs by performing further iterations
of the developed DT method.

##  III Collision dynamics of double-valley dark solitons

The collision dynamics of DVDSs mainly include two cases: (i) a collision
between a DVDS and an SVDS, (ii) a collision between two DVDSs. We investigate
these two types of collisions based on two dark soliton solutions derived by
the three-fold and four-fold DT (see details in Appendix A). A collision
between two DVDSs usually makes each soliton’s profile vary, similar to a
collision between nondegenerate bright solitons BB2 . However, the collision
between a DVDS and an SVDS demonstrates a striking state transition process:
the DVDS transforms into a breather after colliding with the SVDS, and the
latter does not admit any oscillation, just a profile variation.

We first study the interaction between a DVDS and an SVDS based on the exact
solution Eq. (A) by implementing a three-fold DT with the introduced complex
parameters $\xi_{1}=-v_{1}+iw_{1},\xi_{2}=-v_{1}+iw_{2}$ (generating a DVDS in
component $q_{1}$) and $\xi_{3}=-v_{2}+iw_{3}$ (generating an SVDS in
component $q_{1}$) (see Appendix A for the detailed calculation process). A
typical example of the striking state transition process is illustrated in
Fig. 4(a1). The asymmetric DBBS is evidently transformed into a breather,
whose density evolution admits a periodic oscillation behavior. In contrast,
the fundamental DBBS (moving toward the left) maintains a soliton state with
slight profile variation. Such a state transition induced by a collision is
not discussed in the previous literature DD1 ; DS1 ; DS2 ; DB2 ; DB1 . Figs.
4(a2) and (a3) show the density evolution of bright solitons in components
$q_{2}$ and $q_{3}$, respectively, demonstrating that this breathing behavior
also emerges after the interaction.

Figure 4: The collision dynamics of DVDSs. (a1): Collision between one DVDS
and an SVDS. There is a striking state transition process in which a DVDS
transitions to a breather after the collision, whereas the SVDS does not
breath and experiences only slight profile variation. (a2) and (a3) show
similar state transition dynamics in the other bright soliton components.
(b1): Collision between two DVDSs. For this case, there is no state transition
occurrence for the DVDSs. (b2) and (b3) display similar collision dynamics in
the other bright soliton components. The parameters are $v_{1}=0$,
$w_{1}=0.35$, $w_{2}=0.6$, $v_{2}=-0.2$, $w_{3}=0.3$, $\alpha_{1}=1$,
$\beta_{2}=1$, and $\alpha_{3}=\beta_{3}=1$ for (a1-a3). The parameters are
$v_{1}=-0.15$, $w_{1}=0.3$, $w_{2}=0.4$, $v_{2}=0$, $w_{3}=0.35$,
$w_{4}=0.45$, $\alpha_{1}=1$, $\beta_{2}=1$, $\alpha_{3}=1$, and $\beta_{4}=1$
for (b1-b3).

To understand the state transition phenomenon, we further analyze the exact
solution given by Eq. (A) by using the asymptotic analysis technique
asymptotic1 ; asymptotic2 ; asymptotic3 . Our analysis suggests that the state
transition is induced by the mixture of effective energies of the solitons in
the three components during the collision process. Before the collision
(within the limit $t\rightarrow-\infty$, with $v_{1}>v_{2}$), the DVDS-related
vector soliton takes the following asymptotic forms:

$\displaystyle
q_{1}^{i}=\frac{f_{1}e^{-\nu_{1}}+f_{2}e^{\nu_{1}}+f_{3}e^{-\nu_{2}}+f_{4}e^{\nu_{2}}}{m_{1}e^{-\nu_{1}}+m_{2}e^{\nu_{1}}+m_{3}e^{-\nu_{2}}+m_{4}e^{\nu_{2}}}e^{-\textrm{i}t},$
(4a) $\displaystyle
q_{2}^{i}=\frac{g_{1}(g_{2}e^{-\kappa_{2}}+g_{3}e^{\kappa_{2}})}{m_{1}e^{-\nu_{1}}+m_{2}e^{\nu_{1}}+m_{3}e^{-\nu_{2}}+m_{4}e^{\nu_{2}}}e^{\textrm{i}\varphi_{1}},$
(4b) $\displaystyle
q_{3}^{i}=\frac{h_{1}(h_{2}e^{-\kappa_{1}}-h_{3}e^{\kappa_{1}})}{m_{1}e^{-\nu_{1}}+m_{2}e^{\nu_{1}}+m_{3}e^{-\nu_{2}}+m_{4}e^{\nu_{2}}}e^{\textrm{i}\varphi_{2}}.$
(4c)

After the collision (in the limit $t\rightarrow+\infty$), the asymptotic
analysis expressions for the vector soliton take the following forms:

$\displaystyle
q_{1}^{f}\\!=\\!\frac{\delta_{1}\\!\left(\delta_{2}e^{-\textrm{i}\theta}\\!+\\!\delta_{3}e^{\textrm{i}\theta}\right)\\!+\\!\delta_{4}(\\!\delta_{5}e^{\nu_{1}}\\!+\\!\delta_{6}e^{-\nu_{1}}\\!)\\!+\\!\delta_{7}(\\!\delta_{8}e^{\nu_{2}}\\!+\\!\delta_{9}e^{-\nu_{2}}\\!)}{\varrho_{1}(\varrho_{2}e^{-\textrm{i}\theta}\\!+\\!\varrho_{2}^{*}e^{\textrm{i}\theta})\\!+\\!\varrho_{3}(\\!\varrho_{4}e^{-\nu_{1}}\\!+\\!\varrho_{5}e^{\nu_{1}}\\!+\\!\varrho_{6}e^{\nu_{2}}\\!+\\!\varrho_{7}e^{-\nu_{2}}\\!)}e^{-\textrm{i}t},$
(5a) $\displaystyle
q_{2}^{f}\\!=\\!\frac{e^{\textrm{i}\varphi_{1}}\zeta_{1}(\zeta_{2}e^{-\kappa_{2}}+\zeta_{3}e^{\kappa_{2}})+e^{\textrm{i}\varphi_{2}}\zeta_{4}(\zeta_{5}e^{-\kappa_{1}}+\zeta_{6}e^{\kappa_{1}})}{\varrho_{1}(\varrho_{2}e^{-\textrm{i}\theta}\\!+\\!\varrho_{2}^{*}e^{\textrm{i}\theta})\\!+\\!\varrho_{3}(\\!\varrho_{4}e^{-\nu_{1}}\\!+\\!\varrho_{5}e^{\nu_{1}}\\!+\\!\varrho_{6}e^{\nu_{2}}\\!+\\!\varrho_{7}e^{-\nu_{2}}\\!)},$
(5b) $\displaystyle
q_{3}^{f}=\\!\frac{e^{\textrm{i}\varphi_{1}}\varsigma_{1}(\varsigma_{2}e^{-\kappa_{2}}+\varsigma_{3}e^{\kappa_{2}})+e^{\textrm{i}\varphi_{2}}\varsigma_{4}(\varsigma_{5}e^{-\kappa_{1}}+\varsigma_{6}e^{\kappa_{1}})}{\varrho_{1}(\varrho_{2}e^{-\textrm{i}\theta}\\!+\\!\varrho_{2}^{*}e^{\textrm{i}\theta})\\!+\\!\varrho_{3}(\\!\varrho_{4}e^{-\nu_{1}}\\!+\\!\varrho_{5}e^{\nu_{1}}\\!+\\!\varrho_{6}e^{\nu_{2}}\\!+\\!\varrho_{7}e^{-\nu_{2}}\\!)}.$
(5c)

In Eqs. (4) and Eqs. (5),

$\displaystyle\nu_{1}=\kappa+\kappa_{2},\nu_{2}=\kappa_{1}-\kappa_{2},\theta=\varphi_{1}-\varphi_{2},$
$\displaystyle\kappa_{1}=w_{1}(x-v_{1}t),\varphi_{1}=v_{1}x-\frac{1}{2}(v_{1}^{2}-w_{1}^{2}+2)t,$
$\displaystyle\kappa_{2}=w_{2}(x-v_{1}t),\varphi_{2}=v_{1}x-\frac{1}{2}(v_{1}^{2}-w_{2}^{2}+2)t.$

In addition, the expressions for $f_{j}$, $g_{j}$, $h_{j}$, $m_{j}$,
$\delta_{j}$, $\zeta_{j}$, $\varsigma_{j}$, and $\varrho_{j}$ are given in
Appendix B, and they are all complex constants related to $\xi_{j}(j=1,2,3)$.
From Eqs. (4), we can see that each initial soliton admits a certain velocity.
As shown in TVPzhao , a moving soliton solution can be transformed to an
eigenstate with an eigenenergy if we choose the soliton center as a frame of
reference. Therefore, it is reasonable to define an effective energy for the
solitons as $E_{j}^{*}=\frac{d\phi_{j}}{dt}$ (where $\phi_{j}$ is the phase of
the wave function). We know that the effective energy of the soliton in
components $q_{1}$, $q_{2}$ and $q_{3}$ is $E_{1}^{*}=-1$,
$E_{2}^{*}=-\frac{1}{2}(v_{1}^{2}-w_{1}^{2}+2)$ and
$E_{3}^{*}=-\frac{1}{2}(v_{1}^{2}-w_{2}^{2}+2)$, respectively. For Eqs. (5),
we can see that the effective energies of bright soliton components mix and
emerge in the DVDS component after the collision process. The breathing
behavior obviously originates from the energy mixing term in $e^{\pm
i\theta}$. Namely, the effective energy difference of bright solitons
determines the oscillation period, and the period is

$\displaystyle T$ $\displaystyle=$
$\displaystyle\frac{2\pi}{|E_{2}^{*}-E_{3}^{*}|}=\frac{4\pi}{|w_{1}^{2}-w_{2}^{2}|}.$
(6)

The oscillation period is identical among the three components. Based on this
discussion, we also revisit the collision dynamics of nondegenerate bright
solitons BB2  and find that an effective energy mixture can also emerge. This
means that the abovementioned breather behavior is also observable during the
collision between a nondegenerate bright soliton and a degenerate bright
soliton.

Next, we investigate the interaction between two DVDSs based on the exact
solution given in Eq. (33) by implementing a four-fold DT with the parameters
$\xi_{1}=-v_{1}+iw_{1}$, $\xi_{2}=-v_{1}+iw_{2}$ (which generates one DVDS),
$\xi_{3}=-v_{2}+iw_{3}$, and $\xi_{4}=-v_{2}+iw_{4}$ (which generates another
DVDS) (see Appendix A for the detailed calculation process). For this case, a
typical example of the density distribution is displayed in the second panel
of Fig. 4, for which (b1), (b2) and (b3) correspond to component $q_{1}$,
component $q_{2}$ and component $q_{3}$, respectively. As shown in Fig. 4(b1),
the collision between the two DVDSs causes their profiles to change,
accompanied by a phase shift. However, for this case, there is no state
transition occurrence for either DVDS, which is dramatically different from
the collision dynamics between the DVDS and SVDS described in Fig. 4(a1).
Moreover, the effective energy mixture no longer emerges in this case. The
inelastic collision dynamics for the other two bright soliton components (see
Figs. 4(b2) and (b3)) are similar to the those for the two nondegenerate
bright solitons that collide in BB2 .

##  IV Triple-valley dark solitons in four-component repulsive condensates

We now extend our discussion to four-component BECs with repulsive
interactions, which are described by the following repulsive four-component
Manakov model:

$\displaystyle\rm{i}q_{j,t}+\frac{1}{2}q_{j,xx}-(|q_{1}|^{2}+|q_{2}|^{2}+|q_{3}|^{2}+|q_{4}|^{2})q_{j}=0,$
(7)

where $q_{j}(x,t)(j=1,2,3,4)$ denote the four component fields in BECs. By
applying the three-fold DT with $\lambda_{j}=1/2(\xi_{j}+1/\xi_{j})$ and
$\xi_{j}=-v_{1}+iw_{j}$ ($j=1,2,3$), which is similar to the DT presented in
Appendix A, we can obtain a class of dark-bright-bright-bright soliton
solutions, which admits a TVDS in component $q_{1}$ and triple-hump bright
solitons in the other three components. The exact soliton solutions can be
expressed in the following form (we do not provide the explicit solution
process here for brevity):

$\displaystyle
q_{1}=\frac{1}{\xi_{1}\xi_{2}\xi_{3}}\frac{N_{1}}{M_{1}}e^{-\textrm{i}t},$
(8a) $\displaystyle
q_{2}=-\textrm{i}\frac{2w_{1}}{\xi_{1}}\alpha\sqrt{1-v_{1}^{2}-w_{1}^{2}}\frac{N_{2}}{M_{1}}e^{\textrm{i}\left[v_{1}x-\frac{1}{2}(2+v_{1}^{2}-w_{1}^{2})t\right]},$
(8b) $\displaystyle
q_{3}=-\textrm{i}\frac{2w_{2}}{\xi_{2}}\beta\sqrt{1-v_{1}^{2}-w_{2}^{2}}\frac{N_{3}}{M_{1}}e^{\textrm{i}\left[v_{1}x-\frac{1}{2}(2+v_{1}^{2}-w_{2}^{2})t\right]},$
(8c) $\displaystyle
q_{4}=-\textrm{i}\frac{2w_{3}}{\xi_{3}}\gamma\sqrt{1-v_{1}^{2}-w_{3}^{2}}\frac{N_{4}}{M_{1}}e^{\textrm{i}\left[v_{1}x-\frac{1}{2}(2+v_{1}^{2}-w_{3}^{2})t\right]},$
(8d)

with

$\displaystyle\kappa_{1}=$ $\displaystyle
2w_{1}(x-v_{1}t),\kappa_{2}=2w_{2}(x-v_{1}t),\kappa_{3}=2w_{3}(x-v_{1}t),$
$\displaystyle N_{1}\\!=$
$\displaystyle\xi_{1}\xi_{2}^{*}\xi_{3}^{*}\alpha^{2}\eta_{2}e^{\kappa_{1}}+\xi_{1}^{*}\xi_{2}\xi_{3}^{*}\beta^{2}\eta_{3}e^{\kappa_{2}}+\xi_{1}^{*}\xi_{2}^{*}\xi_{3}\gamma^{2}\eta_{4}e^{\kappa_{3}}+$
$\displaystyle\eta_{5}\\!\Big{[}\xi_{1}\xi_{2}\xi_{3}\alpha^{2}\beta^{2}\gamma^{2}e^{\kappa_{1}+\kappa_{2}+\kappa_{3}}\\!+\\!\xi_{1}^{*}\xi_{2}\xi_{3}\beta^{2}\gamma^{2}e^{\kappa_{2}+\kappa_{3}}\\!+\\!$
$\displaystyle\xi_{1}\xi_{2}\xi_{3}^{*}\alpha^{3}\beta^{2}e^{\kappa_{1}+\kappa_{2}}\\!+\\!\xi_{1}\xi_{2}^{*}\xi_{3}\alpha^{2}\gamma^{2}e^{\kappa_{1}+\kappa_{3}}\\!\Big{]}\\!+\\!\xi_{1}^{*}\xi_{2}^{*}\xi_{3}^{*}\eta_{1},$
$\displaystyle N_{2}\\!=$
$\displaystyle\Big{[}\rho_{1}+\rho_{2}\beta^{2}e^{\kappa_{2}}+\rho_{3}\gamma^{2}e^{\kappa_{3}}+\eta_{5}\beta^{2}\gamma^{2}e^{\kappa_{2}+\kappa_{3}}\Big{]}e^{\kappa_{1}/2},$
$\displaystyle N_{3}\\!=$
$\displaystyle\Big{[}\rho_{4}+\rho_{5}\alpha^{2}e^{\kappa_{1}}-\rho_{3}\gamma^{2}e^{\kappa_{3}}+\eta_{5}\alpha^{2}\gamma^{2}e^{\kappa_{1}+\kappa_{3}}\Big{]}e^{\kappa_{2}/2},$
$\displaystyle N_{4}\\!=$
$\displaystyle\Big{[}\rho_{6}-\rho_{5}\alpha^{2}e^{\kappa_{1}}-\rho_{2}\beta^{2}e^{\kappa_{3}}+\eta_{5}\alpha^{2}\beta^{2}e^{\kappa_{1}+\kappa_{2}}\Big{]}e^{\kappa_{3}/2},$
$\displaystyle M_{1}\\!=$
$\displaystyle\Big{[}\alpha^{2}\beta^{2}e^{\kappa_{1}+\kappa_{2}}\\!+\\!\alpha^{2}\gamma^{2}e^{\kappa_{1}+\kappa_{3}}\\!+\\!\beta^{2}\gamma^{2}e^{\kappa_{2}+\kappa_{3}}\Big{]}\eta_{5}\\!+\\!\eta_{1}$
$\displaystyle\\!+\\!\eta_{2}\alpha^{2}e^{\kappa_{1}}\\!+\\!\eta_{3}\beta^{2}e^{\kappa_{2}}\\!+\\!\eta_{4}\gamma^{2}e^{\kappa_{3}}\\!+\\!\eta_{5}\alpha^{2}\beta^{2}\gamma^{2}e^{\kappa_{1}\\!+\\!\kappa_{2}\\!+\\!\kappa_{3}}.$

In the expressions listed above,
$\eta_{1}=(w_{1}\\!-\\!w_{2})^{2}(w_{1}\\!-\\!w_{3})^{2}(w_{2}\\!-\\!w_{3})^{2}$,
$\eta_{2}=(w_{1}\\!+\\!w_{2})^{2}(w_{2}\\!-\\!w_{3})^{2}(w_{1}\\!+\\!w_{3})^{2}$,
$\eta_{3}=(w_{1}+w_{2})^{2}(w_{1}-w_{3})^{2}(w_{2}+w_{3})^{2}$,
$\eta_{4}=(w_{1}-w_{2})^{2}(w_{2}+w_{3})^{2}(w_{2}+w_{3})^{2}$,
$\eta_{5}=(w_{1}+w_{2})^{2}(w_{1}+w_{3})^{2}(w_{2}+w_{3})^{2}$,
$\rho_{1}=(w_{1}^{2}-w_{2}^{2})(w_{1}^{2}-w_{3}^{2})(w_{2}-w_{3})^{2}$,
$\rho_{2}=(w_{1}+w_{2})^{2}(w_{1}^{2}-w_{3}^{2})(w_{2}+w_{3})^{2}$,
$\rho_{3}=(w_{1}^{2}-w_{2}^{2})(w_{1}+w_{3})^{2}(w_{2}+w_{3})^{2}$,
$\rho_{4}=(w_{2}^{2}-w_{1}^{2})(w_{1}-w_{3})^{2}(w_{2}^{2}-w_{3}^{2})$,
$\rho_{5}=(w_{1}+w_{2})^{2}(w_{1}+w_{3})^{2}(w_{2}^{2}-w_{3}^{2})$, and
$\rho_{6}=(w_{1}-w_{2})^{2}(w_{1}^{2}-w_{3}^{2})(w_{2}^{2}-w_{3}^{2})$. The
parameters $v_{1}$, $w_{1}$, $w_{2}$, $w_{3}$, $\alpha$, $\beta$, and $\gamma$
are real constants that co-determine the profile and position of the TVDS. The
parameters $w_{1}$, $w_{2}$, and $w_{3}$ are three soliton width-dependent
parameters ($w_{1}\neq w_{2}\neq w_{3}$). The parameter $v_{1}$ is the soliton
velocity, which should satisfy the constraint $v_{1}^{2}+w_{j}^{2}<1$
$(j=1,2,3)$. The parameters $\alpha$, $\beta$ and $\gamma$ are three free
parameters associated with the center positions and relative values of the
three valleys. The profile of the TVDS is generally asymmetric but becomes
symmetric for certain values of the parameters $\alpha$, $\beta$, and
$\gamma$, which is similar to the above three-component case.

Figure 5: The density profile (a) and phase distribution (b) of an asymmetric
TVSD. Three phase jumps across the three density valleys can be observed. The
parameters are $v_{1}=0.1$, $w_{1}=0.2$, $w_{2}=0.28$, $w_{3}=0.33$,
$\alpha=\frac{1}{6}$, $\beta=\frac{1}{10}$, $\gamma=1$.

With the arbitrary setting of these parameters, the solution given in Eqs. (8)
shows a TVDS in component $q_{1}$, a triple-hump bright soliton with two nodes
in component $q_{2}$, a bright soliton with one node in component $q_{3}$ and
a bright soliton without nodes in component $q_{4}$. The bright solitons in
the last three components are similar to the solitons reported in three-
component attractive BECs BB2 . In this section, we mainly discuss the TVDS in
the first component.

The density profile and phase distribution of an asymmetric TVDS are displayed
in Fig. 5. Fig. 5(a) depicts the density distribution of an asymmetric TVDS.
The density profile of the TVDS is somewhat similar to those of stationary
multisoliton complexes on a background MCDS1 ; MCDS2 : both are formed as the
special nonlinear superpositions of pairs of bright and dark solitons.
However, the derivation method for the TVDS is distinctive. Consequently, our
developed DT method can be used to derive more general MVDS solutions with
moving velocities and investigate their collisions analytically through
further iterations. Fig. 5(b) demonstrates that the TVDS features three phase
jumps across three valleys. The soliton width-dependent parameter $w_{j}$ also
affects the phase jump region of the TVDS. As an example, we demonstrate the
variation in the phase jump value with the changes in the soliton width-
dependent parameters ($w_{2}$ and $w_{3}$) and the velocity in Fig. 6. The
parameter $w_{1}$ is fixed to be $0.2$. For the static TVDS solution, the
phase jump is $3\pi$. The phase jump value achieves the minimum value when the
velocity tends toward the limit of the maximum velocity. With $w_{1}=0.2$, the
phase jump range of the TVDS is $[0.14\pi,3\pi]$, but the phase jump range of
the TVDS can vary within $[0,3\pi]$ by further decreasing the value of
$w_{1}$. The phase jumps of triple-hump bright solitons in component $q_{2}$,
component $q_{3}$, and component $q_{4}$ are always $2\pi$, $\pi$ and $0$,
respectively. Moreover, the parameters $w_{j}$ can also vary the soliton
velocity region for a TVDS. The velocity range is $[0,\sqrt{1-w_{m}^{2}})$,
where $w_{m}$ is the largest of the three width-dependent parameters. These
characteristics are similar to those of the DVDS case, which could hold for
dark solitons with more valleys.

Figure 6: The phase jump (in $\pi$ unit) of a TVDS with varying velocity and
soliton width-dependent parameters $w_{2}$ and $w_{3}$. The width-dependent
parameters significantly affect the phase jump and velocity region. $w_{1}$ is
fixed to be $0.2$. The other parameters are $\alpha=1/6$, $\beta=1/10$,
$\gamma=1$, and $t=0.$

Moreover, the collision dynamics of TVDSs include three main cases: a
collision between a TVDS and an SVDS, a collision between a TVDS and a DVDS,
and a collision between two TVDSs. A breather transition also emerges for the
TVDS in the first and second cases, whereas in the third case, while the
soliton profiles similarly vary after the collision, there is no breather
behavior. The oscillation behavior is much more abundant in these cases than
in the cases for the above DVDS. We will systemically discuss their properties
in the future with the further development of analysis techniques.

##  V Conclusion and discussion

We obtain exact DVDS and TVDS solutions by further developing the DT method.
Their velocity and phase jump characteristics are characterized in detail, and
in particular, we demonstrate that changes in the soliton width-dependent
parameters have considerable influences on the velocity and phase jump ranges.
This finding indicates that the effects of the soliton width should be
considered when studying the motion of MVDSs in external potentials. The
collision dynamics of the DVDS and TVDS are also discussed. The collisions
involving DVDSs are discussed in detail, and we report a striking state
transition process in which a DVDS transitions into a breather after colliding
with an SVDS due to the mixture of the effective energies of the soliton
states in the three components. Furthermore, our analyses suggest that
breather transitions exist widely in the collision processes involving MVDSs.

Our discussion can be extended to $N$-component coupled systems. The
$(N-1)$-valley dark soliton solution can be obtained by applying the
$(N-1)$-fold DT with similar constraint conditions on the eigenfunctions of
the Lax pair. The phase jump of the $(N-1)$-valley dark soliton could vary in
the region $[0,(N-1)\pi]$. This argument is supported by our calculation up to
a five-component case. Nevertheless, further study is needed to learn how to
express the analytical solution for arbitrary $N$-component coupled systems.
The collision properties between MVDSs and degenerate vector solitons are
expected to be much more abundant than those of previously reported vector
soliton collisions Lakshman ; DBB ; DB2 ; BD ; DT6 . In particular, the state
transition between the breather and MVDS could have some important hints for
soliton state manipulation fields.

Recently, three-component vector solitons and their collisions were
experimentally observed in BECs with repulsive interactions DBB . This finding
indicates that soliton dynamics could be quantitatively described by the above
integrable repulsive three-component Manakov model. Therefore, we discuss the
possibilities to observe DVDSs in three-component repulsive BECs in
combination with well-developed quantum engineering techniques DB1 ; DDB ; DBB
; DDE ; BBC . Let us consider quasi-one-dimensional elongated BECs of
${}^{87}\rm{Rb}$. The different components are the magnetic sublevels
$m_{F}=0,\pm 1$ of the $F=1$ hyperfine manifold. In the beginning, all atoms
are prepared in the $m_{F}=0$ state. One can use spatial local control beams
to transfer atoms from the initial state $m_{F}=0$ to $m_{F}=\pm 1$. With
knowledge of the density and phase given by the solution of Eqs. (2), one
could transfer the atom and imprint phase on them simultaneously to approach
the initial state for the DVDS and bright soliton states in the corresponding
components. Our numerical simulations suggest that these soliton states are
robust against small deviations and noise, suggesting it is possible to
observe these solitons experimentally.

## ACKNOWLEDGMENTS

This work is supported by the National Natural Science Foundation of China
(Contract No. 12022513,11775176,12047502), the Major Basic Research Program of
Natural Science of Shaanxi Province (Grant No. 2018KJXX-094), Scientific
research program of Education Department of Shaanxi Provience (18JK0098) and
Scientific Research Foundation of SUST (2017BJ-30).

##  Appendix A Derivation of the double-valley dark soliton solution of
equation (1)

The $N$-component repulsive BEC can be described by the following $N$-coupled
Manakov model lingDS ; DBB :

$\textrm{i}\mathbf{q}_{t}+\frac{1}{2}\mathbf{q}_{xx}-\mathbf{q}^{{\dagger}}\mathbf{q}\mathbf{q}=0,$
(9)

where

$\displaystyle\mathbf{q}=(q_{1},q_{2},...q_{N})^{\mathsf{T}},$

which admits the following Lax pair:

$\displaystyle\psi_{x}=U(\lambda;Q)\psi,$ (10a)
$\displaystyle\psi_{t}=V(\lambda;Q)\psi,$ (10b)

with

$\displaystyle U=(\textrm{i}\lambda\sigma_{3}+\textrm{i}Q),$ (11a)
$\displaystyle V=\left[\textrm{i}\lambda^{2}\sigma_{3}+\textrm{i}\lambda
Q-\frac{1}{2}(\textrm{i}\sigma_{3}Q^{2}-\sigma_{3}Q_{x})\right],$ (11b)

where

$\displaystyle Q=\left[\begin{array}[]{cc}0&-\mathbf{q}^{{\dagger}}\\\
\mathbf{q}&\mathbf{0}_{N\times N}\\\
\end{array}\right],~{}~{}\sigma_{3}=\left[\begin{array}[]{cc}1&\mathbf{0}_{1\times
N}\\\ \mathbf{0}_{N\times 1}&-\mathbb{I}_{N\times N}\\\ \end{array}\right],$

As mentioned in the main text, multivalley dark solitons can be obtained in an
$N(N>2)$ component system. As an example, we will take the DVDS solution
derivation process in the three-component case to introduce the calculation
method for MVDSs.

To obtain the DVDS solutions, we use the following seed solutions:

$\displaystyle q_{01}=e^{-\textrm{i}t},q_{02}=0,q_{01}=0.$ (13)

First, we need to solve the Lax pair equation (10) with the above seed
solutions. We use the following gauge transformation:

$\displaystyle S=\rm{diag}(1,e^{it},1,1),$

Which converts the variable coefficient differential equation into a constant-
coefficient equation. Then, we can obtain

$\displaystyle\psi_{0,x}=U_{0}\psi_{0},$ (14a)
$\displaystyle\psi_{0,t}=\left(\frac{\rm{i}}{2}U_{0}^{2}+\lambda
U_{0}+\frac{\rm{i}}{2}\lambda^{2}\right)\psi_{0},$ (14b)

where

$\displaystyle U_{0}=\rm{i}\begin{bmatrix}\lambda&-1&0&0\\\ 1&-\lambda&0&0\\\
0&0&-\lambda&0\\\ 0&0&0&-\lambda\\\ \end{bmatrix},$

In the following, we consider the property of $U_{0}$. We can obtain the
characteristic equation of matrix $U_{0}$ at
$\lambda=\lambda_{j}=a_{j}+\textrm{i}b_{j}$:

$\begin{split}\textrm{det}(\textrm{i}\tau_{j}-U_{0})=(\lambda_{j}+\tau_{j})^{2}(\tau_{j}^{2}-\lambda_{j}^{2}+1)=0.\end{split}$
(15)

The eigenvalues of (15) are

$\displaystyle\tau_{j1}=-\sqrt{\lambda_{j}^{2}-1},\tau_{j2}=\tau_{j3}=-\lambda_{j},\tau_{j4}=\sqrt{\lambda_{j}^{2}-1}.$
(16)

To obtain the vector solution of (10), we further diagonalize the matrix
$U_{0}$. Then, we obtain

$\displaystyle\phi_{x}=\widetilde{U}_{0}\phi,~{}~{}\phi=H^{-1}S\psi,$ (17a)
$\displaystyle\phi_{t}=\left(\frac{\textrm{i}}{2}\widetilde{U}_{0}^{2}+\lambda_{j}\widetilde{U}_{0}+\frac{\textrm{i}}{2}\lambda_{j}^{2}\right)\phi,$
(17b)

where the transformation matrix $H$ can be expressed as the following form:

$H=\left(\\!\begin{array}[]{cccc}\lambda_{j}+1+\tau_{j1}&0&0&\lambda_{j}+1+\tau_{j4}\\\
\lambda_{j}+1-\tau_{j1}&0&0&\lambda_{j}+1-\tau_{j4}\\\ 0&1&0&0\\\ 0&0&1&0\\\
\end{array}\\!\right).$

Thus, we have the vector solution for (17):

$\phi_{j}\\!=\\!\left(\\!\begin{array}[]{c}\phi_{j1}\\\ \phi_{j2}\\\
\phi_{j3}\\\ \phi_{j4}\\\
\end{array}\\!\right)\\!=\\!\left(\\!\begin{array}[]{r}c_{j1}\exp{\textrm{i}[\tau_{j1}x+\frac{1}{2}(\lambda_{j}^{2}\\!+\\!2\lambda\tau_{j1}-\tau_{j1}^{2})t]}\\\
c_{j2}\exp{\textrm{i}[\tau_{j2}x+\frac{1}{2}(\lambda_{j}^{2}\\!+\\!2\lambda\tau_{j2}-\tau_{j2}^{2})t]}\\\
c_{j3}\exp{\textrm{i}[\tau_{j3}x+\frac{1}{2}(\lambda_{j}^{2}\\!+\\!2\lambda\tau_{j3}-\tau_{j3}^{2})t]}\\\
c_{j4}\exp{\textrm{i}[\tau_{j4}x+\frac{1}{2}(\lambda_{j}^{2}\\!+\\!2\lambda\tau_{j4}-\tau_{j4}^{2})t]}\\\
\end{array}\\!\right).$ (18)

The coefficients $c_{j1}$, $c_{j2}$, $c_{j3}$, and $c_{j4}$ are arbitrary
complex parameters. Then, the special solution of Lax pair (10) at
$\lambda_{j}$ can be obtained based on $\psi_{j}=S^{-1}H\phi_{j}$:

$\psi_{j}\\!=\\!\left(\\!\begin{array}[]{c}\psi_{j1}\\\ \psi_{j2}\\\
\psi_{j3}\\\ \psi_{j4}\\\
\end{array}\\!\right)\\!=\\!\left(\\!\begin{array}[]{c}(\\!1\\!+\\!\lambda\\!+\\!\tau_{j1})\phi_{j1}\\!+\\!(\\!1\\!+\\!\lambda\\!+\\!\tau_{j4})\phi_{j4}\\\
\left[(\\!1\\!+\\!\lambda\\!-\\!\tau_{j1})\phi_{j1}\\!+\\!(\\!1\\!+\\!\lambda\\!-\\!\tau_{j4})\phi_{j4}\right]\\!e^{-it}\\\
\phi_{j2}\\\ \phi_{j3}\\\ \end{array}\\!\right)\\!.\\\ $ (19)

We can see from (16) that the algebraic equation (15) has a pair of opposite
complex roots, i.e., $\tau_{j1}=-\tau_{j4}$. To obtain the DVDS, we pick only
one of them in special solution (19). Namely, we need to let $\phi_{j1}=0$
(i.e., the coefficient $c_{j1}=0$) or $\phi_{j4}=0$ (i.e., the coefficient
$c_{j4}=0$). In this paper, we choose $\phi_{j4}=0$. Then, the above special
solutions $\psi_{j}$ for spectral problem (10) at $\lambda_{j}$ are re-
expressed as follows:

$\psi_{j}=\left(\begin{array}[]{c}\psi_{j1}\\\ \psi_{j2}\\\ \psi_{j3}\\\
\psi_{j4}\\\ \end{array}\right)\\\
=\left(\\!\begin{array}[]{c}(1+\lambda+\tau_{j1})\phi_{j1}\\\
\left[(1+\lambda-\tau_{j1})\phi_{j1}\right]e^{-it}\\\ \phi_{j2}\\\
\phi_{j3}\\\ \end{array}\\!\right).$ (20)

Next, we need to perform a two-fold DT using the special solutions (20) to
derive the DVDS solutions. First, we perform the first-step iteration by
applying the DT in lingDS  with $\lambda_{1}=a_{1}+\textrm{i}b_{1}$ and
constrain the eigenfunction $\psi_{13}=\phi_{12}=0$, i.e., the coefficient
$c_{12}=0$ (or the eigenfunction $\psi_{14}=\phi_{13}=0$, i.e., the
coefficient $c_{13}=0$):

$\displaystyle\psi[1]$ $\displaystyle=$ $\displaystyle
T[1]\psi,\,\,\,T[1]=\mathbb{I}+\frac{\lambda_{1}^{*}-\lambda_{1}}{\lambda-\lambda_{1}^{*}}\frac{\psi_{1}\psi_{1}^{\dagger}\Lambda}{\psi_{1}^{\dagger}\Lambda\psi_{1}},$
$\displaystyle Q[1]$ $\displaystyle=$ $\displaystyle
Q+(\lambda_{1}-\lambda_{1}^{*})\left[\sigma_{3},\frac{\psi_{1}\psi_{1}^{\dagger}\Lambda}{\psi_{1}^{\dagger}\Lambda\psi_{1}}\right],$
(21)

where $\Lambda=\rm{diag}(1,-1,-1,-1)$ and a dagger denotes the matrix
transpose and complex conjugate. For the second-step iteration, we employ
$\psi_{2}$, which is mapped to
$\psi_{2}[1]=T[1]|_{\lambda=\lambda_{2}}\psi_{2}$ with
$\lambda_{2}=a_{2}+\textrm{i}b_{2}$, and we constrain the eigenfunction
$\psi_{24}=\phi_{23}=0$, i.e., the coefficient $c_{23}=0$ (or the
eigenfunction $\psi_{23}=\phi_{22}=0$, i.e., the coefficient $c_{22}=0$):

$\displaystyle\psi[2]=T[2]\psi[1],T[2]=\mathbb{I}\\!+\\!\frac{\lambda_{2}^{*}\\!-\\!\lambda_{2}}{\lambda\\!-\\!\lambda_{2}^{*}}\frac{\psi_{2}[1]\psi_{2}[1]^{\dagger}\Lambda}{\psi_{2}[1]^{\dagger}\Lambda\psi_{2}[1]},$
$\displaystyle
Q[2]=Q[1]+(\lambda_{2}-\lambda_{2}^{*})\left[\sigma_{3},\frac{\psi_{2}[1]\psi_{2}[1]^{\dagger}\Lambda}{\psi_{2}[1]^{\dagger}\Lambda\psi_{2}[1]}\right].$
(22)

Furthermore, we require that the solitons’ velocities in these two iteration
processes be equal. Only if the above iterative processes and velocity
requirements are met can the first component of solution $Q[2]$ be the DVDS
solution. Then, we need to analyze the relationship between the soliton
velocity and spectral parameter $\lambda_{j}$.

We find that the solitons’ velocities can be obtained by calculating the
following:

$\displaystyle\widetilde{v}_{j}=\frac{\left[\textrm{Im}(\tau_{j1})-b_{j}\right]\textrm{Re}(\tau_{j1})-\left[\textrm{Im}(\tau_{j1})+3b_{j}\right]a_{j}}{\textrm{Im}(\tau_{j1})+b_{j}},$
(23)

and the inverse soliton width is

$\displaystyle\widetilde{w}_{j}=-\left[\textrm{Im}(\tau_{j1})+b_{j}\right].$
(24)

Moreover, the spectral parameters should satisfy the constraint condition

$\displaystyle\left[b_{j}-\textrm{Im}(\tau_{j1})\right]^{2}+\left[a_{j}-\textrm{Re}(\tau_{j1})\right]^{2}-1>0.$
(25)

With the spectral parameter $\lambda_{j}=a_{j}+\textrm{i}b_{j}$, the velocity
or width of the soliton cannot be directly determined by the real or imaginary
part of the spectral parameter. This greatly increases the difficulty of
taking the equal velocity in the two-fold DT (A)-(A) with different spectral
parameters. Accordingly, we express the real and imaginary parts of the
spectrum parameter by the velocity $\widetilde{v}_{j}$ and width
$\widetilde{w}_{j}$. By combining (23)-(25) and performing further
calculations, we obtain

$\displaystyle
a_{j}=-\frac{\widetilde{v}_{j}(1+\widetilde{v}_{j}^{2}+\widetilde{w}_{j}^{2})}{2(\widetilde{v}_{j}^{2}+\widetilde{w}_{j}^{2})},b_{j}=\frac{\widetilde{w}_{j}(1-\widetilde{v}_{j}^{2}-\widetilde{w}_{j}^{2})}{2(\widetilde{v}_{j}^{2}+\widetilde{w}_{j}^{2})},$
(26)

and the constrain condition is represented as

$\displaystyle\widetilde{v}_{j}^{2}+\widetilde{w}_{j}^{2}<1.$ (27)

Figure A1: The relation between the velocity and spectral parameter expressed
by different representations. (a): The spectral parameters are written as
$\lambda_{1}=\frac{1}{2}(\xi_{1}+\frac{1}{\xi_{1}})$ and
$\xi_{1}=-v_{1}+\textrm{i}w_{1}$. For this case, there is a linear
relationship between the velocity and the real part of the spectral parameter;
namely, $v_{1}$ is just simply the velocity of the soliton. (b): The spectral
parameter is written as $\lambda_{1}=a_{1}+\textrm{i}b_{1}$. This clearly
shows a nonlinear relation between the velocity and spectral parameter.
Therefore, we introduce the parameters $\xi_{j}$ to simplify our soliton
solution, which is also much more convenient for discussing the soliton’s
physical properties. The white line corresponds to $w_{1}\neq 0$ or $b_{1}\neq
0$.

Under this representation, one can guarantee that the solitons’ velocities in
the first-step iteration (A) and the second-step iteration (A) remain exactly
equal, i.e., $\widetilde{v}_{2}=\widetilde{v}_{1}$. In other words, by
performing the iterations in (A) and (A) and setting the spectral parameters
under the conditions (26) and (27), the first component of $Q[2]$ is the DVDS
solution.

As we have mentioned above, the velocity expression plays a crucial role in
the derivation process of multivalley dark solitons. However, for the spectral
parameter with the form $\lambda_{j}=a_{j}+\textrm{i}b_{j}$, we must take the
square root of a complex number to obtain the physical parameters, namely, the
velocity and width. To facilitate this analysis, we employ a more direct
approach. Interestingly, we note that if we describe the spectral parameter
with the form $\lambda_{j}=\frac{1}{2}(\xi_{j}+\frac{1}{\xi_{j}})$ with an
arbitrary complex parameter $\xi_{j}=-v_{j}+\textrm{i}w_{j}$, the solution can
be simplified greatly. Moreover, the physical meanings of the parameters
$v_{j}$ and $w_{j}$ are much clearer than those of the parameters $a_{j}$ and
$b_{j}$, as shown in Fig. A1. Then the eigenvalues of Eq.(15) can be re-
expressed as

$\displaystyle\tau_{j1}\\!=-\\!\frac{1}{2}\left[\xi_{j}\\!-\\!\frac{1}{\xi_{j}}\right],\tau_{j2}\\!=\\!-\frac{1}{2}\left[\xi_{j}\\!+\\!\frac{1}{\xi_{j}}\right],$
$\displaystyle\tau_{j3}\\!=\\!-\frac{1}{2}\left[\xi_{j}\\!+\\!\frac{1}{\xi_{j}}\right],\tau_{j4}\\!=\\!\frac{1}{2}\left[\frac{1}{\xi_{j}}\\!-\\!\xi_{j}\right].$

To simplify the solution forms, we can set the coefficients $\phi_{j1}$,
$\phi_{j2}$, and $\phi_{j3}$ to

$\displaystyle c_{j1}=1,$ (28) $\displaystyle
c_{j2}=\frac{\xi_{j}+1}{\xi_{j}}\sqrt{1-|\xi_{j}|^{2}}\alpha_{j},$ (29)
$\displaystyle
c_{j3}=\frac{\xi_{j}+1}{\xi_{j}}\sqrt{1-|\xi_{j}|^{2}}\beta_{j},$ (30)

where $\alpha_{j}$ and $\beta_{j}$ are arbitrary constants. Then, the specific
solution (20) of spectrum problem (10) at $\lambda_{j}$ can be re-expressed as
the following form:

$\psi_{j}=\left(\begin{array}[]{c}\psi_{j1}\\\ \psi_{j2}\\\ \psi_{j3}\\\
\psi_{j4}\\\
\end{array}\right)=\left(\\!\begin{array}[]{c}(1+\xi_{j})\phi_{j1}\\\
(1+\frac{1}{\xi_{j}})\phi_{j1}e^{-it}\\\ \phi_{j2}\\\ \phi_{j3}\\\
\end{array}\\!\right).$ (31)

Then, we can obtain the DVDS solution in the first component of $Q[2]$ and the
bright solitons in the other two components by performing the following two-
fold DT processes. For the first-step DT (A), the spectral parameters are
$\lambda_{1}=\frac{1}{2}(\xi_{1}+1/\xi_{1})$ and
$\xi_{1}=-v_{1}+\textrm{i}w_{1}$, and the eigenfunction is constrained by
$\psi_{13}=\phi_{12}=0$, i.e., $\alpha_{1}=0$ (or the eigenfunction
$\psi_{14}=\phi_{13}=0$, i.e., $\beta_{1}=0$). For the second-step DT (A), the
spectral parameters are $\lambda_{2}=\frac{1}{2}(\xi_{2}+1/\xi_{2})$ and
$\xi_{2}=-v_{1}+\textrm{i}w_{2}$, and the eigenfunction is constrained by
$\psi_{24}=\phi_{23}=0$, i.e., $\beta_{2}=0$ (or the eigenfunction
$\psi_{23}=\phi_{23}=0$, i.e., $\alpha_{2}=0$). The simplified solution for
$Q[2]$ (A) is presented in Eq. (2), where $\beta_{1}=0,\alpha_{2}=0$,
$\alpha_{1}$ and $\beta_{2}$ are nonzero constants. This indicates that the
parameter $v_{1}$ is the velocity, and $w_{1}$ and $w_{2}$ are the width-
dependent parameters for DVDSs. Therefore, the physical meanings of the
parameters $v_{j}$ and $w_{j}$ are much clearer than those of the parameters
$a_{j}$ and $b_{j}$ (26), as shown in Fig. A1.

To study the collision dynamics of DVDSs, further iterations are needed. For
example, we can investigate the collision between a DVDS and an SVDS by
performing a third-step DT with $\lambda_{3}=\frac{1}{2}(\xi_{3}+1/\xi_{3})$
and $\xi_{3}=-v_{2}+\textrm{i}w_{3}$. We employ $\psi_{3}$, which is mapped to
$\psi_{3}[2]=(T[2]\psi_{3}[1])|_{\lambda=\lambda_{3}}$ with
$\psi_{3}[1]=(T[1]\psi_{3})|_{\lambda=\lambda_{3}}$:

$\displaystyle\psi[3]=T[3]\psi[2],T[3]=\mathbb{I}\\!+\\!\frac{\lambda_{3}^{*}\\!-\\!\lambda_{3}}{\lambda\\!-\\!\lambda_{3}^{*}}\frac{\psi_{3}[2]\psi_{3}[2]^{\dagger}\Lambda}{\psi_{3}[2]^{\dagger}\Lambda\psi_{3}[2]},$
$\displaystyle
Q[3]=Q[2]+(\lambda_{3}-\lambda_{3}^{*})\left[\sigma_{3},\frac{\psi_{3}[2]\psi_{3}[2]^{\dagger}\Lambda}{\psi_{3}[2]^{\dagger}\Lambda\psi_{3}[2]}\right].$
(32)

For this case, we constrain the eigenfunctions $\psi_{33}\neq 0$ and
$\psi_{34}\neq 0$, i.e., $\alpha_{3}\neq 0$ and $\beta_{3}\neq 0$, and the
eigenfunctions $\psi_{1}$ and $\psi_{2}$ are the same as in the first-step and
second-step DT. A typical example of this case is shown in Figs. 4(a1)-(a3),
demonstrating a striking collision process for which a DVDS is transformed
into a breather after colliding with an SVDS.

Naturally, by performing a fourth-step DT, one can investigate the interaction
between two DVDSs. Before performing this iteration, the eigenfunctions of
$\psi_{3}$ in the third-step DT (A) should be made to satisfy the constraint
conditions where the eigenfunction $\psi_{34}=\phi_{33}=0$, i.e.,
$\beta_{3}=0$ (or the eigenfunction $\psi_{33}=\phi_{32}=0$, i.e.,
$\alpha_{3}=0$). Then, we employ $\psi_{4}$, which is mapped to
$\psi_{4}[3]=(T[3]\psi_{4}[2])|_{\lambda=\lambda_{4}}=(T[3]T[2]T[1]\psi_{4})|_{\lambda=\lambda_{4}}$.
Then, the two DVDS solutions can be obtained as follows with the spectral
parameters $\lambda_{4}=\frac{1}{2}(\xi_{4}+1/\xi_{4})$ and
$\xi_{4}=-v_{2}+\textrm{i}w_{4}$, and we can constrain the eigenfunction
$\psi_{43}=\phi_{42}=0$, i.e., $\alpha_{4}=0$ (or the eigenfunction
$\psi_{44}=\phi_{43}=0$, i.e., $\beta_{4}=0$):

$\displaystyle
Q[4]=Q[3]+(\lambda_{4}-\lambda_{4}^{*})\left[\sigma_{3},\frac{\psi_{4}[3]\psi_{4}[3]^{\dagger}\Lambda}{\psi_{4}[3]^{\dagger}\Lambda\psi_{4}[3]}\right].$
(33)

In the four-component case, the TVDS can also be obtained by performing the
above developed DT with three-fold iterations . We do not show the detailed
derivation for this process herein. The simplified solution for TVDS is
presented in (8).

In general, the $n$-fold Darboux matrix can be constructed in the following
form:

$\begin{split}\mathbf{T}_{n}=\mathbb{I}+\mathbf{Y}_{n}\mathbf{M}_{n}^{-1}(\lambda\mathbb{I}-\mathbf{D}_{n})^{-1}\mathbf{Y}_{n}^{{\dagger}}\Lambda,\end{split}$
(34)

with

$\displaystyle\mathbf{Y}_{n}$
$\displaystyle=\left[|\psi_{1}\rangle,|\psi_{2}\rangle,\cdots,|\psi_{n}\rangle\right]=\left[\begin{array}[]{c}\Psi_{1}\\\
\Psi_{2}\end{array}\right],$ $\displaystyle\mathbf{D}_{n}$ $\displaystyle={\rm
diag}\left(\lambda_{1}^{*},\lambda_{2}^{*},\cdots,\lambda_{n}^{*}\right),$
$\displaystyle\mathbf{M}_{n}$
$\displaystyle=\left(\frac{\langle\psi_{i}|\psi_{j}\rangle}{\lambda_{i}^{*}-\lambda_{j}}\right)_{1\leq
i,j\leq n}.$

$\Psi_{1}$ is a $1\times n$ matrix, and $\Psi_{2}$ is an $N\times n$ matrix.
The Bäcklund transformation between the old potential functions and the new
functions is expressed as follows:

$\begin{split}\mathbf{q}[n]&=\mathbf{q}+2\Psi_{2}\mathbf{M}^{-1}\Psi_{1}^{\dagger}.\\\
\end{split}$ (35)

## Appendix B

The explicit expressions for $f_{j}$, $g_{j}$, $h_{j}$, and $m_{j}$ of the
asymptotic expressions in Eqs. (4) are

$\displaystyle f_{1}$
$\displaystyle=\frac{\xi_{1}^{*}\xi_{2}^{*}\xi_{3}^{*}(w_{1}-w_{2})^{2}}{\xi_{1}\xi_{2}\xi_{3}(w_{1}+w_{2})^{2}}\sigma_{5}\sigma_{6},f_{2}=\frac{\xi_{3}^{*}}{\xi_{3}}\sigma_{7}\sigma_{8}\alpha_{1}^{2}\beta_{2}^{2},$
$\displaystyle f_{3}$
$\displaystyle=\frac{\xi_{1}^{*}\xi_{3}^{*}}{\xi_{1}\xi_{3}}\beta_{2}^{2}\sigma_{5}\sigma_{8},f_{4}=\frac{\xi_{2}^{*}\xi_{3}^{*}}{\xi_{2}\xi_{3}}\alpha_{1}^{2}\sigma_{6}\sigma_{7},$
$\displaystyle g_{1}$
$\displaystyle=\frac{2\alpha_{1}w_{1}p_{1}}{v_{1}-\textrm{i}w_{1}}\left[2(v_{1}-v_{2})w_{3}+\textrm{i}\sigma_{1}\right],$
$\displaystyle g_{2}$
$\displaystyle=\frac{w_{1}-w_{2}}{w_{1}+w_{2}}\sigma_{2},g_{3}=\sigma_{3}\beta^{2},$
$\displaystyle h_{1}$
$\displaystyle=\frac{2\beta_{2}w_{2}p_{2}}{v_{1}-\textrm{i}w_{2}}\left[2(v_{2}-v_{1})w_{3}-\textrm{i}\sigma_{2}\right],$
$\displaystyle h_{2}$
$\displaystyle=\frac{w_{1}-w_{2}}{w_{1}+w_{2}}\sigma_{1},h_{3}=\sigma_{4}\alpha^{2},$
$\displaystyle m_{1}$
$\displaystyle=\frac{(w_{1}-w_{2})^{2}}{(w_{1}+w_{2})^{2}}\sigma_{5}\sigma_{6},m_{2}=\sigma_{7}\sigma_{8}\alpha_{1}^{2}\beta_{2}^{2},$
$\displaystyle m_{3}$
$\displaystyle=\sigma_{6}\sigma_{7}\alpha_{1}^{2},m_{4}=\sigma_{5}\sigma_{8}\beta_{2}^{2},$
$\displaystyle p_{1}$
$\displaystyle=\sqrt{1-|\xi_{1}|^{2}},p_{2}=\sqrt{1-|\xi_{2}|^{2}},$
$\displaystyle\xi_{1}$
$\displaystyle=-v_{1}+\textrm{i}w_{1},\xi_{2}=-v_{1}+\textrm{i}w_{2},\xi_{3}=-v_{2}+\textrm{i}w_{3},$
$\displaystyle\sigma_{1}$
$\displaystyle=(v_{1}\\!-\\!v_{2})^{2}\\!+\\!(w_{1}^{2}\\!-\\!w_{3}^{2}),\sigma_{2}\\!=\\!(v_{1}\\!-\\!v_{2})^{2}\\!+\\!(w_{2}^{2}\\!-\\!w_{3}^{2}),$
$\displaystyle\sigma_{3}$
$\displaystyle=(v_{1}\\!-\\!v_{2})^{2}\\!+\\!(w_{2}^{2}\\!+\\!w_{3}^{2}),\sigma_{4}\\!=\\!(v_{1}\\!-\\!v_{2})^{2}\\!+\\!(w_{1}^{2}\\!+\\!w_{3}^{2}),$
$\displaystyle\sigma_{5}$
$\displaystyle=(v_{1}\\!-\\!v_{2})^{2}\\!+\\!(w_{1}\\!-\\!w_{3})^{2},\sigma_{6}\\!=\\!(v_{1}\\!-\\!v_{2})^{2}\\!+\\!(w_{2}\\!-\\!w_{3})^{2},$
$\displaystyle\sigma_{7}$
$\displaystyle=(v_{1}\\!-\\!v_{2})^{2}\\!+\\!(w_{1}\\!+\\!w_{3})^{2},\sigma_{8}\\!=\\!(v_{1}\\!-\\!v_{2})^{2}\\!+\\!(w_{2}\\!+\\!w_{3})^{2}.$

The explicit expressions for $\delta_{j}$, $\zeta_{j}$, $\varsigma_{j}$, and
$\varrho_{j}$ of the asymptotic expressions in Eqs. (5) are

$\displaystyle R_{j}$
$\displaystyle\\!=\\!|1+\xi_{j}|^{2},S_{j}\\!=\\!1+\xi_{j}^{2},\Lambda_{ij}=\xi_{i}\xi_{j}^{*}-1,\Xi_{ij}=\xi_{i}-\xi_{j}^{*},$
$\displaystyle\Gamma_{j}$
$\displaystyle\\!=\\!(\xi_{j}+\xi_{j}^{*})(1+|\xi_{j}|^{2}),G_{j}=\Xi_{ij}\Lambda_{ij},(i,j=1,2,3)$
$\displaystyle\delta_{1}$
$\displaystyle\\!=\\!p_{1}p_{2}\alpha_{1}\beta_{2}\frac{R_{1}R_{2}|\Lambda_{21}|^{2}G_{11}G_{22}G_{33}}{8|\xi_{1}|^{2}|\xi_{2}|^{2}\xi_{1}\xi_{2}},$
$\displaystyle\delta_{2}$
$\displaystyle\\!=\\!\frac{\Xi_{21}}{|\xi_{1}|^{2}G_{31}G_{23}},\delta_{3}=\frac{\Xi_{12}}{|\xi_{2}|^{2}G_{32}G_{13}},$
$\displaystyle\delta_{4}$
$\displaystyle\\!=\\!\frac{R_{1}R_{2}|\Lambda_{21}|^{2}\Lambda_{11}\Lambda_{22}}{4|\xi_{1}|^{2}|\xi_{2}|^{2}\xi_{1}\xi_{2}},\delta_{5}=\frac{|\xi_{1}-\xi_{2}|^{2}}{\xi_{1}\xi_{2}},$
$\displaystyle\delta_{6}$
$\displaystyle\\!=\\!\frac{-\alpha_{1}^{2}\beta_{2}^{2}|\Xi_{12}|^{2}}{4\xi_{1}^{*}\xi_{2}^{*}|G_{13}|^{2}|G_{32}|^{2}}\big{[}G_{11}G_{22}G_{33}^{2}-(2|\xi_{3}|^{2}\\!|S_{1}|^{2}\\!+$
$\displaystyle\\!2|\xi_{1}|^{2}\\!|S_{3}|^{2}\\!-\\!\Gamma_{1}\\!\Gamma_{3})(2|\xi_{3}|^{2}\\!|S_{2}|^{2}\\!+\\!2|\xi_{2}|^{2}\\!|S_{3}|^{2}\\!-\\!\Gamma_{2}\\!\Gamma_{3}\\!)\big{]}\\!,$
$\displaystyle\delta_{7}$
$\displaystyle\\!=\\!\frac{R_{1}R_{2}\Lambda_{11}\Lambda_{22}|G_{12}|^{2}}{8|\xi_{1}|^{2}|\xi_{2}|^{2}\xi_{1}\xi_{2}},$
$\displaystyle\delta_{8}$
$\displaystyle\\!=\\!\alpha_{1}^{2}\frac{2|\xi_{3}|^{2}|S_{1}|^{2}+2|\xi_{1}|^{2}|S_{3}|^{2}-\Gamma_{1}\Gamma_{3}}{\xi_{1}^{*}\xi_{2}|G_{13}|^{2}},$
$\displaystyle\delta_{9}$
$\displaystyle=\\!\beta_{2}^{2}\frac{2|\xi_{3}|^{2}|S_{2}|^{2}+2|\xi_{2}|^{2}|S_{3}|^{2}-\Gamma_{2}\Gamma_{3}}{\xi_{1}\xi_{2}^{*}|G_{23}|^{2}},$
$\displaystyle\zeta_{1}$
$\displaystyle=\\!\frac{p_{1}\alpha_{1}R_{1}R_{2}G_{12}\Lambda_{22}}{8\xi_{1}\xi_{2}^{*}\xi_{3}|\xi_{2}|^{2}G_{13}},$
$\displaystyle\zeta_{2}$
$\displaystyle=\\!\frac{\xi_{1}^{*}-\xi_{2}^{*}}{|\xi_{1}|^{4}\xi_{2}}G_{11}\Lambda_{21}\left[\xi_{1}\Gamma_{3}-2S_{1}|\xi_{3}|^{2}\right],$
$\displaystyle\zeta_{3}$
$\displaystyle=\\!\frac{8\beta_{2}^{2}}{|G_{23}|^{2}}(\lambda_{1}-\lambda_{1}^{*})(\lambda_{1}^{*}-\lambda_{2})\big{[}\big{(}4\lambda_{1}|\xi_{3}|^{2}-\Gamma_{3}\big{)}$
$\displaystyle\\!\left(2|S_{2}|^{2}|\xi_{3}|^{2}+2|S_{3}|^{2}|\xi_{2}|^{2}-\Gamma_{2}\Gamma_{3}\right)\\!+\\!G_{22}G_{33}^{2}\big{]},$
$\displaystyle\zeta_{4}$
$\displaystyle=\\!\frac{p_{2}\beta_{2}R_{1}R_{2}\Xi_{21}|\Lambda_{21}|^{2}\Lambda_{11}G_{22}G_{33}}{8|\xi_{1}|^{4}|\xi_{2}|^{4}\xi_{3}G_{23}},\zeta_{5}\\!=\\!\xi_{2}^{*}-\xi_{1}^{*},$
$\displaystyle\zeta_{6}$
$\displaystyle=\\!\alpha_{1}^{2}\Xi_{12}[|\xi_{3}|^{2}\xi_{1}^{*}S_{1}\Lambda_{11}+|\xi_{1}|^{2}S_{1}\Gamma_{3}-\xi_{1}|\xi_{1}|^{2}R_{3}$
$\displaystyle-|\xi_{1}|^{2}|\xi_{3}|^{2}S_{1}(\xi_{1}+\xi_{1}^{*})]/[\xi_{1}|G_{13}|^{2}],$
$\displaystyle\varsigma_{1}$
$\displaystyle=\\!p_{1}\alpha_{1}\frac{R_{1}R_{2}G_{12}G_{33}\Lambda_{22}}{8\xi_{1}\xi_{2}^{*}\xi_{3}|\xi_{2}|^{2}G_{13}},\varsigma_{2}\\!=\\!\frac{\xi_{1}^{*}-\xi_{2}^{*}}{|\xi_{1}|^{2}\xi_{1}^{*}\xi_{2}}G_{11}\Lambda_{21},$
$\displaystyle\varsigma_{3}$
$\displaystyle=\\!\frac{16\beta_{2}^{2}}{|G_{23}|^{2}}(\\!\lambda_{1}^{*}\\!-\\!\lambda_{1}\\!)\\!(\\!\lambda_{1}^{*}\\!-\\!\lambda_{2}\\!)\\!\big{[}S_{2}^{2}|\xi_{3}|^{2}\xi_{2}^{*}/\xi_{2}\\!-\\!\xi_{2}^{*}S_{2}\Gamma_{3}\\!+\\!|\xi_{2}|^{2}R_{3}\big{]},$
$\displaystyle\varsigma_{4}$
$\displaystyle=\\!\frac{-p_{2}\beta_{2}(\xi_{2}^{*}-\xi_{2})}{8|\xi_{1}|^{2}|\xi_{2}|^{2}\xi_{2}\xi_{3}G_{23}}R_{1}R_{2}|\Lambda_{21}|^{2}G_{22}\Lambda_{11}\Xi_{21},$
$\displaystyle\varsigma_{5}$
$\displaystyle=\\!(\xi_{2}^{*}-\xi_{1}^{*})(2S_{2}|\xi_{3}|^{2}-\xi_{2}\Gamma_{3}),$
$\displaystyle\varsigma_{6}$
$\displaystyle=\\!\alpha_{1}^{2}\Xi_{12}\big{[}\\!-\xi_{2}G_{11}G_{33}^{2}-(2S_{2}|\xi_{3}|^{2}-\xi_{2}\Gamma_{3})$
$\displaystyle(2|\xi_{3}|^{2}R_{1}+2|\xi_{1}|^{2}R_{3}-\Gamma_{1}\Gamma_{3})\big{]}/(2|G_{13}|^{2}),$
$\displaystyle\varrho_{1}$
$\displaystyle=\\!p_{1}p_{2}\alpha_{1}\beta_{2}\frac{R_{1}R_{2}G_{11}G_{22}|\Lambda_{21}|^{2}}{8|\xi_{1}|^{4}|\xi_{2}|^{4}|G_{13}^{2}G_{32}|^{2}},$
$\displaystyle\varrho_{2}$ $\displaystyle=\\!G_{32}G_{13}G_{33}\Xi_{23},$
$\displaystyle\varrho_{3}$
$\displaystyle=\\!\frac{R_{1}R_{2}\Lambda_{11}\Lambda_{22}|\Lambda_{21}|^{2}}{4|\xi_{1}|^{4}|\xi_{2}|^{4}},\varrho_{4}\\!=\\!|\xi_{1}-\xi_{2}|^{2},$
$\displaystyle\varrho_{5}$
$\displaystyle=\\!\frac{-\alpha_{1}^{2}\beta_{2}^{2}|\Xi_{12}|^{2}}{4|G_{13}|^{2}|G_{23}|^{2}}\big{[}-(2|\xi_{3}|^{2}R_{1}+2|\xi_{1}|^{2}R_{3}-\Gamma_{1}\Gamma_{3})$
$\displaystyle(2|\xi_{3}|^{2}R_{2}+2|\xi_{2}|^{2}R_{3}-\Gamma_{2}\Gamma_{3})+G_{11}G_{22}G_{33}^{2}\big{]},$
$\displaystyle\varrho_{6}$
$\displaystyle=\\!\beta_{2}^{2}|\Xi_{12}|^{2}\big{[}2|\xi_{3}|^{2}R_{2}+2|\xi_{2}|^{2}R_{3}-\Gamma_{2}\Gamma_{3}\big{]}/(2|G_{23}|^{2}),$
$\displaystyle\varrho_{7}$
$\displaystyle=\\!\alpha_{1}^{2}|\Xi_{12}|^{2}\big{[}2|\xi_{3}|^{2}R_{1}+2|\xi_{1}|^{2}R_{3}-\Gamma_{1}\Gamma_{3}\big{]}/(2|G_{13}|^{2}).$

## References

  * (1) P. G. Drazin and R. S. Johnson, _Solitons: an introduction_ (Cambridge University,1989).
  * (2) P. G. Kevrekidis, D. J. Frantzeskakis, and R. Carretero-González, _Emergent Nonlinear Phenomena in Bose-Einstein Condensates: Theory and Experiment_ (Springer, Berlin Heidelberg, 2008).
  * (3) D. J. Frantzeskakis, Dark solitons in atomic Bose-Einstein condensates: from theory to experiments,  J. Phys. A: Math. Theor. 43, 213001 (2010).
  * (4) P. Öhberg and L. Santos, Dark solitons in a two-component Bose-Einstein condensate,  Phys. Rev. Lett. 86, 2918 (2001).
  * (5) T. Busch and J. R. Anglin, Dark-bright solitons in inhomogeneous Bose-Einstein condensates,  Phys. Rev. Lett. 87, 010401 (2001).
  * (6) J. Ieda, T. Miyakawa, M. Wadati, Matter-Wave Solitons in an F=1 Spinor Bose-Einstein Condensate,  J. Phys. Soc. Jpn. 73, 2996-3007 (2004).
  * (7) F Kh. Abdullaev, A. Gammal, A. M. Kamchatnov and L. Tomio, Dynamics of bright matter wave solitons in Bose-Einstein condensate,  Int. J. Mod. Phys. B 19, 3415-3473 (2005).
  * (8) C. Becker, S. Stellmer, P. Soltan-Panahi, S. Dörscher, M. Baumert, E.-M. Richter, J. Kronjäger, K. Bongs, and K. Sengstock, Oscillations and interactions of dark and dark-bright solitons in Bose-Einstein condensates,  Nat. Phys. 4, 496 (2008).
  * (9) X.-F. Zhang, X.-H. Hu, X.-X. Liu and W. M. Liu, Vector solitons in two-component Bose-Einstein condensates with tunable interactions and harmonic potential,  Phys. Rev. A 79, 033630 (2009).
  * (10) M. A. Hoefer, J. J. Chang, C. Hamner, and P. Engels, Dark-dark solitons and modulational instability in miscible two-component Bose-Einstein condensates,  Phys. Rev.A 84, 041605(R) (2011).
  * (11) C. Hamner, J. J. Chang, P. Engels, and M. A. Hoefer, Generation of Dark-Bright Soliton Trains in Superfluid-Superfluid Counterflow,  Phys. Rev. Lett. 106, 065302 (2011).
  * (12) I. Danaila, M. A. Khamehchi, V. Gokhroo, P. Engels, and P. G. Kevrekidis, Vector dark-antidark solitary waves in multicomponent Bose-Einstein condensates,  Phys. Rev. A 94, 053617 (2016).
  * (13) M. O. D. Alotaibi and L. D. Carr, Dynamics of dark-bright vector solitons in Bose-Einstein condensates,  Phys. Rev. A 96, 013601 (2017).
  * (14) I. Morera, A. Muñoz Mateo, A. Polls, and B. Juliá-Díaz, Dark-dark-soliton dynamics in two density-coupled Bose-Einstein condensates,  Phys. Rev. A 97, 043621 (2018).
  * (15) T. M. Bersano, V. Gokhroo, M. A. Khamehchi, J. D’ Ambroise, D. J. Frantzeskakis, Three-Component Soliton States in Spinor F = 1 Bose-Einstein Condensates,  Phys. Rev. Lett. 120, 063202 (2018).
  * (16) S. Lannig, C.-M. Schmied, M. Prüfer, P. Kunkel, R. Strohmaier, H. Strobel, T. Gasenzer, P. G. Kevrekidis and M. K. Oberthaler, Collisions of three-component vector solitons in Bose-Einstein condensates,  Phys. Rev. Lett. 125, 170401 (2020).
  * (17) L.-C. Zhao, W. Wang, Q. Tang, Z.-Y. Yang, W.-L. Yang, and J. Liu, Spin soliton with a negative-positive mass transition,  Phys. Rev. A 101, 043621 (2020).
  * (18) K. Tai, A. Hasegawa, and A. Tomita, Observation of modulational instability in optical fibers,  Phys. Rev. Lett. 56, 135 (1986).
  * (19) Y. S. Kivshar and B. Luther-Davies,  Dark optical solitons: physics and applications,  Phys. Rep. 298, 81 (1998).
  * (20) Y. S. Kivshar and G. P. Agrawal, _Optical Solitons: From Fibers to Photonic Crystals_ (Academic, NewYork, 2003).
  * (21) L.-C. Zhao, Y.-H. Qin, C. Lee and J. Liu, Classification of Dark Solitons via Topological Vector Potentials, arXiv:2006.10234v1 (2020).
  * (22) S. Stalin, R. Ramakrishnan, M. Senthilvelan, and M. Lakshmanan, Nondegenerate solitons in Manakov system,  Phys. Rev. Lett.122, 043901 (2019).
  * (23) Y.-H. Qin, L.-C.Zhao, L. Ling, Nondegenerate bound-state solitons in multicomponent Bose-Einstein condensates,  Phys. Rev. E 100, 022212 (2019).
  * (24) S. Stalin, R. Ramakrishnan, and M. Lakshmanan, Nondegenerate soliton solutions in certain coupled nonlinear Schrödinger systems,  Phys. Lett. A 384 126201 (2020).
  * (25) R. Ramakrishnan, S. Stalin, and M. Lakshmanan, Nondegenerate Solitons and their Collisions in Manakov System,  Phys. Rev. E 102, 042212 (2020).
  * (26) V. E. Zakharov and A. B. Shabat, Interaction between solitons in a stable medium, Sov. Phys.-JETP 37, 823 (1973).
  * (27) K. J. Blow and N. J. Doran, Multiple dark soliton solutions of the nonlinear Schrödinger equation, Phys. Lett. 107, 55 (1985).
  * (28) V. B. Matveev and M. A. Salle, Darboux Transformations and Solitons (Springer-Verlag, Berlin, 1991).
  * (29) C. H. Gu, H. S.Hu, and Z. X. Zhou, Darboux Transformations in Integrable Systems: Theory and Their Applications to Geometry (Springer, Berlin, 2006).
  * (30) R. Radhakrishnan and M. Lakshmanan, Bright and dark soliton solutions to coupled nonlinear Schrödinger equations,  J. Phys. A: Math. Gen. 28, 2683 (1995).
  * (31) A. P. Sheppard, Y. S. Kivshar, Polarized dark solitons in isotropic Kerr media,  Phys. Rev. E 55 4773 (1997).
  * (32) Y. Ohta, D.-S. Wang, and J. Yang, General N-Dark-Dark Solitons in the Coupled Nonlinear Schrödinger Equations,  Studies in Applied Mathematics 127, 345-371 (2001).
  * (33) B.-F. Feng, General N-soliton solution to a vector nonlinear Schrödinger equation,  J. Phys. A: Math. Theor 47, 355203(2014) .
  * (34) Q. H. Park and H. J. Shin, Systematic construction of multicomponent optical solitons,  Phys. Rev. E 61, 3093 (2000).
  * (35) T. Kanna and M. Lakshmanan, Exact soliton solutions, shape changing collisions, and partially coherent solitons in coupled nonlinear Schrödinger equations,  Phys. Rev. Lett. 86, 5043 (2001).
  * (36) E. V. Doktorov and S. B. Leble, A Dressing Method in Mathematical Physics (Springer-Verlag, Berlin, 2007).
  * (37) N. Akhmediev, W. Krolikowski, and A. W. Snyder, Partially coherent solitons of variable shape,  Phys. Rev. Lett. 81, 4632 (1998); A. Ankiewicz, W. Krolikowski, and N. N. Akhmediev, Partially coherent solitons of variable shape in a slow Kerr-like medium: exact solutions,  Phys. Rev. E 59, 6079 (1999); N. Akhmediev and A. Ankiewicz, Multi-soliton complexes,  Chaos 10, 600 (2000).
  * (38) L. C. Zhao, Z. Y. Yang, and W. L. Yang, Solitons in nonlinear systems and eigen-states in quantum wells,  Chin. Phys. B 28, 010501 (2019).
  * (39) B. L. Guo, L. Ling, Q. P. Liu, Nonlinear Schrödinger equation: generalized Darboux transformation and rogue wave solutions,  Phys. Rev. E 85, 026607 (2012).
  * (40) J. S. He, H. R. Zhang, L. H. Wang, K. Porsezian, and A. S. Fokas, Generating mechanism for higher-order rogue waves,  Phys. Rev. E 87, 052914 (2013).
  * (41) L. Ling, B. Guo L.-C. Zhao, High-order rogue waves in vector nonlinear Schrödinger equations,  Phys. Rev. E 89, 041201(R) (2014).
  * (42) C. Liu, Z.-Y. Yang, L.-C. Zhao, and W.-L. Yang, Vector breathers and the inelastic interaction in a three-mode nonlinear optical fiber,  Phys. Rev. A 89, 055803 (2014).
  * (43) L. Ling, L. C. Zhao, B. Guo, Darboux transformation and multi-dark soliton for N-component nonlinear Schrödinger equations,  Nonlinearity 28, 3243-3261 (2015).
  * (44) L. Wang, Y.-J. Zhu, F.-H. Qi, M. Li and R. Guo, Modulational instability, higher-order localized wave structures, and nonlinear wave interactions for a nonautonomous Lenells-Fokas equation in inhomogeneous fibers,  Chaos 25, 063111 (2015).
  * (45) L. D. Landau and E. M. Lifshitz, Quantum Mechanics (Nauka, Moscow, 1989).
  * (46) N. Akhmediev and A. Ankiewicz, Partially Coherent Solitons on a Finite Background,  Phys. Rev. Lett. 82, 2661, (1999).
  * (47) Y.-H. Qin, Y. Wu, L.-C.Zhao, Z. Y. Yang, Interference properties of two-component matter wave solitons,  Chin. Phys. B 29, 020303 (2020).
  * (48) W. Królikowski, N. Akhmediev, B. Luther-Davies, Darker-than-black solitons: Dark solitons with total phase shift greater than $\rm{\pi}$,  Phys. Rev. E 48, 3980 (1993).
  * (49) W. Królikowski, B. Luther-Davies, Dark optical solitons in saturable nonlinear media,  Opt. Lett. 18, 188-190 (1993).
  * (50) Th. Busch and J. R. Anglin, Motion of Dark Solitons in Trapped Bose-Einstein Condensates,  Phys. Rev. Lett. 84, 2298 (2000).
  * (51) A. Weller, J. P. Ronzheimer, C. Gross, J. Esteve, and M. K. Oberthaler, Experimental Observation of Oscillating and Interacting Matter Wave Dark Solitons,  Phys. Rev. Lett. 101, 130401 (2008).
  * (52) J. H. V. Nguyen, P. Dyke, D. Luo, B. A. Malomed and R. G. Hulet, Collisions of matter-wave solitons,  Nature Phys. 10, 918-922 (2014).
  * (53) S. V. Manakov, On the theory of two-dimensional stationary self-focusing of electromagnetic waves, Zh. Eksp. Tecr. Fiz. 65, 505-516 (1973) [Sov. Phys. JETP 38, 248 (1974)]
  * (54) R. Radhakrishnan, M. Lakshmanan and J. Hietarinta, Inelastic collision and switching of coupled bright solitons in optical fibers,  Phys. Rev. E 56, 2213 (1997).
  * (55) A. A. Sukhorukov and N. N. Akhmediev, Multisoliton complexes on a background,  Phys. Rev. E 61, 5893 (2000).
  * (56) A. A. Sukhorukov, A. Ankiewicz and N. N. Akhmediev, Multisoliton complexes in a sea of radition modes,  Opt. Commun. 195, 293-302 (2001).