Domain walls and vector solitons in the coupled nonlinear Schrödinger equation

We outline a program to classify domain walls (DWs) and vector solitons in the 1D two-component coupled nonlinear Schrödinger (CNLS) equation without restricting the signs or magnitudes of any coefficients. The CNLS equation is reduced first to a complex ordinary differential equation (ODE), and then to a real ODE after imposing a restriction. In the real ODE, we identify four possible equilibria including ZZ, ZN, NZ, and NN, with Z(N) denoting a zero (nonzero) value in a component, and analyze their spatial stability. We identify two types of DWs including asymmetric DWs between ZZ and NN and symmetric DWs between ZN and NZ. We identify three codimension-1 mechanisms for generating vector solitons in the real ODE including heteroclinic cycles, local bifurcations, and exact solutions. Heteroclinic cycles are formed by assembling two DWs back-to-back and generate extended bright-bright (BB), dark-dark (DD), and dark-bright (DB) solitons. Local bifurcations include the Turing (Hamiltonian–Hopf) bifurcation that generates Turing solitons with oscillatory tails and the pitchfork bifurcation that generates DB, bright-antidark, DD, and dark-antidark solitons with monotonic tails. Exact solutions include scalar bright and dark solitons with vector amplitudes. Any codimension-1 real vector soliton can be numerically continued into a codimension-0 family. Complex vector solitons have two more parameters: a dark or antidark component can be numerically continued in the wavenumber, while a bright component can be multiplied by a constant phase factor. We introduce a numerical continuation method to find real and complex vector solitons and show that DWs and DB solitons in the immiscible regime can be related by varying bifurcation parameters. We show that collisions between two DB solitons with a nonzero phase difference in their bright components typically feature a mass exchange that changes the frequencies and phases of the two bright components and the two soliton velocities.


Introduction
The nonlinear Schrödinger (NLS) equation is well known to provide a universal description of the envelope dynamics of quasi-monochromatic plane waves in weakly nonlinear dispersive systems [1,2].The most remarkable feature of the NLS equation is complete integrability, which enables analytical solutions to be found via the inverse scattering transform for arbitrary initial conditions [3,4].As a natural progression, coupled NLS (CNLS) equations arise when multiple waves interact nonlinearly, as already noted in early studies [5,6].Subsequently, such equations find prominent applications in diverse settings including birefringent optical fibers [7], photorefractive materials [8,9], topological photonics [10], water waves [11,12], coupled sine-Gordon equations [13], and many others.Much recent interest in CNLS equations relates to Bose-Einstein condensates (BECs), where an N-component CNLS equation, often called the Gross-Pitaevskii equation (GPE), describes BECs with N bosonic components in the mean field limit [14,15].Here we consider a 1D two-component CNLS equation with general coefficients, which is a homogeneous GPE without trapping potential in certain parameter regimes: where A and B are complex fields, and d 1−2 and g 1−4 are real constants.If d 1 d 2 > 0, then Eq. ( 1) is completely integrable when d 1 = d 2 , g 1 = g 4 , g 2 = g 3 [16,17,18].If in addition g 1 g 2 > 0, Eq. ( 1) becomes the celebrated Manakov system [5], which is focusing when d 1 g 1 > 0 and defocusing when d 1 g 1 < 0.
In the binary BEC setting, there is much interest in DB and DD solitons in the defocusing regime of Eq. ( 1) [36].Most of these studies assume d 1 = d 2 , namely equal masses for the two species, with free parameters g 1−4 .The oscillation and collision dynamics of DB solitons in a harmonic trap are revealed in a pioneering theoretical study [37].The formation of DD solitons is also predicted from the interaction of dark solitons in different components in the miscible regime [38].Subsequently, families of vector solitons are produced that include not only DB and DD solitons, but also brightantidark (BAD) and dark-antidark (DAD) solitons, where antidark means a hump rather than a dip in a nonzero background [39].DB solitons are first realized experimentally using a phase imprinting method [40].Subsequently, an alternative method is proposed to generate trains of DB solitons from the counterflow of two superfluids [41], which motivates further theoretical and experimental studies on the interaction dynamics of DB solitons [42,43,44].DB solitons are also found in Eq. ( 1) with unequal masses (dispersion coefficients) in both d 1 d 2 < 0 [45] and d 1 d 2 > 0 [46,47] cases.In both cases, it is found that DB solitons with a one-hump bright component could coexist with those whose bright component has two or more humps.
In addition to various solitons reviewed above, some nonlinear partial differential equations (PDEs) like Eq. ( 1) also admit localized traveling waves connecting two different plane waves, which are known by various names including domain walls (DWs), fronts, and kinks.Such solutions play a central role in dissipative systems in various settings [48,49,50,51].They are also found in some scalar 1D conservative (dispersive) systems, possibly with dissipative (diffusive) perturbations, such as the modified Korteweg-de Vries-Burgers equation [52] and the NLS equations with stimulated Raman scattering [53].In general, the time evolution of an initial step between two plane waves in a dispersive system is the core topic in the field of dispersive hydrodynamics, which reformulates the original PDE as a system of hyperbolic conservation laws [54].In this framework, dispersive shock waves arise when the flux is convex, while rarefaction waves and DWs (kinks) arise when the flux is non-convex [55].
In two-component 1D conservative systems, such as Eq. ( 1) in the defocusing regime, DWs can form between the two components [51].In the optics setting, such DWs in the symmetric CNLS equation, known as polarization domain walls (PDWs), are indeed predicted theoretically [56,57,58,59].While there are early experimental indications of PDWs [60,61], the first direct observation of PDWs is only made recently in optical fibres with potential applications to optical data transmission [62].Applications of optical bistability to the design of photonic memory devices are also explored in recent work [63,64].In the binary BEC setting, DWs are found in the symmetric CNLS equation with a trapping potential [39] and possibly also linear interconversion terms [65].More generally, in the asymmetric CNLS equation, i.e., Eq. ( 1) with d 1 d 2 > 0 and free parameters g 1−4 , DWs are found when the immiscibility condition g 2 g 4 − g 1 g 3 > 0 holds [66,67,68].
The above studies typically make two assumptions on the coefficients of Eq. ( 1) that are relevant in the optics and/or BEC setting, the first being d 1 d 2 > 0, or more often d 1 = d 2 , and the second being g 2 g 4 > 0, or equivalently g 2 = g 4 by rescaling variables.However, even if both are assumed, certain parameter regimes beyond the traditional dichotomy between the focusing and the defocusing regimes are rarely explored.Moreover, in the general problem of two interacting quasi-monochromatic plane waves, g 2 g 4 > 0 often applies, but d 1 d 2 > 0 and d 1 d 2 < 0 are both very common.As we shall see, for DWs and vector solitons, there is a simple relation between solutions found at g 2 g 4 > 0 and those found at g 2 g 4 < 0. Thus, in this paper, we outline a program to classify DWs and vector solitons in Eq. (1) without restricting the signs or magnitudes of any coefficients.
This paper is organized as follows.In Section 2, we reduce Eq. ( 1) first to an ordinary differential equation (ODE) with complex coefficients, and then to an ODE with real coefficients after imposing a restriction.In the real ODE, we identify all possible equilibria and analyze their spatial stability.In Section 3, we identify all possible DWs between two equilibria and construct their parameter space.In Section 4, we identify three codimension-1 mechanisms for generating vector solitons in the real ODE including heteroclinic cycles, local bifurcations, and exact solutions.Heteroclinic cycles are formed by assembling two DWs back-to-back.All possible local bifurcations on the four equilibria are identified.Exact solutions include scalar bright and dark solitons with vector amplitudes.We classify these codimension-1 real vector solitons by their types and parities and discuss numerical continuation of any such solution into a codimension-0 family.We construct the parameter space of complex vector solitons depending on whether the two components are bright, dark, or antidark.In Section 5, we introduce a numerical continuation method to find real and complex vector solitons and apply this method to DWs and DB solitons in the immiscible regime.In Section 6, we study collisions between two DB solitons with a nonzero phase difference in their bright components.The paper concludes in Section 7.

Equilibria and their spatial stability
To find DWs connecting two arbitrary plane waves and vector solitons embedded in an arbitrary plane wave, we apply the following traveling wave reduction to Eq. ( 1): , where φ, ψ are complex envelopes dependent only on the comoving variable ζ = x − C g t, with C g ∈ R the group velocity of the envelope, and k A , k B ∈ R and ω A , ω B ∈ R are the respective wavenumbers and frequencies.This reduces Eq. ( 1) to a coupled system of ODEs in ζ: In this framework, a plane wave in Eq. ( 1) corresponds to an equilibrium in Eq. ( 2), namely a ζ-independent solution to Eq. ( 2).Thus, a DW connecting two plane waves in Eq. ( 1) corresponds to a heteroclinic orbit between two equilibria in Eq. ( 2), and a vector soliton embedded in a plane wave in Eq. ( 1) corresponds to a homoclinic orbit to an equilibrium in Eq. ( 2).We note that prior work, e.g.Ref. [45], often sets C g = 0 since traveling waves can be obtained from stationary waves using the Galilean transformation, but we shall keep C g as a free parameter in this paper to emphasize that both DWs and vector solitons can have nonzero velocities; see e.g.Ref. [69].
Once DWs and vector solitons are found in Eq. ( 2), their temporal stability in Eq. ( 1) can be determined using established methods.The temporal stability of plane waves, namely modulational instability, in Eq. ( 1) is already addressed in an early study [70].However, a localized structure can be temporally unstable even if its plane wave background(s) are temporally stable [71].In this section, we shall focus on the general existence of DWs and vector solitons, and we shall keep such solutions even if they are temporally unstable.Although temporally unstable solutions are less interesting than temporally stable ones, they can still be important for understanding the nonlinear dynamics of modulational instability and other complex spatiotemporal dynamics.Besides, our overall strategy is to first find solutions in a (typically codimension-1) subset of the parameter space that is analytically tractable, and then use these solutions Domain Walls and Vector Solitons in the Coupled Nonlinear Schrödinger Equation 5as starting points for numerical continuations into the full parameter space.Even if a starting point is temporally unstable, it should not be discarded since the numerical continuation could yield temporally stable solutions.

The complex ODE
The key property of Eq. ( 2) that enables the existence of homoclinic orbits is its reversibility under the action where Θ A , Θ B ∈ R and * denotes complex conjugation.This reversibility represents a reflection symmetry in the phase space of Eq. ( 2) and implies the generic existence of homoclinic orbits to an equilibrium [72].We observe that if (φ, ψ) is reversible under R 0,0 , then (e iΘ A φ, e iΘ B ψ) is reversible under R Θ A ,Θ B .Thus, we only need to find homoclinic orbits reversible under R 0,0 , and the full solution family follows from multiplying each component by a constant phase factor.We need Θ A , Θ B ∈ [0, 2π) to obtain all possible solutions, but we only need Θ A , Θ B ∈ [0, π) to obtain all possible reversibility symmetries in Eq. ( 3).The equilibria of Eq. ( 2), denoted as (φ 0 , ψ 0 ) ∈ C 2 , can be classified into four types, denoted respectively as ZZ, ZN, NZ, and NN, with Z (N) corresponding to zero (nonzero) values.There is a circle of equivalent equilibria in the complex plane for each nonzero component, so ZZ represents a single equilibrium, ZN or NZ represents a circle of equivalent equilibria, and NN represents a torus of equivalent equilibria.Thus, under suitable conditions, we expect a discrete, not continuous, family of homoclinic orbits connecting the two equilibria (φ 0 , ψ 0 ) and (φ * 0 , ψ * 0 ) that are equivalent under the reversibility R 0,0 .
Our theoretical framework of seeking solutions to a nonlinear PDE via its reduction to a spatial ODE, known as spatial dynamics [72,50], provides a systematic method to find localized structures.For generic parameter choices, neither the original PDE in Eq. ( 1) nor the spatial ODE in Eq. ( 2) is integrable, and this non-integrable regime is the focus of this paper.However, even in the integrable regime of Eq. ( 1), namely the scalar NLS equation or the Manakov system, spatial dynamics can still be used to reproduce scalar or vector solitons originally derived using integrable systems theory.

The real ODE
Although Eq. (2) provides the general setting for finding DWs and vector solitons, it is rather technical to work with since it is a second-order ODE with two complex fields, or equivalently an 8D real dynamical system.Therefore, we will first seek DWs and vector solitons in Eq. ( 2) with the restriction Domain Walls and Vector Solitons in the Coupled Nonlinear Schrödinger Equation 6and then extend these solutions to Eq. ( 2) without this restriction.After applying Eq. ( 4), Eq. ( 2) becomes the following ODE with real coefficients: where Thus, Eq. ( 6) implies that Eq. ( 5) is invariant under the parameter transformation (C g , ω A , ω B ) → (0, ωA , ωB ), which is an ODE version of the Galilean transformation.
The four types of equilibria of Eq. ( 5) are respectively given by ZZ : φ 0 = 0, ψ 0 = 0; The restriction of Eq. ( 5) to φ, ψ ∈ R yields the real ODE The equilibria of Eq. ( 8) are given by Eq. ( 7) but restricted to the real axis, so the N component can have either sign.The possible reversibility symmetries of Eq. ( 8) must be given by Eq. (3) but restricted to real coefficients.Thus, the only possible symmetries of Eq. ( 8) are up-down symmetries: where s A , s B ∈ {+, −}.Since s A,B = e 2iΘ A,B , s A,B = + corresponds to Θ A,B = 0 and s A,B = − corresponds to Θ A,B = π/2.Thus, an even component is already reversible under R 0,0 in Eq. (3), while an odd component must be divided, or equivalently multiplied, by i = e iπ/2 , to be reversible under R 0,0 in Eq. (3).The coefficients of Eq. ( 8) can be simplified by rescaling variables.This will not be carried out explicitly, but we record two key features.We shall refer to (d 1−2 , g 1−4 ) as the CNLS parameters, and (ω A , ωB ), or equivalently (C g , ω A , ω B ), as the traveling wave parameters.First, in Section 1, we stated that the only possible restriction on the CNLS parameters is g 2 g 4 > 0 on the PDE level of Eq. (1).However, this restriction is optional at the ODE level of Eq. ( 8) since the g 2 g 4 > 0 and g 2 g 4 < 0 cases are related by flipping the sign of one of the two equations in Eq. ( 8).Second, (ω A , ωB ) can be multiplied by a positive constant by rescaling (ζ, φ, ψ), so the traveling wave parameter space can be taken as a unit circle in the (ω A , ωB )-plane rather than the whole plane.Thus, (ω A , ωB ) are effectively a single parameter that can be taken as ω ≡ ωB /ω A with the caveat that each ω corresponds to two antipodal points on the unit circle.
The existence of heteroclinic and homoclinic orbits in Eq. ( 8) depends crucially on the existence and stability of its equilibria, where the stability is defined purely within Eq. ( 8) with respect to the spatial coordinate ζ.An equilibrium must be spatially hyperbolic, namely its spatial eigenvalues are either stable or unstable, not neutrally stable, to ensure the generic existence of heteroclinic or homoclinic orbits from or to it.To obtain the spatial eigenvalues λ of an equilibrium, we first perturb it as where [φ 1 , ψ 1 ] T is the spatial eigenvector.Then, we substitute this into Eq.( 8) and linearize to obtain where I 2 is the 2 × 2 identity matrix and is the Jacobian matrix with We note that Eq. ( 8) as a 4D dynamical system would have a 4 ×4 Jacobian matrix, but Eqs.(11)(12) provide a simpler and equivalent formulation.The characteristic equation of the eigenvalue problem in Eq. ( 11) is which is a quadratic equation in λ 2 .The equilibrium (φ 0 , ψ 0 ) is spatially hyperbolic if and only if Re(λ) = 0 for all four roots λ of Eq. ( 13), or equivalently Eq. ( 13) has no non-positive real roots λ 2 .The Jacobian matrix J in Eq. ( 12) is generally a 2 × 2 real matrix without any symmetry, so the configuration of its two eigenvalues λ 2 given by the two roots λ 2 of Eq. ( 13) is always symmetric with respect to the real axis.There can be two positive real roots, two negative real roots, one positive and one negative real roots, or two complex roots.After taking the square root, the configuration of these four spatial eigenvalues λ in the complex plane is always symmetric with respect to both the real and the imaginary axes, so a spatially hyperbolic equilibrium always has two stable and two unstable spatial eigenvalues; see Ref. [72] for a detailed discussion of spatial eigenvalues in 4D reversible dynamical systems and their ramifications.One can similarly analyze the spatial stability in Eq. ( 5) of any complex equilibrium in Eq. ( 7).This will not be shown explicitly, but the key conclusion is that an equilibrium and its stable or unstable eigenspace always have the same phases in both φ and ψ.Thus, the stable or unstable manifold of a complex equilibrium in Eq. ( 5) is equivalent to that of a real equilibrium in Eq. ( 8) modulo phase rotations, so we can seek any solution tending to an equilibrium in Eq. ( 8) rather than Eq. ( 5) without loss of generality.
One can check that for all localized structures including both DWs and vector solitons reviewed in Section 1, the background equilibria are always spatially hyperbolic in Eq. ( 8).These include the ZZ background of BB solitons in both the focusing Manakov system (integrable) and the focusing CNLS equation (non-integrable), the ZN or NZ background of DB solitons and the NN background of DD solitons in both the defocusing Manakov system (integrable) and the defocusing CNLS equation (nonintegrable), and so on.However, outside these well-known focusing and defocusing regimes, the spatial eigenvalues of all possible equilibria must be calculated using Eq. ( 13) to determine the possible existence of various localized structures.

Domain walls
First, we seek DWs as heteroclinic orbits between two equilibria in Eq. (8).Such an orbit requires a transverse intersection between the 2D invariant manifolds of the two equilibria in the 4D phase space, which generally exists on a codimension-1 subset of the parameter space.Thus, a heteroclinic orbit typically exists at a single value of a bifurcation parameter known as the Maxwell point [50].The Maxwell point of Eq. ( 8) can be determined analytically thanks to the existence of a real-valued conserved quantity that is the Hamiltonian if Eq. ( 8) is considered as a mechanical system.For later discussions, we provide the following expression of the Hamiltonian applicable to Eq. ( 5), namely the complex version of Eq. ( 8): It can be verified that E is constant along any trajectory of Eq. ( 5), namely E ′ = 0 assuming Eq. ( 5).In this Hamiltonian, the two terms containing first derivatives form the kinetic energy, while the remaining terms form the potential energy.The Maxwell point between two equilibria i and j is then determined by where E i represents Eq. ( 14) evaluated at equilibrium i.We note that in Ref. [67], the Maxwell point between two equilibria is determined by the potential energy alone, but the inclusion of the kinetic energy in Eq. ( 14) becomes essential if the Maxwell point needs to be determined between ζ-dependent states such as plane waves.We have checked whether Eq. ( 15) can be satisfied for each pair of equilibria in Eq. ( 7) in any parameter regime, and found that the only possibilities are either between ZZ and NN or between ZN and NZ.
The Maxwell point between ZZ and NN is remarkable since ZZ and NN are not related by any symmetry.We shall refer to such DWs as "asymmetric domain walls" and discuss them qualitatively as part of our classification program.An example time evolution in Eq. ( 1) is shown in Fig. 1 with the initial condition being two asymmetric DWs assembled back-to-back to satisfy the periodic boundary conditions in x.We observe that the asymmetric DW is a valid solution of Eq. ( 1) and persists for a short time before instability develops.We stress that Fig. 1  By contrast, the Maxwell point between ZN and NZ are more natural since ZN and NZ can be related by exchanging φ and ψ in Eq. ( 5), although this exchange symmetry is inexact for generic parameter choices.In this case, Eq. ( 15) yields the following condition for the existence of heteroclinic orbits between ZN and NZ: This heteroclinic condition defines four points on the unit circle in the (ω A , ωB )-plane, but the existence conditions for both ZN and NZ in Eq. ( 7) choose a unique point depending on the signs of g 1 and g 3 .Moreover, the spatially hyperbolic conditions for both ZN and NZ require the defocusing conditions d 1 g 1 < 0 and d 2 g 3 < 0 and the immiscibility condition d 1 d 2 (g 2 g 4 − g 1 g 3 ) > 0. The latter reduces to the immiscibility condition in the BEC setting when d 1 d 2 > 0 [66].In this setting, DWs are well studied both physically [67] and mathematically [68].
To complete the classification of DWs between ZN and NZ, we express the heteroclinic condition in Eq. ( 16) in the original traveling wave parameter space (C g , ω A , ω B ) using the Galilean transformation in Eq. ( 6).Since Eq. ( 16) defines a half-line from the origin in the (ω A , ωB )-plane, Eqs.(6,16) define a 2D angular sector denoted by A in the 3D (C 2 g , ω A , ω B )-space where C 2 g ≥ 0 by definition.For generic choices of the CNLS parameters (d 1−2 , g 1−4 ), the sector A has a 2D projection onto the (ω A , ω B )-plane, so there is a unique Maxwell point with C 2 g as the bifurcation parameter and (ω A , ω B ) fixed: However, for special choices of (d 1−2 , g 1−4 ) satisfying the sector A is perpendicular to the (ω A , ω B )-plane, so (ω A , ω B ) must satisfy Eq. ( 16) with ω replaced by ω, and C 2 g ≥ 0 is arbitrary at fixed (ω A , ω B ).This applies most notably to symmetric CNLS ( [56,57,58,39,65], and more generally to a codimension-1 subset of the CNLS parameter space defined by Eq. ( 18).
We stress that asymmetric DWs between ZZ and NN and symmetric DWs between ZN and NZ cannot be related by continuous deformation of parameters, namely they are homotopically nonequivalent.To see this, we first define the focusing and defocusing regimes for Eq. ( 8) with general coefficients.At fixed parameters, we call Eq. ( 8) focusing if the ZZ equilibrium is spatially hyperbolic, and defocusing if at least one of the two equilibria, ZN and NZ, is spatially hyperbolic.By this definition, Eq. ( 8) may be focusing, defocusing, or neither, but not both.Thus, asymmetric DWs only exist in the focusing regime, while symmetric DWs only exist in the defocusing regime.Note that prior work on DWs between ZN and NZ in binary BECs define asymmetric and symmetric DWs depending on the symmetry of their spatial profiles [69,73].In our language, both are symmetric DWs that are typically homotopically equivalent.
After classifying heteroclinic orbits in Eq. ( 5), the next question is whether these orbits persist into Eq.( 2) with the restriction in Eq. ( 4) removed.Since equilibria of Eq. ( 2) correspond to plane waves in Eq. ( 5), the Maxwell point between two equilibria in Eq. ( 2) can be determined as the Maxwell point between two plane waves in Eq. ( 5) using the Hamiltonian in Eq. ( 14).Thus, one may naively expect the existence of a heteroclinic orbit between the ZN and NZ equilibria in Eq. ( 2) at this Maxwell point.However, the time evolution in Eq. ( 1) of an initial step between arbitrary ZN and NZ plane waves typically yields a DW with two dispersive shock waves respectively to the left and to the right of the DW, and this DW does satisfy the restriction in Eq. ( 4) even if the initial condition does not.A proper treatment of this time evolution requires modulation theory [54,55], but here we tentatively conclude that the restriction in Eq. ( 4) is indispensable for the existence of DWs.

Vector solitons
In this section, we first seek vector solitons as homoclinic orbits to an equilibrium in the real ODE in Eq. ( 8), and then extend this solution family into the complex ODE in Eq. ( 2).Since Eq. ( 8) is reversible under the action R s A ,s B in Eq. ( 9), the 4D phase Domain Walls and Vector Solitons in the Coupled Nonlinear Schrödinger Equation 11space (φ, ψ, φ ′ , ψ ′ ) of Eq. ( 8) is symmetric with respect to four 2D sections fix(R s A ,s B ) with s A,B = ± specified by two conditions: , and ψ = 0, s B = + A homoclinic orbit to an equilibrium requires a transverse intersection between the 2D stable manifold of the equilibrium and a 2D section fix(R s A ,s B ).This intersection generally exists on a codimension-0 subset of the parameter space, so any spatially hyperbolic equilibrium may support a homoclinic orbit [72].In principle, all possible homoclinic orbits to an equilibrium can be obtained by computing the 2D stable manifold of this equilibrium, but invariant manifolds are generally hard to compute due to their complicated geometry.In practice, homoclinic orbits are often found using numerical continuation in a bifurcation parameter starting from an initial guess.The initial guess for numerical continuation could be an analytical ansatz, but a more systematic approach is to identify bifurcation points where a solution branch originates or terminates.Thus, our analysis centers upon three codimension-1 mechanisms to generate vector solitons in Eq. ( 8) from local and global bifurcation theory.

Heteroclinic cycles
Perhaps the most important type of global bifurcation for homoclinic orbits is the Maxwell point for heteroclinic orbits in Section 3, since a heteroclinic orbit and its symmetric partner can be assembled back-to-back to form a special type of homoclinic orbit known as a heteroclinic cycle [50].A heteroclinic orbit from equilibrium i to equilibrium j, denoted by i → j, has an even partner j → i and an odd partner −j → −i, for either component.A heteroclinic orbit can always be assembled with its even partner to form a heteroclinic cycle i → j → i, but it can only be assembled with its odd partner to form a heteroclinic cycle i → ±j → −i when j = Z since Z = −Z, not when j = N since N = −N.
Two asymmetric DWs between ZZ and NN can form a heteroclinic cycle with NN embedded in ZZ, which is an extended BB soliton, or a heteroclinic cycle with ZZ embedded in NN, which is an extended DD soliton.The extended BB soliton has multiplicity 1 and is always even-even; see the initial condition of Fig. 1 for an example.The extended DD soliton has multiplicity 4 and can be even-even, even-odd, odd-even, or odd-odd.
Two symmetric DWs between ZN and NZ can form a heteroclinic cycle with NZ embedded in ZN, or one with ZN embedded in NZ, both of which are extended DB solitons.An extended DB soliton has multiplicity 2 and can be even-even or odd-even.
Overall, an extended vector soliton with N D dark components has multiplicity 2 N D since a dark component can be even or odd.Besides, this multiplicity persists when the extended vector soliton is numerically continued, at least near the Maxwell point.
Hence, we shall refer to a vector soliton thus obtained near the Maxwell point as a quasi-extended vector soliton.Such a vector soliton can inherit certain properties from the underlying DW even if its spatial profile is not fully extended.

Local bifurcations
In addition to global bifurcations that only occur in certain parameter regimes, e.g., the immiscible regime, Eq. ( 8) also exhibits local bifurcations that can occur in virtually all parameter regimes.A codimension-1 local bifurcation occurs when an equilibrium becomes spatially hyperbolic, namely one or two pairs of spatial eigenvalues given by Eq. ( 13) leave the imaginary axis, as a bifurcation parameter, say ω, varies.A codimension-2 bifurcation occurs when two codimension-1 bifurcations collide, or the nature of a codimension-1 bifurcation changes, e.g., from supercritical to subcritical, as two bifurcation parameters, say ω and one of (d 1−2 , g 1−4 ), vary.
There are two types of codimension-1 local bifurcations that can yield homoclinic orbits [72].The first type is the Turing or Hamiltonian-Hopf bifurcation, where four purely imaginary eigenvalues collide pairwise on the imaginary axis to yield a complex quartet of eigenvalues.The second type is the pitchfork bifurcation, where two purely imaginary eigenvalues collide at the origin to yield two real eigenvalues.Note that the pitchfork bifurcation reflects the up-down symmetry of Eq. ( 8), and a saddle-node bifurcation would occur without such symmetries.
The defining feature of the Turing bifurcation is its ability to generate spatially periodic states.As such, it has played a key role for pattern formation in nonlinear dissipative PDEs for decades.More recently, the subcritical Turing bifurcation has been much studied since it can generate spatially localized states with rich bifurcation structures [50].We have found that a subcritical Turing bifurcation can occur in Eq. ( 8), but only on the NN equilibrium.This bifurcation creates two branches of even-even localized states whose spatial profiles relative to NN differ by an overall sign to leading nontrivial order.We shall refer to localized states near this bifurcation as weakly nonlinear "Turing solitons" and discuss them qualitatively as part of our classification program.An example time evolution in Eq. ( 1) is shown in Fig. 2 with the initial condition being a weakly nonlinear Turing soliton.We observe that the Turing soliton is a valid solution of Eq. ( 1) and persists for a short time before instability develops.We stress that Fig. 2 only provides a preliminary example of Turing solitons, and we cannot rule out the possible existence of stable Turing solitons with other parameter choices.A systematic study of Turing solitons and related solutions is left as future work.
Generally, we define a Turing soliton as any vector soliton embedded in NN with an oscillatory tail due to the complex spatial eigenvalues of NN.A component of a Turing soliton can be called dark if the center is a dip or antidark if the center is a hump.For example, the initial condition of Fig. 2 can be called a DAD Turing soliton since the center of the A (φ) component is a hump and the center of the B (ψ) component is a dip.However, apart from the central hump or dip, a DAD Turing soliton has an infinite number of humps and dips in its oscillatory tail.If the DAD Turing soliton is weakly nonlinear, then the oscillatory tail can have a larger L2-norm than the central hump or dip.By contrast, a normal DAD soliton has a monotonic tail due to the real spatial eigenvalues of NN, so the central hump or dip generally dominates.In a dynamical system with an up-down symmetry, there is always an equilibrium at zero, and the pitchfork bifurcation causes two nonzero equilibria related by the updown symmetry to emerge from the zero equilibrium as a parameter varies.A pitchfork bifurcation in Eq. ( 8) can occur on a zero component regardless of whether the other component is zero or nonzero.Thus, ZN or NZ can bifurcate from ZZ, and NN can bifurcate from ZN or NZ.For the latter, the bifurcation parameter can be chosen as either ωA or ωB .However, for the former, the bifurcation parameter can only be chosen as the ωA,B that the ZN or NZ equilibrium depends on.Thus, there are no ambiguities if the ωA,B value where an equilibrium i bifurcates from an equilibrium j is denoted by ωj A,B .For example, the ωA value where NN bifurcates from ZN is denoted by ωZN A .A network can be constructed to visualize these pitchfork bifurcations with the nodes representing the four equilibria and the edges representing the bifurcation points; see Fig. 3.Note that the two types of Maxwell points in Section 3 represented on this network would occupy the two diagonal edges.
Analytical expressions can be derived for weakly nonlinear vector solitons near these pitchfork bifurcation points using asymptotic analysis.The details of these rather standard calculations will be omitted, but we note that asymptotic analyses near different pitchfork bifurcations always yield the same second-order ODE, known as the normal form, only with different coefficients [72].This normal form ODE can admit, in different parameter regimes, either an even homoclinic orbit to the zero equilibrium, or an odd homoclinic orbit connecting the two nonzero equilibria related by the up-down symmetry.
Near ωZZ A = 0 or ωZZ B = 0, one of the two components remains zero past the bifurcation point, so we obtain the pitchfork bifurcation in the scalar NLS equation that yields the well-known branches of exact solutions.The even homoclinic orbit corresponds to scalar bright solitons in the focusing regime with a sech shape, while the odd homoclinic orbit corresponds to scalar dark solitons in the defocusing regime with a tanh shape.These branches of exact solutions are not vector solitons and so will not be shown explicitly.However, one can identify codimension-1 bifurcation points on the branch of scalar dark solitons from the perspective of quantum mechanics and obtain branches of DB solitons whose bright components can have one or more humps [46,47].In Section 4.3, we shall present another application of exact solutions to scalar NLS by "lifting" them to codimension-1 families of exact solutions to Eq. ( 8) that can be numerically continued directly into vector solitons.
The pitchfork bifurcation point on ZN is given by Near ωZN A , we assume that ωA − ωZN A = O(ǫ 2 ), where 0 < ǫ ≪ 1, and use the expansion where [φ 0 , ψ 0 ] T is the perturbed ZN equilibrium, and [ φ, ψ] T is the ζ-dependent part that determines the spatial profile of the vector soliton.The pitchfork bifurcation is subcritical when d 1 g 3 (g 2 g 4 − g 1 g 3 ) < 0. In this case, the leading nontrivial order for φ is O(ǫ) and contains the even homoclinic orbit of the normal form, which is a bright component of the vector soliton: The pitchfork bifurcation is supercritical when d 1 g 3 (g 2 g 4 − g 1 g 3 ) > 0. In this case, the leading nontrivial order for φ is O(ǫ) and contains the odd homoclinic orbit of the normal form, which is a dark component of the vector soliton: Domain Walls and Vector Solitons in the Coupled Nonlinear Schrödinger Equation 15In both the subcritical and the supercritical cases, the leading nontrivial order for ψ is O(ǫ 2 ) and contains an even localized profile "induced" by φ: where ψ 0 = −ω B /g 3 .In the subcritical case, ψ is a dark component when g 3 g 4 > 0, and an antidark component when g 3 g 4 < 0. Thus, the former is a DB soliton, while the latter is a BAD soliton.In the supercritical case, ψ is an antidark component when g 3 g 4 > 0, and a dark component when g 3 g 4 < 0. Thus, the former is a DAD soliton, while the latter is a DD soliton.In the BEC setting, the existence of these four types of vector solitons is established from the perspective of quantum mechanics in Ref. [39], where the antidark concept is first introduced.In the physics community, DAD solitons are better known as magnetic solitons, which are actively studied theoretically [74,75] and experimentally [76].
The pitchfork bifurcation point on NZ is given by ωNZ = g 4 /g 1 .
Near ωNZ A , we can again assume ωA − ωNZ A = O(ǫ 2 ), use the expansion in Eq. ( 21), and obtain the four types of vector solitons similarly to the ωZN A case.These solutions are summarized next but not shown explicitly.The pitchfork bifurcation is subcritical when d 2 g 1 (g 2 g 4 − g 1 g 3 ) < 0. In this case, ψ is an O(ǫ) bright component that is even.The pitchfork bifurcation is supercritical when d 2 g 1 (g 2 g 4 − g 1 g 3 ) > 0. In this case, ψ is an O(ǫ) dark component that is odd.In both cases, φ is an O(ǫ 2 ) dark or antidark component that is even.In the subcritical case, φ is dark if g 1 g 2 > 0 and antidark if g 1 g 2 < 0. These are respectively a DB soliton and a BAD soliton.In the supercritical case, φ is antidark if g 1 g 2 > 0 and dark if g 1 g 2 < 0. These are respectively a DAD soliton and a DD soliton.Notably, these DD or DAD solitons bifurcating from ωNZ A may coexist with the DD or DAD solitons bifurcating from ωZN A , but these two branches are distinct since they have opposite parities in both components.
The pitchfork bifurcation is the fundamental bifurcation connecting the four equilibria as visualized in Fig. 3. Thus, after analyzing the codimension-1 bifurcations, it is also useful to identify the codimension-2 bifurcations.To identify the codimension-2 pitchfork bifurcation point for any equilibrium from Fig. 3, one can first identify the node representing this equilibrium and then equate the two codimension-1 points represented by the two edges meeting at this node.Since a codimension-2 point selects a triangle in Fig. 3 that contains a diagonal edge, it can provide valuable insights into the interplay between DWs and vector solitons.We shall identify the codimension-2 points but leave a detailed study as future work.
For ZZ, there are no nontrivial codimension-2 points.
For NN, Eqs.(20,25) imply that the codimension-2 point is g 2 g 4 −g 1 g 3 = 0.This is the miscible-immiscible transition point, which separates the miscible regime featuring DD solitons and the immiscible regime featuring symmetric DWs between ZN and NZ [38].
For ZN, Eq. ( 20) and ωZZ B = 0 imply that the codimension-2 point is g 3 = 0.For NZ, Eq. ( 25) and ωZZ A = 0 imply that the codimension-2 point is g 1 = 0.These are focusing-defocusing transition points, which separate the focusing regime featuring asymmetric DWs between ZZ and NN and the defocusing regime featuring DB solitons embedded in either ZN or NZ.

Exact solutions
As an application of exact solutions to scalar NLS including scalar bright and dark solitons, we lift them to exact solutions to Eq. ( 8) by requiring a scalar soliton multiplied by a constant amplitude vector to satisfy Eq. ( 8).
To find scalar bright solitons with vector amplitudes, we substitute the following ansatz into Eq.( 8): where A 1 , A 2 , W > 0 without loss of generality.This yields the following consistency condition: while the soliton parameters are To find scalar dark solitons with vector amplitudes, we substitute the following ansatz into Eq.( 8): where A 1 , A 2 , W > 0 without loss of generality.This yields the same consistency condition as the bright case given by Eq. ( 27), while the soliton parameters are Arguably, scalar bright or dark solitons with vector amplitudes are not true BB or DD vector solitons.However, they are reversible homoclinic orbits that happen to have exact analytical expressions when the consistency condition in Eq. ( 27) is satisfied.Once they are numerically continued such that Eq. ( 27) is no longer satisfied, they become true BB or DD vector solitons whose profiles cannot be exactly a sech or a tanh function.Thus, we shall consider the exact solutions themselves as true vector solitons and call them scalar BB and DD solitons.1. Types of vector solitons and their possible parities.Abbreviations: EE = even-even; EO = even-odd; OE = odd-even; OO = odd-odd.

Summary of real vector solitons
After identifying the three codimension-1 mechanisms for generating vector solitons in Eq. ( 8), we can group these vector solitons by their types and possible parities; see Table 1.Overall, a bright or an antidark component is always even, while a dark component can be either even or odd.Any codimension-1 family of vector solitons in Table 1 can be numerically continued into a codimension-0 family of vector solitons, namely a 2parameter family in (ω A , ωB ) or equivalently a 3-parameter family in (C g , ω A , ω B ) using the Galilean transformation in Eq. ( 6).These codimension-0 families of vector solitons may not be distinct since two vector solitons of the same type and parity but generated by different mechanisms could be parametrically related in Eq. ( 8).BAD, DAD, and Turing solitons can only be generated by local bifurcations, not heteroclinic cycles or exact solutions.Nonetheless, since Turing, DD, and DAD solitons share the NN background, a Turing soliton could be parametrically related to a DD or a DAD soliton via the so-called Belyakov-Devaney point, where the spatial eigenvalues given by Eq. ( 13) transition from a complex quartet of eigenvalues to four real eigenvalues by colliding pairwise on the real axis, or vice versa [72].
There are two types of BB solitons including extended and scalar BB solitons.Both are even-even, so they could be parametrically related.
There are two types of DB solitons including extended and weakly nonlinear DB solitons.The former can be even-even or odd-even, while the latter are even-even.Thus, even-even extended DB solitons could be parametrically related to weakly nonlinear DB solitons.
There are three types of DD solitons including extended, weakly nonlinear, and scalar DD solitons.Extended DD solitons can have all four possible parities.Weakly nonlinear DD solitons are either even-odd or odd-even, say even-odd.Scalar DD solitons are odd-odd.Thus, even-odd extended DD solitons could be parametrically related to weakly nonlinear DD solitons, while odd-odd extended DD solitons could be parametrically related to scalar DD solitons.

Complex vector solitons
Once we have a 3-parameter family of real homoclinic orbits in (C g , ω A , ω B ) in Eq. ( 8), we can remove the restriction in Eq. ( 4) and numerically continue these orbits in (k A , k B ) to obtain a 5-parameter family of complex homoclinic orbits in Eq. ( 2).Although these orbits are different solutions to the ODE in Eq. ( 2), they do not always yield different vector solitons in the PDE in Eq. (1).Consider a numerical continuation from k A to kA ; the k B case is similar.To identify which ODE solutions yield the same vector soliton in the PDE, we observe that the A component of a vector soliton and, recalling that Thus, the two solution families, φ at kA and φ at k A , yield the same family of vector solitons.
If φ is homoclinic to the N equilibrium, namely dark or antidark, then φ is homoclinic to a nonzero plane wave.Since the numerical continuation from k A to kA yields a solution family homoclinic to the N equilibrium, this solution family cannot be equivalent to φ, so this continuation must yield a different family of vector solitons.Notably, this continuation might connect two branches of real vector solitons with different parities.
However, if φ is homoclinic to the Z equilibrium, namely bright, then φ is also homoclinic to the Z equilibrium.Since the numerical continuation from k A to kA also yields a solution family homoclinic to the Z equilibrium, this solution family must be equivalent to φ with the equivalence given by Eqs.(29)(30), so this continuation does not yield a different family of vector solitons.
The nontrivial continuation in the wavenumber for a dark or antidark component, but not a bright component, is perhaps best known in the scalar NLS equation, where scalar bright solitons have two parameters including the group velocity and the frequency, while scalar dark solitons have three parameters including the group velocity, the frequency, and the wavenumber [21].The profile of a scalar dark soliton is a real tanh function plus an imaginary constant iν, where ν is the darkness parameter with ν = 0 for a black soliton and ν = 0 for a gray soliton.In our language, a black soliton is Domain Walls and Vector Solitons in the Coupled Nonlinear Schrödinger Equation 19a real scalar dark soliton, a gray soliton is a complex scalar dark soliton, and the former can be continued in the wavenumber to obtain the latter.
The continuation in the wavenumber of a dark or antidark component from real to complex can always be performed numerically, but a complex dark component may also be expressible analytically as exemplified by the scalar NLS equation.In the CNLS equation, one can similarly generalize a black component described by a real tanh function to a gray component by adding an imaginary constant.There are three types of real vector solitons in Table 1 with such a black component including weakly nonlinear DD, weakly nonlinear DAD, and scalar DD solitons.For the first two types, the O(ǫ) component in Eq. ( 23) can be generalized from black to gray similarly to scalar NLS, which will not be shown explicitly.For the third type, both components can be generalized from black to gray by adding an imaginary vector to Eq. ( 28): where A 1 , A 2 , W > 0 without loss of generality, and B 1 , B 2 ∈ R. Substituting this ansatz into Eq.( 2) with the following shorthand notations for its coefficients we obtain the following consistency condition while the soliton parameters are Thus, the 2-parameter family of black-black solitons in Eq. ( 28) can be generalized to the 4-parameter family of gray-gray solitons in Eq. (31).The two darkness parameters B 1 and B 2 become nonzero respectively when K A and K B become nonzero, so this generalization is essentially a continuation in (K A , K B ), or equivalently the two wavenumbers (k A , k B ).The general family of complex DD solitons have five parameters and can be obtained from this family of gray-gray solitons by numerical continuation in Ω B , or equivalently ω B , such that the above consistency condition is violated.Although a bright component cannot be non-trivially continued in the wavenumber like a dark or an antidark component, the latter cannot be non-trivially multiplied by a constant phase factor e iΘ A or e iΘ B like the former for vector soliton collisions.This is because in a collision between two vector soliton respectively located on the left and on the right, the left and right dark or antidark components must have matching phases, while the left and right bright components can have any phases.Moreover, the left and right bright components can also have different frequencies, so one can expect richer dynamics when a soliton collision involves at least one bright component.If the A component of a vector soliton is bright, then we define the polarization of this vector soliton as the phase difference ρ = Θ A − Θ B regardless of whether the B component is bright, dark, or antidark.

Numerical continuation
Since vector solitons are reversible homoclinic orbits, only half of the soliton profile needs to be numerically continued, say the right half.Thus, we can choose a domain ζ ∈ [0, L] with L large, such that ζ = 0 is the center of the soliton and ζ → L is the tail of the soliton.At ζ = L, we can impose Neumann boundary conditions such that the tail decays to the background equilibrium.At ζ = 0, we must impose boundary conditions required by the reversibility symmetry of the soliton in question.
To numerically continue a real vector soliton with the reversibility in Eq. ( 9), we must impose at ζ = 0 the boundary conditions in Eq. ( 19), which depend on the parities of the two components.As stated earlier, to convert a real vector soliton to a complex vector soliton with the reversibility R 0,0 in Eq. ( 3), any odd component must be multiplied by i.To numerically continue a complex vector soliton with the reversibility R 0,0 , we must impose at ζ = 0 the boundary conditions corresponding to the symmetric section fix(R 0,0 ) given explicitly by As an application of our classification program, we consider an example with the CNLS parameters d 1 = −1, d 2 = −3, g 1 = 2, g 2 = 4, g 3 = 5, and g 4 = 8.This parameter combination is in the defocusing regime relevant to BECs, and the immiscibility condition g 2 g 4 − g 1 g 3 > 0 is satisfied.Thus, the ZN and NZ equilibria can be both spatially hyperbolic, and a DW may form at the Maxwell point in Eq. ( 17).An example combination of the traveling wave parameters satisfying both conditions is (C g , ω A , ω B ) = (5.093,−9.9, −9.8).To compute a DW profile at the Maxwell point, we make Eq.( 8) the steady state of a dissipative PDE and evolve an initial step between ZN and NZ in this PDE to reach the steady state, but there are other numerical methods available.
Once a DW profile between ZN and NZ is computed, it can be treated as half of the profile of an extended DB soliton and numerically continued into a codimension-0 family of real vector solitons.As listed in Table 1, an extended DB soliton has two possible parities including even-even and odd-even.We shall focus on the former and leave the latter as future work.Thus, we impose even-even boundary conditions at ζ = 0, namely Eq. ( 19 continue the DW in one of the traveling wave parameters (C g , ω A , ω B ), say ω A .We use the software AUTO-07p [77] for all numerical continuations in this paper, but other software can also be used.The bifurcation diagram is shown in Fig. 4(a) with the profiles of the accompanying solutions shown in the insets.Note that although the full soliton profile is plotted, only half of the profile is numerically continued.The starting point is the extended DB soliton at A 1 whose profile is shown as φ A and ψ A and consists of two DWs.As ω A varies, the two DWs "merge" into a generic real DB soliton at B 1 whose profile is shown as φ B and ψ B .Although not shown explicitly in Fig. 4(a), this branch of even-even DB solitons terminates at the pitchfork bifurcation point ω ZN A on the ZN equilibrium.
As discussed in Section 4.5, a 3-parameter family of real DB solitons in (C g , ω A , ω B ) can be numerically continued in the wavenumber of the dark component k B into a 4parameter family of complex DB solitons.A key motivation to compute complex vector solitons is to study interactions (collisions) involving such solitons, during which the interacting solitons must share the same plane wave background.Thus, we impose reversible boundary conditions in Eq. ( 32) at ζ = 0, Neumann boundary conditions at ζ = L, and numerically continue complex DB solitons in C g with k B fixed rather than in k B with C g fixed.The bifurcation diagram is shown in Fig. 4(b) with the profiles of the accompanying solutions shown in the insets.The starting point is the real DB soliton at B 1 in Fig. 4(a).As C g increases, one obtains a complex DB soliton at A 2 whose profile is shown as φ A and ψ A .As C g decreases, one obtains a complex DB soliton at B 2 whose profile is shown as φ B and ψ B .The bifurcation diagram is symmetric with respect to B 1 and the solution profiles at A 2 and B 2 are related by complex conjugation.This is expected since Eq. ( 2) to the left and to the right of B 1 differ only in the sign of the coefficient of ψ ′ .
Extended DB solitons formed by two DWs are well studied in the optics setting where they are often called domain wall solitons [56,57,58,59], and in the BEC setting where they are sometimes called bubble droplets [69].The above example shows that DWs and DB solitons can be parametrically related in the immiscible regime of the asymmetric CNLS equation relevant to BECs.An interesting prospect is to deduce analytically the properties of quasi-extended DB solitons near the Maxwell point from the properties of the underlying DW.

Collisions between two polarized dark-bright solitons
In Section 5, we have obtained a 4-parameter family of complex DB solitons in (C , ω A , ω B , k B ) with the B component being dark and the A component being bright.As discussed in Section 4.5, the bright component A can be multiplied by a constant phase factor e iΘ A to yield a 5-parameter family of polarized DB solitons in (C g , ω A , ω B , Θ A , k B ) for soliton collisions.In such collisions, the 2 parameters (ω B , k B ) must be the same between the two colliding solitons such that the solitons share the same plane wave background, but the other 3 parameters (C g , ω A , Θ A ) can be different.Among these, C g must be different for the collision to occur, ω A will be taken as the same for simplicity, and we shall focus on the effect of different Θ A .We shall denote Θ A as the polarization ρ recalling that ρ = Θ A − Θ B and Θ B = 0 for the dark component B. In the BEC setting, it is well known that collisions between two polarized DB solitons can yield a mass exchange between the two bright components [37].We shall examine this mass exchange in the immiscible regime through a numerical example.
To reduce radiation for soliton collisions, we consider an example with the CNLS parameters d 1 = d 2 = −1, g 1 = g 2 = g 3 = 2, and g 4 = 3, which differs from the defocusing Manakov system only by adding 1 to g 4 .Apart from near integrability, this example and the example in Section 5 have similar bifurcation structures since both belong to the immiscible regime.A DW exists at the Maxwell point (C g , ω A , ω B ) = (5.395,−9.5, −10).A real DB soliton can be found by numerical continuation of the DW to ω A = −9.7 similarly to Fig. 4(a).A family of complex DB solitons sharing the same plane wave background can be obtained by numerical continuation of the real DB soliton in C g similarly to Fig. 4(b).These two bifurcation diagrams will not be shown explicitly since they are similar to Fig. 4(a,b).We shall call the group velocity of the real DB soliton C g ) = (5.789,5).Thus, these two solitons have opposite velocities in the frame comoving with velocity C (0) g , and the two soliton profiles are related by complex conjugation.
The initial conditions to be time-evolved in Eq. (1) consist of these two complex g , C DB solitons respectively located on the left and on the right.The left bright component is multiplied by e iρ L , and the right bright component is multiplied by e iρ R , where ρ L,R ∈ [0, 2π) are the polarizations of the two colliding solitons.We can set ρ R = 0 without loss of generality and explore the effect of varying ρ L .The collision dynamics for two choices of ρ L are shown in Fig. 5.As expected, the ρ L = 0 case yields no mass exchange as shown in Fig. 5(a).To explore ρ L = 0, we choose ρ L = π/2, and observe a mass exchange as shown in Fig. 5(b).Note that both collisions produce negligible radiation compared with the solitons, which is consistent with the near integrability of Eq. (1).
We shall examine this mass exchange in terms of both soliton parameters and soliton profiles.
In terms of soliton parameters, the 2 parameters of the dark component, (ω B , k B ), cannot change since the incoming and outgoing solitons must share the same plane wave background.However, the soliton velocity C g , which is a shared parameter between both components, can change as shown in Fig. 5(b).The 2 parameters of the bright component, (ω A , ρ), can also change.
In terms of soliton profiles, the profiles of both components can change.As shown in Fig. 5(c), there is a visible difference in the amplitudes of both components after the polarized collision in Fig. 5(b) (blue) compared with either after the unpolarized collision in Fig. 5(a) (red) or before collisions (black-dashed).Thus, in terms of soliton profiles alone, it is tempting to describe the mass exchange as an energy transfer, but the latter usually applies between two bright components, not between a bright component and a dark one.
Overall, the change in soliton profiles in both components results from the change in the 2 soliton parameters (C g , ω A ), as exemplified by Fig. 4. The polarizations ρ L,R of the two solitons can also change, but they do not affect the soliton profiles.Thus, the mass exchange is better described by the change in soliton parameters than the change in soliton profiles, although the soliton profiles are more visible than the soliton parameters.
We stress that Fig. 5(b) shows a collision between two generic DB solitons whose dark and bright components have comparable amplitudes.In such cases, the relative changes of the amplitudes of both components produced by the mass exchange are moderate as shown in Fig. 5(c).However, if either or both of the two colliding solitons are near bifurcation points, then their dark and bright components could have very different amplitudes.In such cases, the relative changes of the amplitudes of either or both components produced by the mass exchange could be very large or very small.The eigenfunctions in the spectra of vector solitons can yield valuable insights into the collision dynamics of these solitons in non-integrable regimes of Eq. ( 1); see e.g.Ref. [33] for the collision dynamics of BB solitons explained from this perspective.

Conclusion
In this paper, we have outlined a program to classify DWs and vector solitons in the 1D two-component CNLS equation without restricting the signs or magnitudes of any coefficients.The CNLS equation is reduced first to a complex ODE, and then to a real ODE after imposing a restriction.In the real ODE, we have identified four possible equilibria including ZZ, ZN, NZ, and NN, and analyzed their spatial stability.We have identified two types of DWs including asymmetric DWs between ZZ and NN and symmetric DWs between ZN and NZ.We have identified three mechanisms for generating vector solitons in the real ODE including heteroclinic cycles formed by assembling two DWs back-to-back, local bifurcations on the four equilibria, and exact solutions that are scalar solitons with vector amplitudes.Local bifurcations are further divided into the Turing bifurcation that generate Turing solitons and the pitchfork bifurcation that generates DB, BAD, DD, and DAD solitons.We have identified a 3D parameter space for real vector solitons and a 5D parameter space for complex vector solitons, and introduced a numerical continuation method to find these vector solitons.Our classification program applies not only to optics and BECs, but also to the general problem of two interacting quasi-monochromatic plane waves.Further solutions that can be incorporated into our classification program include those originating from secondary bifurcations or codimension-2 bifurcations.
As an application of our classification program, we have shown that real and complex DB solitons can be parametrically related to DWs in the immiscible regime of the CNLS equation relevant to BECs.We have also shown that in this regime, collisions between two polarized DB solitons typically feature a mass exchange that changes the frequencies and polarizations of the two bright components and the two soliton velocities.
Our classification program can be generalized to the N-components CNLS equation for N interacting quasi-monochromatic plane waves.This equation overlaps with the N-component GPE for bosonic BECs [36] in certain parameter regimes.This equation has 2 N possible equilibria, so its complexity increases exponentially with N.
Our classification program is also relevant to the 2D two-component CNLS equation.In the BEC setting, the 2D GPE can exhibit radial solutions like vortices due to the radial symmetry of the Laplacian [36].However, for two interacting quasimonochromatic plane waves in 2D, the dispersion relation of either plane wave can be either elliptic or hyperbolic, so the linear part of the 2D CNLS equation generally lacks radial symmetry.Nonetheless, the 2D CNLS equation can be reduced along any line to the 1D CNLS equation that we studied, so 1D solitons in the latter correspond to 2D line solitons in the former.
We have applied our classification program to DWs and DB solitons in the homogeneous GPE, but the effects of a trapping potential on such solutions can also be remarkable.Stable pinning of DWs near maxima of the potential is shown both physically [65] and mathematically [68].A plethora of DB solitons are found with a harmonic trap using numerical continuation in recent studies [78,79].The GPE in the immiscible regime with a trapping potential continues to provide a great venue for further explorations of DWs and DB solitons.
only provides a preliminary example of asymmetric DWs, and we cannot rule out the possible existence of stable asymmetric DWs with other parameter choices.A systematic study of asymmetric DWs and related solutions is left as future work.

21 AFigure 4 .
Figure 4. (Color online) Bifurcation diagrams and the profiles of the accompanying solutions for numerical continuations in ω A shown in panel (a), followed by C g shown in panel (b).The solution profiles are shown in terms of their real (blue-solid) and imaginary (red-dashed) parts.The CNLS parameters are d 1 = −1, d 2 = −3, g 1 = 2, g 2 = 4, g 3 = 5, and g 4 = 8, and the traveling wave parameters at the starting point A 1 in panel (a) is (C g , ω A , ω B ) = (5.093,−9.9, −9.8).(a) The profile at A 1 is shown as φ A and ψ A , and the profile at B 1 is shown as φ B and ψ B .(b) The profile at A 2 is shown as φ A and ψ A , and the profile at B 2 is shown as φ B and ψ B .

= 5 .
395.For soliton collisions, we shall choose two complex DB solitons with C g values symmetric with respect to C (0) g similarly to Fig. 4(b), which are respectively called (C (1) g , C

g ) = ( 5 .
789, 5), and the polarization of the right soliton is ρ R = 0. (a, b) Space-time plots in the frame comoving with velocity C (0) g = 5.395.The polarization of the left soliton is ρ L = 0 in panel (a) and ρ L = π/2 in panel (b).(c) Soliton profiles after the unpolarized collision in panel (a) (red), after the polarized collision in panel (b) (blue), and before collisions (black-dashed).The left soliton is always shown in the left half-plane, x < 0, and the right soliton is always shown in the right half-plane, x > 0.

Table
Domain Walls and Vector Solitons in the Coupled Nonlinear Schrödinger Equation17