The multi-component Bose–Einstein condensate (BEC) is an ideal platform for studying the dynamics of vector solitons due to the advantages of macroscopic visibility and tunable atomic interactions [1, 2]. Many different vector solitons have been obtained in multi-component coupled BEC systems, such as bright-bright solitons [3–6], bright-dark solitons [7–9], dark-bright solitons [10–14], dark-dark solitons [15–17], bright-bright-dark solitons [9, 18] and dark-dark-bright solitons [18]. These studies on solitons in multi-component BECs are usually based on the integrable systems, such as the Manakov model [19] with equal inter- and intra-species interactions, and the coupled model with spin-dependent interactions [20–22]. Furthermore, in the non-integrable systems, in which interaction strengths are unequal, some vector soliton solutions can be obtained theoretically [23–27]. All the above-mentioned solitons are mainly mass solitons, due to the fact that their total mass densities have localized profiles. Multi-component BECs can be considered as pseudo-spin systems with different spin values [28–30]. In such systems, a localized spin polarization (the magnetic soliton) was the focus of both theoretical [24, 25] and experimental [31, 32] studies, due to its advantages in investigating solitary waves in non-Manakov cases and observing pure spin currents. Very recently, a spin soliton for the spin-1/2 system has been studied. The spin soliton's total mass density is also uniform and a soliton profile is shown on the spin density [33]. The spin solitons can admit spin-balanced or spin-imbalanced density backgrounds, thus having potential in the study of spin currents. Particularly, this type of soliton exhibits different dynamics under the effect of external forces compared with well-known vector solitons for the integrable cases. Spin soliton solutions for the spin-1 system have also been discovered [34]. These solutions are obtained based on the vector soliton solutions for the integrable Manakov model, and thus there are many constraints on the nonlinearity coefficients. These constraints can be relaxed and the variational method can provide much more exact analytical soliton solutions for two-component models with different nonlinearity coefficients [26]. Although the dark-bright soliton solutions for the three-component model are reported in [26], we employ a decoupling method to obtain spin soliton solutions more directly. We expect to reduce the constraints by using these soliton solutions for non-integrable models and obtain more general spin solitons in spin-1 systems.

In this paper, we utilize the soliton solutions in [26] to reduce the constraints in [34], and therefore extend the existence region of spin soliton solutions in nonlinearity coefficient space. We also perform spectral stability analysis and numerical simulations, and show that the spin solitons that admit spectral stability with discrete spectra are stable, while those with continuous spectra are metastable against weak noise. The existence region of spin solitons further enriches the family of solitons in BECs and provides a reference for investigations on non-integrable phenomena and experiments.

The structure of this paper is as follows. In section 2, we give six spin soliton solutions as well as their constraints on nonlinear interaction strengths, and we carefully depict their properties and differences. In section 3, we calculate the Bogoliubov–de Gennes excitation spectra of these spin solitons and test their stability using numerical simulations. In section 4, we investigate the collisions between two spin solitons. When the speed is much lower than the upper limit, the soliton exhibits shape oscillation and slight acceleration after the collision. Finally, in section 5, we summarize and discuss the results.

Under the mean-field approximation, the dynamics of a quasi-one-dimensional ⁸⁷Rb BEC can be described by the following coupled Gross–Pitaevskii equations [20–22]: where ψ_i (i = 1, 2, 3) represents the wave function of the ith component and the symbol ‘*' indicates the complex conjugate, x and t respectively represent the spatial coordinate and time evolution, and V_i(x) denotes the weak axial trapping potential. The coefficients g_ii (g_ij) are the intra- (inter-) species interaction strengths given by the s-wave scattering length. Since g_s/∣g_ij∣ is very small for ⁸⁷Rb [35, 36], the spin-dependent terms can be safely ignored.

Besides, the weak trapping potentials can be set to V_i = 0(i = 1, 2, 3) and the coupled Gross–Pitaevskii equations are reduced to coupled nonlinear Schrödinger equations (NLSEs). If g_ij are all equal, the equations will become the well-known Manakov model [19], which is completely integrable and can be solved exactly by the Darboux transformation [37, 38], Hirota bilinear method [39] and inverse-scattering method [40]. The Manakov case was the focus of previous theoretical studies; therefore, many different vector solitons have been obtained [3–18]. The realization of Feshbach resonance techniques enables the manipulation of the scattering length between three hyperfine states of atoms. In this case, the model is no longer integrable and its exact solutions cannot be obtained by the above-mentioned methods analytically. Very recently, a decoupling method was suggested and exact soliton solutions for the spin-1/2 and spin-1 BECs were derived. These solitons are not mass solitons due to their uniform mass densities (∑∣ψ_i∣²=1 for simplicity), and this character enables the reduced three-component system to still be seen as a pseudo-spin system. In particular, the pseudo-spin densities $F_z={\mathrm{\Psi }}^{†}{\hat{F}}_z\mathrm{\Psi }=|{\psi }_1{|}^2-|{\psi }_3{|}^2$⁴(⁴The spin density $F_z={\mathrm{\Psi }}^{†}{\hat{F}}_z\mathrm{\Psi }=|{\psi }_1{|}^2-|{\psi }_3{|}^2$, where $\mathrm{\Psi }={\left({\psi }_1{\psi }_2{\psi }_3\right)}^{\mathrm{T}}$ and the symbol ‘†' denotes the conjugate transpose. The spin-1 matrices, ${\hat{F}}_x=\frac{1}{\sqrt{2}}\left(\begin{array}{c}010\\ 101\\ 010\end{array}\right)$, ${\hat{F}}_y=\frac{\mathrm{i}}{\sqrt{2}}\left(\begin{array}{c}0-10\\ 10-1\\ 010\end{array}\right)$, and ${\hat{F}}_z=\frac{1}{\sqrt{2}}\left(\begin{array}{c}100\\ 000\\ 00-1\end{array}\right)$.) (where $\mathrm{\Psi }={\left({\psi }_1{\psi }_2{\psi }_3\right)}^T$ and the symbol ‘†' denotes the conjugate transpose) have soliton profiles, and therefore these solitons are called spin solitons.

Three-component coupled NLSEs can be decoupled into a two-component model with one scalar NLSE. The previously reported spin solitons were obtained when the two-component model was integrable, for which the exact soliton solutions can be derived analytically [34]. Recent research indicates that the exact analytical soliton solutions can be obtained for the non-integrable model [26]. Based on this, we can extend the coefficient space for spin soliton solutions, since the constraints on nonlinearity coefficients can be greatly relaxed. Now, one can obtain exact spin soliton solutions as long as the nonlinearity coefficients satisfy G₁ = g₁₂ = g₂₁ and G₂ = g₁₃ = g₃₁ = g₂₃ = g₃₂. According to the restriction on the uniform total mass density, we find that g₃₃ is uniquely determined by the other coefficients:

For the non-integrable two-component model, there are four types of exact analytical soliton solutions, namely, bright-bright solitons, bright-dark solitons, dark-bright solitons and dark-dark solitons. Each of them exist in two separated regions with different constraint conditions [26]. When we go back to the original three-component equations, six types of spin soliton solutions can be obtained and they also have two different parts in the coefficient space. There is no essential difference between the same type of solitons in the two parts, and they are labeled as cases ‘−1' and ‘−2' in figure 1. In the subsequent discussions, we will mainly focus on case-1 solitons. The distribution map for the six spin soliton solutions are shown in figure 1 according to different constraint conditions on nonlinearity coefficients.

It should be emphasized that the existence of spin solitons directly depends on the difference between the nonlinearity coefficients F_i (i = 1, 2, 3, 4) rather than themselves. Here, F_i is the difference between G₂ and other nonlinearity coefficients, where F₁ = G₂ − g₁₁, F₂ = G₂ − g₂₂, F₃ = G₂ − g₃₃ and F₄ = G₂ − G₁. Hence, spin soliton solutions obviously differ from solitons for the integrable model. If F_i=0, the model becomes the Manakov model, and spin solitons cannot exist. Specifically, spin soliton solutions cannot exist on any surface in figure 1. We present the specific classification of spin solitons based on F_i in table 1. Next, we give the generalised six spin soliton solutions, and carefully analyze their characteristics and differences.

In general, the spin solitons can be classified based on the soliton structures in the three components. Firstly, we show spin soliton solutions with two bright soliton components and one dark soliton component, which are a bright–bright–dark spin soliton (BBD), a bright–dark–bright spin soliton (BDB) and a dark–bright–bright spin soliton (DBB). The BBD solution can be written as: where v and $w_1=\sqrt{F_3-v^2}$ are the velocity and inverse width of the soliton. It can be seen that the velocity v must be less than $\sqrt{F_3}$, where F₃ is the difference between intra-species and inter-species interaction strengths of the ψ₃ component. This indicates that the spin soliton differs from the vector soliton for the Manakov model, for which the soliton's maximum speed depends on the nonlinear interaction strength. Since the total mass density is uniform for spin solitons, the width and amplitude of the bright soliton directly depend on its velocity, which is inconsistent with that of the vector solitons [3–18, 41–47]. As an example, the mass density and spin density profiles for a static BBD are shown in figure 2(a). One can find that the background of the spin density is −1, while the maximum spin density can be varied within the region (−1, 1] by changing v and F_i.

Next, the BDB solution is given as follows: where $w_2=\sqrt{(F_1{F}_2-F_4^2)/(F_1-F_4)-v^2}$ is the inverse width of this type of spin soliton. It is easy to get a speed limit of $v<\sqrt{(F_1{F}_2-F_4^2)/(F_1-F_4)}$, which is evidently greater than that of BBD. In addition, BDB has a concave structure in the spin density distribution. Its spin density background is zero and the maximum spin density can be changed within (−1, 0], similar to the magnetic soliton [24, 25]. Figure 2(b) shows the mass density and spin density of a BDB with v = 0.

The exact expression for DBB can be written as follows: where $w_3=\sqrt{(F_1{F}_2-F_4^2)/(F_2-F_4)-v^2}$. For DBB, the speed limit can be derived as $v<\sqrt{(F_1{F}_2-F_4^2)/(F_2-F_4)}$. As an example, we show one DBB with zero velocity in figure 2(c); its spin density shows a concave structure, where the background is 1 and the minimum spin density can be changed within [1, − 1). We emphasize that the above three types of spin solitons cannot be obtained by gauge transformations from each other, due to their different maximum speeds and spin distributions.

When the spin soliton admits two dark soliton components together with one bright soliton component, another three types can be obtained. They are dark–dark–bright spin soliton (DDB), bright–dark–dark spin soliton (BDD) and dark–bright–dark spin soliton (DBD). The exact solution for DDB can be given as: where $w_4=\sqrt{-F_3-v^2}$, and the maximum speed for the DDB depends on the difference between the nonlinearity coefficients F₃, namely, $v<\sqrt{-F_3}$. The mass density and spin density of a static DDB is shown in figure 2(d). Similarly, one can find that for the DDB, the spin density background can be varied in the region of [0, 1], and the minimum can be changed within [−1, 1).

The exact solution for BDD can be given as: where $w_5=\sqrt{(F_4^2-F_1{F}_2)/(F_2-F_4)-v^2}$. Here, v still denotes the soliton velocity, with the maximum speed being determined by the difference between the nonlinearity coefficients, i.e. $v<\sqrt{(F_4^2-F_1{F}_2)/(F_2-F_4)}$. In contrast to the fact that the spin solitons possess two bright soliton components, the spin density background for BDD can be changed within [1, 0) and the maximum value can be varied within (0, 1] by changing v and F_i. As an example, we show the mass density and spin density of a static BDD in figure 2(e).

Lastly, the exact solution for DBD can be given as: where $w_6=\sqrt{(F_4^2-F_1{F}_2)/(F_1-F_4)-v^2}$. It can be directly solved that the maximum spin soliton speed limit is $v<\sqrt{(F_4^2-F_1{F}_2)/(F_1-F_4)}$. In particular, DBD's spin density is a convex or concave structure with different nonlinearity coefficients. If G₁ > G₂, the spin density is a concave structure. If G₁ < G₂, the spin density is a convex structure. In particular, when G₁ = G₂, one can find that ∣ψ₁∣² = ∣ψ₃∣², and thus the spin distribution of DBD is zero. This unique soliton structure is displayed in figure 2(f). Since the soliton possesses a uniform spin density only under this particular condition, we still consider DBD as a spin soliton in this paper. Although BDD, DBD and DDB have similar expressions, their maximum velocity limit and spin density are different; therefore, they are completely different spin solitons and cannot be converted through any trivial gauge transformations. The phases of dark soliton components for spin solitons are similar to that of the scalar dark soliton. They also admit a π phase jump when they are static, and the phase difference tends to zero with the soliton velocity approaching their different maximum speeds.

It is worth noting that the methodology employed in this study differs from [26], which utilizes a variational method. By appropriately choosing trial wave functions, [26] obtained exact soliton solutions in dark-bright form. Obtaining spin solitons using the decoupling method is much more direct and convenient since the preassigned trial wave functions or integrations are abandoned. Therefore, the decoupling method provides a more direct way to obtain spin soliton solutions and can be conveniently extended to higher spin cases. Besides, existence regions for different structures can be classified more clearly based on the decoupling method. The spin soliton solutions in [34] exist in the surface that admits F₂ = F₄. Last, but not least, it is worth mentioning that the integrable cases have been widely studied and various exact soliton solutions have been obtained, based on methods such as the inverse-scattering method [48], the Hirota bilinear method [49] and the Darboux transformation method [50, 51]. A series of multi-soliton solutions can be provided by multiple iterations, but there were strong constraints on the coefficients. In contrast, the decoupling method can provide exact solutions (for example, the spin soliton) in considerably extended coefficient space. However, it is quite difficult to obtain multi-soliton solutions in this way.

In the previous section, we gave exact solutions for six types of spin solitons. The stability of solitons is crucial for their observation in experiments. Hence, we need to check the stability of the solitons we obtained. In this section, we will test the stability of these spin solitons to find spin solitons that can be prepared experimentally. The Bogoliubov–de Gennes excitation spectrum is one of the commonly used methods for testing the stability of solitons. We introduce a small perturbation on the spin soliton, $\mathrm{\Psi }$_P = $\mathrm{\Psi }$ + εe^−iμt, where $ϵ=\mathit{P}{\mathrm{e}}^{\mathrm{i}\lambda t}+{\mathit{Q}}^{\ast }{\mathrm{e}}^{-\mathrm{i}{\lambda }^{\ast }t}$ represents a weak perturbation, λ = λ_R + iλ_I is the eigenvalue, and P=(f₁₊, f₂₊, f₃₊)^T, Q^*=(f₁₋^*, f₂₋^*, f₃₋^*)^T. By substituting ψ_P into equation (1) and ignoring the high-order terms of the perturbation ε, we can solve an eigenvalue problem for the eigenvector (f₁₊, f₂₊, f₃₊, f₁₋, f₂₋, f₃₋)^T. Generally, if the eigenvalue λ is purely real, it indicates that the spin soliton is stable against small perturbations. Otherwise, if the eigenvalue of λ has a non-zero imaginary part, it indicates that the soliton is unstable and will be destroyed by small perturbations that grow exponentially. The calculated excitation spectra are shown in figure 3. One can find that the BBD, BDB and DBB only exist as discrete spectra without imaginary parts; therefore, they enable spectrum stability and can be stable over a long time evolution. By contrast, the DDB possesses non-zero imaginary parts, indicating that it will be unstable under the effects of small perturbations. We also notice that if G₁ ≠ G₂, the excitation spectrum of the BDD and DBD also admit non-zero imaginary parts. Nevertheless, when we take G₁ = G₂, the BDD and DBD will admit continuous spectra without imaginary parts; thus, it seems that these solitons are still stable.

The results of spectral stability analysis indicate that there are three types of spectra for the spin solitons, i.e. the discrete spectrum with a zero imaginary part, the discrete spectrum with a non-zero imaginary part and the continuous spectrum with a zero imaginary part. Next, we choose BBD, DDB and DBD as representatives to perform numerical simulations of equation (1) (with g_s = 0.02) and observe the soliton dynamics. All the parameters are the same as we used in figures 2 and 3. The initial states are given by the exact solution we presented, and 1% random noise is added to each component. The results are shown in figure 4. It can be seen that BBD remains stable over a long evolutionary time (over 600 dimensionless units). However, not all spin solitons are stable during evolution. The Bogoliubov–de Gennes excitation spectrum indicates that DDB is dynamically unstable, since it admits non-zero imaginary parts. These soliton structures will be destroyed rapidly, even in the absence of random noise. If G₁ ≠ G₂, the spectra for DBD (or BDD) also exhibit non-zero imaginary parts and they are dynamically unstable. Besides, when G₁ = G₂, we notice that they have a continuous excitation spectrum without an imaginary part but sudden unstable behavior occurs. As shown in figure 4(c), the soliton is stable for about 100 dimensionless units and then the system abruptly becomes unstable. This sudden unstable behavior is quite different from that with imaginary parts on the excitation spectrum. Because DBD can remain stable for a long time but finally becomes unstable, we refer to these special dynamics of the spin soliton as metastable. This result may indicate that spectral stability is not a sufficient condition for soliton stability during time evolution.

In this section, we investigate the collision between two spin solitons. Due to the non-integrability of the model, exact two-spin soliton solutions cannot be obtained using the analytical methods [37–40]. Therefore, we study the collision process using numerical simulations. As an example, we demonstrate the collision between two BBDs. The initial state can be given by a linear superposition of two spin soliton solutions when they are placed far from each other,

The parameters ${\kappa }_l=\sqrt{F_3-v_l^2}$ (l = 1,2) represent the inverse widths of two spin solitons, and x_l0 determine the initial positions of two solitons. The Δφ denotes the relative phase between bright solitons. By adjusting v_l, one can investigate collisions between spin solitons with different velocities.

In figures 5(a1)–(a3), we illustrate one case with v₁ = − v₂ = 0.1 and Δφ = 0. When the soliton velocity is much lower than the maximum speed, mass transport between bright solitons leads to changes in soliton profiles after the collision. In figures 5(b1)–(b3), we present the case where Δφ = π/2, while the velocities are retained. The relative phase between two bright solitons significantly changes the density redistribution during the collision. This phenomenon is similar to that reported in [34]. Interestingly, spin solitons exhibit weak amplitude and positional oscillations and accelerating behavior after collisions in both cases. The shape oscillation arises due to the presence of pairs of internal modes in the spectrum of spin solitons in the non-integrable model, and it causes resonant energy radiation [52, 53]. As a consequence, the soliton energy decreases, leading to the acceleration of spin solitons. Since the oscillation frequency and energy radiation are much lower, this phenomenon was not observed in [34].

When the speeds of two spin solitons approach the maximum speed, the collisions between them are almost elastic and the relative phase only changes the interference patterns in the collision process, which is similar to that reported in [34]. Numerical simulations show that the collisions for the other two types of spin solitons (BDB and DBB) are similar to those of BBD. During the collisions, the spin density undergoes redistribution and the non-integrability of the system induces weak oscillations and an acceleration phenomenon after collisions.

We obtain six types of exact analytical spin soliton solutions that exist in the extended region of nonlinearity coefficient space. These results further enrich the soliton solutions for the spin-1 BEC systems. By analysing the Bogoliubov–de Gennes excitation spectra and performing numerical simulations of these spin solitons, three types of spin solitons are confirmed to be stable against small deviations caused by spin-dependent terms and weak noise. In particular, two types of spin solitons are shown to be spectrally stable but become unstable over a rather long time period. We therefore refer to this novel case as metastable. The collisions between spin solitons directly exhibit the non-integrability of the model. Spin solitons start oscillation after the collision, and the energy radiation caused by the oscillation results in the acceleration of solitons. These results further enrich the soliton family for non-integrable models and can provide theoretical references for the experimental observation of spin solitons.

The coupled Gross–Pitaevskii equations (equation (1)) can be reduced to the following equations when the external potentials are V_i(x) = 0 and the weak spin-dependent effects are ignored,

By introducing the constraints on total density for spin solitons (∑∣ψ_i∣² = 1) and the nonlinearity coefficients (g₁₃ = g₂₃), equation (A1) can be written as follows: where G₁ = g₁₂ = g₂₁ and G₂ = g₁₃ = g₃₁ = g₂₃ = g₃₂. By taking ${\psi }_i=q_i{\mathrm{e}}^{\mathrm{i}{G}_{2}t}$, equation (A2) can be further simplified to: where F₁ = G₂ − g₁₁, F₂ = G₂ − g₂₂, F₃ = G₂ − g₃₃ and F₄ = G₂ − G₁. By employing this approach, the original model is simplified into the form that consists of two-component coupled nonlinear Schrödinger equations (NLSEs) and a scalar NLSE. Since the exact soliton solutions for the non-integrable two-component NLSE (depart from the integrable limit F₁ = F₂ = F₄) were obtained recently [26], we can greatly extend the existence region for the spin solitons [34]. When F₃ takes a positive(negative) value, the third component admits a dark(bright) soliton solution. Considering the constraint on the total mass density ∑∣q_i∣² = 1, the mass density of the two-component coupled NLSE should admit a density hump (dip). Based on the soliton solutions in [26], the spin soliton solutions can be constructed.

For example, we present how to obtain BBD. If F₃ > 0, ∣q₃∣² is a dark soliton:To keep the total mass density uniform (∑∣q_i∣²=1), the two coupled components (∣q₁∣² + ∣q₂∣²) should be a ‘bright soliton'. Therefore, we set F₁ > F₄, F₂ > F₄, and the coupled q₁ and q₂ admit a bright-bright soliton [26]: where $w=\sqrt{(F_2-F_4)/(F_4^2-F_1{F}_2)}{f}_1$, and $f_2=\sqrt{(F_4-F_1)/(F_4-F_2)}{f}_1$. According to the restriction on the total mass density, one can obtain the conditions of $\sqrt{F_3-v^2}=\sqrt{(F_2-F_4)/(F_4^2-F_1{F}_2)}{f}_1$ and $f_1^2+f_2^2=(F_3-v^2)/F_3$ and then get the expression of f₁ and f₂. Therefore, we can obtain the solution for equation (A3): where $w=\sqrt{F_3-v^2}$. Finally, by taking ${\psi }_i=q_i{\mathrm{e}}^{\mathrm{i}{G}_{2}t}$, the BBD solution (equations (3)–(5)) of equation (A1) is obtained. Other types of spin solitons can also be derived by this method.